Index: /issm/trunk/src/c/modules/Misfitx/Misfitx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Misfitx/Misfitx.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/Misfitx/Misfitx.cpp	(revision 3905)
@@ -0,0 +1,46 @@
+/*!\file Misfitx
+ * \brief: compute misfit between observations and model
+ */
+
+#include "./Misfitx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../SurfaceAreax/SurfaceAreax.h"
+
+void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,Parameters* parameters,
+			int analysis_type,int sub_analysis_type){
+	
+	/*Intermediary*/
+	int fit;
+	double S;
+
+	/*output: */
+	double J;
+	double J_sum;
+	
+	/*First, get elements and loads configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements,loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Compute surface area: */
+	SurfaceAreax(&S,elements,nodes,vertices, loads,materials,parameters,analysis_type,sub_analysis_type);
+
+	/*add surface area to elements :*/
+	elements->UpdateInputsFromVector(&S,SurfaceAreaEnum,ConstantEnum);
+	
+
+	/*Compute Misfit: */
+	elements->Misfit(&J,analysis_type,sub_analysis_type);
+
+	/*Sum all J from all cpus of the cluster:*/
+	MPI_Reduce (&J,&J_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD );
+	MPI_Bcast(&J_sum,1,MPI_DOUBLE,0,MPI_COMM_WORLD); 
+	J=J_sum;
+
+	/*Assign output pointers: */
+	*pJ=J;
+}
Index: /issm/trunk/src/c/modules/Misfitx/Misfitx.h
===================================================================
--- /issm/trunk/src/c/modules/Misfitx/Misfitx.h	(revision 3905)
+++ /issm/trunk/src/c/modules/Misfitx/Misfitx.h	(revision 3905)
@@ -0,0 +1,16 @@
+/*!\file:  Misfitx.h
+ * \brief header file for inverse methods misfit computation
+ */ 
+
+#ifndef _MISFITX_H
+#define _MISFITX_H
+
+#include "../DataSet/DataSet.h"
+#include "../objects/objects.h"
+
+/* local prototypes: */
+void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, 
+			int analysis_type,int sub_analysis_type);
+
+#endif  /* _MISFITX_H */
+
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp	(revision 3905)
@@ -0,0 +1,49 @@
+/*
+ * CreateConstraintsBalancedthickness.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsBalancedthickness(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count=0;
+	DataSet* constraints = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
+
+	count=1; //matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		if(iomodel->my_vertices[i]){
+			
+			if ((int)iomodel->spcthickness[2*i]){
+		
+				constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));//we enforce first translation degree of freedom, for temperature
+				count++;
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->spcthickness);
+
+	/*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();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp	(revision 3905)
@@ -0,0 +1,186 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedthickness.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+		}//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);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->dhdt);
+
+	}
+	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");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
+		IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+		IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+		for (i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,6));
+			}
+		}//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);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->dhdt);
+		xfree((void**)&iomodel->upperelements);
+		xfree((void**)&iomodel->lowerelements);
+
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));
+
+	/*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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+
+	/*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:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp	(revision 3905)
@@ -0,0 +1,30 @@
+/*! \file CreateLoadsBalancedthickness.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.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/modules/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateConstraintsBalancedthickness2.cpp	(revision 3905)
@@ -0,0 +1,22 @@
+/*
+ * CreateConstraintsBalancedthickness2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../io/io.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsBalancedthickness2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	DataSet* constraints = NULL;
+	
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp	(revision 3905)
@@ -0,0 +1,163 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedthickness2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j;
+	int vertex_index;
+	int node_index;
+
+	/*DataSets: */
+	DataSet* elements  = NULL;
+	DataSet* nodes = NULL;
+	DataSet* vertices = NULL;
+	DataSet* materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+		}//for (i=0;i<numberofelements;i++)
+	
+		/*Free data : */
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->dhdt);
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+		ISSMERROR("not implemented yet");
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));
+
+	/*Create nodes and vertices: */
+	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->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+
+	/*Build Vertices dataset*/
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+		}
+	}
+
+	/*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
+	for (i=0;i<iomodel->numberofelements;i++){
+		for (j=0;j<3;j++){
+
+			if(iomodel->my_nodes[3*i+j]){ 
+
+				//Get index of the vertex on which the current node is located
+				vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
+				ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);
+
+				//Compute Node index (id-1)
+				node_index=3*i+j;
+
+				/*Add node to nodes dataset: */
+				nodes->AddObject(new Node(i+1,vertex_index,node_index,iomodel));
+
+			}
+		}
+	}
+
+	/*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->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	xfree((void**)&iomodel->gridonicesheet);
+	xfree((void**)&iomodel->gridoniceshelf);
+	xfree((void**)&iomodel->elements);
+
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateLoadsBalancedthickness2.cpp	(revision 3905)
@@ -0,0 +1,46 @@
+/*! \file CreateLoadsBalancedthickness2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i;
+
+	/*Output*/
+	DataSet* loads=NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Get edges and elements*/
+	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+
+	/*First load data:*/
+	for (i=0;i<iomodel->numberofedges;i++){
+		loads->AddObject(new Numericalflux(i+1,i,iomodel));
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->edges);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->thickness);
+
+	/*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();
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp	(revision 3905)
@@ -0,0 +1,54 @@
+/*
+ * CreateConstraintsBalancedvelocities.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsBalancedvelocities(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count=1;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
+
+	count=1; //matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+
+			if ((int)iomodel->spcvelocity[6*i+0] && (int)iomodel->spcvelocity[6*i+1]){ //spc if vx and vy are constrained
+
+				/*This grid needs to be spc'd: */
+				constraints->AddObject(new Spc(count,i+1,1,pow( pow(*(iomodel->spcvelocity+6*i+4),2.0) + pow(*(iomodel->spcvelocity+6*i+5),2.0) ,0.5)));
+				count++;
+			}
+
+		} //if((my_vertices[i]==1))
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->spcvelocity);
+	
+	/*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();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp	(revision 3905)
@@ -0,0 +1,186 @@
+/*
+ * CreateElementsNodesAndMaterialsBalancedvelocities.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"accumulation_rate");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+
+		}//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);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+
+	}
+	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");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+		IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+
+		for (i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,6));
+
+			}//if(my_elements[i])
+
+		}//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);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->upperelements);
+		xfree((void**)&iomodel->lowerelements);
+
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials
+		
+	/*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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+			
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+	
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp	(revision 3905)
@@ -0,0 +1,30 @@
+/*! \file CreateLoadsBalancedvelocities.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+
+void	CreateLoadsBalancedvelocities(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/modules/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 3905)
@@ -0,0 +1,75 @@
+/*!\file: CreateParametersControl.cpp
+ * \brief driver for creating parameters dataset, for control analysis.
+ */ 
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	int i;
+	double* control_parameter=NULL;
+	
+	Parameters* parameters=NULL;
+
+	/*Get parameters: */
+	parameters=*pparameters;
+	
+	parameters->AddObject(new   BoolParam(ControlAnalysisEnum,iomodel->control_analysis));
+	if(iomodel->control_analysis){
+
+		if (strcmp(iomodel->control_type,"drag")==0){
+			parameters->AddObject(new BoolParam(ExtrudeParamEnum,false));
+			parameters->AddObject(new IntParam(ControlTypeEnum,DragCoefficientEnum));
+		}
+		else  if (strcmp(iomodel->control_type,"B")==0){
+			parameters->AddObject(new BoolParam(ExtrudeParamEnum,true));
+			parameters->AddObject(new IntParam(ControlTypeEnum,RheologyBEnum));
+		}
+		else
+		 ISSMERROR("control_type %s not supported yet!",iomodel->control_type);
+
+		parameters->AddObject(new BoolParam(ControlSteadyEnum,false));
+		parameters->AddObject(new IntParam(NStepsEnum,iomodel->nsteps));
+		parameters->AddObject(new DoubleParam(TolXEnum,iomodel->tolx));
+		parameters->AddObject(new DoubleParam(EpsCmEnum,iomodel->eps_cm));
+		parameters->AddObject(new DoubleParam(MeanVelEnum,iomodel->meanvel));
+		parameters->AddObject(new DoubleParam(CmNoiseDmpEnum,iomodel->cm_noisedmp));
+		parameters->AddObject(new DoubleParam(CmMinDmpValueEnum,iomodel->cm_mindmp_value));
+		parameters->AddObject(new DoubleParam(CmMinDmpSlopeEnum,iomodel->cm_mindmp_slope));
+		parameters->AddObject(new DoubleParam(CmMaxDmpValueEnum,iomodel->cm_maxdmp_value));
+		parameters->AddObject(new DoubleParam(CmMaxDmpSlopeEnum,iomodel->cm_maxdmp_slope));
+		parameters->AddObject(new DoubleParam(CmMinEnum,iomodel->cm_min));
+		parameters->AddObject(new DoubleParam(CmMaxEnum,iomodel->cm_max));
+		parameters->AddObject(new BoolParam(CmGradientEnum,iomodel->cm_gradient));
+
+		/*Now, recover fit, optscal and maxiter as vectors: */
+		IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit");
+		IoModelFetchData(&iomodel->cm_jump,NULL,NULL,iomodel_handle,"cm_jump");
+		IoModelFetchData(&iomodel->optscal,NULL,NULL,iomodel_handle,"optscal");
+		IoModelFetchData(&iomodel->maxiter,NULL,NULL,iomodel_handle,"maxiter");
+
+		parameters->AddObject(new DoubleVecParam(FitEnum,iomodel->fit,iomodel->nsteps));
+		parameters->AddObject(new DoubleVecParam(CmJumpEnum,iomodel->cm_jump,iomodel->nsteps));
+		parameters->AddObject(new DoubleVecParam(OptScalEnum,iomodel->optscal,iomodel->nsteps));
+		parameters->AddObject(new DoubleVecParam(MaxIterEnum,iomodel->maxiter,iomodel->nsteps));
+
+		xfree((void**)&iomodel->fit);
+		xfree((void**)&iomodel->cm_jump);
+		xfree((void**)&iomodel->optscal);
+		xfree((void**)&iomodel->maxiter);
+
+		IoModelFetchData(&control_parameter,NULL,NULL,iomodel_handle,iomodel->control_type);
+		parameters->AddObject(new DoubleVecParam(ControlParameterEnum,control_parameter,iomodel->nsteps));
+		xfree((void**)&control_parameter);
+	}
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 3905)
@@ -0,0 +1,108 @@
+/*!\file: CreateDataSets
+ * \brief general driver for creating all datasets that make a finite element iomodel
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./ModelProcessorx/ModelProcessorx.h"
+
+
+void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*create parameters common to all solutions: */
+	CreateParameters(pparameters,iomodel,iomodel_handle);
+	CreateParametersControl(pparameters,iomodel,iomodel_handle);
+
+	/*This is just a high level driver: */
+	if (iomodel->analysis_type==DiagnosticAnalysisEnum){
+
+		if (iomodel->sub_analysis_type==HorizAnalysisEnum){
+
+			CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
+			CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
+			CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
+		}
+		else if (iomodel->sub_analysis_type==VertAnalysisEnum){
+		
+			CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
+			CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle);
+			CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle);
+		}
+		else if (iomodel->sub_analysis_type==StokesAnalysisEnum){
+
+			CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle);
+			CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle);
+			CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle);
+		}
+		else if (iomodel->sub_analysis_type==HutterAnalysisEnum){
+
+			CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+			CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle);
+			CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle);
+		}
+	}
+	else if (iomodel->analysis_type==SlopecomputeAnalysisEnum){
+
+		CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes, pvertices,pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==ThermalAnalysisEnum){
+
+		CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsThermal(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsThermal(ploads,iomodel,iomodel_handle);
+					
+	}
+	else if (iomodel->analysis_type==MeltingAnalysisEnum){
+			
+		CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsMelting(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsMelting(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==PrognosticAnalysisEnum){
+
+		CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsPrognostic(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsPrognostic(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==Prognostic2AnalysisEnum){
+
+		CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum){
+
+		CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==Balancedthickness2AnalysisEnum){
+
+		CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle);
+	}
+	else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum){
+
+		CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pvertices,pmaterials, iomodel,iomodel_handle);
+		CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle);
+		CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle);
+	}
+	else{
+
+		ISSMERROR("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!");
+
+	}
+	CreateParametersQmu(pparameters,iomodel,iomodel_handle);
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 3905)
@@ -0,0 +1,119 @@
+/*!\file: CreateParameters.cpp
+ * \brief general driver for creating parameters dataset
+ */ 
+
+#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 "../io/io.h"
+
+void CreateParameters(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	int i;
+	
+	Parameters* parameters = NULL;
+	int      numberofdofspernode;
+	char**   parameteroutput=NULL;
+	char*    descriptor=NULL;
+	char*    tag=NULL;
+	#ifdef _SERIAL_
+		mxArray* pfield=NULL;
+		mxArray* pfield2=NULL;
+	#endif
+
+	/*Initialize dataset: */
+	parameters = new Parameters(ParametersEnum);
+
+	if (strcmp(iomodel->meshtype,"2d")==0) parameters->AddObject(new IntParam(DimEnum,2));
+	else parameters->AddObject(new IntParam(DimEnum,3));
+	parameters->AddObject(new    IntParam(AnalysisTypeEnum,iomodel->analysis_type));
+	parameters->AddObject(new    IntParam(SubAnalysisTypeEnum,iomodel->sub_analysis_type));
+	parameters->AddObject(new StringParam(OutputFileNameEnum,iomodel->outputfilename));
+	parameters->AddObject(new   BoolParam(IsHutterEnum,iomodel->ishutter));
+	parameters->AddObject(new   BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn));
+	parameters->AddObject(new   BoolParam(IsStokesEnum,iomodel->isstokes));
+	parameters->AddObject(new    IntParam(VerboseEnum,iomodel->verbose));
+	parameters->AddObject(new DoubleParam(EpsResEnum,iomodel->eps_res));
+	parameters->AddObject(new DoubleParam(EpsRelEnum,iomodel->eps_rel));
+	parameters->AddObject(new DoubleParam(EpsAbsEnum,iomodel->eps_abs));
+	parameters->AddObject(new    IntParam(MaxNonlinearIterationsEnum,iomodel->max_nonlinear_iterations));
+	parameters->AddObject(new DoubleParam(EpsVelEnum,iomodel->epsvel));
+	parameters->AddObject(new DoubleParam(YtsEnum,iomodel->yts));
+	parameters->AddObject(new DoubleParam(DtEnum,iomodel->dt*iomodel->yts)); //Dt is in years in iomodel
+	parameters->AddObject(new DoubleParam(NdtEnum,iomodel->ndt*iomodel->yts)); //Ndt is in yeats in iomodel
+	parameters->AddObject(new DoubleParam(PenaltyOffsetEnum,iomodel->penalty_offset));
+	parameters->AddObject(new DoubleParam(SparsityEnum,iomodel->sparsity));
+	parameters->AddObject(new   BoolParam(LowmemEnum,iomodel->lowmem));
+	parameters->AddObject(new    IntParam(ConnectivityEnum,iomodel->connectivity));
+	parameters->AddObject(new DoubleParam(BetaEnum,iomodel->beta));
+	parameters->AddObject(new DoubleParam(MeltingPointEnum,iomodel->meltingpoint));
+	parameters->AddObject(new DoubleParam(LatentHeatEnum,iomodel->latentheat));
+	parameters->AddObject(new DoubleParam(HeatCapacityEnum,iomodel->heatcapacity));
+	parameters->AddObject(new   BoolParam(ArtDiffEnum,iomodel->artdiff));
+	parameters->AddObject(new DoubleParam(PenaltyMeltingEnum,iomodel->penalty_melting));
+	parameters->AddObject(new    IntParam(MinThermalConstraintsEnum,iomodel->min_thermal_constraints));
+	parameters->AddObject(new    IntParam(StabilizeConstraintsEnum,iomodel->stabilize_constraints));
+	parameters->AddObject(new DoubleParam(StokesReconditioningEnum,iomodel->stokesreconditioning));
+	parameters->AddObject(new DoubleParam(ViscosityOvershootEnum,iomodel->viscosity_overshoot));
+	parameters->AddObject(new   BoolParam(WaitOnLockEnum,iomodel->waitonlock));
+	parameters->AddObject(new StringParam(SolverStringEnum,iomodel->solverstring));
+	parameters->AddObject(new    IntParam(NumberOfVerticesEnum,iomodel->numberofvertices));
+	parameters->AddObject(new    IntParam(NumberOfElementsEnum,iomodel->numberofelements));
+
+	/*Deal with more complex parameters*/
+	if (
+				iomodel->analysis_type==Prognostic2AnalysisEnum ||
+				iomodel->analysis_type==Balancedthickness2AnalysisEnum
+				)
+		parameters->AddObject(new    IntParam(NumberOfNodesEnum,3*iomodel->numberofelements));
+	else parameters->AddObject(new    IntParam(NumberOfNodesEnum,iomodel->numberofvertices));
+
+	DistributeNumDofs(&numberofdofspernode,iomodel->analysis_type,iomodel->sub_analysis_type);
+	parameters->AddObject(new    IntParam(NumberOfDofsPerNodeEnum,numberofdofspernode));
+
+	IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo");
+	parameters->AddObject(new    IntParam(NumRiftsEnum,iomodel->numrifts));
+	xfree((void**)&iomodel->riftinfo); 
+
+	parameters->AddObject(new    IntParam(NumOutputEnum,iomodel->numoutput));
+	if(iomodel->numoutput){
+		parameteroutput=(char**)xmalloc(iomodel->numoutput*sizeof(char*));
+
+		#ifdef _SERIAL_
+		pfield=mxGetField(iomodel_handle,0,"parameteroutput");
+		for(i=0;i<iomodel->numoutput;i++){
+			pfield2=mxGetCell(pfield,i);
+			FetchData(&descriptor,pfield2);
+			parameteroutput[i]=descriptor;
+		}
+		#else
+		tag=(char*)xmalloc((strlen("parameteroutput_i")+1)*sizeof(char));
+		for(i=0;i<iomodel->numoutput;i++){
+			sprintf(tag,"%s%i","parameteroutput_",i);
+			IoModelFetchData(&descriptor,iomodel_handle,tag);
+			parameteroutput[i]=descriptor;
+		}
+		#endif
+
+		/*Ok, we have all the parameter output  descriptors. Build a parameter with it: */
+		parameters->AddObject(new StringArrayParam(ParameterOutputEnum,parameteroutput,iomodel->numoutput));
+	}
+
+	/*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	parameters->Presort();
+
+	/*Free data: */
+	xfree((void**)&tag);
+	for(i=0;i<iomodel->numoutput;i++){
+		char* descriptor=parameteroutput[i];
+		xfree((void**)&descriptor);
+	}
+	xfree((void**)&parameteroutput);
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp	(revision 3905)
@@ -0,0 +1,59 @@
+/*!\file:  CreateSingleNodeToElementConnectivity.cpp
+ * \brief: create connectivity table
+ */ 
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../shared/shared.h"
+#include "../io/io.h"
+#include "../include/include.h"
+#include "./ModelProcessorx.h"
+
+void CreateSingleNodeToElementConnectivity(IoModel* iomodel){
+
+	/*Intermediary*/
+	int i,j;
+	int vertexid;
+	int elementswidth;
+
+	/*output*/
+	int* connectivity=NULL;
+
+	/*Some checks if debugging*/
+	ISSMASSERT(iomodel->numberofvertices);
+	ISSMASSERT(iomodel->numberofelements);
+	ISSMASSERT(iomodel->my_elements);
+	ISSMASSERT(iomodel->elements);
+
+	/*Allocate ouput*/
+	connectivity=(int*)xcalloc(iomodel->numberofvertices,sizeof(int));
+
+	/*Get element width (3 or 6)*/
+	if (strcmp(iomodel->meshtype,"2d")==0){
+		elementswidth=3;
+	}
+	else{
+		elementswidth=6;
+	}
+
+	/*Create connectivity table*/
+	for (i=0;i<iomodel->numberofelements;i++){
+		/*!! in parallel we do not want the vertex to be connected to an element that is not in its partition!!*/
+		if(iomodel->my_elements[i]){
+			for (j=0;j<elementswidth;j++){
+				vertexid=(int)iomodel->elements[elementswidth*i+j];
+				ISSMASSERT(vertexid>0 && vertexid-1<iomodel->numberofvertices);
+				connectivity[vertexid-1]=i+1;
+			}
+		}
+	}
+
+	/*Assign to iomodel*/
+	iomodel->singlenodetoelementconnectivity=connectivity;
+
+	/*Do not free connectivity!!!*/
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 3905)
@@ -0,0 +1,91 @@
+/*
+ * CreateConstraintsDiagnosticHoriz.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+
+void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+
+	int i,j;
+	int count=0;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+	Rgb*    rgb  = NULL;
+	int     node1,node2;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
+	
+	/*Spcs: fetch data: */
+	IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
+	IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+
+	count=1; //matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		if(iomodel->my_vertices[i]){
+
+			if ((int)iomodel->spcvelocity[6*i+0] | (int)iomodel->gridonhutter[i]) 
+				constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcvelocity+6*i+3)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 1 to vx.
+
+			count++;
+			
+			if ((int)iomodel->spcvelocity[6*i+1] | (int)iomodel->gridonhutter[i]) 
+				constraints->AddObject(new Spc(count,i+1,2,*(iomodel->spcvelocity+6*i+4)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 2 to vy
+
+			count++;
+		}
+	}
+	  
+	/*Free data: */
+	xfree((void**)&iomodel->spcvelocity);
+	xfree((void**)&iomodel->gridonhutter);
+
+	
+	/*penalty loads: */
+	IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");
+	
+	count=1; //matlab indexing
+	if(iomodel->numpenalties){
+
+		for (i=0;i<iomodel->numpenalties;i++){
+
+			for(j=1;j<iomodel->numlayers;j++){
+	
+				/*We are pairing nodes along a vertical profile.*/
+				node1=(int)*(iomodel->penalties+iomodel->numlayers*i);
+				node2=(int)*(iomodel->penalties+iomodel->numlayers*i+j);
+
+				constraints->AddObject(new Rgb(count,node1,node2,1));  //add count'th Rgb on dof 1 between node1 and node2
+				
+				count++;
+				
+				constraints->AddObject(new Rgb(count,node1,node2,2));  //add count'th Rgb on dof 1 between node1 and node2
+
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->penalties);
+
+	/*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();
+	
+	cleanup_and_return:
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 3905)
@@ -0,0 +1,234 @@
+/*
+ * CreateElementsNodesAndMaterialsDiagnosticHoriz.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+	
+	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+		IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p");
+		IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q");
+		IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+		IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+		IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		if (iomodel->control_analysis){
+			IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs");
+			IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
+			IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
+		}
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+				
+				if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+
+					/*Create and add tria element to elements dataset: */
+					elements->AddObject(new Tria(i+1,i,iomodel));
+					
+					/*Create and add material property to materials dataset: */
+					materials->AddObject(new Matice(i+1,i,iomodel,3));
+				}
+			}
+		}//for (i=0;i<numberofelements;i++)
+	
+		/*Free data : */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->elements_type);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->drag_coefficient);
+		xfree((void**)&iomodel->drag_p);
+		xfree((void**)&iomodel->drag_q);
+		xfree((void**)&iomodel->rheology_B);
+		xfree((void**)&iomodel->rheology_n);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		if (iomodel->control_analysis){
+			xfree((void**)&iomodel->vx_obs);
+			xfree((void**)&iomodel->vy_obs);
+			xfree((void**)&iomodel->weights);
+		}
+		
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+		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->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+		IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p");
+		IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q");
+		IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+		IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		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->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz");
+		if (iomodel->control_analysis){
+			IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs");
+			IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
+			IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
+		}
+		IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+		IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+
+	
+		for (i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum || *(iomodel->elements_type+2*i+0)==PattynFormulationEnum){ //elements of type 1 are Hutter type Tria. Don't create this elements.
+					/*Create and add penta element to elements dataset: */
+					elements->AddObject(new Penta(i+1,i,iomodel));
+					
+					/*Create and add material property to materials dataset: */
+					materials->AddObject(new Matice(i+1,i,iomodel,6));
+					
+				}
+			}//if(my_elements[i])
+
+		}//for (i=0;i<numberofelements;i++)
+
+		/*Free data: */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->elements_type);
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->drag_coefficient);
+		xfree((void**)&iomodel->drag_p);
+		xfree((void**)&iomodel->drag_q);
+		xfree((void**)&iomodel->rheology_n);
+		xfree((void**)&iomodel->rheology_B);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonbed);
+		xfree((void**)&iomodel->elementonsurface);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->vz);
+		if (iomodel->control_analysis){
+			xfree((void**)&iomodel->vx_obs);
+			xfree((void**)&iomodel->vy_obs);
+			xfree((void**)&iomodel->weights);
+		}
+		xfree((void**)&iomodel->upperelements);
+		xfree((void**)&iomodel->lowerelements);
+
+
+
+	} //if (strcmp(meshtype,"2d")==0)
+	
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials
+	
+	/*Create nodes and vertices: */
+	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->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+	
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+			
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	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: */
+	elements->Presort();
+	nodes->Presort();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp	(revision 3905)
@@ -0,0 +1,97 @@
+/*! \file CreateLoadsDiagnosticHoriz.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+	Icefront*   icefront = NULL;
+	Riftfront*  riftfront= NULL;
+
+	int segment_width;
+	int element;
+	int i;
+	int counter=1;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
+	if (!iomodel->ismacayealpattyn)goto cleanup_and_return;
+	
+	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
+	 * referenced by a certain load must belong to the cluster node): */
+	IoModelFetchData(&iomodel->pressureload,&iomodel->numberofpressureloads,NULL,iomodel_handle,"pressureload");
+	IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+	IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+
+	/*First load data:*/
+	for (i=0;i<iomodel->numberofpressureloads;i++){
+		
+		/*Retrieve element to which this icefront belongs: */
+		if (strcmp(iomodel->meshtype,"2d")==0) segment_width=4; 
+		else segment_width=6;
+		element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill)
+
+		/*Now, if this element is not in the partition, pass: */
+		if(!iomodel->my_elements[element]) continue;
+		
+		/*Do not create ice front if Hutter or Stokes elements*/
+		if ((int)*(iomodel->elements_type+2*element+0)==(HutterFormulationEnum || StokesFormulationEnum)) continue;
+
+		/*Create and  add load: */
+		loads->AddObject(new Icefront(counter,i,iomodel));
+		counter++;
+
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->pressureload);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+
+	/*Create Riffront loads for rifts: */
+	IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+	IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+	IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+	
+	if(iomodel->numrifts){
+		for(i=0;i<iomodel->numrifts;i++){
+				
+			if(iomodel->my_elements[(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)]){
+
+				loads->AddObject(new Riftfront(counter,i,iomodel));
+				counter++;
+			}
+		}
+	}
+				
+	/*free ressources: */
+	xfree((void**)&iomodel->riftinfo);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->surface);
+	xfree((void**)&iomodel->gridoniceshelf);
+
+	cleanup_and_return:
+
+	/*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();
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 3905)
@@ -0,0 +1,57 @@
+/*
+ * CreateConstraintsDiagnosticHutter.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsDiagnosticHutter(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count=0;
+
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Now, is the flag ishutter on? otherwise, do nothing: */
+	if (!iomodel->ishutter) goto cleanup_and_return;
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+
+	count=1; //matlab indexing
+	/*vx and vy are spc'd if we are not on gridonhutter: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i])){
+			if (!(int)iomodel->gridonhutter[i]){
+
+				constraints->AddObject(new Spc(count,i+1,1,0));
+				count++;
+
+				constraints->AddObject(new Spc(count,i+1,2,0));
+				count++;
+			}
+		} //if((my_vertices[i]))
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->gridonhutter);
+
+	/*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();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 3905)
@@ -0,0 +1,170 @@
+/*
+ * CreateElementsNodesAndMaterialsDiagnosticHutter.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Now, is the flag ishutter on? otherwise, do nothing: */
+	if (!iomodel->ishutter)goto cleanup_and_return;
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*Hutter elements can be partitioned using epart, even if
+	 * each hutter elements either lies on a node (in 2d), or a pair of vertically juxtaposed nodes (in 3d): */
+
+	/*Fetch data temporarily needed: */
+	IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	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->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface");
+	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
+	IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+	IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+	IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+
+	/*2d mesh: */
+	if (strcmp(iomodel->meshtype,"2d")==0){
+
+		for (i=0;i<iomodel->numberofvertices;i++){
+
+			if(iomodel->my_vertices[i]){
+
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Sing(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,1));
+
+			}
+
+		} //for (i=0;i<iomodel->numberofvertices;i++)
+	} //if (strcmp(iomodel->meshtype,"2d")==0)
+	else{
+
+		for (i=0;i<iomodel->numberofvertices;i++){
+
+			if(iomodel->my_vertices[i]){
+				if(iomodel->gridonhutter[i]){
+					if(!iomodel->gridonsurface[i]){ 
+
+						/*Create and add penta element to elements dataset: */
+						elements->AddObject(new Beam(i+1,i,iomodel));
+
+						/*Create and add material property to materials dataset: */
+						materials->AddObject(new Matice(i+1,i,iomodel,2));
+
+					}
+				}
+			}
+
+		} //for (i=0;i<iomodel->numberofvertices;i++)
+
+	} //if (strcmp(iomodel->meshtype,"2d")==0)
+	
+	/*Free data: */
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->gridonhutter);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->surface);
+	xfree((void**)&iomodel->bed);
+	xfree((void**)&iomodel->gridonsurface);
+	xfree((void**)&iomodel->uppernodes);
+	xfree((void**)&iomodel->drag_coefficient);
+	xfree((void**)&iomodel->rheology_B);
+	xfree((void**)&iomodel->rheology_n);
+
+	/*Add new constrant material property to materials, at the end: */
+	if (strcmp(iomodel->meshtype,"2d")==0){
+		materials->AddObject(new Matpar(iomodel->numberofvertices+1,iomodel));                          //put it at the end of the materials
+	}
+	else{
+		materials->AddObject(new Matpar(iomodel->numberofvertices2d*(iomodel->numlayers-1)+1,iomodel)); //put it at the end of the materials
+	}
+		
+	/*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->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	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: */
+	elements->Presort();
+	nodes->Presort();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp	(revision 3905)
@@ -0,0 +1,28 @@
+/*! \file CreateLoadsDiagnosticHutter.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Now, is the flag ishutter on? otherwise, do nothing: */
+	if (!iomodel->ishutter)goto cleanup_and_return;
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp	(revision 3905)
@@ -0,0 +1,61 @@
+/*
+ * CreateConstraintsDiagnosticStokes.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsDiagnosticStokes(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+	int count=0;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Now, is the flag ishutter on? otherwise, do nothing: */
+	if (!iomodel->isstokes)goto cleanup_and_return;
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes");
+
+	count=1; //matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+
+			if ((int)!iomodel->gridonstokes[i]){
+
+				/*This grid will see its vx,vy and vz dofs spc'd to pattyn velocities: */
+				constraints->AddObject(new Spc(count,i+1,1,0));
+				count++;
+
+				constraints->AddObject(new Spc(count,i+1,2,0));
+				count++;
+
+				constraints->AddObject(new Spc(count,i+1,3,0));
+				count++;
+			}
+		} //if((my_vertices[i]==1))
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->gridonstokes);
+
+	/*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();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}	
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 3905)
@@ -0,0 +1,175 @@
+/*
+ * CreateElementsNodesAndMaterialsDiagnosticStokes.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Now, do we have Stokes elements?*/
+	if (strcmp(iomodel->meshtype,"2d")==0) goto cleanup_and_return;
+	if (!iomodel->isstokes)                goto cleanup_and_return;
+
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+	IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p");
+	IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q");
+	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->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+	IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+	IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+	IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+	IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz");
+	IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+	IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+	if (iomodel->control_analysis){
+		IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs");
+		IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
+		IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights");
+	}
+
+	for (i=0;i<iomodel->numberofelements;i++){
+
+		if(iomodel->my_elements[i]){
+			 if (*(iomodel->elements_type+2*i+1)==StokesFormulationEnum){
+
+				 /*Create and add penta element to elements dataset: */
+				 elements->AddObject(new Penta(i+1,i,iomodel));
+
+				 /*Create and add material property to materials dataset: */
+				 materials->AddObject(new Matice(i+1,i,iomodel,6));
+			 }
+
+		}
+
+	}//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->drag_coefficient);
+	xfree((void**)&iomodel->drag_p);
+	xfree((void**)&iomodel->drag_q);
+	xfree((void**)&iomodel->elementoniceshelf);
+	xfree((void**)&iomodel->elementonbed);
+	xfree((void**)&iomodel->elementonsurface);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->rheology_n);
+	xfree((void**)&iomodel->rheology_B);
+	xfree((void**)&iomodel->accumulation_rate);
+	xfree((void**)&iomodel->melting_rate);
+	xfree((void**)&iomodel->vx);
+	xfree((void**)&iomodel->vy);
+	xfree((void**)&iomodel->vz);
+	if (iomodel->control_analysis){
+		xfree((void**)&iomodel->vx_obs);
+		xfree((void**)&iomodel->vy_obs);
+		xfree((void**)&iomodel->weights);
+	}
+	xfree((void**)&iomodel->upperelements);
+	xfree((void**)&iomodel->lowerelements);
+
+
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));
+	
+	/*First fetch data: */
+	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");
+	IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes");
+	IoModelFetchData(&iomodel->borderstokes,NULL,NULL,iomodel_handle,"borderstokes");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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->gridonstokes);
+	xfree((void**)&iomodel->borderstokes);
+	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: */
+	elements->Presort();
+	nodes->Presort();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 3905)
@@ -0,0 +1,101 @@
+/*! \file CreateLoadsDiagnosticStokes.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	
+	DataSet*    loads    = NULL;
+	Pengrid*    pengrid  = NULL;
+
+	int segment_width;
+	int element;
+
+	/*pengrid intermediary data: */
+	int numberofpressureloads;
+	int counter=1;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Now, is the flag ishutter on? otherwise, do nothing: */
+	if (!iomodel->isstokes)goto cleanup_and_return;
+
+	/*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids
+	 * referenced by a certain load must belong to the cluster node): */
+	IoModelFetchData(&iomodel->pressureload,&numberofpressureloads,NULL,iomodel_handle,"pressureload");
+	IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+	IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+
+	/*First load data:*/
+	for (i=0;i<numberofpressureloads;i++){
+		
+		segment_width=6;
+
+		element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill)
+
+		/*Now, if this element is not in the partition, pass: */
+		if(!iomodel->my_elements[element]) continue;
+	
+		/*Do not create ice front if it is not a Stokes element*/
+		if ((int)*(iomodel->elements_type+2*element+1)!=StokesFormulationEnum){ 
+			continue;
+		}
+
+		/*Create and  add load: */
+		loads->AddObject(new Icefront(counter,i,iomodel));
+		counter++;
+
+	}
+	/*Free data: */
+	xfree((void**)&iomodel->pressureload);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->thickness);
+	xfree((void**)&iomodel->bed);
+
+	//create penalties for grids on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy
+
+	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridonstokes,NULL,NULL,iomodel_handle,"gridonstokes");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	CreateSingleNodeToElementConnectivity(iomodel);
+	
+	for (i=0;i<iomodel->numberofvertices;i++){
+		
+		if((iomodel->my_vertices[i]==1)){
+
+			if ((iomodel->gridonbed[i]) && (iomodel->gridonicesheet[i]) && (iomodel->gridonstokes[i])){
+				
+				loads->AddObject(new Pengrid(counter,i,iomodel));
+				counter++;
+			}
+		} //if((iomodel->my_vertices[i]==1))
+	}
+
+	xfree((void**)&iomodel->gridonbed);
+	xfree((void**)&iomodel->gridonstokes);
+	xfree((void**)&iomodel->gridonicesheet);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->singlenodetoelementconnectivity);
+
+	/*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/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp	(revision 3905)
@@ -0,0 +1,53 @@
+/*
+ * CreateConstraintsDiagnosticHoriz.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsDiagnosticVert(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	DataSet* constraints = NULL;
+	int count=0;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*return if 2d mesh*/
+	if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity");
+
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	count=1;
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*keep only this partition's nodes:*/
+		if(iomodel->my_vertices[i]){
+
+			if ((int)iomodel->spcvelocity[6*i+2]){
+				constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcvelocity+6*i+5)/iomodel->yts)); //add count'th spc, on node i+1, setting dof 1 to vx.
+				count++;
+
+			}
+		} 
+	}
+
+	/*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();
+
+	cleanup_and_return:
+	/*Free data: */
+	xfree((void**)&iomodel->spcvelocity);
+
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 3905)
@@ -0,0 +1,145 @@
+/*
+ * CreateElementsNodesAndMaterialsDiagnosticVert.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Now, is the flag macayaealpattyn on? otherwise, do nothing: */
+	if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*Create 3d elements: */
+	/*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->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+	IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+	IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+	IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+	for (i=0;i<iomodel->numberofelements;i++){
+
+		if(iomodel->my_elements[i]){
+
+			/*Create and add penta element to elements dataset: */
+			elements->AddObject(new Penta(i+1,i,iomodel));
+
+			/*Create and add material property to materials dataset: */
+			materials->AddObject(new Matice(i+1,i,iomodel,6));
+
+		}
+
+	}//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->elements_type);
+	xfree((void**)&iomodel->melting_rate);
+	xfree((void**)&iomodel->accumulation_rate);
+	xfree((void**)&iomodel->elementonwater);
+	xfree((void**)&iomodel->vx);
+	xfree((void**)&iomodel->vy);
+	xfree((void**)&iomodel->upperelements);
+	xfree((void**)&iomodel->lowerelements);
+
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials
+	
+	/*First fetch data: */
+	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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+	
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp	(revision 3905)
@@ -0,0 +1,23 @@
+/*! \file CreateLoadsDiagnosticVert.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateConstraintsMelting.cpp	(revision 3905)
@@ -0,0 +1,29 @@
+/*
+ * CreateConstraintsMelting.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsMelting(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count;
+	
+	DataSet* constraints = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*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();
+
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp	(revision 3905)
@@ -0,0 +1,149 @@
+/*
+ * CreateElementsNodesAndMaterialsMelting.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../include/include.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+	
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+	IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p");
+	IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q");
+	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->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+	IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+	IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+	IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+	IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure");
+	IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+	IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+	for (i=0;i<iomodel->numberofelements;i++){
+		
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Penta(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,6));
+			}
+
+	}//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->drag_coefficient);
+	xfree((void**)&iomodel->drag_p);
+	xfree((void**)&iomodel->drag_q);
+	xfree((void**)&iomodel->elementoniceshelf);
+	xfree((void**)&iomodel->elementonbed);
+	xfree((void**)&iomodel->elementonsurface);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->rheology_n);
+	xfree((void**)&iomodel->rheology_B);
+	xfree((void**)&iomodel->accumulation_rate);
+	xfree((void**)&iomodel->melting_rate);
+	xfree((void**)&iomodel->elementonwater);
+	xfree((void**)&iomodel->pressure);
+	xfree((void**)&iomodel->upperelements);
+	xfree((void**)&iomodel->lowerelements);
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials
+
+	/*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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+			
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp	(revision 3905)
@@ -0,0 +1,53 @@
+/*! \file CreateLoadsMelting.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsMelting(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+
+	DataSet*    loads    = NULL;
+
+	/*if 2d: return*/
+	if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	//create penalties for grids: no grid can have a temperature over the melting point
+	IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	CreateSingleNodeToElementConnectivity(iomodel);
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+		
+		if((iomodel->my_vertices[i]==1)){
+
+			if (iomodel->gridonbed[i]){ 
+				
+				loads->AddObject(new Pengrid(i+1,i,iomodel));
+			}
+		}
+	}
+	xfree((void**)&iomodel->gridonbed);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->singlenodetoelementconnectivity);
+
+	/*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/modules/ModelProcessorx/ModelProcessorx.h
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h	(revision 3905)
@@ -0,0 +1,94 @@
+/* \file ModelProcessorx.h
+ * \brief  Header file for model processor
+ */
+
+#ifndef _MODEL_PROCESSORX_H_
+#define _MODEL_PROCESSORX_H_
+
+#define RIFTINFOSIZE 11
+
+class DataSet;
+class IoModel;
+class Parameters;
+#include "../io/io.h"
+
+
+/*Creation of fem datasets: general drivers*/
+void  CreateDataSets(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateParameters(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+
+/*Create of fem datasets: specialised drivers: */
+
+/*diagnostic horizontal*/
+void	CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*diagnostic vertical*/
+void	CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*diagnostic hutter*/
+void	CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*diagnostic stokes*/
+void	CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*slope compute*/
+void	CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*control:*/
+void  CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+
+/*thermal:*/
+void	CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*melting:*/
+void	CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*prognostic:*/
+void	CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*prognostic2:*/
+void	CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*balancedthickness:*/
+void	CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+void	CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+
+/*balancedvelocities:*/
+void	CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle);
+void	CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
+void  CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*qmu: */
+void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle);
+
+/*partitioning: */
+void  Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+/*Conectivity*/
+void CreateSingleNodeToElementConnectivity(IoModel* iomodel);
+
+#endif
Index: /issm/trunk/src/c/modules/ModelProcessorx/Partitioning.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Partitioning.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Partitioning.cpp	(revision 3905)
@@ -0,0 +1,274 @@
+/*!\file:  Partitioning.cpp
+ * \brief: partition elements and nodes and vertices
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include <string.h>
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../io/io.h"
+#include "../include/include.h"
+#include "../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx/ModelProcessorx.h"
+
+void  DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle);
+void  ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle);
+
+void  Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){
+	
+	if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum)
+		DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle);
+	else
+		ContinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle);
+}
+
+void  ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){
+
+	/*as many nodes as there are vertices */
+
+	int i;
+
+	extern int my_rank;
+	extern int num_procs;
+
+	/*output: */
+	bool* my_elements=NULL;
+	bool* my_vertices=NULL;
+	bool* my_nodes=NULL;
+	bool* my_bordervertices=NULL;
+
+	/*intermediary: */
+	int* epart=NULL; //element partitioning.
+	int* npart=NULL; //node partitioning.
+	int  elements_width; //number of columns in elements (2d->3, 3d->6)
+	Vec  bordervertices=NULL;
+	double* serial_bordervertices=NULL;
+	int  el1,el2;
+
+	/*Number of vertices per elements, needed to correctly retrieve data: */
+	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->numberofvertices,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofvertices2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs);
+
+	/*Free elements and elements2d: */
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->elements2d);
+
+	#else
+	/*In serial mode, epart is full of 0: all elements belong to cpu 0: */
+	epart=(int*)xcalloc(iomodel->numberofelements,sizeof(int));
+	#endif
+
+	/*Deal with rifts, they have to be included into one partition only, not several: */
+	IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo");
+
+	for(i=0;i<iomodel->numrifts;i++){
+		el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing
+		el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing
+		epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids;
+	}
+
+	/*Free rifts: */
+	xfree((void**)&iomodel->riftinfo); 
+
+	/*Used later on: */
+	my_vertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool));
+	my_elements=(bool*)xcalloc(iomodel->numberofelements,sizeof(bool));
+
+	/*Start figuring out, out of the partition, which elements belong to this cpu: */
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	for (i=0;i<iomodel->numberofelements;i++){
+
+		/*!All elements have been partitioned above, only deal with elements for this cpu: */
+		if(my_rank==epart[i]){ 
+
+			my_elements[i]=1;
+			
+			/*Now that we are here, we can also start building the list of vertices belonging to this cpu 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 vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
+			 will hold which vertices belong to this partition*/
+			my_vertices[(int)*(iomodel->elements+elements_width*i+0)-1]=1;
+			my_vertices[(int)*(iomodel->elements+elements_width*i+1)-1]=1;
+			my_vertices[(int)*(iomodel->elements+elements_width*i+2)-1]=1;
+			
+			if(elements_width==6){
+				my_vertices[(int)*(iomodel->elements+elements_width*i+3)-1]=1;
+				my_vertices[(int)*(iomodel->elements+elements_width*i+4)-1]=1;
+				my_vertices[(int)*(iomodel->elements+elements_width*i+5)-1]=1;
+			}
+		}
+	}//for (i=0;i<numberofelements;i++)
+	/*Free data : */
+	xfree((void**)&iomodel->elements);
+
+	#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:*/
+		bordervertices=NewVec(iomodel->numberofvertices);
+
+		for (i=0;i<iomodel->numberofvertices;i++){
+			if(my_vertices[i])VecSetValue(bordervertices,i,1,ADD_VALUES);
+		}
+		VecAssemblyBegin(bordervertices);
+		VecAssemblyEnd(bordervertices);
+
+		VecToMPISerial(&serial_bordervertices,bordervertices);
+
+		/*now go through serial_bordervertices, and booleanize it: */
+		my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool));
+		for(i=0;i<iomodel->numberofvertices;i++){
+			if(serial_bordervertices[i]>1)my_bordervertices[i]=1;
+		}
+	#else
+		/*No border vertices: */
+		my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool));
+	#endif
+
+	/*Deal with my_nodes: */
+	my_nodes=(bool*)xmalloc(iomodel->numberofvertices*sizeof(bool));
+	memcpy(my_nodes,my_vertices,iomodel->numberofvertices*sizeof(bool));
+
+	/*Free ressources:*/
+	xfree((void**)&npart);
+	xfree((void**)&epart);
+	xfree((void**)&serial_bordervertices);
+	VecFree(&bordervertices);
+
+	/*Assign output pointers:*/
+	*pmy_elements=my_elements;
+	*pmy_vertices=my_vertices;
+	*pmy_nodes=my_nodes;
+	*pmy_bordervertices=my_bordervertices;
+}
+
+
+void  DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){
+
+	/*each element has it own nodes (as many as vertices) + additional nodes from neighbouring elements for each edge. This yields to a very different partition for 
+	 * the nodes and the vertices. The vertices are similar to continuous galerkin, but the nodes partitioning involves edges, which mess up sorting of 
+	 * ids. */
+	
+	int i,j;
+
+	/*output: */
+	bool*   my_elements=NULL;
+	bool*   my_vertices=NULL;
+	bool*   my_nodes=NULL;
+	bool*   my_nodescontinuous=NULL;
+	bool*   my_bordervertices=NULL;
+
+	int     i1,i2;
+	int     cols;
+	double  e1,e2;
+	int     pos;
+
+	/*First: get element and vertices partitioning from Continuous Galerkin: only the nodes are partitioned differently*/
+	ContinuousGalerkinPartitioning(&my_elements,&my_vertices,&my_nodescontinuous,&my_bordervertices,iomodel,iomodel_handle);
+	xfree((void**)&my_nodescontinuous);
+
+	/*Now we must build node partitioning
+	 *  - there are three nodes per element (discontinous)
+	 *  - for each element present of each partition, its three nodes will be in this partition
+	 *  - the edges require the dofs of the 2 nodes of each elements sharing the edge.
+	 *    if the 2 elements sharing the edge are on 2 different cpus, we must duplicate
+	 *    the two nodes that are not on the cpus so that the edge can access the dofs of
+	 *    all its 4 nodes
+	 */
+
+	/*Allocate*/
+	my_nodes=(bool*)xcalloc(3*iomodel->numberofelements,sizeof(int));
+
+	/*First: add all the nodes of all the elements belonging to this cpu*/
+	if (strcmp(iomodel->meshtype,"2d")==0){
+		for (i=0;i<iomodel->numberofelements;i++){
+			if (my_elements[i]){
+				my_nodes[3*i+0]=1;
+				my_nodes[3*i+1]=1;
+				my_nodes[3*i+2]=1;
+			}
+		}
+	}
+	else{
+		ISSMERROR("not implemented yet");
+	}
+
+	/*Second: add all missing nodes*/
+
+	/*Get edges and elements*/
+	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,&cols,iomodel_handle,"edges");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	if (cols!=4) ISSMERROR("field edges should have 4 columns");
+
+	/*!All elements have been partitioned above, only create elements for this CPU: */
+	for (i=0;i<iomodel->numberofedges;i++){
+
+		/*Get left and right elements*/
+		e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
+		e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
+
+		/* 1) If the element e1 is in the current partition
+		 * 2) and if the edge of the element is shared by another element (internal edge)
+		 * 3) and if this element is not in the same partition:
+		 * we must clone the nodes on this partition so that the loads (Numericalflux)
+		 * will have access to their properties (dofs,...)*/
+		if(my_elements[(int)e1] && !isnan(e2) && !my_elements[(int)e2]){ 
+
+			/*1: Get vertices ids*/
+			i1=(int)iomodel->edges[4*i+0];
+			i2=(int)iomodel->edges[4*i+1];
+
+			/*2: Get the column where these ids are located in the index*/
+			pos=UNDEF;
+			for(j=0;j<3;j++){
+				if ((int)iomodel->elements[3*(int)e2+j]==i1) pos=j;
+			}
+
+			/*3: We have the id of the elements and the position of the vertices in the index
+			 * we can now create the corresponding nodes:*/
+			if (pos==0){
+				my_nodes[(int)e2*3+0]=1;
+				my_nodes[(int)e2*3+2]=1;
+			}
+			else if(pos==1){
+				my_nodes[(int)e2*3+1]=1;
+				my_nodes[(int)e2*3+0]=1;
+			}
+			else if (pos==2){
+				my_nodes[(int)e2*3+2]=1;
+				my_nodes[(int)e2*3+1]=1;
+			}
+			else{
+				ISSMERROR("Problem in edges creation");
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->edges);
+
+	/*Assign output pointers:*/
+	*pmy_elements=my_elements;
+	*pmy_vertices=my_vertices;
+	*pmy_nodes=my_nodes;
+	*pmy_bordervertices=my_bordervertices;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp	(revision 3905)
@@ -0,0 +1,52 @@
+/*
+ * CreateConstraintsPrognostic.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../io/io.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsPrognostic(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count=0;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spcthickness,NULL,NULL,iomodel_handle,"spcthickness");
+
+	count=1;//matlab indexing
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+
+			if ((int)iomodel->spcthickness[2*i]){
+
+				constraints->AddObject(new Spc(count,i+1,1,*(iomodel->spcthickness+2*i+1)));
+				count++;
+			}
+		} //if((my_vertices[i]==1))
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->spcthickness);
+	
+	/*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();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp	(revision 3905)
@@ -0,0 +1,187 @@
+/*
+ * CreateElementsNodesAndMaterialsPrognostic.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../io/io.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+	
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+		}//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);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+
+	}
+	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");
+		IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure");
+		IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+		IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+
+		for (i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,6));
+			}
+		}//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);
+		xfree((void**)&iomodel->pressure);
+		xfree((void**)&iomodel->temperature);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+		xfree((void**)&iomodel->upperelements);
+		xfree((void**)&iomodel->lowerelements);
+
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property to materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials
+
+	/*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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp	(revision 3905)
@@ -0,0 +1,29 @@
+/*! \file CreateLoadsPrognostic.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsPrognostic(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/modules/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp	(revision 3905)
@@ -0,0 +1,22 @@
+/*
+ * CreateConstraintsPrognostic2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsPrognostic2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	DataSet* constraints = NULL;
+	
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp	(revision 3905)
@@ -0,0 +1,160 @@
+/*
+ * CreateElementsNodesAndMaterialsPrognostic2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../io/io.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j;
+	int vertex_index;
+	int node_index;
+
+	/*DataSets: */
+	DataSet* elements  = NULL;
+	DataSet* nodes = NULL;
+	DataSet* vertices = NULL;
+	DataSet* materials = NULL;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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");
+		IoModelFetchData(&iomodel->melting_rate,NULL,NULL,iomodel_handle,"melting_rate");
+		IoModelFetchData(&iomodel->accumulation_rate,NULL,NULL,iomodel_handle,"accumulation_rate");
+		IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+		IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+		}//for (i=0;i<numberofelements;i++)
+	
+		/*Free data : */
+		xfree((void**)&iomodel->thickness);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementoniceshelf);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->melting_rate);
+		xfree((void**)&iomodel->accumulation_rate);
+		xfree((void**)&iomodel->vx);
+		xfree((void**)&iomodel->vy);
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+		ISSMERROR("not implemented yet");
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));
+
+	/*Create nodes and vertices: */
+	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->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter");
+	IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet");
+	IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf");
+	if (strcmp(iomodel->meshtype,"3d")==0){
+		IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids");
+		IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids");
+	}
+
+	/*Build Vertices dataset*/
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+		}
+	}
+
+	/*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/
+	for (i=0;i<iomodel->numberofelements;i++){
+		for (j=0;j<3;j++){
+
+			if(iomodel->my_nodes[3*i+j]){ 
+
+				//Get index of the vertex on which the current node is located
+				vertex_index=(int)*(iomodel->elements+3*i+j)-1; //(Matlab to C indexing)
+				ISSMASSERT(vertex_index>=0 && vertex_index<iomodel->numberofvertices);
+
+				//Compute Node index (id-1)
+				node_index=3*i+j;
+
+				/*Add node to nodes dataset: */
+				nodes->AddObject(new Node(node_index+1,vertex_index,node_index,iomodel));
+
+			}
+		}
+	}
+
+	/*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->gridonhutter);
+	xfree((void**)&iomodel->uppernodes);
+	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: */
+	elements->Presort();
+	nodes->Presort();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp	(revision 3905)
@@ -0,0 +1,55 @@
+/*! \file CreateLoadsPrognostic2.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../io/io.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i;
+	int element;
+
+	/*Output*/
+	DataSet* loads=NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Get edges and elements*/
+	IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
+
+	/*First load data:*/
+	for (i=0;i<iomodel->numberofedges;i++){
+
+		/*Get left and right elements*/
+		element=(int)iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
+
+		/*Now, if this element is not in the partition, pass: */
+		if(!iomodel->my_elements[element]) continue;
+
+		/* Add load */
+		loads->AddObject(new Numericalflux(i+1,i,iomodel));
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->edges);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->thickness);
+
+	/*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();
+
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp	(revision 3905)
@@ -0,0 +1,16 @@
+/*!\file: CreateParametersPrognostic.cpp
+ * \brief driver for creating parameters dataset, for prognostic analysis.
+ */ 
+
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void CreateParametersPrognostic2(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Qmu/CreateParametersQmu.cpp	(revision 3905)
@@ -0,0 +1,242 @@
+/*!\file: CreateParametersQmu.cpp
+ * \brief general driver for creating parameters dataset
+ */ 
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../io/io.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
+	
+	int i,j,k;
+	
+	Parameters* parameters = NULL;
+	int      second_count;
+	
+	int*     part=NULL;
+	double*  dpart=NULL;
+
+	char**  responsedescriptors=NULL;
+	char**  variabledescriptors=NULL;
+	char*   descriptor=NULL;
+	char*   tag=NULL;
+	double* dakota_parameter=NULL;
+
+	//qmu files
+	char* qmuinname=NULL;
+	char* qmuerrname=NULL;
+	char* qmuoutname=NULL;
+	extern int my_rank;
+				
+	
+	/*parameters for mass flux: */
+	double* qmu_mass_flux_segments=NULL;
+	double* my_qmu_mass_flux_segments=NULL;
+	int num_qmu_mass_flux_segments=0;
+	int my_num_qmu_mass_flux_segments=0;
+
+	#ifdef _SERIAL_
+		mxArray* pfield=NULL;
+		mxArray* pfield2=NULL;
+	#endif
+	
+	/*recover parameters : */
+	parameters=*pparameters;
+
+	parameters->AddObject(new   BoolParam(QmuAnalysisEnum,iomodel->qmu_analysis));
+	if(iomodel->qmu_analysis){
+
+		//name of qmu input, error and output files
+		qmuinname=(char*)xmalloc((strlen(iomodel->name)+strlen(".qmu.in")+1)*sizeof(char));
+		sprintf(qmuinname,"%s%s",iomodel->name,".qmu.in");
+		parameters->AddObject(new   StringParam(QmuInNameEnum,qmuinname));
+
+		qmuoutname=(char*)xmalloc((strlen(iomodel->name)+strlen(".qmu.out")+1)*sizeof(char));
+		sprintf(qmuoutname,"%s%s",iomodel->name,".qmu.out");
+		parameters->AddObject(new   StringParam(QmuOutNameEnum,qmuoutname));
+
+		qmuerrname=(char*)xmalloc((strlen(iomodel->name)+strlen(".qmu.err")+1)*sizeof(char));
+		sprintf(qmuerrname,"%s%s",iomodel->name,".qmu.err");
+		parameters->AddObject(new   StringParam(QmuErrNameEnum,qmuerrname));
+		
+		parameters->AddObject(new    IntParam(QmuNPartEnum,iomodel->qmu_npart));
+
+		/*Deal with variables for qmu iomodeling: */
+		variabledescriptors=(char**)xmalloc(iomodel->numberofvariables*sizeof(char*));
+
+
+		/*Fetch descriptors: logic varies if we are running parallel or serial. In parallel, qmumarshall 
+		 * took care of marshalling all the variable descriptors, so it's easy. In serial mode, 
+		 * the variables are in md.variables(md.ivar), as a strucuture: */
+
+		#ifdef _SERIAL_
+		pfield=mxGetField(iomodel_handle,0,"variabledescriptors");
+		for(i=0;i<iomodel->numberofvariables;i++){
+			pfield2=mxGetCell(pfield,i);
+			FetchData(&descriptor,pfield2);
+			variabledescriptors[i]=descriptor;
+		}
+		#else
+		tag=(char*)xmalloc((strlen("variabledescriptori")+1)*sizeof(char));
+		for(i=0;i<iomodel->numberofvariables;i++){
+			sprintf(tag,"%s%i","variabledescriptor",i);
+			IoModelFetchData(&descriptor,iomodel_handle,tag);
+			variabledescriptors[i]=descriptor;
+		}
+		#endif
+
+		/*Ok, we have all the variable descriptors. Build a parameter with it: */
+		parameters->AddObject(new StringArrayParam(VariableDescriptorsEnum,variabledescriptors,iomodel->numberofvariables));
+
+		/*Deal with responses and partition for qmu iomodeling: */
+		responsedescriptors=(char**)xmalloc(iomodel->numberofresponses*sizeof(char*));
+
+		/*Fetch descriptors: */
+		#ifdef _SERIAL_
+		pfield=mxGetField(iomodel_handle,0,"responsedescriptors");
+		for(i=0;i<iomodel->numberofresponses;i++){
+			pfield2=mxGetCell(pfield,i);
+			FetchData(&descriptor,pfield2);
+			responsedescriptors[i]=descriptor;
+		}
+		#else
+		xfree((void**)&tag);
+		tag=(char*)xmalloc((strlen("responsedescriptori")+1)*sizeof(char));
+
+		for(i=0;i<iomodel->numberofresponses;i++){
+			sprintf(tag,"%s%i","responsedescriptor",i);
+			IoModelFetchData(&descriptor,iomodel_handle,tag);
+			responsedescriptors[i]=descriptor;
+		}
+		#endif
+
+		/*Ok, we have all the response descriptors. Build a parameter with it: */
+		parameters->AddObject(new StringArrayParam(ResponseDescriptorsEnum,responsedescriptors,iomodel->numberofresponses));
+
+		/*partition vertices in iomodel->qmu_npart parts, unless a partition is already present: */
+		IoModelFetchData(&dpart,NULL,NULL,iomodel_handle,"part");
+
+		if(!dpart){
+
+			/*Partition elements and vertices and nodes: */
+			Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+			dpart=(double*)xmalloc(iomodel->numberofvertices*sizeof(double));
+			for(i=0;i<iomodel->numberofvertices;i++)dpart[i]=iomodel->my_vertices[i];
+		}
+		parameters->AddObject(new DoubleVecParam(QmuPartEnum,dpart,iomodel->numberofvertices));
+
+		/*Ok, now if any of the variables input from Dakota are distributed, we are going to need the parameters: */
+		for(i=0;i<iomodel->numberofvariables;i++){
+
+			descriptor=variabledescriptors[i];
+
+			if ((strcmp(descriptor,"thickness")==0) ||
+				(strcmp(descriptor,"drag")     ==0)
+				){
+
+				IoModelFetchData(&dakota_parameter,NULL,NULL,iomodel_handle,descriptor);
+				parameters->AddObject(new DoubleVecParam(DakotaParameterEnum,dakota_parameter,iomodel->numberofvertices));
+				xfree((void**)&dakota_parameter);
+
+			}
+		}
+
+		/*Deal with data needed for some responses: */
+		for(i=0;i<iomodel->numberofresponses;i++){
+			char* descriptor=responsedescriptors[i];
+			if (strcmp(descriptor,"mass_flux")==0){
+
+				/*We need the qmu_mass_flux_segments to be able to compute the mass flux: */
+				IoModelFetchData(&qmu_mass_flux_segments,&num_qmu_mass_flux_segments,NULL,iomodel_handle,"qmu_mass_flux_segments");
+
+				#ifdef _PARALLEL_
+
+					/*Only if partitioning exist do we care about the segments: */
+					if(iomodel->my_elements){
+
+						/*Use the element partitioning vector from the iomodel to down select qmu_mass_flux_segments to only segments that are relevant 
+						 * to this cpu: */
+						my_num_qmu_mass_flux_segments=0;
+						for(j=0;j<num_qmu_mass_flux_segments;j++){
+							if (  iomodel->my_elements[(int)(*(qmu_mass_flux_segments+5*j+4))-1])my_num_qmu_mass_flux_segments++;
+						}
+
+					
+						if(my_num_qmu_mass_flux_segments){
+							my_qmu_mass_flux_segments=(double*)xcalloc(5*my_num_qmu_mass_flux_segments,sizeof(double));
+							second_count=0;
+							for(j=0;j<num_qmu_mass_flux_segments;j++){
+								if (iomodel->my_elements[(int)*(qmu_mass_flux_segments+5*j+4)-1]){
+									for(k=0;k<5;k++)*(my_qmu_mass_flux_segments+5*second_count+k)=*(qmu_mass_flux_segments+5*j+k);
+									second_count++;
+								}
+							}
+						}
+
+						parameters->AddObject(new DoubleMatParam(QmuMassFluxSegmentsEnum,my_qmu_mass_flux_segments,my_num_qmu_mass_flux_segments,5));
+
+					}
+					
+				#else
+
+					parameters->AddObject(new DoubleMatParam(QmuMassFluxSegmentsEnum,qmu_mass_flux_segments,num_qmu_mass_flux_segments,5));
+
+				#endif
+
+				xfree((void**)&qmu_mass_flux_segments);
+				xfree((void**)&my_qmu_mass_flux_segments);
+			}
+			if (strcmp(descriptor,"misfit")==0){
+
+				/*We need the observed velocity: */
+				IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs");
+				IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs");
+				IoModelFetchData(&iomodel->fit,NULL,NULL,iomodel_handle,"fit");
+
+				if(iomodel->vx_obs)for(i=0;i<iomodel->numberofvertices;i++)iomodel->vx_obs[2*i+0]=iomodel->vx_obs[i]/iomodel->yts;
+				if(iomodel->vy_obs)for(i=0;i<iomodel->numberofvertices;i++)iomodel->vy_obs[2*i+1]=iomodel->vy_obs[i]/iomodel->yts;
+
+				parameters->AddObject(new DoubleVecParam(VxObsEnum,iomodel->vx_obs,iomodel->numberofvertices));
+				parameters->AddObject(new DoubleVecParam(VyObsEnum,iomodel->vy_obs,iomodel->numberofvertices));
+				parameters->AddObject(new DoubleVecParam(FitEnum,iomodel->fit,iomodel->nsteps));
+
+				xfree((void**)&iomodel->vx_obs);
+				xfree((void**)&iomodel->vy_obs);
+				xfree((void**)&iomodel->fit);
+
+			}
+		}
+
+		/*Free data: */
+		xfree((void**)&tag);
+		for(i=0;i<iomodel->numberofresponses;i++){
+			char* descriptor=responsedescriptors[i];
+			xfree((void**)&descriptor);
+		}
+		xfree((void**)&responsedescriptors);
+
+		for(i=0;i<iomodel->numberofvariables;i++){
+			char* descriptor=variabledescriptors[i];
+			xfree((void**)&descriptor);
+		}
+		xfree((void**)&variabledescriptors);
+
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->elements2d);
+		xfree((void**)&part);
+		xfree((void**)&dpart);
+		xfree((void**)&qmuinname);
+		xfree((void**)&qmuerrname);
+		xfree((void**)&qmuoutname);
+	}
+
+	/*Assign output pointer: */
+	*pparameters=parameters;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp	(revision 3905)
@@ -0,0 +1,22 @@
+/*
+ * CreateConstraintsSlopeCompute.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsSlopeCompute(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet* constraints = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 3905)
@@ -0,0 +1,158 @@
+/*
+ * CreateElementsNodesAndMaterialsSlopeCompute.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+		/*output: int* epart, int* my_grids, double* my_bordergrids*/
+
+	int i,j,k,n;
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+	
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*2d mesh: */
+	if (strcmp(iomodel->meshtype,"2d")==0){
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		
+		for (i=0;i<iomodel->numberofelements;i++){
+
+			if(iomodel->my_elements[i]){
+				
+				/*Create and add tria element to elements dataset: */
+				elements->AddObject(new Tria(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,3));
+			}
+
+		}//for (i=0;i<numberofelements;i++)
+
+	
+		/*Free data : */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementonwater);
+
+	}
+	else{ //	if (strcmp(meshtype,"2d")==0)
+
+		/*Fetch data needed: */
+		IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+		IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");
+		IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");
+		IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");
+		IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+		IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+		IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+		for (i=0;i<iomodel->numberofelements;i++){
+			if(iomodel->my_elements[i]){
+				/*Create and add penta element to elements dataset: */
+				elements->AddObject(new Penta(i+1,i,iomodel));
+
+				/*Create and add material property to materials dataset: */
+				materials->AddObject(new Matice(i+1,i,iomodel,6));
+
+			}//if(my_elements[i])
+		}//for (i=0;i<numberofelements;i++)
+
+		/*Free data: */
+		xfree((void**)&iomodel->elements);
+		xfree((void**)&iomodel->surface);
+		xfree((void**)&iomodel->bed);
+		xfree((void**)&iomodel->elementonbed);
+		xfree((void**)&iomodel->elementonwater);
+		xfree((void**)&iomodel->upperelements);
+		xfree((void**)&iomodel->lowerelements);
+
+	} //if (strcmp(meshtype,"2d")==0)
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials
+	
+	/*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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+			
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+	
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp	(revision 3905)
@@ -0,0 +1,30 @@
+/*! \file CreateLoadsSlopeCompute.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	DataSet*    loads    = NULL;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	/*Now, is the flag isstokes on? otherwise, do nothing: */
+	if (!iomodel->isstokes)goto cleanup_and_return;
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*ploads=loads;
+
+}
+
+
Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp	(revision 3905)
@@ -0,0 +1,58 @@
+/*
+ * CreateConstraintsThermal.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../io/io.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../ModelProcessorx.h"
+
+void	CreateConstraintsThermal(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+	int count=0;
+	
+	DataSet* constraints = NULL;
+	Spc*    spc  = NULL;
+
+	/*Create constraints: */
+	constraints = new DataSet(ConstraintsEnum);
+
+	/*return if 2d mesh*/
+	if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
+
+	/*Fetch data: */
+	IoModelFetchData(&iomodel->spctemperature,NULL,NULL,iomodel_handle,"spctemperature");
+
+	count=1;
+
+	/*Create spcs from x,y,z, as well as the spc values on those spcs: */
+	for (i=0;i<iomodel->numberofvertices;i++){
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+
+			if ((int)iomodel->spctemperature[2*i]){
+
+				constraints->AddObject(new Spc(count,i+1,1,iomodel->spctemperature[2*i+1]));
+				count++;
+
+			}
+
+		} //if((my_grids[i]==1))
+	}
+
+	/*Free data: */
+	xfree((void**)&iomodel->spctemperature);
+
+	/*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();
+
+	cleanup_and_return:
+	
+	/*Assign output pointer: */
+	*pconstraints=constraints;
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp	(revision 3905)
@@ -0,0 +1,152 @@
+/*
+ * CreateElementsNodesAndMaterialsThermal.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../include/include.h"
+#include "../../shared/shared.h"
+#include "../../MeshPartitionx/MeshPartitionx.h"
+#include "../ModelProcessorx.h"
+
+void	CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	/*Intermediary*/
+	int i,j,k,n;
+
+	/*DataSets: */
+	DataSet*    elements  = NULL;
+	DataSet*    nodes = NULL;
+	DataSet*    vertices = NULL;
+	DataSet*    materials = NULL;
+	
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum);
+	nodes     = new DataSet(NodesEnum);
+	vertices  = new DataSet(VerticesEnum);
+	materials = new DataSet(MaterialsEnum);
+
+	/*Partition elements and vertices and nodes: */
+	Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle);
+
+	/*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->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient");
+	IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p");
+	IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q");
+	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->elements_type,NULL,NULL,iomodel_handle,"elements_type");
+	IoModelFetchData(&iomodel->geothermalflux,NULL,NULL,iomodel_handle,"geothermalflux");
+	IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B");
+	IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n");
+	IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");
+	IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure");
+	IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature");
+	IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx");
+	IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy");
+	IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements");
+	IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements");
+	
+	for (i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+
+			/*Create and add tria element to elements dataset: */
+			elements->AddObject(new Penta(i+1,i,iomodel));
+
+			/*Create and add material property to materials dataset: */
+			materials->AddObject(new Matice(i+1,i,iomodel,6));
+		}
+	}//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->drag_coefficient);
+	xfree((void**)&iomodel->drag_p);
+	xfree((void**)&iomodel->drag_q);
+	xfree((void**)&iomodel->elementoniceshelf);
+	xfree((void**)&iomodel->elementonbed);
+	xfree((void**)&iomodel->elementonsurface);
+	xfree((void**)&iomodel->elements_type);
+	xfree((void**)&iomodel->geothermalflux);
+	xfree((void**)&iomodel->rheology_n);
+	xfree((void**)&iomodel->rheology_B);
+	xfree((void**)&iomodel->elementonwater);
+	xfree((void**)&iomodel->pressure);
+	xfree((void**)&iomodel->temperature);
+	xfree((void**)&iomodel->vx);
+	xfree((void**)&iomodel->vy);
+	xfree((void**)&iomodel->upperelements);
+	xfree((void**)&iomodel->lowerelements);
+
+	/*Add new constrant material property tgo materials, at the end: */
+	materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials
+	
+	/*Create nodes and vertices: */
+	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");
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+
+		/*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+		if(iomodel->my_vertices[i]){
+			
+			/*Add vertex to vertices dataset: */
+			vertices->AddObject(new Vertex(i+1,i,iomodel));
+
+			/*Add node to nodes dataset: */
+			nodes->AddObject(new Node(i+1,i,iomodel));
+
+		}
+	}
+
+	/*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);
+	
+	/*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();
+	vertices->Presort();
+	materials->Presort();
+
+	cleanup_and_return:
+
+	/*Assign output pointer: */
+	*pelements=elements;
+	*pnodes=nodes;
+	*pvertices=vertices;
+	*pmaterials=materials;
+
+}
Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp
===================================================================
--- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp	(revision 3905)
@@ -0,0 +1,56 @@
+/*! \file CreateLoadsThermal.c:
+ */
+
+#include "../../DataSet/DataSet.h"
+#include "../../toolkits/toolkits.h"
+#include "../../io/io.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+#include "../../objects/objects.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../ModelProcessorx.h"
+
+void	CreateLoadsThermal(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
+
+	int i;
+
+	DataSet*    loads    = NULL;
+	Pengrid*    pengrid  = NULL;
+
+	/*return if 2d mesh*/
+	if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return;
+
+	/*Create loads: */
+	loads   = new DataSet(LoadsEnum);
+
+	//create penalties for grids: no grid can have a temperature over the melting point
+	IoModelFetchData(&iomodel->spctemperature,NULL,NULL,iomodel_handle,"spctemperature");
+	IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
+	CreateSingleNodeToElementConnectivity(iomodel);
+
+	for (i=0;i<iomodel->numberofvertices;i++){
+	
+		/*keep only this partition's nodes:*/
+		if((iomodel->my_vertices[i]==1)){
+
+			if (!iomodel->spctemperature[2*i]){ //No penalty applied on spc grids!
+
+				loads->AddObject(new Pengrid(i+1,i,iomodel));
+						
+			}
+		}
+	}
+	xfree((void**)&iomodel->spctemperature);
+	xfree((void**)&iomodel->elements);
+	xfree((void**)&iomodel->singlenodetoelementconnectivity);
+
+	/*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/modules/MpcNodesx/MpcNodesx.cpp
===================================================================
--- /issm/trunk/src/c/modules/MpcNodesx/MpcNodesx.cpp	(revision 3905)
+++ /issm/trunk/src/c/modules/MpcNodesx/MpcNodesx.cpp	(revision 3905)
@@ -0,0 +1,55 @@
+/*!\file MpcNodesx
+ * \brief: establish matrix for wringing out rigid body modes. (multi-point constraints)
+ */
+
+#include "./MpcNodesx.h"
+
+#include "../shared/shared.h"
+#include "../include/include.h"
+#include "../toolkits/toolkits.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+
+int MpcNodesx( Mat* pRmg, DataSet* nodes,DataSet* constraints){
+
+	int i;
+	int numberofdofs;
+	int numberofequations;
+
+	/*output: */
+	Mat Rmg=NULL;
+
+	
+	/*fake values: */
+	int Rmg_connectivity;
+	int Rmg_numberofdofspernode;
+
+	/*output: */
+	Vec yg=NULL;
+
+	/*First, recover number of dofs from nodes: */
+	numberofdofs=nodes->NumberOfDofs();
+
+	/*Get number of equations from number of constraints: */
+	numberofequations=constraints->NumberOfRgbs();
+
+	if(numberofequations){
+
+		/*Allocate Rmg: */
+		Rmg_connectivity=2; Rmg_numberofdofspernode=1; 
+		
+		Rmg=NewMat(numberofequations,numberofdofs,NULL,&Rmg_connectivity,&Rmg_numberofdofspernode);
+
+		/*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
+		constraints->SetupMpcs(Rmg,nodes);
+
+		/*Assemble matrix: */
+		MatAssemblyBegin(Rmg,MAT_FINAL_ASSEMBLY);
+		MatAssemblyEnd(Rmg,MAT_FINAL_ASSEMBLY);
+	}
+
+	/*Assign output pointers: */
+	*pRmg=Rmg;
+	
+	return 1;
+	
+}
Index: /issm/trunk/src/c/modules/MpcNodesx/MpcNodesx.h
===================================================================
--- /issm/trunk/src/c/modules/MpcNodesx/MpcNodesx.h	(revision 3905)
+++ /issm/trunk/src/c/modules/MpcNodesx/MpcNodesx.h	(revision 3905)
@@ -0,0 +1,15 @@
+/*!\file:  MpcNodesx.h
+ * \brief header file for rigid body motion
+ */ 
+
+#ifndef _MPCNODESX_H
+#define _MPCNODESX_H
+
+#include "../DataSet/DataSet.h"
+#include "../toolkits/toolkits.h"
+
+/* local prototypes: */
+int		MpcNodesx( Mat* pRmg, DataSet* nodes,DataSet* constraints);
+
+#endif  /* _MPCNODESX_H */
+