Index: /issm/trunk/src/c/BuildNodeSetsx/BuildNodeSetsx.cpp
===================================================================
--- /issm/trunk/src/c/BuildNodeSetsx/BuildNodeSetsx.cpp	(revision 303)
+++ /issm/trunk/src/c/BuildNodeSetsx/BuildNodeSetsx.cpp	(revision 304)
@@ -40,45 +40,48 @@
 	gsize=nodes->NumberOfDofs();
 
-	/*Now, allocate parallel vectors for all the sets: */
-	flag_pv_g=NewVec(gsize);
-	flag_pv_m=NewVec(gsize);
-	flag_pv_n=NewVec(gsize);
-	flag_pv_f=NewVec(gsize);
-	flag_pv_s=NewVec(gsize);
+	if(gsize){
 
-	/*Now, go through all nodes and have every node plug into 
-	 * pv_m, pv_n, pv_f and pv_s, their own set flags, for each dof: */
-	nodes->FlagNodeSets(flag_pv_g,flag_pv_m,flag_pv_n,flag_pv_f,flag_pv_s);
+		/*Now, allocate parallel vectors for all the sets: */
+		flag_pv_g=NewVec(gsize);
+		flag_pv_m=NewVec(gsize);
+		flag_pv_n=NewVec(gsize);
+		flag_pv_f=NewVec(gsize);
+		flag_pv_s=NewVec(gsize);
 
-	/*Now, every cpu has 4 flag vectors, of size gsize. Create partition vectors (like a pos=find(flag_pv_g) in matlab*/
-	PartitionSets(&vec_pv_m,&vec_pv_n,flag_pv_g,flag_pv_m,flag_pv_n,gsize); /*! split g set into m and n sets*/
-	PartitionSets(&vec_pv_f,&vec_pv_s,flag_pv_n,flag_pv_f,flag_pv_s,gsize); /*! split n set into f and s sets*/
+		/*Now, go through all nodes and have every node plug into 
+		 * pv_m, pv_n, pv_f and pv_s, their own set flags, for each dof: */
+		nodes->FlagNodeSets(flag_pv_g,flag_pv_m,flag_pv_n,flag_pv_f,flag_pv_s);
 
-	/*Free ressources:*/
-	VecFree(&flag_pv_g);
-	VecFree(&flag_pv_m);
-	VecFree(&flag_pv_n);
-	VecFree(&flag_pv_f);
-	VecFree(&flag_pv_s);
+		/*Now, every cpu has 4 flag vectors, of size gsize. Create partition vectors (like a pos=find(flag_pv_g) in matlab*/
+		PartitionSets(&vec_pv_m,&vec_pv_n,flag_pv_g,flag_pv_m,flag_pv_n,gsize); /*! split g set into m and n sets*/
+		PartitionSets(&vec_pv_f,&vec_pv_s,flag_pv_n,flag_pv_f,flag_pv_s,gsize); /*! split n set into f and s sets*/
 
-	/*Serialize vectors: */
-	VecGetSize(vec_pv_m,&msize);
-	VecGetSize(vec_pv_n,&nsize);
-	VecGetSize(vec_pv_f,&fsize);
-	VecGetSize(vec_pv_s,&ssize);
-	VecToMPISerial(&pv_m,vec_pv_m);
-	VecToMPISerial(&pv_n,vec_pv_n);
-	VecToMPISerial(&pv_f,vec_pv_f);
-	VecToMPISerial(&pv_s,vec_pv_s);
+		/*Free ressources:*/
+		VecFree(&flag_pv_g);
+		VecFree(&flag_pv_m);
+		VecFree(&flag_pv_n);
+		VecFree(&flag_pv_f);
+		VecFree(&flag_pv_s);
 
-	/*Free ressources:*/
-	VecFree(&vec_pv_m);
-	VecFree(&vec_pv_n);
-	VecFree(&vec_pv_f);
-	VecFree(&vec_pv_s);
+		/*Serialize vectors: */
+		VecGetSize(vec_pv_m,&msize);
+		VecGetSize(vec_pv_n,&nsize);
+		VecGetSize(vec_pv_f,&fsize);
+		VecGetSize(vec_pv_s,&ssize);
+		VecToMPISerial(&pv_m,vec_pv_m);
+		VecToMPISerial(&pv_n,vec_pv_n);
+		VecToMPISerial(&pv_f,vec_pv_f);
+		VecToMPISerial(&pv_s,vec_pv_s);
+
+		/*Free ressources:*/
+		VecFree(&vec_pv_m);
+		VecFree(&vec_pv_n);
+		VecFree(&vec_pv_f);
+		VecFree(&vec_pv_s);
 
 
-	/*Create NodeSets* object: */
-	nodesets=new NodeSets(pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
+		/*Create NodeSets* object: */
+		nodesets=new NodeSets(pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
+	}
 
 	/*Assign output pointers: */
Index: /issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 303)
+++ /issm/trunk/src/c/DataSet/DataSet.cpp	(revision 304)
@@ -150,7 +150,8 @@
 	vector<Object*>::iterator iterator;
 
-	iterator = find(objects.begin(), objects.end(),object);
-
-	objects.erase(iterator);
+	if(object){
+		iterator = find(objects.begin(), objects.end(),object);
+		objects.erase(iterator);
+	}
 
 }
@@ -1139,2 +1140,21 @@
 
 }
+
+
+void DataSet::ComputePressure(Vec p_g){
+
+	vector<Object*>::iterator object;
+	Element* element=NULL;
+
+	for ( object=objects.begin() ; object < objects.end(); object++ ){
+		
+		if(EnumIsElement((*object)->Enum())){
+
+			element=(Element*)(*object);
+			element->ComputePressure(p_g);
+		}
+	}
+
+}
+
+
Index: /issm/trunk/src/c/DataSet/DataSet.h
===================================================================
--- /issm/trunk/src/c/DataSet/DataSet.h	(revision 303)
+++ /issm/trunk/src/c/DataSet/DataSet.h	(revision 304)
@@ -76,4 +76,5 @@
 		void  VelocityExtrude(Vec ug,double* ug_serial);
 		int   DeleteObject(Object* object);
+		void  ComputePressure(Vec p_g);
 
 };
Index: /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 303)
+++ /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp	(revision 304)
@@ -39,4 +39,13 @@
 		return DiagnosticStokesAnalysisEnum();
 	}
+	else if (strcmp(analysis_type,"diagnostic_hutter")==0){
+		return DiagnosticHutterAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"surface_slope_compute")==0){
+		return SurfaceSlopeComputeAnalysisEnum();
+	}
+	else if (strcmp(analysis_type,"bed_slope_compute")==0){
+		return BedSlopeComputeAnalysisEnum();
+	}
 	else throw ErrorException(__FUNCT__," could not recognized analysis type!");
 
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 303)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp	(revision 304)
@@ -25,4 +25,7 @@
 int MeltingAnalysisEnum(void){          return          29; }
 int DiagnosticStokesAnalysisEnum(void){ return          30; }
+int DiagnosticHutterAnalysisEnum(void){ return          31; }
+int SurfaceSlopeComputeAnalysisEnum(void){ return       32; }
+int BedSlopeComputeAnalysisEnum(void){ return           33; }
 
 /*datasets: */
Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 303)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 304)
@@ -39,4 +39,7 @@
 int MeltingAnalysisEnum(void);
 int DiagnosticStokesAnalysisEnum(void);
+int DiagnosticHutterAnalysisEnum(void);
+int SurfaceSlopeComputeAnalysisEnum(void);
+int BedSlopeComputeAnalysisEnum(void);
 
 
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 303)
+++ /issm/trunk/src/c/Makefile.am	(revision 304)
@@ -165,7 +165,20 @@
 					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
 					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
 					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
 					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
 					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
+					./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
+					./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
+					./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
+					./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
+					./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
+					./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
+					./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
+					./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
 					./ModelProcessorx/Control/CreateParametersControl.cpp\
 					./Dofx/Dofx.h\
@@ -187,4 +200,6 @@
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
+					./ComputePressurex/ComputePressurex.h\
+					./ComputePressurex/ComputePressurex.cpp\
 					./BuildNodeSetsx/BuildNodeSetsx.h\
 					./BuildNodeSetsx/BuildNodeSetsx.cpp\
@@ -387,7 +402,20 @@
 					./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
 					./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
+					./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
 					./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
 					./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
 					./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
+					./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
+					./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
+					./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
+					./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
+					./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
+					./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
+					./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
+					./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
+					./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
 					./ModelProcessorx/Control/CreateParametersControl.cpp\
 					./Dofx/Dofx.h\
@@ -409,4 +437,6 @@
 					./ConfigureObjectsx/ConfigureObjectsx.h\
 					./ConfigureObjectsx/ConfigureObjectsx.cpp\
+					./ComputePressurex/ComputePressurex.h\
+					./ComputePressurex/ComputePressurex.cpp\
 					./BuildNodeSetsx/BuildNodeSetsx.h\
 					./BuildNodeSetsx/BuildNodeSetsx.cpp\
Index: /issm/trunk/src/c/NormalizeConstraintsx/NormalizeConstraintsx.cpp
===================================================================
--- /issm/trunk/src/c/NormalizeConstraintsx/NormalizeConstraintsx.cpp	(revision 303)
+++ /issm/trunk/src/c/NormalizeConstraintsx/NormalizeConstraintsx.cpp	(revision 304)
@@ -26,38 +26,40 @@
 	int     null=0;
 
-	/*Recover data: */
-	msize=nodesets->GetMSize();
-	nsize=nodesets->GetNSize();
-	pv_m=nodesets->GetPV_M();
-	pv_n=nodesets->GetPV_N();
+	if(nodesets){
+		/*Recover data: */
+		msize=nodesets->GetMSize();
+		nsize=nodesets->GetNSize();
+		pv_m=nodesets->GetPV_M();
+		pv_n=nodesets->GetPV_N();
 
 
-	if (msize){
+		if (msize){
 
-		/*Build row partitioning vector for Rmm and Rmn: */
-		row_m=(double*)xmalloc(msize*sizeof(double));
-		for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing
+			/*Build row partitioning vector for Rmm and Rmn: */
+			row_m=(double*)xmalloc(msize*sizeof(double));
+			for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing
 
-		/*Partition Rmg into Rmm and Rmn: */
-		MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
-		MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);
+			/*Partition Rmg into Rmm and Rmn: */
+			MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
+			MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);
 
-		/*Invert Rmm: */
-		MatInvert(&InvRmm,Rmm);
+			/*Invert Rmm: */
+			MatInvert(&InvRmm,Rmm);
 
-		/*Do Gmn=-(Rmm-1)*Rmn :*/
-		MatScale(InvRmm,-1.0);
-		MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);
+			/*Do Gmn=-(Rmm-1)*Rmn :*/
+			MatScale(InvRmm,-1.0);
+			MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);
 
+		}
+		else{
+			Gmn=NULL;
+		}
+
+		/*Free ressources:*/
+		xfree((void**)&row_m);
+		MatFree(&Rmm);
+		MatFree(&Rmn);
+		MatFree(&InvRmm);
 	}
-	else{
-		Gmn=NULL;
-	}
-
-	/*Free ressources:*/
-	xfree((void**)&row_m);
-	MatFree(&Rmm);
-	MatFree(&Rmn);
-	MatFree(&InvRmm);
 	
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp
===================================================================
--- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 303)
+++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp	(revision 304)
@@ -51,6 +51,6 @@
 
 	/*First serialize partition vector: */
-	VecToMPISerial(&partition,part);
-	
+	if(part)VecToMPISerial(&partition,part);
+
 	
 	if (   (analysis_type==ControlAnalysisEnum()) || 
@@ -64,40 +64,41 @@
 		parameters->FindParam((void*)&vz,"vz");
 
-		u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
+		if(numberofnodes){
+			u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
 
-		if(vx){
-			for(i=0;i<numberofnodes;i++){
-				u_g[(int)(3*partition[i]+0)]=vx[i];  
+			if(vx){
+				for(i=0;i<numberofnodes;i++){
+					u_g[(int)(3*partition[i]+0)]=vx[i];  
+				}
 			}
-		}
 
-		if(vy){
-			for(i=0;i<numberofnodes;i++){
-				u_g[(int)(3*partition[i]+1)]=vy[i];  
+			if(vy){
+				for(i=0;i<numberofnodes;i++){
+					u_g[(int)(3*partition[i]+1)]=vy[i];  
+				}
 			}
-		}
-
-		if(vz){
-			for(i=0;i<numberofnodes;i++){
-				u_g[(int)(3*partition[i]+2)]=vz[i];  
+			if(vz){
+				for(i=0;i<numberofnodes;i++){
+					u_g[(int)(3*partition[i]+2)]=vz[i];  
+				}
 			}
-		}
 
 
-		/*Now, create new parameters: */
-		count++;
-		param= new Param(count,"u_g",DOUBLEVEC);
-		param->SetDoubleVec(u_g,3*numberofnodes);
-		parameters->AddObject(param);
+			/*Now, create new parameters: */
+			count++;
+			param= new Param(count,"u_g",DOUBLEVEC);
+			param->SetDoubleVec(u_g,3*numberofnodes);
+			parameters->AddObject(param);
 
-		/*erase old parameters: */
-		param=(Param*)parameters->FindParamObject("vx");
-		parameters->DeleteObject((Object*)param);
+			/*erase old parameters: */
+			param=(Param*)parameters->FindParamObject("vx");
+			parameters->DeleteObject((Object*)param);
 
-		param=(Param*)parameters->FindParamObject("vy");
-		parameters->DeleteObject((Object*)param);
+			param=(Param*)parameters->FindParamObject("vy");
+			parameters->DeleteObject((Object*)param);
 
-		param=(Param*)parameters->FindParamObject("vz");
-		parameters->DeleteObject((Object*)param);
+			param=(Param*)parameters->FindParamObject("vz");
+			parameters->DeleteObject((Object*)param);
+		}
 
 	}
Index: /issm/trunk/src/c/Reducevectorgtosx/Reducevectorgtosx.cpp
===================================================================
--- /issm/trunk/src/c/Reducevectorgtosx/Reducevectorgtosx.cpp	(revision 303)
+++ /issm/trunk/src/c/Reducevectorgtosx/Reducevectorgtosx.cpp	(revision 304)
@@ -18,22 +18,26 @@
 	Vec yn=NULL;
 
-	if (nodesets->GetNSize() && nodesets->GetSSize()){
+	if(nodesets){
 
-		VecPartition(&yn,yg,nodesets->GetPV_N(),nodesets->GetNSize());
-	
-		VecPartition(&ys,yn,nodesets->GetPV_S(),nodesets->GetSSize());
-	
+
+		if (nodesets->GetNSize() && nodesets->GetSSize()){
+
+			VecPartition(&yn,yg,nodesets->GetPV_N(),nodesets->GetNSize());
+		
+			VecPartition(&ys,yn,nodesets->GetPV_S(),nodesets->GetSSize());
+		
+		}
+
+		/*Create ys0, full of 0: */
+		if(ys){
+			VecDuplicate(ys,&ys0);
+			VecSet(ys0,0.0);
+			VecAssemblyBegin(ys0);
+			VecAssemblyEnd(ys0);
+		}
+
+		/*Free ressources:*/
+		VecFree(&yn);
 	}
-
-	/*Create ys0, full of 0: */
-	if(ys){
-		VecDuplicate(ys,&ys0);
-		VecSet(ys0,0.0);
-		VecAssemblyBegin(ys0);
-		VecAssemblyEnd(ys0);
-	}
-
-	/*Free ressources:*/
-	VecFree(&yn);
 	
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/SpcNodesx/SpcNodesx.cpp
===================================================================
--- /issm/trunk/src/c/SpcNodesx/SpcNodesx.cpp	(revision 303)
+++ /issm/trunk/src/c/SpcNodesx/SpcNodesx.cpp	(revision 304)
@@ -24,9 +24,11 @@
 	numberofdofs=nodes->NumberOfDofs();
 
-	/*Allocate yg: */
-	yg=NewVec(numberofdofs);
+	if(numberofdofs){
+		/*Allocate yg: */
+		yg=NewVec(numberofdofs);
 
-	/*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
-	constraints->SetupSpcs(nodes,yg);
+		/*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
+		constraints->SetupSpcs(nodes,yg);
+	}
 
 	/*Assign output pointers: */
Index: /issm/trunk/src/c/io/FetchNodeSets.cpp
===================================================================
--- /issm/trunk/src/c/io/FetchNodeSets.cpp	(revision 303)
+++ /issm/trunk/src/c/io/FetchNodeSets.cpp	(revision 304)
@@ -32,16 +32,22 @@
 	int ssize;
 
-	FetchData((void**)&pv_m,NULL,NULL,mxGetField(dataref,0,"pv_m"),"Vector","Vec");
-	FetchData((void**)&pv_n,NULL,NULL,mxGetField(dataref,0,"pv_n"),"Vector","Vec");
-	FetchData((void**)&pv_f,NULL,NULL,mxGetField(dataref,0,"pv_f"),"Vector","Vec");
-	FetchData((void**)&pv_s,NULL,NULL,mxGetField(dataref,0,"pv_s"),"Vector","Vec");
-	
-	gsize=(int)mxGetScalar(mxGetField(dataref,0,"gsize"));
-	msize=(int)mxGetScalar(mxGetField(dataref,0,"msize"));
-	nsize=(int)mxGetScalar(mxGetField(dataref,0,"nsize"));
-	fsize=(int)mxGetScalar(mxGetField(dataref,0,"fsize"));
-	ssize=(int)mxGetScalar(mxGetField(dataref,0,"ssize"));
+	if(mxIsEmpty(dataref)){
+		nodesets=NULL;
+	}
+	else{
 
-	nodesets=new NodeSets( pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
+		FetchData((void**)&pv_m,NULL,NULL,mxGetField(dataref,0,"pv_m"),"Vector","Vec");
+		FetchData((void**)&pv_n,NULL,NULL,mxGetField(dataref,0,"pv_n"),"Vector","Vec");
+		FetchData((void**)&pv_f,NULL,NULL,mxGetField(dataref,0,"pv_f"),"Vector","Vec");
+		FetchData((void**)&pv_s,NULL,NULL,mxGetField(dataref,0,"pv_s"),"Vector","Vec");
+		
+		gsize=(int)mxGetScalar(mxGetField(dataref,0,"gsize"));
+		msize=(int)mxGetScalar(mxGetField(dataref,0,"msize"));
+		nsize=(int)mxGetScalar(mxGetField(dataref,0,"nsize"));
+		fsize=(int)mxGetScalar(mxGetField(dataref,0,"fsize"));
+		ssize=(int)mxGetScalar(mxGetField(dataref,0,"ssize"));
+
+		nodesets=new NodeSets( pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
+	}
 
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/io/WriteNodeSets.cpp
===================================================================
--- /issm/trunk/src/c/io/WriteNodeSets.cpp	(revision 303)
+++ /issm/trunk/src/c/io/WriteNodeSets.cpp	(revision 304)
@@ -41,83 +41,88 @@
 	mxArray*    field=NULL;
 
-	/*Recover data from the nodesets class: */
-	gsize=nodesets->GetGSize();
-	msize=nodesets->GetMSize();
-	nsize=nodesets->GetNSize();
-	fsize=nodesets->GetFSize();
-	ssize=nodesets->GetSSize();
+	if(nodesets){
+		/*Recover data from the nodesets class: */
+		gsize=nodesets->GetGSize();
+		msize=nodesets->GetMSize();
+		nsize=nodesets->GetNSize();
+		fsize=nodesets->GetFSize();
+		ssize=nodesets->GetSSize();
 
-	if(msize){
-		pv_m=(double*)xmalloc(msize*sizeof(double));
-		memcpy(pv_m,nodesets->GetPV_M(),msize*sizeof(double));
-	}
-	if(nsize){
-		pv_n=(double*)xmalloc(nsize*sizeof(double));
-		memcpy(pv_n,nodesets->GetPV_N(),nsize*sizeof(double));
-	}
-	if(fsize){
-		pv_f=(double*)xmalloc(fsize*sizeof(double));
-		memcpy(pv_f,nodesets->GetPV_F(),fsize*sizeof(double));
-	}
-	if(ssize){
-		pv_s=(double*)xmalloc(ssize*sizeof(double));
-		memcpy(pv_s,nodesets->GetPV_S(),ssize*sizeof(double));
-	}
+		if(msize){
+			pv_m=(double*)xmalloc(msize*sizeof(double));
+			memcpy(pv_m,nodesets->GetPV_M(),msize*sizeof(double));
+		}
+		if(nsize){
+			pv_n=(double*)xmalloc(nsize*sizeof(double));
+			memcpy(pv_n,nodesets->GetPV_N(),nsize*sizeof(double));
+		}
+		if(fsize){
+			pv_f=(double*)xmalloc(fsize*sizeof(double));
+			memcpy(pv_f,nodesets->GetPV_F(),fsize*sizeof(double));
+		}
+		if(ssize){
+			pv_s=(double*)xmalloc(ssize*sizeof(double));
+			memcpy(pv_s,nodesets->GetPV_S(),ssize*sizeof(double));
+		}
 
-	/*Build structure in matlab workspace with all these fields: */
-	fnames[0] = "gsize";
-	fnames[1] = "msize";
-	fnames[2] = "nsize";
-	fnames[3] = "fsize";
-	fnames[4] = "ssize";
-	fnames[5] = "pv_m";
-	fnames[6] = "pv_n";
-	fnames[7] = "pv_f";
-	fnames[8] = "pv_s";
+		/*Build structure in matlab workspace with all these fields: */
+		fnames[0] = "gsize";
+		fnames[1] = "msize";
+		fnames[2] = "nsize";
+		fnames[3] = "fsize";
+		fnames[4] = "ssize";
+		fnames[5] = "pv_m";
+		fnames[6] = "pv_n";
+		fnames[7] = "pv_f";
+		fnames[8] = "pv_s";
 
-	dataref=mxCreateStructArray( ndim,onebyone,nfields,fnames);
-	
-	mxSetField( dataref, 0, "gsize",mxCreateDoubleScalar((double)gsize));
-	mxSetField( dataref, 0, "msize",mxCreateDoubleScalar((double)msize));
-	mxSetField( dataref, 0, "nsize",mxCreateDoubleScalar((double)nsize));
-	mxSetField( dataref, 0, "fsize",mxCreateDoubleScalar((double)fsize));
-	mxSetField( dataref, 0, "ssize",mxCreateDoubleScalar((double)ssize));
+		dataref=mxCreateStructArray( ndim,onebyone,nfields,fnames);
+		
+		mxSetField( dataref, 0, "gsize",mxCreateDoubleScalar((double)gsize));
+		mxSetField( dataref, 0, "msize",mxCreateDoubleScalar((double)msize));
+		mxSetField( dataref, 0, "nsize",mxCreateDoubleScalar((double)nsize));
+		mxSetField( dataref, 0, "fsize",mxCreateDoubleScalar((double)fsize));
+		mxSetField( dataref, 0, "ssize",mxCreateDoubleScalar((double)ssize));
 
-	if(msize){
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-		mxSetM(field,msize); mxSetN(field,1); mxSetPr(field,pv_m);
+		if(msize){
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+			mxSetM(field,msize); mxSetN(field,1); mxSetPr(field,pv_m);
+		}
+		else{
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+		}
+		mxSetField( dataref, 0, "pv_m",field);
+
+		if(nsize){
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+			mxSetM(field,nsize); mxSetN(field,1); mxSetPr(field,pv_n);
+		}
+		else{
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+		}
+		mxSetField( dataref, 0, "pv_n",field);
+
+		
+		if(fsize){
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+			mxSetM(field,fsize); mxSetN(field,1); mxSetPr(field,pv_f);
+		}
+		else{
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+		}
+		mxSetField( dataref, 0, "pv_f",field);
+		
+		if(ssize){
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+			mxSetM(field,ssize); mxSetN(field,1); mxSetPr(field,pv_s);
+		}
+		else{
+			field = mxCreateDoubleMatrix(0,0,mxREAL);
+		}
+		mxSetField( dataref, 0, "pv_s",field);
 	}
 	else{
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
+		dataref=mxCreateDoubleMatrix(0,0,mxREAL);
 	}
-	mxSetField( dataref, 0, "pv_m",field);
-
-	if(nsize){
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-		mxSetM(field,nsize); mxSetN(field,1); mxSetPr(field,pv_n);
-	}
-	else{
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-	}
-	mxSetField( dataref, 0, "pv_n",field);
-
-	
-	if(fsize){
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-		mxSetM(field,fsize); mxSetN(field,1); mxSetPr(field,pv_f);
-	}
-	else{
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-	}
-	mxSetField( dataref, 0, "pv_f",field);
-	
-	if(ssize){
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-		mxSetM(field,ssize); mxSetN(field,1); mxSetPr(field,pv_s);
-	}
-	else{
-		field = mxCreateDoubleMatrix(0,0,mxREAL);
-	}
-	mxSetField( dataref, 0, "pv_s",field);
 
 	/*Assign output pointers:*/
Index: /issm/trunk/src/c/issm.h
===================================================================
--- /issm/trunk/src/c/issm.h	(revision 303)
+++ /issm/trunk/src/c/issm.h	(revision 304)
@@ -46,4 +46,5 @@
 #include "./VelocityExtrudex/VelocityExtrudex.h"
 #include "./GriddataMeshToGridx/GriddataMeshToGridx.h"
+#include "./ComputePressurex/ComputePressurex.h"
 
 
Index: /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp
===================================================================
--- /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp	(revision 303)
+++ /issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp	(revision 304)
@@ -22,4 +22,16 @@
 		numdofs=1;
 	}
+	else if (analysis_type==DiagnosticStokesAnalysisEnum()){
+		numdofs=4;
+	}
+	else if (analysis_type==DiagnosticHutterAnalysisEnum()){
+		numdofs=2;
+	}
+	else if (analysis_type==SurfaceSlopeComputeAnalysisEnum()){
+		numdofs=1;
+	}
+	else if (analysis_type==BedSlopeComputeAnalysisEnum()){
+		numdofs=1;
+	}
 	else if (analysis_type==ControlAnalysisEnum()){
 		numdofs=2;
Index: /issm/trunk/src/c/toolkits/petsc/patches/PetscVectorToMatlabVector.cpp
===================================================================
--- /issm/trunk/src/c/toolkits/petsc/patches/PetscVectorToMatlabVector.cpp	(revision 303)
+++ /issm/trunk/src/c/toolkits/petsc/patches/PetscVectorToMatlabVector.cpp	(revision 304)
@@ -39,10 +39,11 @@
 	VecGetSize(vector,&rows);
 
-	idxm=(int*)xmalloc(rows*sizeof(int));
-	values=(double*)xmalloc(rows*sizeof(double));
-	for(i=0;i<rows;i++)idxm[i]=i;
-	 
-	VecGetValues(vector,rows,idxm,values);
-
+	if(rows){
+		idxm=(int*)xmalloc(rows*sizeof(int));
+		values=(double*)xmalloc(rows*sizeof(double));
+		for(i=0;i<rows;i++)idxm[i]=i;
+		 
+		VecGetValues(vector,rows,idxm,values);
+	}
 
 	/*Using values, build a matlab vector: */
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 303)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 304)
@@ -121,7 +121,7 @@
 
 %figure out solution types
-md.ishutter=any(md.elements_type(:,1)==hutterenum());
-md.ismacayealpattyn=any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() );
-md.isstokes=any(md.elements_type(:,2));
+md.ishutter=double(any(md.elements_type(:,1)==hutterenum()));
+md.ismacayealpattyn=double(any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() ));
+md.isstokes=double(any(md.elements_type(:,2)));
 
 end
Index: /issm/trunk/src/m/classes/public/solve.m
===================================================================
--- /issm/trunk/src/m/classes/public/solve.m	(revision 303)
+++ /issm/trunk/src/m/classes/public/solve.m	(revision 304)
@@ -54,5 +54,5 @@
 end
 
-disp(sprintf('\nlaunching solution sequence\n'));
+displaystring(md.debug,'%s','\nlaunching solution sequence\n');
 
 %If running in parallel, we have a different way of launching the solution
@@ -97,5 +97,5 @@
 
 %Check result is consistent
-disp(sprintf('%s\n','checking result consistency'));
+displaystring(md.debug,'%s\n','checking result consistency');
 if ~isresultconsistent(md,solutiontype),
 	disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)
Index: /issm/trunk/src/m/utils/String/displaystring.m
===================================================================
--- /issm/trunk/src/m/utils/String/displaystring.m	(revision 304)
+++ /issm/trunk/src/m/utils/String/displaystring.m	(revision 304)
@@ -0,0 +1,11 @@
+function displaystring(debug,format,varargin)
+%DISPLAY -  display string in solution sequences. wrapper to disp and sprintf.  
+%
+%   Usage:
+%      displaystring(debug,format,string)
+%      ex:   display(1,'%s\n','string to display');
+%      debug can be used to switch display on and off
+	
+if debug,
+	disp(sprintf(format,varargin{:}));
+end
Index: /issm/trunk/src/mex/Makefile.am
===================================================================
--- /issm/trunk/src/mex/Makefile.am	(revision 303)
+++ /issm/trunk/src/mex/Makefile.am	(revision 304)
@@ -6,4 +6,5 @@
 else
 bin_PROGRAMS =  BuildNodeSets\
+				ComputePressure\
 				ConfigureObjects \
 				ControlOptimization\
@@ -68,4 +69,7 @@
 BuildNodeSets_SOURCES = BuildNodeSets/BuildNodeSets.cpp\
 			  BuildNodeSets/BuildNodeSets.h
+
+ComputePressure_SOURCES = ComputePressure/ComputePressure.cpp\
+			  ComputePressure/ComputePressure.h
 
 ConfigureObjects_SOURCES = ConfigureObjects/ConfigureObjects.cpp\
Index: /issm/trunk/src/mex/ProcessParams/ProcessParams.cpp
===================================================================
--- /issm/trunk/src/mex/ProcessParams/ProcessParams.cpp	(revision 303)
+++ /issm/trunk/src/mex/ProcessParams/ProcessParams.cpp	(revision 304)
@@ -25,5 +25,5 @@
 
 	/*Shit partition: */
-	VecShift(partition,-1.0);
+	if(partition)VecShift(partition,-1.0);
 
 	/*!Generate internal degree of freedom numbers: */
Index: /issm/trunk/todo
===================================================================
--- /issm/trunk/todo	(revision 303)
+++ /issm/trunk/todo	(revision 304)
@@ -6,2 +6,3 @@
 Update control to get new inputs, as well as parallel diagnostic
 
+Finish diagnostic_core.m