Index: issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16782)
+++ issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16803)
@@ -100,5 +100,130 @@
 }/*}}}*/
 ElementVector* StressbalanceVerticalAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorVolume(element);
+	ElementVector* pe2=CreatePVectorBase(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceVerticalAnalysis::CreatePVectorVolume(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         approximation;
+	IssmDouble  Jdet;
+	IssmDouble  dudx,dvdy,dwdz;
+	IssmDouble  du[3],dv[3],dw[3];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and basis functions*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetInputValue(&approximation,ApproximationEnum);
+	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
+	Input* vzFS_input=NULL;
+	if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+		vzFS_input=element->GetInput(VzFSEnum); _assert_(vzFS_input);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		vx_input->GetInputDerivativeValue(&du[0],xyz_list,gauss);
+		vy_input->GetInputDerivativeValue(&dv[0],xyz_list,gauss);
+		if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+			vzFS_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
+			dwdz=dw[2];
+		}
+		else dwdz=0;
+		dudx=du[0];
+		dvdy=dv[1];
+
+		for(int i=0;i<numnodes;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceVerticalAnalysis::CreatePVectorBase(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         i,j,approximation;
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_base = NULL;
+	IssmDouble  Jdet;
+	IssmDouble  vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
+	IssmDouble  slope[3];
+
+	if(!element->IsOnBed()) return NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->GetInputValue(&approximation,ApproximationEnum);
+	Input* bed_input=element->GetInput(BedEnum);                                _assert_(bed_input);
+	Input* basal_melting_input=element->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
+	Input* vx_input=element->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);                                  _assert_(vy_input);
+	Input* vzFS_input=NULL;
+	if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+		vzFS_input=element->GetInput(VzFSEnum);       _assert_(vzFS_input);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basal_melting_input->GetInputValue(&basalmeltingvalue,gauss);
+		bed_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+			vzFS_input->GetInputValue(&vz,gauss);
+		}
+		else vz=0.;
+
+		dbdx=slope[0];
+		dbdy=slope[1];
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		for(i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
+
 }/*}}}*/
 void StressbalanceVerticalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
