Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 2715)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 2715)
@@ -0,0 +1,78 @@
+/*
+ * CreateConstraintsBalancedthickness.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateConstraintsBalancedthickness"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+
+void	CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+
+	int i;
+	int count;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*spc intermediary data: */
+	int spc_sid;
+	int spc_node;
+	int spc_dof;
+	double spc_value;
+	
+	double* spcthickness=NULL;
+	
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum());
+
+	/*Fetch data: */
+	IoModelFetchData(&spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
+
+	count=0;
+
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((iomodel->my_grids[i]==1)){
+	#endif
+
+		if ((int)spcthickness[2*i]){
+	
+			/*This grid needs to be spc'd: */
+
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=1; //we enforce first translation degree of freedom, for temperature
+			spc_value=*(spcthickness+2*i+1);
+
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+		}
+
+	#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: */
+	constraints->Presort();
+
+	/*Free data: */
+	xfree((void**)&spcthickness);
+	
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 2715)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 2715)
@@ -0,0 +1,541 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedthickness.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateElementsNodesAndMaterialsBalancedthickness"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../IoModel.h"
+
+
+void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+
+	/*output: int* epart, int* my_grids, double* my_bordergrids*/
+
+
+	int i,j,k,n;
+	extern int my_rank;
+	extern int num_procs;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    materials = NULL;
+	
+	/*Objects: */
+	Node*       node   = NULL;
+	Tria*       tria = NULL;
+	Penta*      penta = NULL;
+	Matice*     matice  = NULL;
+	Matpar*     matpar  = 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_mid;
+	int tria_mparid;
+	int tria_numparid;
+	int tria_g[3];
+	double tria_h[3];
+	double tria_s[3];
+	double tria_b[3];
+	double tria_k[3];
+	double tria_melting[3];
+	double tria_accumulation[3];
+	double tria_geothermalflux[3];
+	int    tria_friction_type;
+	double tria_p;
+	double tria_q;
+	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_mid;
+	int penta_mparid;
+	int penta_numparid;
+	int penta_g[6];
+	double penta_h[6];
+	double penta_s[6];
+	double penta_b[6];
+	double penta_k[6];
+	int penta_friction_type;
+	double penta_p;
+	double penta_q;
+	int penta_shelf;
+	int penta_onbed;
+	int penta_onsurface;
+	int penta_collapse;
+	double penta_melting[6];
+	double penta_accumulation[6];
+	double penta_geothermalflux[6];
+	int penta_thermal_steadystate;
+	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 node_partitionborder=0;
+	double node_x[3];
+	double node_sigma;
+	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: */
+	elements  = new DataSet(ElementsEnum());
+	nodes     = new DataSet(NodesEnum());
+	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: */
+
+	/*2d mesh: */
+	if (strcmp(iomodel->meshtype,"2d")==0){
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		
+		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_mid=-1; //no need for materials
+			tria_mparid=-1; //no need for materials
+			tria_numparid=1;
+
+			/*vertices offsets: */
+			tria_g[0]=(int)*(iomodel->elements+elements_width*i+0);
+			tria_g[1]=(int)*(iomodel->elements+elements_width*i+1);
+			tria_g[2]=(int)*(iomodel->elements+elements_width*i+2);
+
+			/*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 tria element using its constructor:*/
+			tria=new Tria(tria_id, tria_mid, tria_mparid, tria_numparid,tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_onwater);
+
+			/*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
+
+		}//for (i=0;i<numberofelements;i++)
+
+	
+		/*Free data : */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
+		IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
+		IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	
+		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
+
+			
+			/*name and id: */
+			penta_id=i+1; //matlab indexing.
+			penta_mid=-1; 
+			penta_mparid=-1; //no need for materials
+			penta_numparid=1;
+
+			/*vertices,thickness,surface,bed and drag: */
+			for(j=0;j<6;j++){
+				penta_g[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)); 
+			}
+
+			/*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 Penta using its constructor:*/
+			penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type,
+					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,
+					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,
+					penta_thermal_steadystate,penta_onwater); 
+
+			/*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++)
+
+		/*Free data: */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonbed);
+		xfree((void**)&iomodel->elementonsurface);
+		xfree((void**)&iomodel->elementonwater);
+
+	} //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);
+
+		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){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+	IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x");
+	IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y");
+	IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+	IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
+	IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+
+
+	/*Get number of dofs per node: */
+	DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);
+
+	for (i=0;i<iomodel->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((my_grids[i]==1)){
+	#endif
+
+		node_id=i+1; //matlab indexing
+			
+		
+		
+		#ifdef _PARALLEL_
+		if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border
+			node_partitionborder=1;
+		}
+		else{
+			node_partitionborder=0;
+		}
+		#else
+			node_partitionborder=0;
+		#endif
+
+		node_x[0]=iomodel->x[i];
+		node_x[1]=iomodel->y[i];
+		node_x[2]=iomodel->z[i];
+		node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]);
+		
+		node_onbed=(int)iomodel->gridonbed[i];
+		node_onsurface=(int)iomodel->gridonsurface[i];	
+		node_onshelf=(int)iomodel->gridoniceshelf[i];	
+		node_onsheet=(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,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet);
+
+		/*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();
+	nodes->Presort();
+	materials->Presort();
+
+	/*Clean fetched data: */
+	xfree((void**)&iomodel->deadgrids);
+	xfree((void**)&iomodel->x);
+	xfree((void**)&iomodel->y);
+	xfree((void**)&iomodel->z);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->gridonbed);
+	xfree((void**)&iomodel->gridonsurface);
+	xfree((void**)&iomodel->uppernodes);
+	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
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 2715)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 2715)
@@ -0,0 +1,35 @@
+/*! \file CreateLoadsBalancedthickness.c:
+ */
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateLoadsBalancedthickness"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/macros.h"
+#include "../IoModel.h"
+
+
+void	CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum());
+	
+	
+	/*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: */
+	loads->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
+
+
Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp	(revision 2715)
+++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp	(revision 2715)
@@ -0,0 +1,90 @@
+/*!\file: CreateParametersBalancedthickness.cpp
+ * \brief driver for creating parameters dataset, for prognostic analysis.
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "CreateParametersBalancedthickness"
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../IoModel.h"
+
+void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	Param*   param = NULL;
+	DataSet* parameters=NULL;
+	int      count;
+	int      i;
+	int      dim;
+
+	double* vx=NULL;
+	double* vy=NULL;
+	double* vz=NULL;
+	double* u_g=NULL;
+	double* pressure=NULL;
+	double* temperature=NULL;
+	double* thickness=NULL;
+	double* surface=NULL;
+	double* bed=NULL;
+	double* accumulation=NULL;
+	double* melting=NULL;
+
+	/*recover parameters : */
+	parameters=*pparameters;
+
+	count=parameters->Size();
+
+	/*Get vx and vy: */
+	IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz");
+
+	u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double));
+
+	if(vx) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts;
+	if(vy) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts;
+	if(vz) for(i=0;i<iomodel->numberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts;
+
+	count++;
+	param= new Param(count,"u_g",DOUBLEVEC);
+	param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3);
+	parameters->AddObject(param);
+
+
+	xfree((void**)&vx);
+	xfree((void**)&vy);
+	xfree((void**)&vz);
+	xfree((void**)&u_g);
+
+	/*Get melting: */
+	IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting");
+	if(melting) for(i=0;i<iomodel->numberofnodes;i++)melting[i]=melting[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"m_g",DOUBLEVEC);
+	if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1);
+	else param->SetDoubleVec(melting,0,1);
+	parameters->AddObject(param);
+
+	/*Free melting: */
+	xfree((void**)&melting);
+
+	/*Get accumulation: */
+	IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation");
+	if(accumulation) for(i=0;i<iomodel->numberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts;
+	
+	count++;
+	param= new Param(count,"a_g",DOUBLEVEC);
+	if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1);
+	else param->SetDoubleVec(accumulation,0,0);
+	parameters->AddObject(param);
+
+	/*Free accumulation: */
+	xfree((void**)&accumulation);
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: /issm/trunk/src/c/parallel/balancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 2715)
+++ /issm/trunk/src/c/parallel/balancedthickness.cpp	(revision 2715)
@@ -0,0 +1,160 @@
+/*!\file:  balancedthickness.cpp
+ * \brief: balancedthickness solution
+ */ 
+
+#include "../issm.h"
+#include "./parallel.h"
+
+#undef __FUNCT__ 
+#define __FUNCT__ "balancedthickness"
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+
+int main(int argc,char* *argv){
+	
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	int   numberofnodes;
+	double waitonlock=0;
+
+	Model* model=NULL;
+
+	Vec     u_g=NULL;
+	double* u_g_serial=NULL;
+	double* melting_g=NULL;
+	double* accumulation_g=NULL;
+	double  dt;
+	double  yts;
+	int     qmu_analysis;
+
+	/*Results: */
+	DataSet* results=NULL;
+	DataSet* processedresults=NULL;
+	Result*  result=NULL;
+
+	ParameterInputs* inputs=NULL;
+	Param*   param=NULL;
+
+	/*time*/
+	double   start, finish;
+	double   start_core, finish_core;
+	double   start_init, finish_init;
+
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	/*Initialize Petsc and get start time*/
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Initialize model structure: */
+	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
+	model=new Model();
+
+	/*Open handle to data on disk: */
+	fid=pfopen(inputfilename,"rb");
+
+	/*Initialize model structure: */
+	model=new Model();
+
+	_printf_("read and create finite element model:\n");
+	model->AddFormulation(fid,BalancedthicknessAnalysisEnum());
+
+	/*recover parameters: */
+	model->FindParam(&waitonlock,"waitonlock");
+	model->FindParam(&qmu_analysis,"qmu_analysis");
+
+	_printf_("initialize inputs:\n");
+	
+	model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedthicknessAnalysisEnum());
+	model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedthicknessAnalysisEnum());
+	model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedthicknessAnalysisEnum());
+	model->FindParam(&numberofnodes,"numberofnodes");
+	
+	inputs=new ParameterInputs;
+	inputs->Add("velocity",u_g_serial,3,numberofnodes);
+	inputs->Add("melting",melting_g,1,numberofnodes);
+	inputs->Add("accumulation",accumulation_g,1,numberofnodes);
+
+	_printf_("initialize results:\n");
+	results=new DataSet(ResultsEnum());
+	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
+
+	/*are we running the solutoin sequence, or a qmu wrapper around it? : */
+	if(!qmu_analysis){
+
+		/*run balancedthickness analysis: */
+		_printf_("call computational core:\n");
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		balancedthickness_core(results,model,inputs);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+	}
+	else{
+
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on balancedthickness core:\n");
+	
+		#ifdef _HAVE_DAKOTA_ 
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		Qmux(model,inputs,BalancedthicknessAnalysisEnum(),NoneAnalysisEnum());
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	 	#else
+		throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	/*Add analysis_type to results: */
+	result=new Result(results->Size()+1,0,1,"analysis_type","balancedthickness");
+	results->AddObject(result);
+	
+	_printf_("process results:\n");
+	ProcessResults(&processedresults,results,model,BalancedthicknessAnalysisEnum());
+	
+	_printf_("write results to disk:\n");
+	OutputResults(processedresults,outputfilename);
+
+	if (waitonlock>0){
+		_printf_("write lock file:\n");
+		WriteLockFile(lockname);
+	}
+
+	/*Free ressources:*/
+	delete processedresults;
+	delete results;
+	delete model;
+	delete inputs;
+
+	/*Get finish time and close*/
+	MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
+	_printf_("\n   %-34s %f seconds  \n","Model initialization elapsed time:",finish_init-start_init);
+	_printf_("   %-34s %f seconds  \n","Core solution elapsed time:",finish_core-start_core);
+	_printf_("   %-34s %f seconds\n\n","Total elapsed time:",finish-start);
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+
+	/*end module: */
+	MODULEEND();
+
+	return 0; //unix success return;
+}
Index: /issm/trunk/src/c/parallel/balancedthickness_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/balancedthickness_core.cpp	(revision 2715)
+++ /issm/trunk/src/c/parallel/balancedthickness_core.cpp	(revision 2715)
@@ -0,0 +1,64 @@
+/*!\file: balancedthickness_core.cpp
+ * \brief: core of the balancedthickness solution 
+ */ 
+
+#undef __FUNCT__ 
+#define __FUNCT__ "balancedthickness_core"
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./parallel.h"
+#include "../issm.h"
+
+void balancedthickness_core(DataSet* results,Model* model,ParameterInputs* inputs){
+
+	extern int my_rank;
+
+	/*output: */
+	Result* result=NULL;
+
+	/*intermediary: */
+	Vec u_g=NULL;
+
+	/*solutions: */
+	Vec h_g=NULL;
+
+	/*flags: */
+	int verbose=0;
+	int numberofdofspernode;
+	int numberofnodes;
+	int dofs[2]={1,1};
+
+	/*fem balancedthickness model: */
+	FemModel* fem_p=NULL;
+
+
+	/*recover fem model: */
+	fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum());
+
+	//first recover parameters common to all solutions
+	model->FindParam(&verbose,"verbose");
+	model->FindParam(&numberofnodes,"numberofnodes");
+	model->FindParam(&numberofdofspernode,"numberofdofspernode");
+
+	_printf_("depth averaging velocity...\n");
+	u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
+	FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
+	inputs->Add("velocity_average",u_g,2,numberofnodes);
+	
+	_printf_("call computational core:\n");
+	diagnostic_core_linear(&h_g,fem_p,inputs,BalancedthicknessAnalysisEnum(),NoneAnalysisEnum());
+
+	_printf_("extrude computed thickness on all layers:\n");
+	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+
+	/*Plug results into output dataset: */
+	result=new Result(results->Size()+1,0,1,"h_g",h_g);
+	results->AddObject(result);
+
+	/*Free ressources:*/
+	VecFree(&u_g);
+	VecFree(&h_g);
+}
Index: /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m
===================================================================
--- /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m	(revision 2715)
+++ /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m	(revision 2715)
@@ -0,0 +1,9 @@
+function macro=BalancedthicknessAnalysisEnum()
+%BALANCEDTHICKNESSANALYSISENUM - Enum of BalancedthicknessAnalysis
+%
+%   file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum
+%
+%   Usage:
+%      macro=BalancedthicknessAnalysisEnum()
+
+macro=251;
Index: /issm/trunk/src/m/solutions/jpl/balancedthickness.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/balancedthickness.m	(revision 2715)
+++ /issm/trunk/src/m/solutions/jpl/balancedthickness.m	(revision 2715)
@@ -0,0 +1,35 @@
+function md=balancedthickness(md)
+%BALANCEDTHICKNESS - balancedthickness solution sequence.
+%
+%   Usage:
+%      md=balancedthickness(md)
+
+	%timing
+	t1=clock;
+
+	%Build all models requested for diagnostic simulation
+	models.analysis_type=BalancedthicknessAnalysisEnum; %needed for processresults
+	
+	displaystring(md.verbose,'%s',['reading balancedthickness model data']);
+	md.analysis_type=BalancedthicknessAnalysisEnum; models.bt=CreateFemModel(md);
+
+	% figure out number of dof: just for information purposes.
+	md.dof=modelsize(models);
+
+	%initialize inputs
+	displaystring(md.verbose,'\n%s',['setup inputs...']);
+	inputs=inputlist;
+	inputs=add(inputs,'velocity',models.bt.parameters.u_g,'doublevec',3,models.bt.parameters.numberofnodes);
+	inputs=add(inputs,'melting',models.bt.parameters.m_g,'doublevec',1,models.bt.parameters.numberofnodes);
+	inputs=add(inputs,'accumulation',models.bt.parameters.a_g,'doublevec',1,models.bt.parameters.numberofnodes);
+
+	displaystring(md.verbose,'\n%s',['call computational core:']);
+	results=balancedthickness_core(models,inputs,BalancedthicknessAnalysisEnum(),NoneAnalysisEnum());
+
+	displaystring(md.verbose,'\n%s',['load results...']);
+	if ~isstruct(md.results), md.results=struct(); end
+	md.results.balancedthickness=processresults(models,results);
+
+	%stop timing
+	t2=clock;
+	displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);	
Index: /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m
===================================================================
--- /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m	(revision 2715)
+++ /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m	(revision 2715)
@@ -0,0 +1,24 @@
+function results=balancedthickness_core(models,inputs,analysis_type,sub_analysis_type)
+%BALANCEDTHICKNESS_CORE - linear solution sequence
+%
+%   Usage:
+%      h_g=balancedthickness_core(m,inputs,analysis_type,sub_analysis_type)
+
+	%get FE model
+	m=models.bt;
+	results.time=0;
+	results.step=1;
+
+	displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']);
+	%Take only the first two dofs of m.parameters.u_g
+	u_g=get(inputs,'velocity',[1 1 0 0]);
+	u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
+	inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes);
+
+	displaystring(m.parameters.verbose,'\n%s',['call computational core:']);
+	results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type);
+
+	displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']);
+	results.h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0);
+
+end %end function