Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21503)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 21504)
@@ -83,4 +83,10 @@
 	#ifdef _HAVE_NEOPZ_
 	this->InitializeAdaptiveRefinement();
+	FemModel *Test=this->ReMesh();//itapopo: just to test!;
+	printf("   AFTER REMESH!!!\n");
+	//Test->CleanUp();
+	printf("   AFTER CLEANUP!!!\n");
+	//delete Test;
+	printf("   AFTER DELETE!!!\n");
 	#endif
 
@@ -2936,5 +2942,5 @@
 	/*Get initial mesh*/
 	/*Get total number of elements and total numbers of elements in the initial mesh*/
-	long numberofvertices, numberofelements, numberofsegments;
+	int numberofvertices, numberofelements, numberofsegments;
 	numberofvertices = this->vertices->NumberOfVertices();
 	numberofelements = this->elements->NumberOfElements();
@@ -2973,16 +2979,16 @@
 	IssmDouble *id3 = vid3->ToMPISerial();
 	
-	long **elementsptr;
-	elementsptr = new long*[numberofelements];
+	int **elementsptr;
+	elementsptr = new int*[numberofelements];
    for(int i=0;i<numberofelements;i++){
-		elementsptr[i] = new long[elementswidth];
-      elementsptr[i][0] = (long)id1[i];
-		elementsptr[i][1] = (long)id2[i];
-      elementsptr[i][2] = (long)id3[i];
+		elementsptr[i] = new int[elementswidth];
+      elementsptr[i][0] = (int)id1[i];
+		elementsptr[i][1] = (int)id2[i];
+      elementsptr[i][2] = (int)id3[i];
 	}
 
 	/*Create initial mesh (coarse mesh) in neopz data structure*/ 
    if(my_rank==0){ 
-		long nsegments=0;
+		int nsegments=0;
 		this->amr->CreateInitialMesh(numberofvertices, numberofelements, nsegments, elementswidth, x, y, z, elementsptr, NULL);
 	}
@@ -3007,74 +3013,303 @@
 	
 	/*All indexing here is in C type: 0..n-1*/
-
-	//////// to test! using the current femmodel
-	int numberofvertices2, numberofelements2;
-	numberofvertices2 = this->vertices->NumberOfVertices();
-	numberofelements2 = this->elements->NumberOfElements();
-
-	/*Get vertices coordinates*/
-	IssmDouble *x = NULL;
-	IssmDouble *y = NULL;
-	IssmDouble *z = NULL;
-	VertexCoordinatesx(&x, &y, &z, this->vertices,false) ;
-
-	/*Get element vertices*/
-	int elementswidth2 = 3; //just 2D mesh in this version (just tria elements)
-	int* elem_vertices2=xNew<int>(elementswidth2);
-	Vector<IssmDouble>* vid1= new Vector<IssmDouble>(numberofelements2);
-	Vector<IssmDouble>* vid2= new Vector<IssmDouble>(numberofelements2);
-	Vector<IssmDouble>* vid3= new Vector<IssmDouble>(numberofelements2);
-
-	/*Go through elements, and for each element, object, report it cpu:*/
-    for(int i=0;i<this->elements->Size();i++){
-    	Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-    	element->GetVerticesSidList(elem_vertices2);
-    	vid1->SetValue(element->sid,elem_vertices2[0],INS_VAL);
-    	vid2->SetValue(element->sid,elem_vertices2[1],INS_VAL);
-    	vid3->SetValue(element->sid,elem_vertices2[2],INS_VAL);
-    }
-		
+	int my_rank=IssmComm::GetRank();
+   const int elementswidth=3;//just 2D mesh, tria elements
+	
+	/*Solutions which will be used to refine the elements*/
+	IssmDouble *vx=NULL; //This will be used in constraints
+	IssmDouble *vy=NULL; //This will be used in constraints
+	IssmDouble *masklevelset=NULL;
+ 
+	this->GetMaskLevelSet(&masklevelset);
+	// #1: TEST MASK LEVEL SET
+	_printf_("   PRINTING MASKLEVELSET... \n");
+	int numberofvertices=this->vertices->NumberOfVertices();	    
+	if(my_rank==0) for(int i=0;i<numberofvertices;i++) _printf_("vertex: " << i << "\t" << masklevelset[i] << "\n");
+
+	/*Refine the mesh and get the new mesh*/
+	IssmDouble *newx;
+	IssmDouble *newy;
+	IssmDouble *newz;
+	int **newelements;
+	int **newsegments;
+	int newnumberofvertices, newnumberofelements, newnumberofsegments;
+	int type_process=1; //1: it refines father mesh. See AdaptiveMeshRefinement.h (.cpp)
+	if(my_rank==0) this->amr->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,&newx,&newy,&newz,&newelements,&newsegments);
+	
+	if(newnumberofvertices<=0 || newnumberofelements<=0 /*|| newnumberofsegments<=0*/) _error_("Error in the refinement process.");
+	
+	/*Cleanup masklevetset*/
+	xDelete<IssmDouble>(masklevelset);
+
+	// #2: TEST NEW MESH
+	if(my_rank==0){
+		_printf_("   PRINTING COORDINATES... \n");	    
+ 		for(int i=0;i<newnumberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << newx[i] << "\t" << "Y: " << newy[i] << "\n");
+    	_printf_("   PRINTING ELEMENTS... \n");	   
+   	for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelements[i][0] << "\t" << newelements[i][1] << "\t" << newelements[i][2] << "\n");
+	}
+	
+	/*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
+	/*Fill the element list to partitioning*/	
+	int* newelementslist=NULL;
+	newelementslist=xNew<int>(newnumberofelements*elementswidth);
+	for(int i=0;i<newnumberofelements;i++){ //itapopo os elementos poderão sair do ExecuteRefinement nesse formato
+		for(int j=0;j<elementswidth;j++){
+			newelementslist[elementswidth*i+j]=newelements[i][j]; //C indexing
+		}
+	}	
+	bool* my_elements=NULL; 
+	int* my_vertices=NULL;
+
+	this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
+	
+	/*Creating new femmodel with new mesh*/
+	FemModel* output=NULL;
+	output=new FemModel(*this);
+
+	/*Copy basic attributes:*/
+	output->nummodels=this->nummodels;
+	output->solution_type=this->solution_type;
+	output->analysis_counter=this->analysis_counter;
+
+	/*Now, deep copy arrays:*/
+	output->analysis_type_list=xNew<int>(this->nummodels);
+	xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
+
+	output->profiler=static_cast<Profiler*>(this->profiler->copy());
+	output->parameters=static_cast<Parameters*>(this->parameters->Copy());	    
+
+	/*Create vertices*/
+	output->vertices=new Vertices();
+	this->CreateVertices(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,my_vertices,newx,newy,newz,output->vertices);
+ 
+	/* #4: TEST THE VERTICES!*/
+  	//_printf_("     Old vertices deep echo: \n");
+   //this->vertices->DeepEcho();
+   //_printf_("     New vertices deep echo: \n");
+   //output->vertices->DeepEcho();
+
+	/*Creating elements*/
+	/*Just Tria in this version*/
+	output->elements=new Elements();
+	this->CreateElements(newnumberofelements,elementswidth,newelementslist,my_elements,output->elements);
+
+	/*Cleanup*/
+	xDelete<int>(newelementslist);	    
+
+	/*Creating materials*/
+	output->materials=new Materials();
+	this->CreateMaterials(newnumberofelements,my_elements,output->materials);
+	/*Creating nodes*/
+	/*Just SSA (2D) and P1 in this version*/
+	output->nodes=new Nodes();
+
+	int nodecounter=0;
+	int lid=0;
+	for(int i=0;i<output->nummodels;i++){
+		int analysis_enum = output->analysis_type_list[i];
+		//itapopo as the domain is 2D, it is not necessary to create nodes for this analysis
+		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;	    
+		for(int j=0;j<newnumberofvertices;j++){
+			if(my_vertices[j]){				
+				Node* newnode=new Node();	
+				/*id: */
+				newnode->id=nodecounter+j+1;
+				newnode->sid=j;
+				newnode->lid=lid++;
+				newnode->analysis_enum=analysis_enum;
+				/*Initialize coord_system: Identity matrix by default*/
+				for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
+				for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
+				/*indexing:*/
+				newnode->indexingupdate=true;
+				Analysis* analysis=EnumToAnalysis(analysis_enum);
+				int *doftypes=NULL;
+				int numdofs=analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum);
+				newnode->indexing.Init(numdofs,doftypes);
+				xDelete<int>(doftypes);
+				delete analysis;
+				if(analysis_enum==StressbalanceAnalysisEnum)
+					newnode->SetApproximation(SSAApproximationEnum);
+				else
+					newnode->SetApproximation(0);
+
+				/*Stressbalance Horiz*/
+				if(analysis_enum==StressbalanceAnalysisEnum){
+					// itapopo
+					/*Coordinate system provided, convert to coord_system matrix*/
+					//XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
+					//_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);
+
+				}
+				output->nodes->AddObject(newnode);
+			}
+		}
+
+		if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
+	}
+	
+	/* #4: TEST THE NODES!*/
+	_printf_("		Old nodes deep echo: \n");
+	this->nodes->DeepEcho();
+	_printf_("		New nodes deep echo: \n");
+	output->nodes->DeepEcho();
+
+	printf(" I arrived here!!!!!!!\n");
+#ifdef _CONTINUE_NEOPZ_
+
+	/*Create constraints*/
+	IssmDouble *spcvx=NULL;
+	IssmDouble *spcvy=NULL;
+	int numberofnodes_analysistype=this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
+	Vector<IssmDouble>* vspcvx=new Vector<IssmDouble>(numberofnodes_analysistype);
+	Vector<IssmDouble>* vspcvy=new Vector<IssmDouble>(numberofnodes_analysistype);
+
+	IssmDouble BigNumber=1.e8;
+
+	for(int i=0;i<numberofnodes_analysistype;i++){
+		vspcvx->SetValue(i,BigNumber,INS_VAL);
+		vspcvy->SetValue(i,BigNumber,INS_VAL);
+	}
+
+	IssmDouble absmaxspcvx=0;
+	IssmDouble absmaxspcvy=0;
+
+	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(); 
+		int nodeindex=node-1;
+		if(dof==0) {
+			vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
+			if(fabs(spcvalue)>absmaxspcvx) absmaxspcvx=fabs(spcvalue);
+		}
+		else {
+			vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
+			if(fabs(spcvalue)>absmaxspcvy) absmaxspcvy=fabs(spcvalue);
+		}
+	}
+
 	/*Assemble*/
-    vid1->Assemble();
-    vid2->Assemble();
-    vid3->Assemble();
-
-    /*Serialize*/
-	IssmDouble *id1 = vid1->ToMPISerial();
-    IssmDouble *id2 = vid2->ToMPISerial();
-	IssmDouble *id3 = vid3->ToMPISerial();
-
-	delete vid1;
-	delete vid2;
-	delete vid3;
-	
-	int **elementsptr;
-	elementsptr = new int*[numberofelements2];
-    for(int i=0;i<numberofelements2;i++){
-        elementsptr[i] = new int[elementswidth2];
-        elementsptr[i][0] = (int)id1[i];
-        elementsptr[i][1] = (int)id2[i];
-        elementsptr[i][2] = (int)id3[i];
-    }
-
-   delete id1;
-	delete id2;
-	delete id3;
-	//////////////////////////////////////////////
-
-	int my_rank=IssmComm::GetRank();
-
-	/*Solutions which will be used to refine the elements*/
-	double *vx = NULL;
-	double *vy = NULL;
-	double *masklevelset = NULL;
-
-	int elementswidth = 3;//just 2D mesh, tria elements
-	int numberofelements = this->elements->NumberOfElements();
-	int numberofvertices = this->vertices->NumberOfVertices();
-
-	IssmDouble* elementlevelset = xNew<IssmDouble>(elementswidth);
-	int* elem_vertices = xNew<int>(elementswidth);
-	Vector<IssmDouble>* vmasklevelset = new Vector<IssmDouble>(numberofvertices);
+	vspcvx->Assemble();
+	vspcvy->Assemble();
+
+	/*Serialize*/
+	spcvx=vspcvx->ToMPISerial();
+	spcvy=vspcvy->ToMPISerial();
+
+	/*Free the data*/
+	delete vspcvx;
+	delete vspcvy;
+	
+	IssmDouble *newspcvx=NULL;
+	IssmDouble *newspcvy=NULL;
+	int *oldelements=newelementslist; //itapopo
+	int nods_data=numberofnodes_analysistype;
+	int nels_data=newnumberofelements;
+	int M_data=numberofnodes_analysistype;
+	int N_data=1;
+	int N_interp=newnumberofvertices;//itapopo
+	Options *options=NULL;
+
+	InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+	InterpFromMeshToMesh2dx(&newspcvx,oldelements,vx,vy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
+// TEST: New SPCVX and VY must be tested!
+
+	output->constraints = new Constraints();
+
+	nodecounter 			= 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
+	int count 				= 0;
+	int constraintcounter 	= 0; //itapopo
+	IssmDouble eps			= 1.e-2;
+
+	for(int i=0;i<newnumberofvertices;i++){
+		if(my_vertices[i])
+		/*spcvx*/
+		if(fabs(spcvx[i]) < absmaxspcvx+eps){//itapopo
+			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,spcvx[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(fabs(spcvy[i]) < absmaxspcvy+eps){//itapopo
+			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
+			count++;
+		}
+
+	}
+
+	/* #5: TEST CONSTRAINTS*/	    
+	this->constraints->DeepEcho();
+	output->constraints->DeepEcho();
+
+	// output->loads=static_cast<Loads*>(this->loads->Copy());
+	// output->constraints=static_cast<Constraints*>(this->constraints->Copy());
+	// output->results=static_cast<Results*>(this->results->Copy());
+
+	// /*reset hooks for elements, loads and nodes: */
+	// output->elements->ResetHooks();
+	// output->loads->ResetHooks();
+	// output->materials->ResetHooks();
+
+	// /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
+	// for(i=0;i<nummodels;i++){
+	// 	analysis_type=output->analysis_type_list[i];
+	// 	output->SetCurrentConfiguration(analysis_type);
+	// 	if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
+	// 	SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
+	// 	NodesDofx(output->nodes,output->parameters,analysis_type);
+	// 	ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
+	// }
+
+	// /*Reset current configuration: */
+	// analysis_type=output->analysis_type_list[analysis_counter];
+	// output->SetCurrentConfiguration(analysis_type);
+
+	/** TODO
+
+	- generate the required input objects for NeoPZ
+		AMR has fathermesh and previousmesh. On first call, these meshes are generated. 
+
+	- call NeoPZ
+		This creates a newmesh which will be refined. 
+
+	- get the new mesh (index,x,y) from NeoPZ
+		Is is doing by GetNewMesh method.
+
+	- Create a new FemModel* that will be consistent with the new mesh
+		It can be done by a method. (attention with CPU #)
+
+	- Initialize new FemModel based on new mesh (and copy the old FemModel parameters)
+		It can be done by a method. (attention with CPU #)
+
+	- The last and most difficult thing to do is to update the inputs of the elements of the new mesh:
+		It can be done in just one method, which calls GatherInputs, InterpFromMeshToMesh and InputUpdateFromVector
+
+		- CPU #0 will gather all inputs from FemModel
+        	- call InterpFromMeshToMesh to interpolate them onto the new Mesh
+        	- broadcast the new fields to all cpus
+		- call InputUpdateFromVector so that all inputs are updated
+	
+	- return new FemModel
+	
+	*/
+#endif
+	return output;
+}
+/*}}}*/
+void FemModel::GetMaskLevelSet(IssmDouble **pmasklevelset){/*{{{*/
+
+	int elementswidth=3;//just 2D mesh, tria elements
+	int numberofelements=this->elements->NumberOfElements();
+	int numberofvertices=this->vertices->NumberOfVertices();
+
+	IssmDouble* elementlevelset=xNew<IssmDouble>(elementswidth);
+	int* elem_vertices=xNew<int>(elementswidth);
+	Vector<IssmDouble>* vmasklevelset=new Vector<IssmDouble>(numberofvertices);
 
 	for(int i=0;i<this->elements->Size();i++){
@@ -3083,79 +3318,124 @@
 		element->GetVerticesSidList(elem_vertices);
 		vmasklevelset->SetValue(elem_vertices[0],elementlevelset[0],INS_VAL);
-    	vmasklevelset->SetValue(elem_vertices[1],elementlevelset[1],INS_VAL);
-    	vmasklevelset->SetValue(elem_vertices[2],elementlevelset[2],INS_VAL);
-	}
-
-	/*Assemble*/
+      vmasklevelset->SetValue(elem_vertices[1],elementlevelset[1],INS_VAL);
+      vmasklevelset->SetValue(elem_vertices[2],elementlevelset[2],INS_VAL);
+	}
+
+   /*Assemble*/
 	vmasklevelset->Assemble();
-
-	/*Serialize*/
-	masklevelset = vmasklevelset->ToMPISerial();
-
+	
+	/*Serialize and set output*/
+	(*pmasklevelset)=vmasklevelset->ToMPISerial();
+
+	/*Cleanup*/
 	xDelete<IssmDouble>(elementlevelset);
 	xDelete<int>(elem_vertices);
 	delete vmasklevelset;
 
-	//_printf_("   PRINTING MASKLEVELSET... \n");	    
-	//if(my_rank==0){
-	//	for(int i=0;i<numberofvertices;i++) _printf_("vertex: " << i << "\t" << masklevelset[i] << "\n");
-	//}
-
-	/*Refine the mesh*/
-	double *newx;
-	double *newy;
-	double *newz;
-	int **newelements;
-	int **newsegments;
-	int newnumberofvertices, newnumberofelements, newnumberofsegments;
-	int type_process=1;
-	//AMR->ExecuteRefinement(type_process,vx,vy,masklevelset,newnumberofvertices,newnumberofelements,newnumberofsegments,&newx,&newy,&newz,&newelements,&newsegments);
-	
-	//if(newnumberofvertices<=0 || newnumberofelements<=0 || newnumberofsegments=<0) _error_("Error in the mesh refinement process.");
-
-	delete masklevelset;
-
-	///////////////////////// to test using the current femmodel
-	newx = x;
-	newy = y;
-	newz = z;
-	newelements = elementsptr;
-	newnumberofvertices = numberofvertices2;
-	newnumberofelements = numberofelements2;
-
-	//if(my_rank==0){
-	//	_printf_("   PRINTING COORDINATES... \n");	    
-   // 		for(int i=0;i<newnumberofvertices;i++) _printf_("ID: " << i << "\t" << "X: " << newx[i] << "\t" << "Y: " << newy[i] << "\n");
-   // 	_printf_("   PRINTING ELEMENTS... \n");	   
-   // 	for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelements[i][0] << "\t" << newelements[i][1] << "\t" << newelements[i][2] << "\n");
-	//}
-	/////////////////////////
-
-	/*Partitioning the new mesh*/
-	int        *epart          		= NULL; //element partitioning.
-	int        *npart          		= NULL; //node partitioning.
-	int        *newelementslist     = NULL;
-	int        	edgecut				= 1;
-	int 		numflag				= 0;
-	int 		etype 	 			= 1;
-	int 		numprocs 			= IssmComm::GetSize();
-	bool 		*my_elements 		= NULL;
-	int  		*my_vertices 		= NULL;
-	
-	epart 							= xNew<int>(newnumberofelements);
-	npart 							= xNew<int>(newnumberofvertices);
-	newelementslist 				= xNew<int>(newnumberofelements*elementswidth);
-
-	/*Fill the element list to partitioning*/
-	for(int i=0;i<newnumberofelements;i++){ //itapopo os elementos poderão sair do ExecuteRefinement nesse formato
-		for(int j=0;j<elementswidth;j++){
-			newelementslist[elementswidth*i+j] = newelements[i][j]; //C indexing
-		}
-	}
-
-	//if(my_rank==0){
-	//	_printf_("   PRINTING ELEMENTS in ELEMENTLIST... \n");	   
-   // 	for(int i=0;i<newnumberofelements;i++) _printf_("El: " << i << "\t" << newelementslist[i*elementswidth+0] << "\t" << newelementslist[i*elementswidth+1] << "\t" << newelementslist[i*elementswidth+2] << "\n");
-    //}
+}
+/*}}}*/
+void FemModel::CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices){/*{{{*/
+
+	/*Creating connectivity table*/
+	int* connectivity=xNew<int>(newnumberofvertices);
+	for(int i=0;i<newnumberofvertices;i++) connectivity[i]=0;
+
+	for (int i=0;i<newnumberofelements;i++){
+		for (int j=0;j<elementswidth;j++){
+			int vertexid = newelementslist[elementswidth*i+j];
+			_assert_(vertexid>-1 && vertexid<newnumberofvertices);
+			connectivity[vertexid]+=1;
+		}
+	}	
+
+	/*Create vertex and insert in vertices*/
+	for(int i=0;i<newnumberofvertices;i++){
+		if(my_vertices[i]){
+			Vertex *newvertex=new Vertex();	
+			newvertex->id=i+1;
+			newvertex->sid=i;
+			newvertex->pid=UNDEF;
+			newvertex->x=newx[i];
+			newvertex->y=newy[i];
+			newvertex->z=newz[i];
+			newvertex->domaintype=Domain2DhorizontalEnum;
+			newvertex->sigma=0.;
+			newvertex->connectivity=connectivity[i];
+			vertices->AddObject(newvertex);	
+		} 
+	}
+
+	xDelete<int>(connectivity);
+}
+/*}}}*/
+void FemModel::CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* elements){/*{{{*/
+
+	for(int i=0;i<newnumberofelements;i++){
+		if(my_elements[i]){
+			Tria *newtria=new Tria();
+			newtria->id=i+1;
+			newtria->sid=i;
+			newtria->parameters=NULL;
+			newtria->inputs=new Inputs();
+			newtria->nodes=NULL;
+			newtria->vertices=NULL;
+			newtria->material=NULL;
+			newtria->matpar=NULL;
+			if(this->nummodels>0){
+				newtria->element_type_list=xNew<int>(this->nummodels);
+				for(int j=0;j<nummodels;j++) newtria->element_type_list[j]=0;
+			}
+			else newtria->element_type_list=NULL;
+			/*Element hook*/
+			int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
+			int material_id=i+1; // retrieve material_id = i+1;
+			/*retrieve vertices ids*/
+			int* vertex_ids=xNew<int>(elementswidth);
+			for(int j=0;j<elementswidth;j++)	vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);	
+			/*Setting the hooks*/
+			newtria->numanalyses =this->nummodels;
+			newtria->hnodes		=new Hook*[this->nummodels];
+			newtria->hvertices   =new Hook(&vertex_ids[0],elementswidth);
+			newtria->hmaterial   =new Hook(&material_id,1);
+			newtria->hmatpar     =new Hook(&matpar_id,1);
+			newtria->hneighbors  =NULL;
+			/*Initialize hnodes as NULL*/
+			for(int j=0;j<this->nummodels;j++) newtria->hnodes[j]=NULL;
+			/*Clean up*/
+			xDelete<int>(vertex_ids);
+			elements->AddObject(newtria);	
+		} 
+	}
+}
+/*}}}*/
+void FemModel::CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials){/*{{{*/
+
+	/*Just Matice in this version*/
+	for(int i=0;i<newnumberofelements;i++){
+		if(my_elements[i]){
+			materials->AddObject(new Matice(i+1,i,MaticeEnum));	
+		} 
+	}
+	
+	/*Add new constant material property to materials, at the end: */
+	Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
+	materials->AddObject(newmatpar);//put it at the end of the materials	    
+
+}
+/*}}}*/
+void FemModel::ElementsAndVerticesPartitioning(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices){/*{{{*/
+
+	int *epart=NULL; //element partitioning.
+	int *npart=NULL; //node partitioning.
+	int edgecut=1;
+	int numflag=0;
+	int etype=1;
+	int my_rank = IssmComm::GetRank();
+	int numprocs=IssmComm::GetSize();
+	bool *my_elements=NULL;
+	int *my_vertices=NULL;
+	
+	epart=xNew<int>(newnumberofelements);
+	npart=xNew<int>(newnumberofvertices);
 
 	/*Partition using Metis:*/
@@ -3169,19 +3449,16 @@
 	else if (numprocs==1){
 		/*METIS does not know how to deal with one cpu only!*/
-		for (int i=0;i<newnumberofelements;i++)	epart[i]=0;
-		for (int i=0;i<newnumberofvertices;i++)	npart[i]=0;
-	}
-	else _error_("At least one processor is required");
-
-	//_printf_("   PRINTING AFTER MESTIS... \n");	    
+		for (int i=0;i<newnumberofelements;i++) epart[i]=0;
+		for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
+	}
+	else _error_("At least one processor is required");	    
 
 	my_vertices=xNew<int>(newnumberofvertices);
 	my_elements=xNew<bool>(newnumberofelements);
-	for(int i=0;i<newnumberofvertices;i++) my_vertices[i] = 0;
-	for(int i=0;i<newnumberofelements;i++) my_elements[i] = false;
+	for(int i=0;i<newnumberofvertices;i++) my_vertices[i]=0;
+	for(int i=0;i<newnumberofelements;i++) my_elements[i]=false;
 
 	/*Start figuring out, out of the partition, which elements belong to this cpu: */
 	for(int i=0;i<newnumberofelements;i++){
-
 		/*!All elements have been partitioned above, only deal with elements for this cpu: */
 		if(my_rank==epart[i]){ 
@@ -3197,411 +3474,12 @@
 	}
 
+	/*Free ressources:*/
 	xDelete<int>(epart);
-	xDelete<int>(npart);
-
-	//_printf_("   PRINTING AFTER MY_VERTICES... \n");	    
-
-	/*Creating new femmodel with new mesh*/
-	FemModel* output = NULL;
-	int       analysis_type;
-
-	output = new FemModel(*this);
-
-	//_printf_("   PRINTING AFTER FEMMODEL... \n");	    
-
-	/*Copy basic attributes:*/
-	output->nummodels = this->nummodels;
-	output->solution_type = this->solution_type;
-	output->analysis_counter = this->analysis_counter;
-
-	//_printf_("   PRINTING AFTER BASIC COPIES... \n");	    
-
-	/*Now, deep copy arrays:*/
-	output->analysis_type_list=xNew<int>(this->nummodels);
-	xMemCpy<int>(output->analysis_type_list,this->analysis_type_list,this->nummodels);
-
-	output->profiler=static_cast<Profiler*>(this->profiler->copy());
-	output->parameters=static_cast<Parameters*>(this->parameters->Copy());
-
-	//_printf_("   PRINTING AFTER PARAMETERS... \n");	    
-
-	/*Creating connectivity table*/
-	int* connectivity = xNew<int>(newnumberofvertices);
-	for(int i=0;i<newnumberofvertices;i++) connectivity[i] = 0;
-
-	for (int i=0;i<newnumberofelements;i++){
-		for (int j=0;j<elementswidth;j++){
-			int vertexid = newelementslist[elementswidth*i+j];
-			_assert_(vertexid>-1 && vertexid<newnumberofvertices);
-			connectivity[vertexid]+=1;
-		}
-	}	
-
-	//_printf_("   PRINTING AFTER CONNECTIVITY... \n");	    
-
-	/*Creating vertices*/
-	output->vertices = new Vertices();
-
-	for(int i=0;i<newnumberofvertices;i++){
-		if(my_vertices[i]){
-			Vertex *newvertex = new Vertex();
-			
-			newvertex->id 		= i+1;
-			newvertex->sid 		= i;
-			newvertex->pid 		= UNDEF;
-
-			newvertex->x         	= newx[i];
-			newvertex->y         	= newy[i];
-			newvertex->z         	= newz[i];
-			newvertex->domaintype	= Domain2DhorizontalEnum;
-			newvertex->sigma		= 0.;
-
-			newvertex->connectivity = connectivity[i];
-
-			output->vertices->AddObject(newvertex);	
-		} 
-	}
-
-	xDelete<int>(connectivity);
-
-	//_printf_("   PRINTING AFTER VERTICES... \n");	    
-
-	/*Creating elements*/
-	/*Just Tria in this version*/
-	output->elements = new Elements();
-
-	for(int i=0;i<newnumberofelements;i++){
-		if(my_elements[i]){
-
-			Tria *newtria = new Tria();
-
-			newtria->id  = i+1;
-			newtria->sid = i;
-			newtria->parameters = NULL;
-
-			newtria->inputs  = new Inputs();
-
-			newtria->nodes    = NULL;
-			newtria->vertices = NULL;
-			newtria->material = NULL;
-			newtria->matpar   = NULL;
-
-			if(this->nummodels>0){
-				newtria->element_type_list=xNew<int>(this->nummodels);
-				for(int j=0;j<nummodels;j++) newtria->element_type_list[j] = 0;
-			}
-			else newtria->element_type_list = NULL;
-
-			/*Element hook*/
-			int matpar_id		= newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
-			int material_id 	= i+1; // retrieve material_id = i+1;
-
-			/*retrieve vertices ids*/
-			int* vertex_ids = xNew<int>(elementswidth);
-			
-			for(int j=0;j<elementswidth;j++){ 
-				vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);
-			}
-
-			newtria->numanalyses = this->nummodels;
-			newtria->hnodes      = new Hook*[this->nummodels];
-			newtria->hvertices   = new Hook(&vertex_ids[0],elementswidth);
-			newtria->hmaterial   = new Hook(&material_id,1);
-			newtria->hmatpar     = new Hook(&matpar_id,1);
-			newtria->hneighbors  = NULL;
-
-			/*Initialize hnodes as NULL*/
-			for(int j=0;j<this->nummodels;j++){
-				newtria->hnodes[j]=NULL;
-			}
-
-			/*Clean up*/
-			xDelete<int>(vertex_ids);
-			output->elements->AddObject(newtria);	
-		} 
-	}
-
-	xDelete<int>(newelementslist);
-	//_printf_("   PRINTING AFTER ELEMENTS... \n");	    
-
-	/*Creating materials*/
-	/*Just Matice in this version*/
-	output->materials = new Materials();
-
-	for(int i=0;i<newnumberofelements;i++){
-		if(my_elements[i]){
-			output->materials->AddObject(new Matice(i+1,i,MaticeEnum));	
-		} 
-	}
-	/*This is done to follow CreateElementsVerticesAndMaterials, line 57*/
-	//output->elements->InputDuplicate(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
-
-	/*Add new constant material property to materials, at the end: */
-	Matpar *newmatpar = static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
-	output->materials->AddObject(newmatpar);//put it at the end of the materials
-
-	//_printf_("   PRINTING AFTER MATERIALS... \n");	    
-
-	delete x;
-	delete y;
-	delete z;
-	for(int i=0;i<numberofelements;i++) delete elementsptr[i];
-    if(elementsptr) delete elementsptr;
-
-	//itapopo to test and print the elements and coordinates
-	//output->InitializeAdaptiveRefinement();
-
-	/*Creating nodes*/
-	/*Just SSA (2D) and P1 in this version*/
-	output->nodes = new Nodes();
-
-	int nodecounter = 0;
-	int lid 		= 0;
-	for(int i=0;i<output->nummodels;i++){
-		
-		int analysis_enum = output->analysis_type_list[i];
-
-		//_printf_("   	Analysis type:" << EnumToStringx(analysis_enum) << "\n");
-
-		//itapopo as the domain is 2D, it is not necessary to create nodes for this analysis
-		if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;	    
-
-		for(int j=0;j<newnumberofvertices;j++){
-			if(my_vertices[j]){
-				
-				Node* newnode = new Node();	
-
-				/*id: */
-				newnode->id            = nodecounter+j+1;
-				newnode->sid           = j;
-				newnode->lid           = lid++;
-				newnode->analysis_enum = analysis_enum;
-
-				/*Initialize coord_system: Identity matrix by default*/
-				for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
-				for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
-
-				/*indexing:*/
-				newnode->indexingupdate = true;
-
-				Analysis* analysis = EnumToAnalysis(analysis_enum);
-				int *doftypes = NULL;
-				int numdofs        = analysis->DofsPerNode(&doftypes,Domain2DhorizontalEnum,SSAApproximationEnum);
-				newnode->indexing.Init(numdofs,doftypes);
-				xDelete<int>(doftypes);
-				delete analysis;
-
-				if(analysis_enum==StressbalanceAnalysisEnum)
-				 newnode->SetApproximation(SSAApproximationEnum);
-				else
-				 newnode->SetApproximation(0);
-
-				/*Stressbalance Horiz*/
-				if(analysis_enum==StressbalanceAnalysisEnum){
-					// itapopo
-					/*Coordinate system provided, convert to coord_system matrix*/
-					//XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
-					//_assert_(sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]) >1.e-4);
-
-				}
-
-				output->nodes->AddObject(newnode);
-			}
-		}
-
-		if(output->nodes->Size()) nodecounter = output->nodes->MaximumId();
-	}
-	
-	//_printf_("   Old node size: " << this->nodes->Size() << " \n");	    
-
-	//for(int i=0;i<this->nummodels;i++){
-	//	int analysis_type = this->analysis_type_list[i];
-	//	_printf_("   	Analysis type:" << EnumToStringx(analysis_type) << " size: "<< this->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
-	//}
-	
-	//_printf_("   New number of nodes: " << output->nodes->Size() << " \n");	    
-	//for(int i=0;i<output->nummodels;i++){
-	//	int analysis_type = output->analysis_type_list[i];
-	//	_printf_("   	Analysis type:" << EnumToStringx(analysis_type) << " size: "<< output->nodes->NumberOfNodes(analysis_type) << " NDofs: " <<  this->nodes->NumberOfDofs(analysis_type, GsetEnum) << " \n");
-	//}
-	
-	//_printf_("		Old nodes deep echo: \n");
-	//this->nodes->DeepEcho();
-
-	//_printf_("		New nodes deep echo: \n");
-	//output->nodes->DeepEcho();
-
-	//_printf_("   PRINTING AFTER NODES... \n");	    
-
-	/*Create constraints*/
-	IssmDouble *spcvx = NULL;
-	IssmDouble *spcvy = NULL;
-
-	int numberofnodes_analysistype 	= this->nodes->NumberOfNodes(StressbalanceAnalysisEnum);
-	
-	//_printf_("   Number of nodes: " << numberofnodes_analysistype << " \n");	    
-
-	Vector<IssmDouble>* vspcvx 		= new Vector<IssmDouble>(numberofnodes_analysistype);
-	Vector<IssmDouble>* vspcvy		= new Vector<IssmDouble>(numberofnodes_analysistype);
-
-	IssmDouble BigNumber 			= 1.e8;
-
-	for(int i=0;i<numberofnodes_analysistype;i++){
-		vspcvx->SetValue(i,BigNumber,INS_VAL);
-		vspcvy->SetValue(i,BigNumber,INS_VAL);
-	}
-
-	IssmDouble absmaxspcvx = 0;
-	IssmDouble absmaxspcvy = 0;
-
-	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(); 
-		int nodeindex		= node-1;
-
-		if(dof==0) {
-			vspcvx->SetValue(nodeindex,spcvalue,INS_VAL);
-			if(fabs(spcvalue)>absmaxspcvx) absmaxspcvx = fabs(spcvalue);
-		}
-		else {
-			vspcvy->SetValue(nodeindex,spcvalue,INS_VAL);
-			if(fabs(spcvalue)>absmaxspcvy) absmaxspcvy = fabs(spcvalue);
-		}
-	}
-
-	/*Assemble*/
-	vspcvx->Assemble();
-	vspcvy->Assemble();
-
-	/*Serialize*/
-	spcvx = vspcvx->ToMPISerial();
-	spcvy = vspcvy->ToMPISerial();
-
-	/*Free the data*/
-	delete vspcvx;
-	delete vspcvy;
-	
-	double *newspcvx 	= NULL;
-	double *newspcvy 	= NULL;
-	int *oldelements 	= newelementslist; //itapopo
-	double *oldx		= x; //itapopo
-	double *oldy		= y; //itapopo
-	int nods_data		= numberofnodes_analysistype;
-	int nels_data 		= newnumberofelements;
-	int M_data 			= numberofnodes_analysistype;
-	int N_data  		= 1;
-	int N_interp 		= newnumberofvertices;//itapopo
-	Options *options   	= NULL;
-
-	//itapopo voltar aqui
-	//InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
-	//InterpFromMeshToMesh2dx(&newspcvx,oldelements,oldx,oldy,nods_data,nels_data,spcvx,M_data,N_data,newx,newy,N_interp,options);
-
-	output->constraints = new Constraints();
-
-	nodecounter 			= 0; //itapopo deve começar pelo primeiro nó do StressbalanceAnalysis
-	int count 				= 0;
-	int constraintcounter 	= 0; //itapopo
-	IssmDouble eps			= 1.e-2;
-
-	for(int i=0;i<newnumberofvertices;i++){
-		if(my_vertices[i])
-		/*spcvx*/
-		if(fabs(spcvx[i]) < absmaxspcvx+eps){//itapopo
-			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,0,spcvx[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(fabs(spcvy[i]) < absmaxspcvy+eps){//itapopo
-			output->constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+i+1,1,spcvy[i],StressbalanceAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
-			count++;
-		}
-
-	}
-
-	//_printf_("   Old constraints size: " << this->constraints->Size() << " \n");	    
-	//this->constraints->DeepEcho();
-
-	//_printf_("\n");
-
-	//_printf_("   New constraints size: " << output->constraints->Size() << " \n");
-	
-	//_printf_("SPCVX\n");
-	//for(int i = 0;i<newnumberofvertices;i++){
-	//	_printf_("value: " << spcvx[i] << "\n");
-	//}	    
-	//_printf_("SPCVY\n");
-	//for(int i = 0;i<newnumberofvertices;i++){
-	//	_printf_("value: " << spcvy[i] << "\n");
-	//}
-
-	//output->constraints->DeepEcho();
-
-	//_printf_("   PRINTING AFTER CONSTRAINTS... \n");	    
-
-	//_printf_("   Old loads size: " << this->loads->Size() << " \n");	    
-	//this->loads->DeepEcho();
-
-	//_printf_("   PRINTING AFTER LOADS... \n");	    
-
-	// output->loads=static_cast<Loads*>(this->loads->Copy());
-	// output->constraints=static_cast<Constraints*>(this->constraints->Copy());
-	// output->results=static_cast<Results*>(this->results->Copy());
-
-	// /*reset hooks for elements, loads and nodes: */
-	// output->elements->ResetHooks();
-	// output->loads->ResetHooks();
-	// output->materials->ResetHooks();
-
-	// /*do the post-processing of the datasets to get an FemModel that can actually run analyses: */
-	// for(i=0;i<nummodels;i++){
-	// 	analysis_type=output->analysis_type_list[i];
-	// 	output->SetCurrentConfiguration(analysis_type);
-	// 	if(i==0) VerticesDofx(output->vertices,output->parameters); //only call once, we only have one set of vertices
-	// 	SpcNodesx(output->nodes,output->constraints,output->parameters,analysis_type);
-	// 	NodesDofx(output->nodes,output->parameters,analysis_type);
-	// 	ConfigureObjectsx(output->elements,output->loads,output->nodes,output->vertices,output->materials,output->parameters);
-	// }
-
-	// /*Reset current configuration: */
-	// analysis_type=output->analysis_type_list[analysis_counter];
-	// output->SetCurrentConfiguration(analysis_type);
-
-	/** TODO
-
-	- generate the required input objects for NeoPZ
-		AMR has fathermesh and previousmesh. On first call, these meshes are generated. 
-
-	- call NeoPZ
-		This creates a newmesh which will be refined. 
-
-	- get the new mesh (index,x,y) from NeoPZ
-		Is is doing by GetNewMesh method.
-
-	- Create a new FemModel* that will be consistent with the new mesh
-		It can be done by a method. (attention with CPU #)
-
-	- Initialize new FemModel based on new mesh (and copy the old FemModel parameters)
-		It can be done by a method. (attention with CPU #)
-
-	- The last and most difficult thing to do is to update the inputs of the elements of the new mesh:
-		It can be done in just one method, which calls GatherInputs, InterpFromMeshToMesh and InputUpdateFromVector
-
-		- CPU #0 will gather all inputs from FemModel
-        	- call InterpFromMeshToMesh to interpolate them onto the new Mesh
-        	- broadcast the new fields to all cpus
-		- call InputUpdateFromVector so that all inputs are updated
-	
-	- return new FemModel
-	
-	*/
-
-	return output;
+	xDelete<int>(npart);	    
+
+	/*Assign output pointers:*/
+	*pmy_elements=my_elements;
+	*pmy_vertices=my_vertices;
+
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21503)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 21504)
@@ -149,4 +149,9 @@
 		void InitializeAdaptiveRefinement(void);
 		FemModel* ReMesh(void);
+		void GetMaskLevelSet(IssmDouble** pmasklevelset);
+		void CreateVertices(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,int* my_vertices,IssmDouble* newx,IssmDouble* newy,IssmDouble* newz,Vertices* vertices);
+		void CreateElements(int newnumberofelements,int elementswidth,int* newelementslist,bool* my_elements,Elements* elements);
+		void CreateMaterials(int newnumberofelements,bool* my_elements,Materials* materials);
+		void ElementsAndVerticesPartitioning(int newnumberofvertices,int newnumberofelements,int elementswidth,int* newelementslist,bool** pmy_elements,int** pmy_vertices);
 		#endif
 };