Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4923)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4924)
@@ -3318,10 +3318,47 @@
 void  Penta::CreatePVectorDiagnosticHutter(Vec pg){
 
-	/*Collapsed formulation: */
-	Beam*  beam=NULL;
-	int    i;
+	int i,j,k;
+	
+	const int numgrids=6;
+	const int NDOF2=2;
+	const int numdofs=NDOF2*numgrids;
+	int       doflist[numdofs];
+	int       numberofdofspernode;
+	double    pe_g[numdofs]={0.0};
+	double    xyz_list[numgrids][3];
+	double    z_list[numgrids];
+	double    z_segment[2];
+	double    slope2,constant_part;
+	double    gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	int       node0,node1;
+	BeamRef*  beam=NULL;
+
+	/*gaussian points: */
+	int      num_gauss;
+	double*  segment_gauss_coord=NULL;
+	double   gauss_coord;
+	double*  gauss_weights=NULL;
+	double   gauss_weight;
+	int      ig;
+	double   slope[2];
+
+	double   z_g;
+	double   rho_ice,gravity,n,B;
+	double   Jdet;
+	double   ub,vb;
+	double   surface,thickness;
+	double   slopex,slopey;
 
 	/*flags: */
 	bool onwater;
+	bool onbed;
+	bool onsurface;
+	int  connectivity[2];
+
+	/*Inputs*/
+	Input* thickness_input;
+	Input* surface_input;
+	Input* slopex_input;
+	Input* slopey_input;
 
 	/*recover some inputs: */
@@ -3331,13 +3368,98 @@
 	if(onwater)return;
 
-	/*Spawn 3 beam elements: */
+	/*recover doflist: */
+	GetDofList(&doflist[0],&numberofdofspernode);
+
+	/*recover some inputs: */
+	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+	/*recover parameters: */
+	rho_ice=matpar->GetRhoIce();
+	gravity=matpar->GetG();
+	n=matice->GetN();
+	B=matice->GetB();
+
+	//Initilaize beam
+	beam=new BeamRef();
+
+	//recover extra inputs
+	thickness_input=inputs->GetInput(ThicknessEnum);
+	surface_input=inputs->GetInput(SurfaceEnum);
+	slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
+	slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
+
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	for(i=0;i<numgrids;i++)z_list[i]=xyz_list[i][2];
+
+	//Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
+	num_gauss=3;
+	GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+
+	/*Loop on the three segments*/
 	for(i=0;i<3;i++){
-		beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta
-		beam->CreatePVector(pg);
-	}
-
-	/*clean up*/
+		slopex_input->GetParameterValue(&slopex, &gauss[i][0]);
+		slopey_input->GetParameterValue(&slopey, &gauss[i][0]);
+		surface_input->GetParameterValue(&surface, &gauss[i][0]);
+		thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
+
+		//compute slope2 slopex and slopey
+		slope2=pow(slopex,2)+pow(slopey,2);
+
+		//%compute constant_part
+		constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
+
+		z_segment[0]=z_list[i];
+		z_segment[1]=z_list[i+3];
+
+		node0=i;
+		node1=i+3;
+
+		connectivity[0]=nodes[node0]->GetConnectivity();
+		connectivity[1]=nodes[node1]->GetConnectivity();
+
+		/*Loop on the Gauss points: */
+		for(ig=0;ig<num_gauss;ig++){
+
+			//Pick up the gaussian point and its weight:
+			gauss_weight=gauss_weights[ig];
+			gauss_coord=segment_gauss_coord[ig];
+
+			beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
+
+			//Get Jacobian determinant:
+			beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
+
+			if (onsurface){
+				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];
+				}
+			}
+			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];
+				}
+			}
+		}
+
+		//Deal with lower surface
+		if (onbed){
+
+			//compute ub
+			constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
+			ub=constant_part*slopex;
+			vb=constant_part*slopey;
+
+			//Add to pe: 
+			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 */
 	delete beam;
-
 }
 /*}}}*/