Index: issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3433)
+++ issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3434)
@@ -16,10 +16,6 @@
 
 
-	/*output: int* epart, int* my_grids, double* my_bordergrids*/
-
-
+	/*Intermediary*/
 	int i,j,k,n;
-	extern int my_rank;
-	extern int num_procs;
 
 	/*DataSets: */
@@ -28,95 +24,4 @@
 	DataSet*    vertices = NULL;
 	DataSet*    materials = NULL;
-	
-	/*Objects: */
-	Node*       node   = NULL;
-	NodeProperties node_properties;
-	Vertex*     vertex = NULL;
-	Tria*       tria = NULL;
-	Penta*      penta = NULL;
-	Matice*     matice  = NULL;
-	Matpar*     matpar  = NULL;
-	ElementProperties* tria_properties=NULL;
-	ElementProperties* penta_properties=NULL;
-
-	/*output: */
-	int* epart=NULL; //element partitioning.
-	int* npart=NULL; //node partitioning.
-	int* my_grids=NULL;
-	double* my_bordergrids=NULL;
-
-
-	/*intermediary: */
-	int elements_width; //size of elements
-	double B_avg;
-			
-	/*tria constructor input: */
-	int tria_id;
-	int tria_matice_id;
-	int tria_matpar_id;
-	int tria_numpar_id;
-	int tria_node_ids[3];
-	double tria_h[3];
-	double tria_s[3];
-	double tria_b[3];
-	int    tria_shelf;
-	bool   tria_onwater; 
-
-	/*matice constructor input: */
-	int    matice_mid;
-	double matice_B;
-	double matice_n;
-	
-	/*penta constructor input: */
-
-	int penta_id;
-	int penta_matice_id;
-	int penta_matpar_id;
-	int penta_numpar_id;
-	int penta_node_ids[6];
-	double penta_h[6];
-	double penta_s[6];
-	double penta_b[6];
-	int penta_shelf;
-	int penta_onbed;
-	int penta_onsurface;
-	int penta_collapse;
-	bool   penta_onwater;
-
-	/*matpar constructor input: */
-	int	matpar_mid;
-	double matpar_rho_ice; 
-	double matpar_rho_water;
-	double matpar_heatcapacity;
-	double matpar_thermalconductivity;
-	double matpar_latentheat;
-	double matpar_beta;
-	double matpar_meltingpoint;
-	double matpar_mixed_layer_capacity;
-	double matpar_thermal_exchange_velocity;
-	double matpar_g;
-
-	/* node constructor input: */
-	int node_id;
-	int vertex_id;
-	int partitionborder=0;
-	int node_onbed;
-	int node_onsurface;
-	int node_onshelf;
-	int node_onsheet;
-	int node_upper_node_id;
-	int node_numdofs;
-
-		
-	#ifdef _PARALLEL_
-	/*Metis partitioning: */
-	int  range;
-	Vec  gridborder=NULL;
-	int  my_numgrids;
-	int* all_numgrids=NULL;
-	int  gridcount;
-	int  count;
-	#endif
-	int  first_grid_index;
 
 	/*First create the elements, nodes and material properties: */
@@ -126,37 +31,6 @@
 	materials = new DataSet(MaterialsEnum());
 
-	/*Width of elements: */
-	if(strcmp(iomodel->meshtype,"2d")==0){
-		elements_width=3; //tria elements
-	}
-	else{
-		elements_width=6; //penta elements
-	}
-
-	#ifdef _PARALLEL_
-	/*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
-	if(strcmp(iomodel->meshtype,"2d")==0){
-		/*load elements: */
-		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
-	}
-	else{
-		/*load elements2d: */
-		IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d");
-	}
-
-
-	MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
-
-	/*Free elements and elements2d: */
-	xfree((void**)&iomodel->elements);
-	xfree((void**)&iomodel->elements2d);
-
-	/*Used later on: */
-	my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));
-	#endif
-
-
-
-	/*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
 
 	/*2d mesh: */
@@ -173,67 +47,14 @@
 		for (i=0;i<iomodel->numberofelements;i++){
 
-		#ifdef _PARALLEL_
-		/*!All elements have been partitioned above, only create elements for this CPU: */
-		if(my_rank==epart[i]){ 
-		#endif
-			
-			
-			/*ids: */
-			tria_id=i+1; //matlab indexing.
-			tria_matice_id=-1; //no need for materials
-			tria_matpar_id=-1; //no need for materials
-			tria_numpar_id=1;
+			if(iomodel->my_elements[i]){
 
-			/*vertices offsets: */
-			tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0);
-			tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1);
-			tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2);
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i,iomodel));
 
-			/*thickness,surface and bed:*/
-			tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing.
-			tria_h[1]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-			tria_h[2]=*(iomodel->thickness+  ((int)*(iomodel->elements+elements_width*i+2)-1)) ;
-
-			tria_s[0]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+0)-1)); 
-			tria_s[1]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-			tria_s[2]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+2)-1)); 
-
-			tria_b[0]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+0)-1)); 
-			tria_b[1]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+1)-1)); 
-			tria_b[2]=*(iomodel->bed+        ((int)*(iomodel->elements+elements_width*i+2)-1)); 
-
-			/*element on iceshelf?:*/
-			tria_shelf=(int)*(iomodel->elementoniceshelf+i);
-			tria_onwater=(bool)*(iomodel->elementonwater+i);
-
-			/*Create properties: */
-			tria_properties=new ElementProperties(3,tria_h, tria_s, tria_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, tria_shelf, UNDEF, tria_onwater, UNDEF, UNDEF, UNDEF);
-
-			/*Create tria element using its constructor:*/
-			tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties);
-
-			/*Delete properties: */
-			delete tria_properties;
-
-			/*Add tria element to elements dataset: */
-			elements->AddObject(tria);
-
-			#ifdef _PARALLEL_
-			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
-			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
-			 will hold which grids belong to this partition*/
-			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
-			#endif
-
-		#ifdef _PARALLEL_
-		}//if(my_rank==epart[i])
-		#endif
-
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i,iomodel,3));
+			}
 		}//for (i=0;i<numberofelements;i++)
 
-	
 		/*Free data : */
 		xfree((void**)&iomodel->elements);
@@ -258,61 +79,11 @@
 	
 		for (i=0;i<iomodel->numberofelements;i++){
-		#ifdef _PARALLEL_
-		/*We are using our element partition to decide which elements will be created on this node: */
-		if(my_rank==epart[i]){
-		#endif
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i,iomodel));
 
-			
-			/*name and id: */
-			penta_id=i+1; //matlab indexing.
-			penta_matice_id=-1; 
-			penta_matpar_id=-1; //no need for materials
-			penta_numpar_id=1;
-
-			/*vertices,thickness,surface,bed and drag: */
-			for(j=0;j<6;j++){
-				penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);
-				penta_h[j]=*(iomodel->thickness+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-				penta_s[j]=*(iomodel->surface+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
-				penta_b[j]=*(iomodel->bed+    ((int)*(iomodel->elements+elements_width*i+j)-1)); 
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i,iomodel,6));
 			}
-
-			/*diverse: */
-			penta_shelf=(int)*(iomodel->elementoniceshelf+i);
-			penta_onbed=(int)*(iomodel->elementonbed+i);
-			penta_onsurface=(int)*(iomodel->elementonsurface+i);
-			penta_collapse=1;
-			penta_onwater=(bool)*(iomodel->elementonwater+i);
-	
-
-			/*Create properties: */
-			penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);
-
-			/*Create Penta using its constructor:*/
-			penta=new Penta(penta_id, penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);
-
-			/*Delete properties: */
-			delete penta_properties;
-
-			/*Add penta element to elements dataset: */
-			elements->AddObject(penta);
-	
-			#ifdef _PARALLEL_
-			/*Now that we are here, we can also start building the list of grids belonging to this node partition: we use 
-			 *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
-			 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids 
-			 will hold which grids belong to this partition*/
-			my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
-			my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
-			#endif
-
-		#ifdef _PARALLEL_
-		}//if(my_rank==epart[i])
-		#endif
-
 		}//for (i=0;i<numberofelements;i++)
 
@@ -329,89 +100,7 @@
 	} //if (strcmp(meshtype,"2d")==0)
 
-	#ifdef _PARALLEL_
-		/*From the element partitioning, we can determine which grids are on the inside of this cpu's 
-		 *element partition, and which are on its border with other nodes:*/
-		gridborder=NewVec(iomodel->numberofnodes);
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel));
 
-		for (i=0;i<iomodel->numberofnodes;i++){
-			if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);
-		}
-		VecAssemblyBegin(gridborder);
-		VecAssemblyEnd(gridborder);
-
-		#ifdef _ISSM_DEBUG_
-		VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);
-		#endif
-		
-		VecToMPISerial(&my_bordergrids,gridborder);
-
-		#ifdef _ISSM_DEBUG_
-		if(my_rank==0){
-			for (i=0;i<iomodel->numberofnodes;i++){
-				printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);
-			}
-		}
-		#endif
-	#endif
-
-	/*Add one constant material property to materials: */
-	matpar_mid=1; //put it at the end of the materials
-	matpar_g=iomodel->g; 
-	matpar_rho_ice=iomodel->rho_ice; 
-	matpar_rho_water=iomodel->rho_water; 
-	matpar_thermalconductivity=iomodel->thermalconductivity; 
-	matpar_heatcapacity=iomodel->heatcapacity; 
-	matpar_latentheat=iomodel->latentheat; 
-	matpar_beta=iomodel->beta; 
-	matpar_meltingpoint=iomodel->meltingpoint; 
-	matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; 
-	matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; 
-
-	/*Create matpar object using its constructor: */
-	matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,
-			matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,
-			matpar_thermal_exchange_velocity,matpar_g);
-		
-	/*Add to materials datset: */
-	materials->AddObject(matpar);
-
-	/*Partition penalties in 3d: */
-	if(strcmp(iomodel->meshtype,"3d")==0){
-	
-		/*Get penalties: */
-		IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
-
-		if(iomodel->numpenalties){
-
-			iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));
-			#ifdef _SERIAL_
-			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;
-			#else
-			for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;
-
-			for(i=0;i<iomodel->numpenalties;i++){
-				first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);
-				if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition  grids
-					/*All grids that are being penalised belong to this node's internal grid partition.:*/
-					iomodel->penaltypartitioning[i]=1;
-				}
-				if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border
-					iomodel->penaltypartitioning[i]=0;
-				}
-			}
-			#endif
-		}
-
-		/*Free penalties: */
-		xfree((void**)&iomodel->penalties);
-	}
-
-	/*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 
-	 We can therefore determine  which grids are internal to this node's partition 
-	 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 
-	 that, go and create the grids*/
-
-	/*Create nodes from x,y,z, as well as the spc values on those grids: */
-		
 	/*First fetch data: */
 	if (strcmp(iomodel->meshtype,"3d")==0){
@@ -429,80 +118,17 @@
 	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
 
+	for (i=0;i<iomodel->numberofvertices;i++){
 
-	/*Get number of dofs per node: */
-	DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
 
-	for (i=0;i<iomodel->numberofnodes;i++){
-	#ifdef _PARALLEL_
-	/*keep only this partition's nodes:*/
-	if((my_grids[i]==1)){
-	#endif
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i,iomodel));
 
-		#ifdef _PARALLEL_
-		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
-			partitionborder=1;
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i,iomodel));
+
 		}
-		else{
-			partitionborder=0;
-		}
-		#else
-			partitionborder=0;
-		#endif
-
-		/*create vertex: */
-		vertex_id=i+1;
-		vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder);
-		vertices->AddObject(vertex);
-
-		/*create node: */
-		node_id=i+1; //matlab indexing
-			
-		
-		node_properties.SetProperties(
-				(int)iomodel->gridonbed[i],
-				(int)iomodel->gridonsurface[i],
-				(int)iomodel->gridoniceshelf[i],
-				(int)iomodel->gridonicesheet[i]);
-
-		if (strcmp(iomodel->meshtype,"3d")==0){
-			if (isnan(iomodel->uppernodes[i])){
-				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
-			}
-			else{
-				node_upper_node_id=(int)iomodel->uppernodes[i];
-			}
-		}
-		else{
-			/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
-			node_upper_node_id=node_id;
-		}
-
-		/*Create node using its constructor: */
-		node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);   
-
-
-		/*set single point constraints.: */
-		if (strcmp(iomodel->meshtype,"3d")==0){
-			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
-			if (!iomodel->gridonbed[i]){
-				for(k=1;k<=node_numdofs;k++){
-					node->FreezeDof(k);
-				}
-			}
-		}
-		/*Add node to nodes dataset: */
-		nodes->AddObject(node);
-
-	#ifdef _PARALLEL_
-	} //if((my_grids[i]==1))
-	#endif
 	}
-
-	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
-	 * datasets, it will not be redone: */
-	elements->Presort();
-	vertices->Presort();
-	nodes->Presort();
-	materials->Presort();
 
 	/*Clean fetched data: */
@@ -518,17 +144,11 @@
 	xfree((void**)&iomodel->gridonicesheet);
 	xfree((void**)&iomodel->gridoniceshelf);
-	
 
-	/*Keep partitioning information into iomodel*/
-	iomodel->epart=epart;
-	iomodel->my_grids=my_grids;
-	iomodel->my_bordergrids=my_bordergrids;
-
-	/*Free ressources:*/
-	#ifdef _PARALLEL_
-	xfree((void**)&all_numgrids);
-	xfree((void**)&npart);
-	VecFree(&gridborder);
-	#endif
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	elements->Presort();
+	vertices->Presort();
+	nodes->Presort();
+	materials->Presort();
 
 	cleanup_and_return:
Index: issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 3433)
+++ issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 3434)
@@ -102,5 +102,5 @@
 	} //if (strcmp(meshtype,"2d")==0)
 
-	/*Add new constrant material property tgo materials, at the end: */
+	/*Add new constrant material property to materials, at the end: */
 	materials->AddObject(new Matpar(iomodel));
 		
Index: issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp
===================================================================
--- issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 3433)
+++ issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 3434)
@@ -14,6 +14,5 @@
 void	CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
 
-	/*output: int* epart, int* my_grids, double* my_bordergrids*/
-
+	/*Intermediary*/
 	int i,j,k,n;
 
@@ -43,5 +42,5 @@
 		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
 		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
-		
+
 		for (i=0;i<iomodel->numberofelements;i++){
 
@@ -56,5 +55,5 @@
 		}//for (i=0;i<numberofelements;i++)
 
-	
+
 		/*Free data : */
 		xfree((void**)&iomodel->elements);
@@ -77,5 +76,5 @@
 		IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
 		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
-	
+
 		for (i=0;i<iomodel->numberofelements;i++){
 			if(iomodel->my_elements[i]){
@@ -102,5 +101,5 @@
 	/*Add new constrant material property to materials, at the end: */
 	materials->AddObject(new Matpar(iomodel));
-		
+
 	/*First fetch data: */
 	if (strcmp(iomodel->meshtype,"3d")==0){
@@ -144,5 +143,5 @@
 	xfree((void**)&iomodel->gridonicesheet);
 	xfree((void**)&iomodel->gridoniceshelf);
-	
+
 	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
 	 * datasets, it will not be redone: */
@@ -152,5 +151,5 @@
 	materials->Presort();
 
-	cleanup_and_return:
+cleanup_and_return:
 
 	/*Assign output pointer: */
