Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21523)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21524)
@@ -3103,9 +3103,4 @@
 	output->materials->ResetHooks();
 
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
-	/*Finally: interpolate all inputs*/
-	this->InterpolateInputs(output);
-
-	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
 	/*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
 	int analysis_type;
@@ -3122,4 +3117,8 @@
 	}
 
+	/*Finally: interpolate all inputs*/
+	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	this->InterpolateInputs(output);
+	
 	/*Reset current configuration: */
 	analysis_type=output->analysis_type_list[this->analysis_counter];
@@ -3214,23 +3213,61 @@
 	printarray(P0inputsold,numelementsold,numP0inputs);
 
-	_error_("stop");
+	/*========== Deal with P1 inputs ==========*/
+	IssmDouble* P1inputsold = xNew<IssmDouble>(numverticesold*numP1inputs);
+	IssmDouble* P1inputsnew = NULL;
+	vector      = NULL;
+
+	for(int i=0;i<numP1inputs;i++){
+		printf("Dealing with %s (type: %s)\n",EnumToStringx(P1input_enums[i]),EnumToStringx(P1input_interp[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);
+	}
+	printf("Printing array to make sure it looks alright\n");
+	printarray(P1inputsold,numverticesold,numP1inputs);
+	
 	/*Old mesh coordinates*/
 	IssmDouble *Xold     = NULL;
 	IssmDouble *Yold     = NULL;
+	IssmDouble *Zold		= NULL;
 	int        *Indexold = NULL;
+	int        *Indexnew = NULL;
 	IssmDouble *Xnew     = NULL;
 	IssmDouble *Ynew     = NULL;
-
-
+	IssmDouble *Znew		= NULL;
+	
+	/*Get the old mesh*/
 	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	this->GetMesh(this,&Xold,&Yold,&Zold,&Indexold);
+
+	/*Get the new mesh*/
+	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	this->GetMesh(newfemmodel,&Xnew,&Ynew,&Znew,&Indexnew);
+	
+	printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+	/*Interplate P0 inputs in the new mesh*/
 	InterpFromMeshToMesh2dx(&P0inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
 				P0inputsold,numelementsold,numP0inputs,
 				Xnew,Ynew,numelementsnew,NULL);
-
+	
+	printf("Printing array AFTER INTERP - P0 INPUTS - to make sure it looks alright\n");
+	printarray(P0inputsnew,numelementsnew,numP0inputs);
+	
+	/*Interpolate P1 inputs in the new mesh*/
+	InterpFromMeshToMesh2dx(&P1inputsnew,Indexold,Xold,Yold,numverticesold,numelementsold,
+				P1inputsold,numverticesold,numP1inputs,
+				Xnew,Ynew,numverticesnew,NULL);
+
+	printf("Printing array AFTER INTERP - P1 INPUTS - to make sure it looks alright\n");
+	printarray(P1inputsnew,numverticesnew,numP1inputs);
+	
 	_error_("STOP");
 
 	xDelete<IssmDouble>(P0inputsold);
 	xDelete<IssmDouble>(P0inputsnew);
-
+	xDelete<IssmDouble>(P1inputsold);
+	xDelete<IssmDouble>(P1inputsnew);
 
 	_error_("STOP");
@@ -3411,5 +3448,5 @@
 
 	int lid=0;
-
+//itapopo use the GetMesh!!
 	for(int j=0;j<newnumberofvertices;j++){
 		if(my_vertices[j]){				
@@ -3455,4 +3492,63 @@
 }
 /*}}}*/
+void FemModel::GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist){/*{{{*/
+
+	if(!femmodel) _error_("GetMesh: NULL femmodel.");
+	
+	int numberofvertices, numberofelements;
+	int elementswidth = 3; //itapopo just 2D mesh in this version (just tria elements)
+	IssmDouble *x		= NULL;
+	IssmDouble *y		= NULL;
+	IssmDouble *z		= NULL;	
+	int* elementslist = NULL;
+	
+	/*Get vertices coordinates*/
+	VertexCoordinatesx(&x, &y, &z, femmodel->vertices,false) ;
+
+	numberofvertices = femmodel->vertices->NumberOfVertices();
+	numberofelements = femmodel->elements->NumberOfElements();
+	
+	/*Get element vertices*/
+	int* elem_vertices=xNew<int>(elementswidth);
+	Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements);
+	Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements);
+	Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements);
+
+	/*Go through elements, and for each element, get vertices*/
+   for(int i=0;i<femmodel->elements->Size();i++){
+    	Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+    	element->GetVerticesSidList(elem_vertices);
+    	vid1->SetValue(element->sid,elem_vertices[0],INS_VAL);
+    	vid2->SetValue(element->sid,elem_vertices[1],INS_VAL);
+    	vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
+   }
+		
+	/*Assemble*/
+   vid1->Assemble();
+   vid2->Assemble();
+   vid3->Assemble();
+
+   /*Serialize*/
+	IssmDouble *id1 = vid1->ToMPISerial();
+   IssmDouble *id2 = vid2->ToMPISerial();
+	IssmDouble *id3 = vid3->ToMPISerial();
+	
+	/*Construct elements list*/
+	elementslist=xNew<int>(numberofelements*elementswidth);
+	if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
+	for(int i=0;i<numberofelements;i++){
+		std::cout << (int)id1[i]+1 << " " << (int)id2[i]+1 <<" " << (int)id3[i]+1 <<  std::endl;
+		elementslist[elementswidth*i+0] = (int)id1[i]+1; //InterpMesh wants Matlab indexing
+		elementslist[elementswidth*i+1] = (int)id2[i]+1; //InterpMesh wants Matlab indexing
+		elementslist[elementswidth*i+2] = (int)id3[i]+1; //InterpMesh wants Matlab indexinf
+	}
+	
+	/*Assign pointers*/
+	*px				= x;
+	*py				= y;
+	*pz				= z;
+	*pelementslist = elementslist;
+}
+/*}}}*/
 void FemModel::CreateConstraints(int newnumberofvertices,int newnumberofelements,int* newelementslist,IssmDouble* newx,IssmDouble* newy,int* my_vertices,Constraints* constraints){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21523)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21524)
@@ -149,4 +149,5 @@
 		void InitializeAdaptiveRefinement(void);
 		FemModel* ReMesh(void);
+		void GetMesh(FemModel* femmodel,IssmDouble** px, IssmDouble** py, IssmDouble** pz, int** pelementslist);
 		void ExecuteRefinement(int &newnumberofvertices,int &newnumberofelements,IssmDouble** newx,IssmDouble** newy,IssmDouble** newz,int** newelementslist);
 		void GetMaskLevelSet(IssmDouble** pmasklevelset);