Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5928)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5929)
@@ -753,13 +753,15 @@
 			pe=CreatePVectorDiagnosticHoriz();
 			if(pe) pe->AddToGlobal(pg,pf);
-			delete pe;
+			if(pe) delete pe;
 			break;
 		case AdjointHorizAnalysisEnum:
 			pe=CreatePVectorAdjointHoriz();
 			if(pe) pe->AddToGlobal(pg,pf);
-			delete pe;
+			if(pe) delete pe;
 			break;
 		case DiagnosticHutterAnalysisEnum:
-			CreatePVectorDiagnosticHutter( pg);
+			pe=CreatePVectorDiagnosticHutter();
+			if(pe) pe->AddToGlobal(pg,pf);
+			if(pe) delete pe;
 			break;
 		case DiagnosticVertAnalysisEnum:
@@ -767,5 +769,7 @@
 			break;
 		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
-			CreatePVectorSlope( pg);
+			pe=CreatePVectorSlope();
+			if(pe) pe->AddToGlobal(pg,pf);
+			if(pe) delete pe;
 			break;
 		case PrognosticAnalysisEnum:
@@ -3271,11 +3275,9 @@
 /*}}}*/
 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{1*/
-void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
+ElementVector* Penta::CreatePVectorDiagnosticHutter(void){
 
 	int i,j,k;
-	
+	int      ig;
 	const int numdofs=NDOF2*NUMVERTICES;
-	int*      doflist=NULL;
-	double    pe_g[numdofs]={0.0};
 	double    xyz_list[NUMVERTICES][3];
 	double    xyz_list_segment[2][3];
@@ -3285,7 +3287,5 @@
 	int       node0,node1;
 	GaussPenta* gauss=NULL;
-	int      ig;
 	double   slope[2];
-
 	double   z_g;
 	double   rho_ice,gravity,n,B;
@@ -3293,27 +3293,20 @@
 	double   ub,vb;
 	double   surface,thickness;
-
-	/*flags: */
 	int  connectivity[2];
 
-	/*If on water, skip: */
-	if(IsOnWater())return;
-
-	/*recover doflist: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*recover parameters: */
+	/*Initialize Element vector and return if necessary*/
+	if(IsOnWater()) return NULL;
+	ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	rho_ice=matpar->GetRhoIce();
 	gravity=matpar->GetG();
 	n=matice->GetN();
 	B=matice->GetB();
-
-	//recover extra inputs
 	Input* thickness_input=inputs->GetInput(ThicknessEnum);  ISSMASSERT(thickness_input);
 	Input* surface_input=inputs->GetInput(SurfaceEnum);      ISSMASSERT(surface_input);
 	Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); ISSMASSERT(slopex_input);
 	Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); ISSMASSERT(slopey_input);
-
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
@@ -3349,8 +3342,8 @@
 
 			if (IsOnSurface()){
-				for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1];
+				for(j=0;j<NDOF2;j++) pe->values[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1];
 			}
 			else{//connectivity is too large, should take only half on it
-				for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1];
+				for(j=0;j<NDOF2;j++) pe->values[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1];
 			}
 		}
@@ -3363,14 +3356,11 @@
 			vb=constant_part*slope[1];
 
-			pe_g[2*node0]+=ub/(double)connectivity[0];
-			pe_g[2*node0+1]+=vb/(double)connectivity[0];
-		}
-	}
-
-	/*Add pe_g to global vector pg: */
-	VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
-
-	/*Clean up */
-	xfree((void**)&doflist);
+			pe->values[2*node0]+=ub/(double)connectivity[0];
+			pe->values[2*node0+1]+=vb/(double)connectivity[0];
+		}
+	}
+
+	/*Clean up and return*/
+	return pe;
 }
 /*}}}*/
@@ -3773,7 +3763,7 @@
 /*}}}*/
 /*FUNCTION Penta::CreatePVectorSlope {{{1*/
-void Penta::CreatePVectorSlope( Vec pg){
-
-	if (!IsOnBed() || IsOnWater()) return;
+ElementVector* Penta::CreatePVectorSlope(void){
+
+	if (!IsOnBed() || IsOnWater()) return NULL;
 
 	/*Call Tria function*/
@@ -3781,9 +3771,7 @@
 	ElementVector* pe=tria->CreatePVectorSlope();
 	delete tria->matice; delete tria;
-	pe->AddToGlobal(pg,NULL);
-	delete pe;
 
 	/*clean up and return*/
-	return;
+	return pe;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5928)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5929)
@@ -159,5 +159,5 @@
 		ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
 		ElementVector* CreatePVectorDiagnosticHoriz(void);
-		void	  CreatePVectorDiagnosticHutter( Vec pg);
+		ElementVector* CreatePVectorDiagnosticHutter(void);
 		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
 		ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
@@ -171,5 +171,5 @@
 		void	  CreatePVectorMelting( Vec pg);
 		void	  CreatePVectorPrognostic( Vec pg);
-		void	  CreatePVectorSlope( Vec pg);
+		ElementVector* CreatePVectorSlope(void);
 		void	  CreatePVectorThermal( Vec pg);
 		void	  GetDofList(int** pdoflist,int approximation_enum,int setenum);