Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16802)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16803)
@@ -818,10 +818,28 @@
 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
 
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
 	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorSSADrivingStress(element);
-	ElementVector* pe2=CreatePVectorSSAFront(element);
+	ElementVector* pe1=CreatePVectorSSADrivingStress(basalelement);
+	ElementVector* pe2=CreatePVectorSSAFront(basalelement);
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
 	/*clean-up and return*/
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	delete pe1;
 	delete pe2;
@@ -834,6 +852,6 @@
 	IssmDouble  driving_stress_baseline,thickness;
 	IssmDouble  Jdet;
-	IssmDouble* xyz_list=NULL;
 	IssmDouble  slope[2];
+	IssmDouble* xyz_list     = NULL;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -881,4 +899,5 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
+
 	return NULL;
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp	(revision 16802)
+++ /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){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h	(revision 16802)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.h	(revision 16803)
@@ -23,4 +23,6 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorVolume(Element* element);
+		ElementVector* CreatePVectorBase(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16802)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16803)
@@ -67,4 +67,5 @@
 		virtual void	GetDofListPressure(int** pdoflist,int setenum)=0;
 		virtual void   JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
+		virtual void   JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
 		virtual void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum)=0;
 		virtual int    GetNodeIndex(Node* node)=0;
@@ -90,4 +91,5 @@
 		virtual void   GetInputValue(IssmDouble* pvalue,int enum_type)=0;
 		virtual void   GetVerticesCoordinates(IssmDouble** xyz_list)=0;
+		virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
 		virtual void   GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
 
@@ -110,5 +112,6 @@
 		virtual Gauss* NewGauss(void)=0;
 		virtual Gauss* NewGauss(int order)=0;
-		virtual ElementVector*  NewElementVector(int approximation_enum)=0;
+		virtual Gauss* NewGaussBase(int order)=0;
+		virtual ElementVector*  NewElementVector(int approximation_enum=NoneApproximationEnum)=0;
 		virtual void   InputScale(int enum_type,IssmDouble scale_factor)=0;
 		virtual void   GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16802)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16803)
@@ -1446,4 +1446,14 @@
 
 }/*}}}*/
+/*FUNCTION Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/
+void Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){
+
+	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
+	::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
+
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list;
+
+}/*}}}*/
 /*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
 void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
@@ -2432,4 +2442,12 @@
 }
 /*}}}*/
+/*FUNCTION Penta::JacobianDeterminantBase{{{*/
+void Penta::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetTriaJacobianDeterminant(pJdet,xyz_list_base,(GaussPenta*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Penta::NoIceInElement {{{*/
 bool   Penta::NoIceInElement(){
@@ -2508,4 +2526,9 @@
 Gauss* Penta::NewGauss(int order){
 	return new GaussPenta(order,order);
+}
+/*}}}*/
+/*FUNCTION Penta::NewGaussBase(int order){{{*/
+Gauss* Penta::NewGaussBase(int order){
+	return new GaussPenta(0,1,2,order);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16802)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16803)
@@ -99,4 +99,5 @@
 		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int name_enum);
 		void   GetVerticesCoordinates(IssmDouble** pxyz_list);
+		void   GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
 
 		int    Sid();
@@ -239,7 +240,9 @@
 		bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
 		void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
+		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
 		bool           NoIceInElement(void); 
 		Gauss*         NewGauss(void);
 		Gauss*         NewGauss(int order);
+		Gauss*         NewGaussBase(int order);
 		ElementVector* NewElementVector(int approximation_enum);
 		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16802)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16803)
@@ -102,4 +102,5 @@
 		int         GetNumberOfVertices(void){_error_("not implemented yet");};
 		void        GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");};
+		void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
 		int         Sid(){_error_("not implemented yet");};
 		void        InputChangeName(int input_enum, int enum_type_old){_error_("not implemented yet");};
@@ -108,4 +109,5 @@
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
 		void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
+		void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
 		void        NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		bool        NoIceInElement(){_error_("not implemented yet");};
@@ -125,4 +127,5 @@
 		Gauss*      NewGauss(void){_error_("not implemented yet");};
 		Gauss*      NewGauss(int order){_error_("not implemented yet");};
+		Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
 		ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
 		#ifdef _HAVE_THERMAL_
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16802)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16803)
@@ -110,4 +110,5 @@
 		void        GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);
 		void        GetVerticesCoordinates(IssmDouble** pxyz_list);
+		void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");};
 		void        InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code);
 		void        InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=MeshElementsEnum);
@@ -267,6 +268,8 @@
 		bool	         IsInput(int name);
 		void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
+		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
 		Gauss*         NewGauss(void);
 		Gauss*         NewGauss(int order);
+		Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
 		ElementVector* NewElementVector(int approximation_enum);
 		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);