Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 305)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 306)
@@ -142,13 +142,4 @@
 	int el1,el2;
 	 
-	/*Penalty partitioning: */
-	int num_grids3d_collapsed;
-	double* double_penalties_grids3d_collapsed=NULL;
-	double* double_penalties_grids3d_noncollapsed=NULL;
-	int     grid_ids[6];
-	int     num_grid_ids;
-	int     grid_id;
-
-
 	/*First create the elements, nodes and material properties: */
 	elements  = new DataSet(ElementsEnum());
@@ -336,5 +327,5 @@
 	
 		/*Free data : */
-		/*xfree((void**)&model->elements);
+		xfree((void**)&model->elements);
 		xfree((void**)&model->thickness);
 		xfree((void**)&model->surface);
@@ -345,5 +336,5 @@
 		xfree((void**)&model->elementoniceshelf);
 		xfree((void**)&model->B);
-		xfree((void**)&model->n);*/
+		xfree((void**)&model->n);
 	}
 	else{ //	if (strcmp(meshtype,"2d")==0)
@@ -668,4 +659,5 @@
 	#ifdef _PARALLEL_
 	xfree((void**)&all_numgrids);
+	xfree((void**)&npart);
 	VecFree(&gridborder);
 	#endif
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 305)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp	(revision 306)
@@ -16,7 +16,17 @@
 void	CreateConstraintsDiagnosticHutter(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){
 
+	int i;
+	int count;
 
 	DataSet* constraints = NULL;
 
+	Spc*    spc  = NULL;
+
+	/*spc intermediary data: */
+	int spc_sid;
+	int spc_node;
+	int spc_dof;
+	double spc_value;
+	
 	/*Create constraints: */
 	constraints = new DataSet(ConstraintsEnum());
@@ -24,6 +34,90 @@
 	/*Now, is the flag ishutter on? otherwise, do nothing: */
 	if (!model->ishutter)goto cleanup_and_return;
+
+	count=0;
 	
-	cleanup_and_return:	
+	/*Fetch data: */
+	ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
+
+	/*vx and vy are spc'd if we are not on gridonhutter: */
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((model->npart[i]==1)){
+	#endif
+
+		if ((int)model->gridonhutter[i]){
+	
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=1; //we enforce first x translation degree of freedom
+			spc_value=0;
+
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+
+			spc_sid=count;
+			spc_node=i+1;
+			spc_dof=2; //we enforce first y translation degree of freedom
+			spc_value=0;
+			
+			spc = new Spc(spc_sid,spc_node,spc_dof,spc_value);
+			constraints->AddObject(spc);
+			count++;
+		}
+
+	#ifdef _PARALLEL_
+	} //if((npart[i]==1))
+	#endif
+	}
+
+	/*Free data: */
+	xfree((void**)&model->gridonhutter);
+
+	//deal with mpcs for 2d-3d mesh transitions
+	
+	/*Fetch data: */
+	ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
+
+	if (strcmp(model->meshtype,"3d")==0){
+	
+		if (model->numpenalties){
+			for (i=0;i<model->numpenalties;i++){
+				for (j=0;j<model->numlayers-1;j++){
+
+					//constrain first dof
+					rgb_id=count;
+					rgb_dof=1;
+					rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j);
+					rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j+1);
+					
+					rgb = new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
+					constraints->AddObject(rgb);
+					count++;
+
+					//constrain second dof
+					rgb_id=count;
+					rgb_dof=2;
+					rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j);
+					rgb_nodeid1=*(model->penalties+(model->numlayers-1)*i+j+1);
+					
+					rgb = new Rgb(rgb_id,rgb_nodeid1,rgb_nodeid2,rgb_dof);
+					constraints->AddObject(rgb);
+					count++;
+	
+				}
+			}
+		}
+	}
+
+	/*Free data: */
+	xfree((void**)&model->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: */
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 305)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp	(revision 306)
@@ -28,7 +28,325 @@
 	materials = new DataSet(MaterialsEnum());
 
+	/*Objects: */
+	Node*       node   = NULL;
+	Matice*     matice = NULL;
+	Matpar*     matpar = NULL;
+	Beam*       beam   = NULL;
+	Sing*       sing   = NULL;
+
+	int         analysis_type;
+	
+	/*output: */
+	int* epart=NULL; //element partitioning.
+	int* npart=NULL; //node partitioning.
+	int* my_grids=NULL;
+	double* my_bordergrids=NULL;
+
+	/*intermediary: */
+	int elements_width; //size of elements
+	double B_avg;
+			
+	/*matice constructor input: */
+	int    matice_mid;
+	double matice_B;
+	double matice_n;
+
+	/*sing constructor input: */
+	int   sing_id;
+	int   sing_mid;
+	int   sing_mparid;
+	int   sing_g;
+	double sing_h,sing_k;
+
+	/*beam constructor input: */
+	int   beam_id;
+	int   beam_mid;
+	int   beam_mparid;
+	int   beam_g[2];
+	double beam_h[2];
+	double beam_s[2];
+	double beam_b[2];
+	double beam_k[2];
+					
+	/*matpar constructor input: */
+	int	matpar_mid;
+	double matpar_rho_ice; 
+	double matpar_rho_water;
+	double matpar_heatcapacity;
+	double matpar_thermalconductivity;
+	double matpar_latentheat;
+	double matpar_beta;
+	double matpar_meltingpoint;
+	double matpar_mixed_layer_capacity;
+	double matpar_thermal_exchange_velocity;
+	double matpar_g;
+
+	/* node constructor input: */
+	int node_id;
+	int node_partitionborder=0;
+	double node_x[3];
+	int node_onbed;
+	int node_onsurface;
+	int node_upper_node_id;
+	int node_numdofs;
+
+	/*First create the elements, nodes and material properties: */
+	elements  = new DataSet(ElementsEnum());
+	nodes     = new DataSet(NodesEnum());
+	materials = new DataSet(MaterialsEnum());
+
 	/*Now, is the flag ishutter on? otherwise, do nothing: */
 	if (!model->ishutter)goto cleanup_and_return;
 
+	/*Get analysis_type: */
+	analysis_type=AnalysisTypeAsEnum(model->analysis_type);
+
+	/*Width of elements: */
+	if(strcmp(model->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(model->meshtype,"2d")==0){
+		/*load elements: */
+		ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
+	}
+	else{
+		/*load elements2d: */
+		ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
+	}
+
+	MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs);
+
+	/*Free elements and elements2d: */
+	xfree((void**)&model->elements);
+	xfree((void**)&model->elements2d);
+
+			
+	#endif
+
+	/*Hutter elements can be partitioned using npart, instead of epart for more classic tria and penta elements. The reason is 
+	 * that each hutter elements either lies on a node (in 2d), or a pair of vertically juxtaposed nodes (in 3d): */
+
+	/*Fetch data temporarily needed: */
+	ModelFetchData((void**)&model->gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat");
+	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
+	ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+	ModelFetchData((void**)&model->uppergrids,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
+	ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat");
+	ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat");
+	ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat");
+
+	/*2d mesh: */
+	if (strcmp(model->meshtype,"2d")==0){
+
+		for (i=0;i<model->numberofnodes;i++){
+		#ifdef _PARALLEL_
+		/*keep only this partition's nodes:*/
+		if(my_rank==npart[i]){
+		#endif
+
+			if(model->gridonhutter[i]){
+				
+				/*Deal with sing element: */
+				sing_id=i+1; 
+				sing_mid=i+1; //refers to the corresponding material property card
+				sing_mparid=model->numberofnodes+1;//refers to the corresponding matpar property card
+				sing_g=i+1;
+				sing_h=model->thickness[i];
+				sing_k=model->drag[i];
+
+				/*Create sing element using its constructor:*/
+				sing=new Sing(sing_id, sing_mid, sing_mparid, sing_g, sing_h, sing_k);
+
+				/*Add tria element to elements dataset: */
+				elements->AddObject(sing);
+
+				/*Deal with material property card: */
+				matice_mid=i+1; //same as the material id from the geom2 elements.
+				matice_B=model->B[i];	
+				matice_n=(double)model->n[i];
+			
+				/*Create matice using its constructor:*/
+				matice= new Matice(matice_mid,matice_B,matice_n);
+	
+				/*Add matice element to materials dataset: */
+				materials->AddObject(matice);
+			}
+
+		#ifdef _PARALLEL_
+		} //if(my_rank==npart[i])
+		#endif
+
+		} //for (i=0;i<model->numberofnodes;i++)
+	} //if (strcmp(model->meshtype,"2d")==0)
+	else{
+
+		for (i=0;i<model->numberofnodes;i++){
+		#ifdef _PARALLEL_
+		/*keep only this partition's nodes:*/
+		if(my_rank==npart[i]){
+		#endif
+
+			if(model->gridonhutter[i]){
+
+				if(!model->gridonsurface[i]){ 
+					
+					/*Deal with sing element: */
+					beam_id=i+1; 
+					beam_mid=i+1; //refers to the corresponding material property card
+					beam_mparid=model->numberofnodes+1;//refers to the corresponding matpar property card
+					beam_g[0]=i+1;
+					beam_g[1]=model->uppergrids[i]; //grid that lays right on top
+					beam_h[0]=model->thickness[i];
+					beam_h[1]=model->thickness[model->uppergrids[i]-1];
+					beam_s[0]=model->surface[i];
+					beam_s[1]=model->surface[model->uppergrids[i]-1];
+					beam_b[0]=model->bed[i];
+					beam_b[1]=model->bed[model->uppergrids[i]-1];
+					beam_k[0]=model->drag[i];
+					beam_k[1]=model->drag[model->uppergrids[i]-1];
+
+					/*Create beam element ubeam its constructor:*/
+					beam=new Beam(beam_id, beam_mid, beam_mparid, beam_g, beam_h, beam_s,beam_b,beam_k);
+
+					/*Add tria element to elements dataset: */
+					elements->AddObject(beam);
+
+					/*Deal with material property card: */
+					matice_mid=i+1; //same as the material id from the geom2 elements.
+					matice_B=model->B[i];	
+					matice_n=(double)model->n[i];
+				
+					/*Create matice ubeam its constructor:*/
+					matice= new Matice(matice_mid,matice_B,matice_n);
+		
+					/*Add matice element to materials dataset: */
+					materials->AddObject(matice);
+				}
+			}
+
+		#ifdef _PARALLEL_
+		} //if(my_rank==npart[i])
+		#endif
+
+		} //for (i=0;i<model->numberofnodes;i++)
+
+	} //if (strcmp(model->meshtype,"2d")==0)
+	
+
+	/*Free data: */
+	xfree((void**)&model->gridonhutter);
+	xfree((void**)&model->thickness);
+	xfree((void**)&model->surface);
+	xfree((void**)&model->gridonsurface);
+	xfree((void**)&model->uppergrids);
+	xfree((void**)&model->drag);
+	xfree((void**)&model->B);
+	xfree((void**)&model->n);
+	
+
+	/*Add one constant material property to materials: */
+	matpar_mid=model->numberofnodes+1; //put it at the end of the materials
+	matpar_g=model->g; 
+	matpar_rho_ice=model->rho_ice; 
+	matpar_rho_water=model->rho_water; 
+	matpar_thermalconductivity=model->thermalconductivity; 
+	matpar_heatcapacity=model->heatcapacity; 
+	matpar_latentheat=model->latentheat; 
+	matpar_beta=model->beta; 
+	matpar_meltingpoint=model->meltingpoint; 
+	matpar_mixed_layer_capacity=model->mixed_layer_capacity; 
+	matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; 
+
+	
+	/*Create nodes from x,y,z, as well as the spc values on those grids: */
+		
+	/*First fetch data: */
+	if (strcmp(model->meshtype,"3d")==0){
+		ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat");
+		ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat");
+	}
+	ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat");
+	ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat");
+	ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat");
+
+	
+	/*Get number of dofs per node: */
+	DistributeNumDofs(&node_numdofs,analysis_type);
+
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if(npart[i]==my_rank){
+	#endif
+
+		node_id=i+1; //matlab indexing
+			
+		node_partitionborder=0;
+
+		node_x[0]=model->x[i];
+		node_x[1]=model->y[i];
+		node_x[2]=model->z[i];
+
+		
+		node_onbed=(int)model->gridonbed[i];
+		node_onsurface=(int)model->gridonsurface[i];	
+		if (strcmp(model->meshtype,"3d")==0){
+			if (isnan(model->uppernodes[i])){
+				node_upper_node_id=node_id;  //nodes on surface do not have upper nodes, only themselves.
+			}
+			else{
+				node_upper_node_id=(int)model->uppernodes[i];
+			}
+		}
+		else{
+			/*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/
+			node_upper_node_id=node_id;
+		}
+
+		/*Create node using its constructor: */
+		node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id);
+
+		/*set single point constraints.: */
+		if (strcmp(model->meshtype,"3d")==0){
+			/*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */
+			if (model->deadgrids[i]){
+				for(k=1;k<=node_numdofs;k++){
+					node->FreezeDof(k);
+				}
+			}
+		}
+		/*Add node to nodes dataset: */
+		nodes->AddObject(node);
+
+	#ifdef _PARALLEL_
+	} //if(npart[i]==my_rank)
+	#endif
+	}
+
+	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
+	 * datasets, it will not be redone: */
+	elements->Presort();
+	nodes->Presort();
+	materials->Presort();
+
+	/*Clean fetched data: */
+	xfree((void**)&model->deadgrids);
+	xfree((void**)&model->x);
+	xfree((void**)&model->y);
+	xfree((void**)&model->z);
+	xfree((void**)&model->gridonbed);
+	xfree((void**)&model->gridonsurface);
+	xfree((void**)&model->uppernodes);
+		
 	cleanup_and_return:
 
@@ -37,3 +355,7 @@
 	*pnodes=nodes;
 	*pmaterials=materials;
+
+	/*Keep partitioning information into model*/
+	xfree((void**)&epart);
+	model->npart=npart;
 }
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 305)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 306)
@@ -432,4 +432,5 @@
 	#ifdef _PARALLEL_
 	xfree((void**)&all_numgrids);
+	xfree((void**)&npart);
 	VecFree(&gridborder);
 	#endif
Index: /issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 305)
+++ /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 306)
@@ -147,4 +147,5 @@
 	/*exterior data: */
 	int* epart; /*!element partitioning.*/
+	int* npart; /*!node partitioning.*/
 	int* my_grids; /*! grids that belong to this cpu*/
 	double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/