Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16799)
@@ -806,5 +806,80 @@
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	switch(approximation){
+		case SSAApproximationEnum: 
+			return CreatePVectorSSA(element);
+		default:
+			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
+	}
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorSSADrivingStress(element);
+	ElementVector* pe2=CreatePVectorSSAFront(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         i;
+	IssmDouble  driving_stress_baseline,thickness;
+	IssmDouble  Jdet;
+	IssmDouble* xyz_list=NULL;
+	IssmDouble  slope[2];
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	IssmDouble* icefrontlevel= xNew<IssmDouble>(numnodes);
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe=element->NewElementVector(SSAApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* thickness_input=element->GetInput(ThicknessEnum);          _assert_(thickness_input); 
+	Input* surface_input  =element->GetInput(SurfaceEnum);            _assert_(surface_input);
+	Input* drag_input     =element->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
+	element->GetInputListOnVertices(icefrontlevel,BedEnum);
+
+	/* 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);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		driving_stress_baseline=element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum)*thickness;
+
+		/*Build load vector: */
+		for(i=0;i<numnodes;i++){
+			pe->values[i*NDOF2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*NDOF2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
+
 }/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -952,5 +1027,5 @@
 			return;
 		case L1L2ApproximationEnum:
-			InputUpdateFromSolutionSSA(solution,element);
+			InputUpdateFromSolutionL1L2(solution,element);
 			return;
 		case SSAHOApproximationEnum:
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16799)
@@ -23,4 +23,7 @@
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorSSA(Element* element);
+		ElementVector* CreatePVectorSSADrivingStress(Element* element);
+		ElementVector* CreatePVectorSSAFront(Element* element);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16799)
@@ -24,4 +24,5 @@
 class Input;
 class Gauss;
+class ElementVector;
 template <class doublematrix> class Matrix;
 template <class doubletype> class Vector;
@@ -54,4 +55,9 @@
 		virtual int    FiniteElement(void)=0;
 		virtual IssmDouble GetMaterialParameter(int enum_in)=0;
+		virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
+		virtual void   TransformLoadVectorCoord(ElementVector* pe,int transformenum)=0;
+		virtual void   TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list)=0;
+		virtual void   TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum)=0;
+		virtual void   TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list)=0;
 		virtual void   TransformSolutionCoord(IssmDouble* values,int transformenum)=0;
 		virtual void   TransformSolutionCoord(IssmDouble* values,int* transformenum_list)=0;
@@ -62,4 +68,5 @@
 		virtual void	GetDofListVelocity(int** pdoflist,int setenum)=0;
 		virtual void	GetDofListPressure(int** pdoflist,int setenum)=0;
+		virtual void   JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
 		virtual void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum)=0;
 		virtual int    GetNodeIndex(Node* node)=0;
@@ -104,4 +111,6 @@
 		virtual int    NumberofNodesPressure(void)=0;
 		virtual Gauss* NewGauss(void)=0;
+		virtual Gauss* NewGauss(int order)=0;
+		virtual ElementVector*  NewElementVector(int approximation_enum)=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 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16799)
@@ -2519,4 +2519,12 @@
 }
 /*}}}*/
+/*FUNCTION Penta::JacobianDeterminant{{{*/
+void Penta::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetJacobianDeterminant(pJdet,xyz_list,(GaussPenta*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Penta::NoIceInElement {{{*/
 bool   Penta::NoIceInElement(){
@@ -2587,7 +2595,12 @@
 }
 /*}}}*/
-/*FUNCTION Penta::NewGauss{{{*/
+/*FUNCTION Penta::NewGauss(){{{*/
 Gauss* Penta::NewGauss(void){
 	return new GaussPenta();
+}
+/*}}}*/
+/*FUNCTION Penta::NewGauss(int order){{{*/
+Gauss* Penta::NewGauss(int order){
+	return new GaussPenta(order,order);
 }
 /*}}}*/
@@ -8021,5 +8034,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+	::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
@@ -8099,5 +8112,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
 
 	/*Clean up and return*/
@@ -8183,5 +8196,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+	::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
@@ -8264,5 +8277,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+	::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
@@ -8542,5 +8555,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
+	::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
@@ -8625,5 +8638,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
+	::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
@@ -8733,5 +8746,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
 
 	/*Clean up and return*/
@@ -8887,5 +8900,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+	::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
@@ -8972,5 +8985,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
 
 	/*Clean up and return*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16799)
@@ -240,10 +240,18 @@
 		bool           IsFloating(void); 
 		bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
+		void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		bool           NoIceInElement(void); 
 		Gauss*         NewGauss(void);
+		Gauss*         NewGauss(int order);
+		ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
+		void           NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		IssmDouble     MinEdgeLength(IssmDouble xyz_list[6][3]);
 		void	         SetClone(int* minranks);
 		Tria*	         SpawnTria(int location);
 		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
+		void           TransformLoadVectorCoord(ElementVector* pe,int transformenum){_error_("not implemented yet");};
+		void           TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list){_error_("not implemented yet");};
+		void           TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};      /*Tiling only*/
+		void           TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/
 		void           TransformSolutionCoord(IssmDouble* values,int transformenum);
 		void           TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16799)
@@ -109,4 +109,6 @@
 		bool        IsFloating(){_error_("not implemented yet");};
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
+		void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
+		void        NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		bool        NoIceInElement(){_error_("not implemented yet");};
 		int         NumberofNodesVelocity(void){_error_("not implemented yet");};
@@ -124,4 +126,6 @@
 		void        GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
 		Gauss*      NewGauss(void){_error_("not implemented yet");};
+		Gauss*      NewGauss(int order){_error_("not implemented yet");};
+		ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");};
 		#ifdef _HAVE_THERMAL_
 		void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
@@ -155,4 +159,8 @@
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
 		IssmDouble  TimeAdapt(){_error_("not implemented yet");};
+		void        TransformLoadVectorCoord(ElementVector* pe,int transformenum){_error_("not implemented yet");};
+		void        TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list){_error_("not implemented yet");};
+		void        TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};
+		void        TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};
 		void        TransformSolutionCoord(IssmDouble* values,int transformenum){_error_("not implemented yet");};
 		void        TransformSolutionCoord(IssmDouble* values,int* transformenum_list){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16799)
@@ -2044,4 +2044,12 @@
 }
 /*}}}*/
+/*FUNCTION Tria::JacobianDeterminant{{{*/
+void Tria::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	this->JacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::HasEdgeOnBed {{{*/
 bool Tria::HasEdgeOnBed(){
@@ -2229,7 +2237,17 @@
 }
 /*}}}*/
-/*FUNCTION Tria::NewGauss{{{*/
+/*FUNCTION Tria::NewGauss(){{{*/
 Gauss* Tria::NewGauss(void){
 	return new GaussTria();
+}
+/*}}}*/
+/*FUNCTION Tria::NewGauss(int order){{{*/
+Gauss* Tria::NewGauss(int order){
+	return new GaussTria(order);
+}
+/*}}}*/
+/*FUNCTION Tria::NewElementVector{{{*/
+ElementVector* Tria::NewElementVector(int approximation_enum){
+	return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum);
 }
 /*}}}*/
@@ -2246,4 +2264,12 @@
 	/*If the level set is awlays <=0, there is no ice here*/
 	return true;
+}
+/*}}}*/
+/*FUNCTION Tria::NodalFunctions{{{*/
+void Tria::NodalFunctions(IssmDouble* basis, Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss);
+
 }
 /*}}}*/
@@ -2731,4 +2757,18 @@
 
 	::TransformSolutionCoord(values,this->nodes,this->NumberofNodes(),transformenum_list);
+
+}
+/*}}}*/
+/*FUNCTION Tria::TransformLoadVectorCoord{{{*/
+void Tria::TransformLoadVectorCoord(ElementVector* pe,int transformenum){
+
+	::TransformLoadVectorCoord(pe,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
+/*FUNCTION Tria::TransformLoadVectorCoord{{{*/
+void Tria::TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list){
+
+	::TransformLoadVectorCoord(pe,this->nodes,this->NumberofNodes(),transformenum_list);
 
 }
@@ -3911,5 +3951,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
 
 	/*Clean up and return*/
@@ -3972,5 +4012,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+	::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
@@ -4042,5 +4082,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+	::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
 
 	/*Clean up and return*/
@@ -4111,5 +4151,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
+	::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
@@ -4191,5 +4231,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
+	::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
 
 	/*Clean up and return*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16799)
@@ -268,8 +268,16 @@
 		void	         InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
 		bool	         IsInput(int name);
+		void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
 		Gauss*         NewGauss(void);
+		Gauss*         NewGauss(int order);
+		ElementVector* NewElementVector(int approximation_enum);
+		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void	         SetClone(int* minranks);
 		Seg*	         SpawnSeg(int index1,int index2);
 		void	         SurfaceNormal(IssmDouble* surface_normal, IssmDouble xyz_list[3][3]);
+		void           TransformLoadVectorCoord(ElementVector* pe,int transformenum);
+		void           TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list);
+		void           TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};
+		void           TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};
 		void           TransformSolutionCoord(IssmDouble* values,int transformenum);
 		void           TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.cpp	(revision 16799)
@@ -115,10 +115,4 @@
 void BoolInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){_error_("not supported yet!");}
 /*}}}*/
-/*FUNCTION BoolInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void BoolInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not supported yet!");}
-/*}}}*/
-/*FUNCTION BoolInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
-void BoolInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not supported yet!");}
-/*}}}*/
 /*FUNCTION BoolInput::ChangeEnum{{{*/
 void BoolInput::ChangeEnum(int newenumtype){
Index: /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/BoolInput.h	(revision 16799)
@@ -50,7 +50,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){_error_("not implemented yet");};
 		void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");};
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.cpp	(revision 16799)
@@ -247,10 +247,6 @@
 	gradient->GetInputValue(pvalue,gauss);
 }/*}}}*/
-/*FUNCTION ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){
-	values->GetInputDerivativeValue(derivativevalues,xyz_list,gauss);
-}/*}}}*/
-/*FUNCTION ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
-void ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){{{*/
+void ControlInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss){
 	values->GetInputDerivativeValue(derivativevalues,xyz_list,gauss);
 }/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h	(revision 16799)
@@ -56,7 +56,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss);
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.h	(revision 16799)
@@ -51,7 +51,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");};
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.cpp	(revision 16799)
@@ -123,10 +123,4 @@
 void DoubleInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){*pvalue=this->value;}
 /*}}}*/
-/*FUNCTION DoubleInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void DoubleInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not supported yet!");}
-/*}}}*/
-/*FUNCTION DoubleInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
-void DoubleInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not supported yet!");}
-/*}}}*/
 /*FUNCTION DoubleInput::GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){{{*/
 void DoubleInput::GetVxStrainRate2d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussTria* gauss){
Index: /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/DoubleInput.h	(revision 16799)
@@ -49,7 +49,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/Input.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/Input.h	(revision 16799)
@@ -32,7 +32,5 @@
 		virtual void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time)=0;
 		virtual void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int index)=0;
-		virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss)=0;
-		virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss)=0;
-		virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss)=0;
+		virtual void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, Gauss* gauss)=0;
 		virtual void GetInputAverage(IssmDouble* pvalue)=0;
 		virtual void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes)=0;
Index: /issm/trunk-jpl/src/c/classes/Inputs/IntInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/IntInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/IntInput.cpp	(revision 16799)
@@ -115,10 +115,4 @@
 void IntInput::GetInputValue(IssmDouble* pvalue,Gauss* gauss){_error_("not supported yet!");}
 /*}}}*/
-/*FUNCTION IntInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void IntInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not supported yet!");}
-/*}}}*/
-/*FUNCTION IntInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
-void IntInput::GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not supported yet!");}
-/*}}}*/
 /*FUNCTION IntInput::ChangeEnum{{{*/
 void IntInput::ChangeEnum(int newenumtype){
Index: /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/IntInput.h	(revision 16799)
@@ -50,7 +50,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void GetInputAverage(IssmDouble* pvalue){_error_("not implemented yet");};
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 16799)
@@ -146,9 +146,10 @@
 }
 /*}}}*/
-/*FUNCTION PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
-void PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){{{*/
+void PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list,Gauss* gauss){
 
 	/*Call PentaRef function*/
-	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	PentaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussPenta*)gauss);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h	(revision 16799)
@@ -49,7 +49,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss);
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.cpp	(revision 16799)
@@ -105,9 +105,10 @@
 }
 /*}}}*/
-/*FUNCTION SegInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussSeg* gauss){{{*/
-void SegInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussSeg* gauss){
+/*FUNCTION SegInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){{{*/
+void SegInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list,Gauss* gauss){
 
 	/*Call SegRef function*/
-	SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
+	_assert_(gauss->Enum()==GaussSegEnum);
+	SegRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussSeg*)gauss);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/SegInput.h	(revision 16799)
@@ -48,10 +48,7 @@
 		void GetInputValue(IssmDouble* pvalue){_error_("not implemented yet");}
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss);
-		void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss);
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp	(revision 16799)
@@ -192,6 +192,6 @@
 }
 /*}}}*/
-/*FUNCTION TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){
+/*FUNCTION TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){{{*/
+void TransientInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){
 
 	IssmDouble time;
@@ -207,5 +207,4 @@
 
 	delete input;
-
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TransientInput.h	(revision 16799)
@@ -55,7 +55,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time);
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss ,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss);
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes);
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp	(revision 16799)
@@ -133,9 +133,10 @@
 }
 /*}}}*/
-/*FUNCTION TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){{{*/
-void TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussTria* gauss){
+/*FUNCTION TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){{{*/
+void TriaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, Gauss* gauss){
 
 	/*Call TriaRef function*/
-	TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,gauss);
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	TriaRef::GetInputDerivativeValue(p,&values[0],xyz_list,(GaussTria*)gauss);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h	(revision 16799)
@@ -50,7 +50,5 @@
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,IssmDouble time){_error_("not implemented yet");};
 		void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int index){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussSeg* gauss){_error_("not implemented yet");};
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list, GaussPenta* gauss){_error_("not implemented yet");};
+		void GetInputDerivativeValue(IssmDouble* derivativevalues, IssmDouble* xyz_list,Gauss* gauss);
 		void GetInputAverage(IssmDouble* pvalue);
 		void GetInputAllTimeAverages(IssmDouble** pvalues,IssmDouble** ptimes, int* pnumtimes);
Index: /issm/trunk-jpl/src/c/classes/gauss/Gauss.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/Gauss.h	(revision 16798)
+++ /issm/trunk-jpl/src/c/classes/gauss/Gauss.h	(revision 16799)
@@ -9,6 +9,7 @@
 
 	public: 
+		IssmDouble   weight;
+
 		virtual        ~Gauss(){};
-
 		virtual int  begin(void)=0;
 		virtual int  end(void)=0;