Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 386)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 387)
@@ -52,4 +52,5 @@
 int IcefrontEnum(void){                 return         201; }
 int PenpairEnum(void){                  return         202; }
+int PengridEnum(void){                  return         203; }
 
 /*Materials: */
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 386)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 387)
@@ -18,4 +18,5 @@
 int SpcEnum(void);
 int PenpairEnum(void);
+int PengridEnum(void);
 int IcefrontEnum(void);
 int ParamEnum(void);
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 386)
+++ /issm/trunk/src/c/Makefile.am	(revision 387)
@@ -56,4 +56,6 @@
 					./objects/Penpair.cpp\
 					./objects/Penpair.h\
+					./objects/Pengrid.cpp\
+					./objects/Pengrid.h\
 					./objects/Icefront.cpp\
 					./objects/Icefront.h\
@@ -297,4 +299,6 @@
 					./objects/Penpair.cpp\
 					./objects/Penpair.h\
+					./objects/Pengrid.cpp\
+					./objects/Pengrid.h\
 					./objects/Icefront.cpp\
 					./objects/Icefront.h\
Index: /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp	(revision 387)
@@ -60,3 +60,4 @@
 		throw ErrorException(__FUNCT__,exprintf("%s%s%s"," analysis_type: ",model->analysis_type," not supported yet!"));
 	}
+
 }
Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp	(revision 387)
@@ -156,8 +156,8 @@
 	parameters->AddObject(param);
 
-	/*reconditioning_number: */
-	count++;
-	param= new Param(count,"reconditioning_number",DOUBLE);
-	param->SetDouble(model->reconditioning_number);
+	/*stokesreconditioning: */
+	count++;
+	param= new Param(count,"stokesreconditioning",DOUBLE);
+	param->SetDouble(model->stokesreconditioning);
 	parameters->AddObject(param);
 
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp	(revision 387)
@@ -99,5 +99,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
-	double penta_reconditioning_number;
+	double penta_stokesreconditioning;
 
 	/*matpar constructor input: */
@@ -414,5 +414,5 @@
 					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-					penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
+					penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); 
 
 			/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp	(revision 387)
@@ -80,5 +80,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
-	double penta_reconditioning_number;
+	double penta_stokesreconditioning;
 
 	/*matpar constructor input: */
@@ -255,6 +255,6 @@
 			penta_viscosity_overshoot=model->viscosity_overshoot;
 
-			/*reconditioning_number: */
-			penta_reconditioning_number=model->reconditioning_number;
+			/*stokesreconditioning: */
+			penta_stokesreconditioning=model->stokesreconditioning;
 
 			
@@ -265,5 +265,5 @@
 						penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 						penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-						penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
+						penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); 
 
 				/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp	(revision 387)
@@ -24,18 +24,11 @@
 	DataSet*    loads    = NULL;
 	Icefront*   icefront = NULL;
-	Penpair*    penpair  = NULL;
+	Pengrid*    pengrid  = NULL;
 
 	int segment_width;
-	int element_type;
 	int i1,i2,i3,i4;
 	int i0,range;
 	int grid1,grid2;
 
-	/*rift penpair: */
-	double  normal[2];
-	double  length;
-	double  friction;
-	int     fill;
-		
 	/*icefront intermediary data: */
 	char icefront_type[ICEFRONTSTRING];
@@ -48,25 +41,15 @@
 	double	icefront_b[MAX_ICEFRONT_GRIDS];
 
-	/*penpair intermediary data: */
-	int penpair_id;
-	int penpair_node_ids[2];
-	double penpair_penalty_offset;
-	int penpair_numdofs;
-	int penpair_dof;
-	int penpair_penalty_lock;
-	int penpair_element_ids[2];
-	double penpair_friction;
-	int penpair_fill;
-	double penpair_normal[2];
-	double penpair_length;
+	/*pengrid intermediary data: */
+	int pengrid_id;
+	int pengrid_node_id;
+	int pengrid_dof;
+	double pengrid_penalty_offset;
+	int pengrid_active=0;
+	int pengrid_thermal_steadystate=1;
 
-	/*Rifts:*/
-	int* riftsnumpenaltypairs=NULL;
-	double** riftspenaltypairs=NULL;
-	int* riftsfill=NULL;
-	int* riftsfriction=NULL;
-	double* riftpenaltypairs=NULL;
-	int el1,el2;
-	
+	int numberofsegs_diag_stokes;
+	int count;
+
 	/*Create loads: */
 	loads   = new DataSet(LoadsEnum());
@@ -77,23 +60,17 @@
 	/*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): */
-	ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat");
+	ModelFetchData((void**)&model->segmentonneumann_diag_stokes,&numberofsegs_diag_stokes,NULL,model_handle,"segmentonneumann_diag_stokes","Matrix","Mat");
 	ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat");
 	ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat");
 	ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat");
 
+	count=0;
+
 	/*First load data:*/
-	for (i=0;i<model->numberofsegs_diag;i++){
+	for (i=0;i<numberofsegs_diag_stokes;i++){
 		
-		if (strcmp(model->meshtype,"2d")==0){
-			segment_width=3;
-			element_type=TriaEnum();
-		}
-		else{
-			segment_width=5;
-			element_type=PentaEnum();
-		}
+		segment_width=5;
 
-
-		element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column
+		element=(int)(*(model->segmentonneumann_diag_stokes+segment_width*i+segment_width-1)-1); //element is in the last column
 
 		#ifdef _PARALLEL_
@@ -106,158 +83,76 @@
 	
 		icefront_mparid=model->numberofelements+1; //matlab indexing
-		icefront_sid=i+1; //matlab indexing
-		icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing
-		icefront_element_type=element_type;
+		icefront_sid=count+1; //matlab indexing
+		icefront_eid=(int)*(model->segmentonneumann_diag_stokes+segment_width*i+segment_width-1); //matlab indexing
+		icefront_element_type=PentaEnum();
 
-		i1=(int)*(model->segmentonneumann_diag+segment_width*i+0);
-		i2=(int)*(model->segmentonneumann_diag+segment_width*i+1);
-			
+		i1=(int)*(model->segmentonneumann_diag_stokes+segment_width*i+0);
+		i2=(int)*(model->segmentonneumann_diag_stokes+segment_width*i+1);
+		i3=(int)*(model->segmentonneumann_diag_stokes+segment_width*i+2);
+		i4=(int)*(model->segmentonneumann_diag_stokes+segment_width*i+3);
+	
 		icefront_node_ids[0]=i1;
 		icefront_node_ids[1]=i2;
-		
+		icefront_node_ids[2]=i3;
+		icefront_node_ids[3]=i4;
+	
 		icefront_h[0]=model->thickness[i1-1];
 		icefront_h[1]=model->thickness[i2-1];
+		icefront_h[2]=model->thickness[i3-1];
+		icefront_h[3]=model->thickness[i4-1];
 
 		icefront_b[0]=model->bed[i1-1];
 		icefront_b[1]=model->bed[i2-1];
-		
-		if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d
-			strcpy(icefront_type,"segment");
-		}
-		else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d.
-			strcpy(icefront_type,"quad");
-			i3=(int)*(model->segmentonneumann_diag+segment_width*i+2);
-			i4=(int)*(model->segmentonneumann_diag+segment_width*i+3);
-			icefront_node_ids[2]=i3;
-			icefront_node_ids[3]=i4;
-			
-			icefront_h[2]=model->thickness[i3-1];
-			icefront_h[3]=model->thickness[i4-1];
-
-			icefront_b[2]=model->bed[i3-1];
-			icefront_b[3]=model->bed[i4-1];
-		}
-		else{
-			throw ErrorException(__FUNCT__," element type not supported yet");
-		}
-
+		icefront_b[2]=model->bed[i3-1];
+		icefront_b[3]=model->bed[i4-1];
+	
 		icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b);
 		
 		loads->AddObject(icefront);
+		count++;
 
 	}
 	/*Free data: */
-	xfree((void**)&model->segmentonneumann_diag);
+	xfree((void**)&model->segmentonneumann_diag_stokes);
 	xfree((void**)&model->elements_type);
 	xfree((void**)&model->thickness);
 	xfree((void**)&model->bed);
 
+
+
+
+	//create penalties for grids on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy
+
+	ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonicesheet,NULL,NULL,model_handle,"gridonicesheet","Matrix","Mat");
+	ModelFetchData((void**)&model->gridonstokes,NULL,NULL,model_handle,"gridonstokes","Matrix","Mat");
+	
+	for (i=0;i<model->numberofnodes;i++){
+	#ifdef _PARALLEL_
+	/*keep only this partition's nodes:*/
+	if((model->my_grids[i]==1)){
+	#endif
+
+		if ((model->gridonbed[i]) && (model->gridonicesheet[i]) && (model->gridonstokes[i]))
+
+		pengrid_id=count+1; //matlab indexing
+		pengrid_node_id=i+1;
+		pengrid_dof=1;
+		pengrid_penalty_offset=model->penalty_offset;
+
+		pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate);
 		
-	/*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */
-
-	/*First fetch penalties: */
-	if (strcmp(model->meshtype,"3d")==0){
-		ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat");
+		loads->AddObject(icefront);
+		count++;
+	#ifdef _PARALLEL_
+	} //if((model->my_grids[i]==1))
+	#endif
 	}
 
-	counter=0;
-	if(model->numpenalties){
-
-		/*First deal with internal grids: */
-		for (i=0;i<model->numpenalties;i++){
-			if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU
-
-				for(j=0;j<model->numlayers-1;j++){
-
-					/*We are pairing grids along a vertical profile.*/
-					grid1=(int)*(model->penalties+model->numlayers*i+j);
-					grid2=(int)*(model->penalties+model->numlayers*i+j+1);
-
-					penpair_id=counter+1; //matlab indexing
-					penpair_dof=1;
-					penpair_node_ids[0]=grid1;
-					penpair_node_ids[1]=grid2;
-					penpair_penalty_offset=model->penalty_offset;
-					penpair_numdofs=1;
-
-					penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
-								penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
-					loads->AddObject(penpair);
 
 
-					counter++;
-
-					penpair_id=counter+1; //matlab indexing
-					penpair_dof=2;
-					penpair_node_ids[0]=grid1;
-					penpair_node_ids[1]=grid2;
-					penpair_penalty_offset=model->penalty_offset;
-					penpair_numdofs=1;
-
-					penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
-								penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
-					loads->AddObject(penpair);
-
-
-					counter++;
-
-				} //for ( i=0; i<numpenalties; i++) 
-			}
-		}
-	}
-
-	/*Free data: */
-	xfree((void**)&model->penalties);
-	
-	/*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */
-	/*First fetch rifts: */
-	if(model->numrifts){
-
-		for(i=0;i<model->numrifts;i++){
-			riftpenaltypairs=model->riftspenaltypairs[i];
-			for(j=0;j<model->riftsnumpenaltypairs[i];j++){
-
-				el1=(int)*(riftpenaltypairs+7*j+2);
-				#ifdef _PARALLEL_
-				if (model->epart[el1-1]!=my_rank){
-					/*This load does not belong to this cluster node, as it references an element 
-					 *that does not belong to this node's partition. Drop this 'j':*/
-					continue;
-				}
-				#endif
-
-				/*Ok, retrieve all the data needed to add a penalty between the two grids: */
-				el2=(int)*(riftpenaltypairs+7*j+3); 
-				grid1=(int)*(riftpenaltypairs+7*j+0); 
-				grid2=(int)*(riftpenaltypairs+7*j+1);
-				normal[0]=*(riftpenaltypairs+7*j+4);
-				normal[1]=*(riftpenaltypairs+7*j+5);
-				length=*(riftpenaltypairs+7*j+6);
-				friction=model->riftsfriction[i];
-				fill=model->riftsfill[i];
-
-				penpair_id=counter+1; //matlab indexing
-				penpair_node_ids[0]=grid1;
-				penpair_node_ids[1]=grid2;
-				penpair_element_ids[0]=el1;
-				penpair_element_ids[1]=el2;
-				penpair_penalty_offset=model->penalty_offset;
-				penpair_penalty_lock=model->penalty_lock;
-				penpair_numdofs=2;
-				penpair_normal[0]=normal[0];
-				penpair_normal[1]=normal[1];
-				penpair_length=length;
-				penpair_friction=friction;
-				penpair_fill=fill;
-
-				penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof,
-						penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length);
-				
-				loads->AddObject(penpair);
-
-				counter++;
-			}
-		}
-	}
+	xfree((void**)&model->gridonbed);
+	xfree((void**)&model->gridonstokes);
+	xfree((void**)&model->gridonicesheet);
 
 	/*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 
@@ -267,15 +162,4 @@
 	cleanup_and_return:
 
-	/*Free ressources:*/
-	xfree((void**)&riftsnumpenaltypairs); 
-	for(i=0;i<model->numrifts;i++){
-		double* temp=riftspenaltypairs[i];
-		xfree((void**)&temp);
-	}
-	xfree((void**)&riftspenaltypairs);
-	xfree((void**)&riftsfill);
-	xfree((void**)&riftsfriction);
-
-
 	/*Assign output pointer: */
 	*ploads=loads;
Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp	(revision 387)
@@ -79,5 +79,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
-	double penta_reconditioning_number;
+	double penta_stokesreconditioning;
 
 	/*matpar constructor input: */
@@ -239,5 +239,5 @@
 				penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 				penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-				penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
+				penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); 
 
 		/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/Model.cpp	(revision 387)
@@ -67,4 +67,5 @@
 	model->bed=NULL;
 	model->elementoniceshelf=NULL;
+	model->gridonicesheet=NULL;
 
 	model->drag_type=0;
@@ -75,4 +76,5 @@
 	model->numberofsegs_diag=0;
 	model->segmentonneumann_diag=NULL;
+	model->segmentonneumann_diag_stokes=NULL;
 	model-> neumannvalues_diag=NULL;
 	model-> gridondirichlet_diag=NULL;
@@ -126,5 +128,5 @@
 	model->yts=0;
 	model->viscosity_overshoot=0;
-	model->reconditioning_number=0;
+	model->stokesreconditioning=0;
 	model->waitonlock=0;
 
@@ -213,5 +215,7 @@
 	xfree((void**)&model->q);
 	xfree((void**)&model->elementoniceshelf);
+	xfree((void**)&model->gridonicesheet);
 	xfree((void**)&model->segmentonneumann_diag);
+	xfree((void**)&model->segmentonneumann_diag_stokes);
 	xfree((void**)&model->neumannvalues_diag);
 	xfree((void**)&model->gridondirichlet_diag);
@@ -255,4 +259,9 @@
 	
 	xfree((void**)&model->control_type);
+	
+	xfree((void**)&model->epart);
+	xfree((void**)&model->npart);
+	xfree((void**)&model->my_grids);
+	xfree((void**)&model->my_bordergrids);
 	
 	/*!Delete entire structure: */
@@ -333,4 +342,5 @@
 	ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
 	ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
+	ModelFetchData((void**)&model->stokesreconditioning,NULL,NULL,model_handle,"stokesreconditioning","Scalar",NULL);
 	ModelFetchData((void**)&model->reconditioning_number,NULL,NULL,model_handle,"reconditioning_number","Scalar",NULL);
 	ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);
Index: /issm/trunk/src/c/ModelProcessorx/Model.h
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/Model.h	(revision 387)
@@ -62,4 +62,5 @@
 	double* bed;
 	double* elementoniceshelf;
+	double* gridonicesheet;
 
 	/*friction: */
@@ -76,4 +77,6 @@
 	double* gridondirichlet_diag;
 	double* dirichletvalues_diag;
+	double* segmentonneumann_diag_stokes;
+
 	//prognostic
 	double* segmentonneumann_prog;
@@ -121,5 +124,5 @@
 	double  yts;
 	double  viscosity_overshoot;
-	double  reconditioning_number;
+	double  stokesreconditioning;
 	int     waitonlock;
 
Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp
===================================================================
--- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 386)
+++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp	(revision 387)
@@ -90,5 +90,5 @@
 	int penta_thermal_steadystate;
 	double penta_viscosity_overshoot;
-	double penta_reconditioning_number;
+	double penta_stokesreconditioning;
 
 	/* node constructor input: */
@@ -287,5 +287,5 @@
 					penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel,
 					penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff,
-					penta_thermal_steadystate,penta_viscosity_overshoot,penta_reconditioning_number); 
+					penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); 
 
 			/*Add penta element to elements dataset: */
Index: /issm/trunk/src/c/objects/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.cpp	(revision 387)
+++ /issm/trunk/src/c/objects/Pengrid.cpp	(revision 387)
@@ -0,0 +1,195 @@
+/*!\file Pengrid.c
+ * \brief: implementation of the Pengrid object
+ */
+
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "stdio.h"
+#include "./Pengrid.h"
+#include <string.h>
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../shared/shared.h"
+#include "../include/typedefs.h"
+
+		
+Pengrid::Pengrid(){
+	return;
+}
+
+Pengrid::Pengrid(int	pengrid_id, int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){
+	
+	id=pengrid_id;
+	dof=pengrid_dof;
+	active=pengrid_active;
+	penalty_offset =pengrid_penalty_offset;
+	thermal_steadystate=pengrid_thermal_steadystate;
+
+	node_id=pengrid_node_id;
+	node_offset=UNDEF;
+	node=NULL;
+
+	return;
+}
+
+Pengrid::~Pengrid(){
+	return;
+}
+		
+void Pengrid::Echo(void){
+
+	printf("Pengrid:\n");
+	printf("   id: %i\n",id);
+	printf("   dof: %i\n",dof);
+	printf("   active: %i\n",active);
+	printf("   penalty_offset: %g\n",penalty_offset);
+	printf("   thermal_steadystate: %i\n",thermal_steadystate);
+	printf("   node_id: [%i]\n",node_id);
+	printf("   node_offset: [%i]\n",node_offset);
+	
+	if(node)node->Echo();
+	return;
+}
+		
+void  Pengrid::Marshall(char** pmarshalled_dataset){
+
+	char* marshalled_dataset=NULL;
+	int   enum_type=0;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*get enum type of Pengrid: */
+	enum_type=PengridEnum();
+	
+	/*marshall enum: */
+	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
+	
+	/*marshall Pengrid data: */
+	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
+	memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
+	memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
+	memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
+	memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
+
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+		
+int   Pengrid::MarshallSize(){
+
+	return sizeof(id)+
+		sizeof(dof)+
+		sizeof(active)+
+		sizeof(penalty_offset)+
+		sizeof(thermal_steadystate)+
+		sizeof(node_id)+
+		sizeof(node_offset)+
+		sizeof(int); //sizeof(int) for enum type
+}
+
+char* Pengrid::GetName(void){
+	return "pengrid";
+}
+		
+
+void  Pengrid::Demarshall(char** pmarshalled_dataset){
+
+	int i;
+	char* marshalled_dataset=NULL;
+
+	/*recover marshalled_dataset: */
+	marshalled_dataset=*pmarshalled_dataset;
+
+	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
+	 *object data (thanks to DataSet::Demarshall):*/
+
+	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
+	memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
+	memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
+	memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
+	memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
+	memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
+	memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
+
+	node=NULL;
+
+	/*return: */
+	*pmarshalled_dataset=marshalled_dataset;
+	return;
+}
+int Pengrid::Enum(void){
+
+	return PengridEnum();
+}
+
+int    Pengrid::GetId(void){ return id; }
+
+int    Pengrid::MyRank(void){ 
+	extern int my_rank;
+	return my_rank; 
+}
+void  Pengrid::DistributeNumDofs(int* numdofspernode,int analysis_type){return;}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::Configure"
+
+void  Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
+
+	DataSet* nodesin=NULL;
+
+	/*Recover pointers :*/
+	nodesin=(DataSet*)pnodesin;
+
+	/*Link this load with its nodes: */
+	ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
+
+}
+
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::CreateKMatrix"
+
+void  Pengrid::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::CreatePVector"
+void  Pengrid::CreatePVector(Vec pg, void* inputs, int analysis_type){
+
+	/*No loads applied, do nothing: */
+	return;
+
+}
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::UpdateFromInputs"
+void  Pengrid::UpdateFromInputs(void* inputs){
+	
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyCreateKMatrix"
+void  Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){
+	throw ErrorException(__FUNCT__," not implemented yet!");
+}
+		
+#undef __FUNCT__ 
+#define __FUNCT__ "Pengrid::PenaltyCreatePVector"
+void  Pengrid::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){
+	throw ErrorException(__FUNCT__," not implemented yet!");
+}
+
+Object* Pengrid::copy() {
+	return new Pengrid(*this); 
+}
+
Index: /issm/trunk/src/c/objects/Pengrid.h
===================================================================
--- /issm/trunk/src/c/objects/Pengrid.h	(revision 387)
+++ /issm/trunk/src/c/objects/Pengrid.h	(revision 387)
@@ -0,0 +1,54 @@
+/*!\file Pengrid.h
+ * \brief: header file for pengrid object */
+
+#ifndef _PENGRID_H_
+#define _PENGRID_H_
+
+#include "./Load.h"
+#include "./Node.h"
+#include "./Element.h"
+
+class Element;
+class Pengrid: public Load{
+
+	private: 
+
+		int		id;
+		int     dof;
+		int     active;
+		double  penalty_offset; 
+		int     thermal_steadystate;
+		
+		/*nodes: */
+		int     node_id;
+		Node*   node;
+		int     node_offset;
+
+	public:
+
+		Pengrid();
+		Pengrid(int	id, int node_id,int dof, int active, double penalty_offset,int thermal_steadystate);
+		~Pengrid();
+
+		void  Echo();
+		void  Marshall(char** pmarshalled_dataset);
+		int   MarshallSize();
+		char* GetName();
+		void  Demarshall(char** pmarshalled_dataset);
+		int   Enum();
+		int   GetId(); 
+		int   MyRank();
+		void  DistributeNumDofs(int* numdofspernode,int analysis_type);
+		void  Configure(void* elements,void* nodes,void* materials);
+		void  CreateKMatrix(Mat Kgg,void* inputs,int analysis_type);
+		void  CreatePVector(Vec pg, void* inputs, int analysis_type);
+		void  UpdateFromInputs(void* inputs);
+		void  PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type);
+		void  PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type);
+		Object* copy();
+
+};
+
+#endif  /* _PENGRID_H_ */
+
+
Index: /issm/trunk/src/c/objects/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Penta.cpp	(revision 386)
+++ /issm/trunk/src/c/objects/Penta.cpp	(revision 387)
@@ -22,5 +22,5 @@
 				double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel, 
 				int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], 
-				int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_reconditioning_number){
+				int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
 	
 	int i;
@@ -59,5 +59,5 @@
 	thermal_steadystate = penta_thermal_steadystate;
 	viscosity_overshoot = penta_viscosity_overshoot;
-	reconditioning_number = penta_reconditioning_number;
+	stokesreconditioning = penta_stokesreconditioning;
 
 	return;
@@ -100,5 +100,5 @@
 	printf("   thermal_steadystate: %i\n",thermal_steadystate);
 	printf("   viscosity_overshoot: %i\n",viscosity_overshoot);
-	printf("   reconditioning_number: %i\n",reconditioning_number);
+	printf("   stokesreconditioning: %i\n",stokesreconditioning);
 	return;
 }
@@ -148,5 +148,5 @@
 	memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
 	memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
-	memcpy(marshalled_dataset,&reconditioning_number,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
+	memcpy(marshalled_dataset,&stokesreconditioning,sizeof(stokesreconditioning));marshalled_dataset+=sizeof(stokesreconditioning);
 	
 	*pmarshalled_dataset=marshalled_dataset;
@@ -185,5 +185,5 @@
 		sizeof(thermal_steadystate) +
 		sizeof(viscosity_overshoot) +
-		sizeof(reconditioning_number) +
+		sizeof(stokesreconditioning) +
 		sizeof(int); //sizeof(int) for enum type
 }
@@ -233,5 +233,5 @@
 	memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
 	memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
-	memcpy(&reconditioning_number,marshalled_dataset,sizeof(reconditioning_number));marshalled_dataset+=sizeof(reconditioning_number);
+	memcpy(&stokesreconditioning,marshalled_dataset,sizeof(stokesreconditioning));marshalled_dataset+=sizeof(stokesreconditioning);
 
 	/*nodes and materials are not pointing to correct objects anymore:*/
Index: /issm/trunk/src/c/objects/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Penta.h	(revision 386)
+++ /issm/trunk/src/c/objects/Penta.h	(revision 387)
@@ -53,5 +53,5 @@
 		int    thermal_steadystate;
 		double viscosity_overshoot;
-		double reconditioning_number;
+		double stokesreconditioning;
 	
 	public:
@@ -61,5 +61,5 @@
 				double p, double q, int    shelf, int    onbed, int    onsurface, double meanvel,double epsvel, 
 				int    collapse, double melting[6], double accumulation[6], double geothermalflux[6], 
-				int    artdiff, int    thermal_steadystate,double viscosity_overshoot,double reconditioning_number);
+				int    artdiff, int    thermal_steadystate,double viscosity_overshoot,double stokesreconditioning);
 		~Penta();
 
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 386)
+++ /issm/trunk/src/c/objects/objects.h	(revision 387)
@@ -22,4 +22,5 @@
 #include "./Icefront.h"
 #include "./Penpair.h"
+#include "./Pengrid.h"
 #include "./Param.h"
 #include "./Element.h" 
Index: /issm/trunk/src/m/classes/@model/model.m
===================================================================
--- /issm/trunk/src/m/classes/@model/model.m	(revision 386)
+++ /issm/trunk/src/m/classes/@model/model.m	(revision 387)
@@ -131,4 +131,7 @@
 	md.gridondirichlet_diag=NaN;
 	md.dirichletvalues_diag=NaN;
+	
+	%Diagnostic stokes
+	md.segmentonneumann_diag_stokes=NaN;
 
 	%Thermal
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 386)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 387)
@@ -126,4 +126,7 @@
 md.borderstokes=borderflags-md.gridonstokes; %remove stokes grids from this list
 
+%icefront for stokes
+md.segmentonneumann_diag_stokes=md.segmentonneumann_diag(find(md.elements_type(md.segmentonneumann_diag(:,end),2)),:);
+
 %figure out solution types
 md.ishutter=double(any(md.elements_type(:,1)==hutterenum()));
Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 386)
+++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m	(revision 387)
@@ -72,5 +72,5 @@
 
 		%"recondition" pressure 
-		p_g=p_g/m_ds.parameters.reconditioning_number;
+		p_g=p_g/m_ds.parameters.stokesreconditioning;
 
 		displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
@@ -95,5 +95,5 @@
 	
 		%"decondition" pressure
-		p_g=u_g(4:4:end)*m_dh.parameters.reconditioning_number;
+		p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
 	end
 end