Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22359)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 22360)
@@ -91,4 +91,19 @@
 	#if !defined(_HAVE_ADOLC_)
 	if(amr_frequency){
+		/*Verifications. AMR supports SSA, P1 and horizontal 2D domain*/
+		bool isSSA;
+		int domaintype,element_type,analysis_counter=-1;
+		this->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
+		this->parameters->FindParam(&domaintype,DomainTypeEnum);
+		for(int i=0;i<this->nummodels;i++) if(this->analysis_type_list[i]==StressbalanceAnalysisEnum){analysis_counter=i;break;}
+		if(analysis_counter==-1) _error_("Could not find alias for analysis_type StressbalanceAnalysisEnum in list of FemModel analyses\n");
+		for(int i=0;i<this->elements->Size();i++){
+			Element* element	= xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+			element_type		= element->element_type_list[analysis_counter];
+			if(element_type!=P1Enum) _error_("Element type "<<EnumToStringx(element_type)<<" not supported with AMR yet!\n");
+		}
+		if(!isSSA) _error_("Flow equation not supported with AMR yet!\n ");
+		if(domaintype!=Domain2DhorizontalEnum) _error_("Domain "<<EnumToStringx(domaintype)<<" not supported with AMR yet!\n");
+
 		this->parameters->FindParam(&amrtype,AmrTypeEnum);
 		switch(amrtype){
@@ -2456,5 +2471,5 @@
 		default: _error_("not implemented yet");
 	}
-
+	
 	/*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
 	this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
@@ -2477,8 +2492,8 @@
 	/*Creating nodes and constraints*/
 	/*Just SSA (2D) and P1 in this version*/
-	Nodes* new_nodes=new Nodes();
-	Constraints* new_constraints=new Constraints();
-
-	int nodecounter=0;
+	Nodes* new_nodes					= new Nodes();
+	Constraints* new_constraints	= new Constraints();
+
+	int nodecounter		=0;
 	int constraintcounter=0;
 	for(int i=0;i<this->nummodels;i++){//create nodes for each analysis in analysis_type_list
@@ -2487,9 +2502,8 @@
 
 		/*As the domain is 2D, it is not necessary to create nodes for this analysis*/
-		/*itapopo must verify if domain is not 3D. Only 2D in this version!*/
-		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;	    
+		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;    
 
 		this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
-		if(analysis_enum==StressbalanceAnalysisEnum) this->CreateConstraints(newnumberofvertices,newnumberofelements,nodecounter,constraintcounter,newx,newy,my_vertices,new_constraints);
+		this->CreateConstraints(new_vertices,nodecounter,constraintcounter,analysis_enum,new_constraints);
 		this->UpdateElements(newnumberofelements,newelementslist,my_elements,nodecounter,i,new_elements);
 
@@ -2660,32 +2674,42 @@
 }
 /*}}}*/
-void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
-
-	int maxinputs = MaximumNumberOfDefinitionsEnum;
+void FemModel::GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp){/*{{{*/
+
+	int maxinputs										= MaximumNumberOfDefinitionsEnum;
+	int numberofvertices								= this->vertices->NumberOfVertices();
+	int numberofelements								= this->elements->NumberOfElements();
+	int elementswidth									= this->GetElementsWidth();
+	int numP0inputs									= -1;
+	IssmDouble* P0inputs								= NULL;
+	Vector<IssmDouble>* vP0inputs					= NULL;
+	int* P0input_enums								= NULL;
+	int* P0input_interp 								= NULL;
+	int numP1inputs									= -1;
+	IssmDouble* P1inputs								= NULL;
+	Vector<IssmDouble>* vP1inputs					= NULL;
+	int* P1input_enums  								= NULL;
+	int* P1input_interp 								= NULL;
+	Vector<IssmDouble>* input_interpolations	= NULL;
+	IssmDouble* input_interpolations_serial	= NULL;
+   int* pos												= NULL; 
+	IssmDouble value									= 0;
 
 	/*Figure out how many inputs we have and their respective interpolation*/
-	Vector<IssmDouble>* input_interpolations=new Vector<IssmDouble>(maxinputs);
+	input_interpolations=new Vector<IssmDouble>(maxinputs);
 	if(this->elements->Size()){
 		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(0));
 		element->GetInputsInterpolations(input_interpolations);
 	}
+
+	/*Assemble and serialize*/
 	input_interpolations->Assemble();
-
-	/*Serialize and set output*/
-	IssmDouble* input_interpolations_serial = input_interpolations->ToMPISerial();
-	delete input_interpolations;
+	input_interpolations_serial = input_interpolations->ToMPISerial();
 
 	/*Count and get enums of all inputs in old mesh*/
-	int  numP0inputs;
-	int  numP1inputs;
-	int *P0input_enums  = NULL;
-	int *P1input_enums  = NULL;
-	int *P0input_interp = NULL;
-	int *P1input_interp = NULL;
 	for(int step=0;step<2;step++){
 		if(step){
 			P0input_enums  = xNew<int>(numP0inputs);
+			P0input_interp = xNew<int>(numP0inputs);
 			P1input_enums  = xNew<int>(numP1inputs);
-			P0input_interp = xNew<int>(numP0inputs);
 			P1input_interp = xNew<int>(numP1inputs);	
 		}
@@ -2720,146 +2744,168 @@
 		}
 	}
-
-	/*========== Deal with P0 inputs ==========*/
-	int numelementsold = this->elements->NumberOfElements();
-	int numelementsnew = newfemmodel_elements->NumberOfElements();
-	int numverticesold = this->vertices->NumberOfVertices();
-	int numverticesnew = newfemmodel_vertices->NumberOfVertices();
-	IssmDouble* P0inputsold = xNew<IssmDouble>(numelementsold*numP0inputs);
-	IssmDouble* P0inputsnew = NULL;
-	IssmDouble* vector      = NULL;
-
-	for(int i=0;i<numP0inputs;i++){
-		GetVectorFromInputsx(&vector,this,P0input_enums[i],ElementSIdEnum);
-
-		/*Copy vector in matrix*/
-		for(int j=0;j<numelementsold;j++) P0inputsold[j*numP0inputs+i] = vector[j];
-		xDelete<IssmDouble>(vector);
-	}
-
-	/*========== Deal with P1 inputs ==========*/
-	IssmDouble* P1inputsold = xNew<IssmDouble>(numverticesold*numP1inputs);
-	IssmDouble* P1inputsnew = NULL;
-	vector      = NULL;
-
-	for(int i=0;i<numP1inputs;i++){
-		GetVectorFromInputsx(&vector,this,P1input_enums[i],VertexSIdEnum);
-
-		/*Copy vector in matrix*/
-		for(int j=0;j<numverticesold;j++) P1inputsold[j*numP1inputs+i] = vector[j];
-		xDelete<IssmDouble>(vector);
-	}
-
-	/*Old mesh coordinates*/
-	IssmDouble *Xold     = NULL;
-	IssmDouble *Yold     = NULL;
-	IssmDouble *Zold		= NULL;
-	int        *Indexold = NULL;
-	IssmDouble *Xnew     = NULL;
-	IssmDouble *Ynew     = NULL;
-	IssmDouble *Znew		= NULL;
-	IssmDouble* XC_new   = NULL;
-	IssmDouble* YC_new   = NULL;
-	int        *Indexnew = NULL;
-
-	/*Get the old mesh*/
-	this->GetMesh(this->vertices,this->elements,&Xold,&Yold,&Zold,&Indexold);
-
-	/*Get the new mesh*/
-	this->GetMesh(newfemmodel_vertices,newfemmodel_elements,&Xnew,&Ynew,&Znew,&Indexnew);
-
-	/*Calculate the center points xc and xy*/
-	int elementswidth = this->GetElementsWidth(); //just tria in this version
-	/*new mesh*/
-	XC_new=xNewZeroInit<IssmDouble>(numelementsnew);
-	YC_new=xNewZeroInit<IssmDouble>(numelementsnew);
-	for(int i=0;i<numelementsnew;i++){
+	
+	/*Get P0 and P1 inputs over the elements*/	
+	pos		= xNew<int>(elementswidth);
+	vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
+	vP1inputs= new Vector<IssmDouble>(numberofvertices*numP1inputs);
+	for(int i=0;i<this->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+		
+		/*Get P0 inputs*/
+		for(int j=0;j<numP0inputs;j++){
+			TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));		
+			input->GetInputAverage(&value);
+			pos[0]=element->Sid()*numP0inputs+j;
+			/*Insert input in the vector*/	
+			vP0inputs->SetValues(1,pos,&value,INS_VAL);
+		}
+		
+		/*Get P1 inputs*/
+		for(int j=0;j<numP1inputs;j++){
+			TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P1input_enums[j]));
+			pos[0]=element->vertices[0]->Sid()*numP1inputs+j;
+			pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
+			pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
+			/*Insert input in the vector*/	
+			vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);	
+		}
+	}
+
+	/*Assemble and serialize*/
+	vP0inputs->Assemble();
+	vP1inputs->Assemble();
+	P0inputs=vP0inputs->ToMPISerial();
+	P1inputs=vP1inputs->ToMPISerial();
+	
+	/*Assign pointers*/
+	*pnumP0inputs		= numP0inputs; 
+	*pP0inputs			= P0inputs; 
+	*pP0input_enums	= P0input_enums;
+	*pP0input_interp	= P0input_interp;
+	*pnumP1inputs		= numP1inputs; 
+	*pP1inputs			= P1inputs; 
+	*pP1input_enums	= P1input_enums;
+	*pP1input_interp	= P1input_interp;	
+
+	/*Cleanup*/
+	delete input_interpolations;
+	delete vP0inputs;
+	delete vP1inputs;
+	xDelete<IssmDouble>(input_interpolations_serial);
+	xDelete<int>(pos);
+}
+/*}}}*/
+void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
+	
+	int numberofelements			= this->elements->NumberOfElements();	//global, entire old mesh
+	int newnumberofelements		= newfemmodel_elements->Size();			//just on the new partition
+	int numberofvertices			= this->vertices->NumberOfVertices();	//global, entire old mesh
+	int newnumberofvertices 	= newfemmodel_vertices->Size();			//just on the new partition
+	int elementswidth				= this->GetElementsWidth(); //just tria in this version
+	int numP0inputs				= -1;
+	IssmDouble* P0inputs			= NULL; //global, entire old mesh
+	IssmDouble* newP0inputs		= NULL; //just on the new partition
+	int* P0input_enums			= NULL;
+	int* P0input_interp 			= NULL;
+	int numP1inputs				= -1;
+	IssmDouble* P1inputs			= NULL; //global, entire old mesh
+	IssmDouble* newP1inputs 	= NULL; //just on the new partition
+	int* P1input_enums  			= NULL;
+	int* P1input_interp 			= NULL;
+	IssmDouble* values			= NULL;	
+   IssmDouble* vector      	= NULL;
+	IssmDouble* x					= NULL;//global, entire old mesh
+	IssmDouble* y					= NULL;//global, entire old mesh
+	IssmDouble* z					= NULL;//global, entire old mesh
+	int* elementslist				= NULL;//global, entire old mesh
+	IssmDouble* newx				= NULL;//just on the new partition
+	IssmDouble* newy				= NULL;//just on the new partition
+	IssmDouble* newz				= NULL;//just on the new partition
+	IssmDouble* newxc				= NULL;//just on the new partition
+	IssmDouble* newyc				= NULL;//just on the new partition
+	int* newelementslist			= NULL;//just on the new partition
+	int* sidtoindex				= NULL;//global vertices sid to partition index 
+
+	/*Get old P0 and P1  inputs (entire mesh)*/
+	this->GetInputs(&numP0inputs,&P0inputs,&P0input_enums,&P0input_interp,&numP1inputs,&P1inputs,&P1input_enums,&P1input_interp);
+
+	/*Get the old mesh (global, entire mesh)*/
+	this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
+
+	/*Get the new mesh (just on the new partition)*/
+	this->GetMeshOnPartition(newfemmodel_vertices,newfemmodel_elements,&newx,&newy,&newz,&newelementslist,&sidtoindex);
+
+	/*Calculate the center points xc and xy (new mesh, new partition)*/
+	newxc=xNewZeroInit<IssmDouble>(newnumberofelements);
+	newyc=xNewZeroInit<IssmDouble>(newnumberofelements);
+	for(int i=0;i<newnumberofelements;i++){
 		for(int j=0;j<elementswidth;j++){
-			int vid = Indexnew[i*elementswidth+j]-1;//Transform to C indexing
-			XC_new[i]+=Xnew[vid];
-			YC_new[i]+=Ynew[vid];
-		}
-		XC_new[i]=XC_new[i]/3.;
-		YC_new[i]=YC_new[i]/3.;
-	}
-
-	/*Interplate P0 inputs in the new mesh*/
-	InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
-				P0inputsold,numelementsold,numP0inputs,
-				XC_new,YC_new,numelementsnew,NULL);
-
-	/*Interpolate P1 inputs in the new mesh*/
-	InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
-				P1inputsold,numverticesold,numP1inputs,
-				Xnew,Ynew,numverticesnew,NULL);
-
-	/*Insert P0 inputs into the new elements.*/
-	vector=NULL;
-	for(int i=0;i<numP0inputs;i++){
-
-		/*Get P0 input vector from the interpolated matrix*/
-		vector=xNew<IssmDouble>(numelementsnew);
-		for(int j=0;j<numelementsnew;j++) vector[j]=P0inputsnew[j*numP0inputs+i];//vector has values for all elements (serial)
-
-		/*Update elements from inputs: */
-		for(int j=0;j<newfemmodel_elements->Size();j++){
-			Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j));
-			switch(P0input_interp[i]){	
+			int vid = newelementslist[i*elementswidth+j]-1;//Transform to C indexing
+			newxc[i]+=newx[vid]/elementswidth;
+			newyc[i]+=newy[vid]/elementswidth;
+		}
+	}
+
+	/*Interplate P0 inputs in the new mesh (just on the new partition)*/
+	InterpFromMeshToMesh2dx(&newP0inputs,elementslist,x,y,numberofvertices,numberofelements,
+				P0inputs,numberofelements,numP0inputs,
+				newxc,newyc,newnumberofelements,NULL);
+
+	/*Interpolate P1 inputs in the new mesh (just on the new partition)*/
+	InterpFromMeshToMesh2dx(&newP1inputs,elementslist,x,y,numberofvertices,numberofelements,
+				P1inputs,numberofvertices,numP1inputs,
+				newx,newy,newnumberofvertices,NULL);
+
+	/*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
+	values=xNew<IssmDouble>(elementswidth);	
+	for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
+		Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
+		/*newP0inputs is just on the new partition*/
+		for(int j=0;j<numP0inputs;j++){
+			switch(P0input_interp[j]){	
 				case P0Enum:
 				case DoubleInputEnum:
-					element->AddInput(new DoubleInput(P0input_enums[i],vector[element->sid]));//sid because newfemmodel has just a partitioning 
+					element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j]));
 					break;
 				case IntInputEnum: 
-					element->AddInput(new IntInput(P0input_enums[i],reCast<int>(vector[element->sid])));//sid because newfemmodel has just a partitioning
+					element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j])));
 					break;
 				case BoolInputEnum:
-					element->AddInput(new BoolInput(P0input_enums[i],reCast<bool>(vector[element->sid])));//sid because newfemmodel has just a partitioning
+					element->AddInput(new BoolInput(P0input_enums[j],reCast<bool>(newP0inputs[i*numP0inputs+j])));
 					break;
 				default:
-					_error_(EnumToStringx(P0input_enums[i])<<" Not supported yet");
+					_error_(EnumToStringx(P0input_enums[j])<<" Not supported yet");
 			}
 		}
-
-		xDelete<IssmDouble>(vector);
-	}
-
-	/*Insert P1 inputs into the new elements.*/
-	vector=NULL;
-	for(int i=0;i<numP1inputs;i++){
-
-		/*Get P1 input vector from the interpolated matrix*/
-		vector=xNew<IssmDouble>(numverticesnew);
-		for(int j=0;j<numverticesnew;j++) vector[j]=P1inputsnew[j*numP1inputs+i];//vector has all vertices	(serial)
-
-		/*Update elements from inputs: */
-		//InputUpdateFromVectorx(newfemmodel,vector,P1input_enums[i],VertexSIdEnum);//VertexSId because vector is serial in SId indexing
-		for(int j=0;j<newfemmodel_elements->Size();j++){
-			Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(j));
-			element->InputUpdateFromVector(vector,P1input_enums[i],VertexSIdEnum);
-		}
-
-		xDelete<IssmDouble>(vector);
-	}
-
+		/*newP1inputs is just on the new partition*/
+		for(int j=0;j<numP1inputs;j++){
+			values[0]=newP1inputs[sidtoindex[element->vertices[0]->Sid()]*numP1inputs+j];
+			values[1]=newP1inputs[sidtoindex[element->vertices[1]->Sid()]*numP1inputs+j];
+			values[2]=newP1inputs[sidtoindex[element->vertices[2]->Sid()]*numP1inputs+j];
+			element->inputs->AddInput(new TriaInput(P1input_enums[j],values,P1Enum));
+		}
+	}
+	
 	/*Cleanup*/
-	xDelete<IssmDouble>(input_interpolations_serial);
-	xDelete<IssmDouble>(P0inputsold);
-	xDelete<IssmDouble>(P0inputsnew);
+	xDelete<IssmDouble>(P0inputs);
+	xDelete<IssmDouble>(newP0inputs);
 	xDelete<int>(P0input_enums);
 	xDelete<int>(P0input_interp);
-	xDelete<IssmDouble>(P1inputsold);
-	xDelete<IssmDouble>(P1inputsnew);
+	xDelete<IssmDouble>(P1inputs);
+	xDelete<IssmDouble>(newP1inputs);
 	xDelete<int>(P1input_enums);
 	xDelete<int>(P1input_interp);
-	xDelete<IssmDouble>(Xold);
-	xDelete<IssmDouble>(Yold);
-	xDelete<IssmDouble>(Zold);
-	xDelete<int>(Indexold);
-	xDelete<IssmDouble>(Xnew);
-	xDelete<IssmDouble>(Ynew);
-	xDelete<IssmDouble>(Znew);
-	xDelete<IssmDouble>(XC_new);
-	xDelete<IssmDouble>(YC_new);
-	xDelete<int>(Indexnew);
+	xDelete<IssmDouble>(x);
+	xDelete<IssmDouble>(y);
+	xDelete<IssmDouble>(z);
+	xDelete<int>(elementslist);
+	xDelete<IssmDouble>(newx);
+	xDelete<IssmDouble>(newy);
+	xDelete<IssmDouble>(newz);
+	xDelete<IssmDouble>(newxc);
+	xDelete<IssmDouble>(newyc);
+	xDelete<int>(newelementslist);
+	xDelete<int>(sidtoindex);
+	xDelete<IssmDouble>(values);
 }
 /*}}}*/
@@ -3109,5 +3155,4 @@
 		}
 	}
-	return;
 }
 /*}}}*/
@@ -3117,19 +3162,21 @@
 	if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
 	
-	int numberofvertices, numberofelements;
-	int elementswidth = this->GetElementsWidth(); // just 2D mesh in this version (just tria elements)
-	IssmDouble *x		= NULL;
-	IssmDouble *y		= NULL;
-	IssmDouble *z		= NULL;	
-	int* elementslist = NULL;
-	
+	int numberofvertices = femmodel_vertices->NumberOfVertices();
+	int numberofelements = femmodel_elements->NumberOfElements();
+	int elementswidth		= this->GetElementsWidth(); // just 2D mesh in this version (just tria elements)
+	IssmDouble* x			= NULL;
+	IssmDouble* y			= NULL;
+	IssmDouble* z			= NULL;	
+	int* elementslist 	= NULL;
+	int* elem_vertices	= NULL;
+	IssmDouble *id1		= NULL;
+   IssmDouble *id2 		= NULL;
+	IssmDouble *id3 		= NULL;
+
 	/*Get vertices coordinates*/
-	VertexCoordinatesx(&x, &y, &z, femmodel_vertices,false) ;
-
-	numberofvertices = femmodel_vertices->NumberOfVertices();
-	numberofelements = femmodel_elements->NumberOfElements();
+	VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
 	
 	/*Get element vertices*/
-	int* elem_vertices=xNew<int>(elementswidth);
+	elem_vertices				= xNew<int>(elementswidth);
 	Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
 	Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
@@ -3151,7 +3198,7 @@
 
    /*Serialize*/
-	IssmDouble *id1 = vid1->ToMPISerial();
-   IssmDouble *id2 = vid2->ToMPISerial();
-	IssmDouble *id3 = vid3->ToMPISerial();
+	id1 = vid1->ToMPISerial();
+   id2 = vid2->ToMPISerial();
+	id3 = vid3->ToMPISerial();
 	
 	/*Construct elements list*/
@@ -3180,99 +3227,150 @@
 }
 /*}}}*/
-void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
-
-	/*itapopo ATTENTION: JUST SPCVX AND SPCVY TO TEST!!!*/
-	/*OTHERS CONSTRAINTS MUST BE IMPLEMENTED!!!*/
-	
-	/*Get x and y of the mesh i-1*/
-	int numberofvertices			= this->vertices->NumberOfVertices();
-	int numberofelements			= this->elements->NumberOfElements();
-	IssmDouble *x					= NULL;
-	IssmDouble *y					= NULL;
-	IssmDouble *z					= NULL;
-	int        *elementslist	= NULL;
-
-	/*elementslist is in Matlab indexing*/
+void FemModel::GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px,IssmDouble** py,IssmDouble** pz,int** pelementslist,int** psidtoindex){/*{{{*/
+
+	if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
+	if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
+	
+	int numberofvertices			= femmodel_vertices->Size();	//number of vertices of this partition
+	int numbertotalofvertices	= femmodel_vertices->NumberOfVertices();	//number total of vertices (entire mesh)
+	int numberofelements			= femmodel_elements->Size();  //number of elements of this partition
+	int elementswidth				= this->GetElementsWidth();	//just 2D mesh in this version (just tria elements)
+	IssmDouble* x					= NULL;
+	IssmDouble* y					= NULL;
+	IssmDouble* z					= NULL;	
+	int* elementslist				= NULL;
+	int* sidtoindex				= NULL;
+	int* elem_vertices			= NULL;
+	
+	/*Get vertices coordinates of this partition*/
+	sidtoindex	= xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
+	x				= xNew<IssmDouble>(numberofvertices);//just this partition
+	y				= xNew<IssmDouble>(numberofvertices);//just this partitio;
+	z				= xNew<IssmDouble>(numberofvertices);//just this partitio;
+	
+	/*Go through in this partition (vertices)*/
+	for(int i=0;i<numberofvertices;i++){//just this partition
+		Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);	
+		/*Attention: no spherical coordinates*/
+		x[i]=vertex->GetX();
+		y[i]=vertex->GetY();
+		z[i]=vertex->GetZ();
+		/*Keep the index and sid pair*/
+		sidtoindex[vertex->Sid()]=i;
+	}
+
+	/*Go through in this partition (elements) and build the element list*/
+	elem_vertices= xNew<int>(elementswidth);
+	elementslist = xNew<int>(numberofelements*elementswidth);
+	if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
+	
+	for(int i=0;i<numberofelements;i++){//just this partition
+    	Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
+    	element->GetVerticesSidList(elem_vertices);
+		elementslist[elementswidth*i+0] = sidtoindex[elem_vertices[0]]+1; //InterpMesh wants Matlab indexing
+		elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
+		elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
+	}	
+		
+	/*Assign pointers*/
+	*px				= x;
+	*py				= y;
+	*pz				= z;
+	*pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
+	*psidtoindex	= sidtoindex;  //it is ncessary to insert inputs 
+
+	/*Cleanup*/
+	xDelete<int>(elem_vertices);
+}
+/*}}}*/
+void FemModel::CreateConstraints(Vertices* newfemmodel_vertices,int nodecounter,int constraintcounter,int analysis_enum,Constraints* newfemmodel_constraints){/*{{{*/
+
+	/*ATTENTION: JUST SPCVX AND SPCVY*/
+	/*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
+	if(analysis_enum!=StressbalanceAnalysisEnum) return;
+	
+	int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum);
+	int dofpernode						= 2;														//vx and vy 
+	int numberofcols					= dofpernode*2;										//to keep dofs and flags in the vspc vector
+	int numberofvertices				= this->vertices->NumberOfVertices();			//global, entire old mesh
+	int numberofelements				= this->elements->NumberOfElements();			//global, entire old mesh
+	int newnumberofvertices			= newfemmodel_vertices->Size();					//local, just the new partition
+	int count							= 0;
+	IssmDouble* x						= NULL;													//global, entire old mesh
+	IssmDouble* y						= NULL;													//global, entire old mesh
+	IssmDouble* z						= NULL;													//global, entire old mesh
+	int*			elementslist		= NULL;													//global, entire old mesh
+	IssmDouble* spc					= NULL;													//global, entire old mesh
+	IssmDouble* newx					= NULL;													//local, just new partition
+	IssmDouble* newy					= NULL;													//local, just new partition
+	IssmDouble* newspc				= NULL;													//local, just new partition
+	IssmDouble eps						= 1.e-8;
+	Vector<IssmDouble>* vspc		= new Vector<IssmDouble>(numberofnodes_analysistype*numberofcols);
+
+	/*Get old mesh (global, entire mesh). Elementslist comes in Matlab indexing*/
 	this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elementslist);
-	
-	/*Get spcvx and spcvy for mesh i-1*/
-	IssmDouble *spcvx						= NULL;
-	IssmDouble *spcvy						= NULL;
-	IssmDouble *spcvxflag				= NULL;
-	IssmDouble *spcvyflag				= NULL;
-	int numberofnodes_analysistype	= this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
-	Vector<IssmDouble>* vspcvx			= new Vector<IssmDouble>(numberofnodes_analysistype);
-	Vector<IssmDouble>* vspcvy			= new Vector<IssmDouble>(numberofnodes_analysistype);
-	Vector<IssmDouble>* vspcvxflag	= new Vector<IssmDouble>(numberofnodes_analysistype);
-	Vector<IssmDouble>* vspcvyflag	= new Vector<IssmDouble>(numberofnodes_analysistype);
-	
+
+	/*Get vertices coordinates of the new partition*/
+	newx=xNew<IssmDouble>(newnumberofvertices);//just the new partition
+	newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
+	for(int i=0;i<newnumberofvertices;i++){//just the new partition
+		Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);	
+		/*Attention: no spherical coordinates*/
+		newx[i]=vertex->GetX();
+		newy[i]=vertex->GetY();
+	}
+
+	/*Get spcvx and spcvy of old mesh*/
 	for(int i=0;i<this->constraints->Size();i++){
-		SpcStatic* spc			= xDynamicCast<SpcStatic*>(this->constraints->GetObjectByOffset(i));
-		int dof					= spc->GetDof();
-		int node					= spc->GetNodeId();
-		IssmDouble spcvalue	= spc->GetValue(); 
+		
+		Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
+		if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n");
+
+		SpcStatic* spcstatic = xDynamicCast<SpcStatic*>(constraint);
+		int dof					= spcstatic->GetDof();
+		int node					= spcstatic->GetNodeId();
+		IssmDouble spcvalue	= spcstatic->GetValue(); 
 		int nodeindex			= node-1;
 		
+		/*vx and vx flag insertion*/
 		if(dof==0) {//vx
-			vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
-			vspcvxflag->SetValue(nodeindex,1,INS_VAL);
+			vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL);    //vx
+			vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
 		} 
+		/*vy and vy flag insertion*/
 		if(dof==1){//vy
-			vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
-			vspcvyflag->SetValue(nodeindex,1,INS_VAL);
-		}
-	}
-
-	/*Assemble*/
-	vspcvx->Assemble();
-	vspcvy->Assemble();
-	vspcvxflag->Assemble();
-	vspcvyflag->Assemble();
-
-	/*Serialize*/
-	spcvx		 = vspcvx->ToMPISerial();
-	spcvy		 = vspcvy->ToMPISerial();
-	spcvxflag = vspcvxflag->ToMPISerial();
-	spcvyflag = vspcvyflag->ToMPISerial();
-	
-	IssmDouble *newspcvx		 = NULL;
-	IssmDouble *newspcvy		 = NULL;
-	IssmDouble *newspcvxflag = NULL;
-	IssmDouble *newspcvyflag = NULL;
-	int nods_data				 = numberofvertices;
-	int nels_data				 = numberofelements;
-	int M_data					 = numberofvertices;
-	int N_data					 = 1;
-	int N_interp				 = newnumberofvertices;
-
-  /*Interpolate spcvx and spcvy in the new mesh*/
-	InterpFromMeshToMesh2dx(&newspcvx,elementslist,x,y,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,NULL);
-	InterpFromMeshToMesh2dx(&newspcvy,elementslist,x,y,nods_data,nels_data,spcvy,M_data,N_data,newx,newy,N_interp,NULL);
-	InterpFromMeshToMesh2dx(&newspcvxflag,elementslist,x,y,nods_data,nels_data,spcvxflag,M_data,N_data,newx,newy,N_interp,NULL);
-	InterpFromMeshToMesh2dx(&newspcvyflag,elementslist,x,y,nods_data,nels_data,spcvyflag,M_data,N_data,newx,newy,N_interp,NULL);
-
-	int count					= 0;
-	IssmDouble eps				= 1.e-8;
-
+			vspc->SetValue(nodeindex*numberofcols+1,spcvalue,INS_VAL);	//vy
+			vspc->SetValue(nodeindex*numberofcols+dofpernode+1,1,INS_VAL);//vyflag
+		}
+	}
+
+	/*Assemble and serialize*/
+	vspc->Assemble();
+	spc=vspc->ToMPISerial();
+
+	/*Interpolate spc values and flags in the new partition*/
+	InterpFromMeshToMesh2dx(&newspc,elementslist,x,y,numberofvertices,numberofelements,
+								spc,numberofvertices,numberofcols,
+								newx,newy,newnumberofvertices,NULL);
+	
 	/*Now, insert the interpolated constraints in the data set (constraints)*/
-	for(int i=0;i<newnumberofvertices;i++){
-		if(my_vertices[i]){ 
-			/*spcvx*/
-			if(!xIsNan<IssmDouble>(newspcvx[i]) && newspcvxflag[i]>(1-eps)){
-				constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,newspcvx[i],StressbalanceAnalysisEnum));
-				//add count'th spc, on node i+1, setting dof 1 to vx.
-				count++;
-			}
-		}
-	}
 	count=0;
-	for(int i=0;i<newnumberofvertices;i++){
-		if(my_vertices[i]){
-			/*spcvy*/
-			if(!xIsNan<IssmDouble>(newspcvy[i]) && newspcvyflag[i]>(1-eps) ){
-				constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,newspcvy[i],StressbalanceAnalysisEnum)); 
-				//add count'th spc, on node i+1, setting dof 1 to vx.
-				count++;
-			}
+	for(int i=0;i<newnumberofvertices;i++){//just in the new partition
+		Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
+		/*spcvx*/
+		if(!xIsNan<IssmDouble>(newspc[i*numberofcols]) && newspc[i*numberofcols+dofpernode]>(1-eps)){
+			newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,0,newspc[i*numberofcols],analysis_enum));
+			//add count'th spc, on node i+1, setting dof 1 to vx.
+			count++;
+		}
+	}
+	count=0;
+	for(int i=0;i<newnumberofvertices;i++){//just in the new partition
+		Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
+		/*spcvy*/
+		if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
+			newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 
+			//add count'th spc, on node i+1, setting dof 1 to vx.
+			count++;
 		}
 	}
@@ -3283,16 +3381,9 @@
 	xDelete<IssmDouble>(z);
 	xDelete<int>(elementslist);
-	xDelete<IssmDouble>(spcvx);
-	xDelete<IssmDouble>(spcvy);
-	xDelete<IssmDouble>(spcvxflag);
-	xDelete<IssmDouble>(spcvyflag);
-	xDelete<IssmDouble>(newspcvx);
-	xDelete<IssmDouble>(newspcvy);
-	xDelete<IssmDouble>(newspcvxflag);
-	xDelete<IssmDouble>(newspcvyflag);
-	delete vspcvx;
-	delete vspcvy;	
-	delete vspcvxflag;
-	delete vspcvyflag;
+	xDelete<IssmDouble>(spc);
+	xDelete<IssmDouble>(newspc);
+	xDelete<IssmDouble>(newx);
+	xDelete<IssmDouble>(newy);
+	delete vspc;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 22359)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 22360)
@@ -167,4 +167,5 @@
 		void AdjustBaseThicknessAndMask(void);
 		void GetMesh(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
+		void GetMeshOnPartition(Vertices* femmodel_vertices,Elements* femmodel_elements,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist,int** psidtoindex);
 		void GetGroundediceLevelSet(IssmDouble** pmasklevelset);
 		void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices);
@@ -172,5 +173,6 @@
 		void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
 		void CreateNodes(int newnumberofvertices,int* my_vertices,int nodecounter,int analysis_enum,Nodes* nodes);
-		void CreateConstraints(int newnumberofvertices,int newnumberofelements,int nodecounter,int constraintcounter,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints);
+		void CreateConstraints(Vertices* newfemmodel_vertices,int nodecounter,int constraintcounter,int analysis_enum,Constraints* newfemmodel_constraints);
+		void GetInputs(int* pnumP0inputs,IssmDouble** pP0inputs,int** pP0input_enums,int** pP0input_interp,int* pnumP1inputs,IssmDouble** pP1inputs,int** pP1input_enums,int** pP1input_interp);
 		void InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements);
 		void UpdateElements(int newnumberofelements,int* newelementslist,bool* my_elements,int nodecounter,int analysis_counter,Elements* newelements);