Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 16808)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp	(revision 16809)
@@ -216,5 +216,109 @@
 }/*}}}*/
 ElementVector* MasstransportAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	if(!element->IsOnBed()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
+
+	switch(element->FiniteElement()){
+		case P1Enum: case P2Enum:
+			return CreatePVectorCG(basalelement);
+		case P1DGEnum:
+			return CreatePVectorDG(basalelement);
+		default:
+			_error_("Element type " << EnumToStringx(element->FiniteElement()) << " not supported yet");
+	}
+
+}/*}}}*/
+ElementVector* MasstransportAnalysis::CreatePVectorCG(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  Jdet,dt;
+	IssmDouble  ms,mb,mb_correction=0.,thickness;
+	IssmDouble* xyz_list;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* mb_correction_input = element->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+	Input* mb_input            = element->GetInput(BasalforcingsMeltingRateEnum);     _assert_(mb_input);
+	Input* ms_input            = element->GetInput(SurfaceforcingsMassBalanceEnum);   _assert_(ms_input);
+	Input* thickness_input     = element->GetInput(ThicknessEnum);                    _assert_(thickness_input);
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* 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);
+
+		ms_input->GetInputValue(&ms,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		thickness_input->GetInputValue(&thickness,gauss);
+		if(mb_correction_input)
+		 mb_correction_input->GetInputValue(&mb_correction,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* MasstransportAnalysis::CreatePVectorDG(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  Jdet,dt;
+	IssmDouble  ms,mb,mb_correction=0.,thickness;
+	IssmDouble* xyz_list;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* mb_correction_input = element->GetInput(BasalforcingsMeltingRateCorrectionEnum);
+	Input* mb_input            = element->GetInput(BasalforcingsMeltingRateEnum);     _assert_(mb_input);
+	Input* ms_input            = element->GetInput(SurfaceforcingsMassBalanceEnum);   _assert_(ms_input);
+	Input* thickness_input     = element->GetInput(ThicknessEnum);                    _assert_(thickness_input);
+
+	/*Initialize mb_correction to 0, do not forget!:*/
+	/* 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);
+
+		ms_input->GetInputValue(&ms,gauss);
+		mb_input->GetInputValue(&mb,gauss);
+		thickness_input->GetInputValue(&thickness,gauss);
+		if(mb_correction_input)
+		 mb_correction_input->GetInputValue(&mb_correction,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 16808)
+++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h	(revision 16809)
@@ -23,4 +23,6 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorCG(Element* element);
+		ElementVector* CreatePVectorDG(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16808)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16809)
@@ -1975,5 +1975,4 @@
 bool Tria::IsOnBed(){
 
-	bool onbed;
 	int meshtype;
 	this->parameters->FindParam(&meshtype,MeshTypeEnum);
@@ -1981,11 +1980,8 @@
 		case Mesh2DverticalEnum:
 			return HasEdgeOnBed();
-			break;
 		case Mesh2DhorizontalEnum:
-			inputs->GetInputValue(&onbed,MeshElementonbedEnum);
-			break;
+			return true;
 		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
 	}
-	return onbed;
 }
 /*}}}*/