Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15730)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15731)
@@ -270,5 +270,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element matrix and vectors*/
@@ -289,7 +288,7 @@
 		D=gauss->weight*Jdet;
 
-		TripleMultiply(basis,1,numdof,1,
+		TripleMultiply(basis,1,numnodes,1,
 					&D,1,1,0,
-					basis,1,numdof,0,
+					basis,1,numnodes,0,
 					&Ke->values[0],1);
 	}
@@ -1583,5 +1582,5 @@
 			break;
 		case AdjointBalancethicknessAnalysisEnum:
-			InputUpdateFromSolutionAdjointBalancethickness(solution);
+			InputUpdateFromSolutionOneDof(solution,AdjointEnum);
 			break;
 		#endif
@@ -1625,13 +1624,16 @@
 void  Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
 
-	const int  numdof         = NDOF1*NUMVERTICES;
-	int*       doflist        = NULL;
-	IssmDouble values[numdof];
-
-	/*Get dof list: */
+	/*Intermediary*/
+	int* doflist = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numdof;i++){
+	for(int i=0;i<numnodes;i++){
 		values[i]=solution[doflist[i]];
 		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
@@ -1642,4 +1644,5 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
 	xDelete<int>(doflist);
 }
@@ -1649,22 +1652,23 @@
 
 	/*Intermediaries*/
-	const int numdof = NDOF1*NUMVERTICES;
-
-	int       i,hydroadjustment;
-	int*      doflist=NULL;
-	IssmDouble    rho_ice,rho_water,minthickness;
-	IssmDouble    newthickness[numdof];
-	IssmDouble    newbed[numdof];
-	IssmDouble    newsurface[numdof];
-	IssmDouble    oldbed[NUMVERTICES];
-	IssmDouble    oldsurface[NUMVERTICES];
-	IssmDouble    oldthickness[NUMVERTICES];
-
-	/*Get dof list: */
+	int        i,hydroadjustment;
+	int*       doflist=NULL;
+	IssmDouble rho_ice,rho_water,minthickness;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
+	IssmDouble* newbed       = xNew<IssmDouble>(numnodes);
+	IssmDouble* newsurface   = xNew<IssmDouble>(numnodes);
+	IssmDouble* oldthickness = xNew<IssmDouble>(numnodes);
+	IssmDouble* oldbed       = xNew<IssmDouble>(numnodes);
+	IssmDouble* oldsurface   = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
 	this->parameters->FindParam(&minthickness,PrognosticMinThicknessEnum);
-	for(i=0;i<numdof;i++){
+	for(i=0;i<numnodes;i++){
 		newthickness[i]=solution[doflist[i]];
 		if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
@@ -1674,7 +1678,7 @@
 
 	/*Get previous bed, thickness and surface*/
-	GetInputListOnVertices(&oldbed[0],BedEnum);
-	GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
-	GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
+	GetInputListOnNodes(&oldbed[0],BedEnum);
+	GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
+	GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
 
 	/*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
@@ -1683,19 +1687,18 @@
 	rho_water=matpar->GetRhoWater();
 
-	for(i=0;i<numdof;i++) {
+	for(i=0;i<numnodes;i++) {
 		/*If shelf: hydrostatic equilibrium*/
 		if (this->nodes[i]->IsGrounded()){
-			newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
-			newbed[i]=oldbed[i];               //same bed: do nothing
+			newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
+			newbed[i]     = oldbed[i];                 //same bed: do nothing
 		}
 		else{ //this is an ice shelf
-
 			if(hydroadjustment==AbsoluteEnum){
-				newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
-				newbed[i]=newthickness[i]*(-rho_ice/rho_water);
+				newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
+				newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
 			}
 			else if(hydroadjustment==IncrementalEnum){
-				newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH 
-				newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
+				newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
+				newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
 			}
 			else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
@@ -1709,4 +1712,10 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(newthickness);
+	xDelete<IssmDouble>(newbed);
+	xDelete<IssmDouble>(newsurface);
+	xDelete<IssmDouble>(oldthickness);
+	xDelete<IssmDouble>(oldbed);
+	xDelete<IssmDouble>(oldsurface);
 	xDelete<int>(doflist);
 }
@@ -1715,7 +1724,7 @@
 void  Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type){
 
-	const int   numdof         = NDOF1 *NUMVERTICES;
-	int        *doflist        = NULL;
-	IssmDouble  values[numdof];
+	const int   numnodes = NUMVERTICES;
+	int        *doflist  = NULL;
+	IssmDouble  values[numnodes];
 
 	/*Check that name is an element input*/
@@ -1750,9 +1759,10 @@
 
 	case NodesEnum:
+		_assert_(this->element_type==P1Enum);
 		/*Get dof list: */
 		GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 
 		/*Use the dof list to index into the vector: */
-		for(int i=0;i<numdof;i++){
+		for(int i=0;i<numnodes;i++){
 			values[i]=vector[doflist[i]];
 			if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
@@ -1766,4 +1776,5 @@
 
 	case NodeSIdEnum:
+		_assert_(this->element_type==P1Enum);
 		for(int i=0;i<NUMVERTICES;i++){
 			values[i]=vector[nodes[i]->Sid()];
@@ -1772,7 +1783,4 @@
 		/*Add input to the element: */
 		this->inputs->AddInput(new TriaInput(name,values,P1Enum));
-
-		/*Free ressources:*/
-		xDelete<int>(doflist);
 		return;
 
@@ -2567,7 +2575,6 @@
 IssmDouble Tria::MassFlux( IssmDouble* segment){
 
-	const int  numdofs=2;
 	int        dim;
-	IssmDouble mass_flux=0;
+	IssmDouble mass_flux=0.;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble normal[2];
@@ -5019,5 +5026,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF2;
 
 	/*Initialize Element vector*/
@@ -5255,8 +5261,4 @@
 ElementMatrix* Tria::CreateKMatrixAdjointSSA(void){
 
-	/*Constants*/
-	const int numnodes = this->NumberofNodes();
-	const int numdof   = NDOF2*numnodes;
-
 	/*Intermediaries */
 	int         i,j;
@@ -5271,4 +5273,7 @@
 	GaussTria  *gauss=NULL;
 
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
 	/*Initialize Jacobian with regular SSA (first part of the Gateau derivative)*/
 	parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
@@ -5328,20 +5333,25 @@
 void  Tria::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){
 
-	const int numdof=NDOF2*NUMVERTICES;
-
-	int       i;
-	int*      doflist=NULL;
-	IssmDouble    values[numdof];
-	IssmDouble    lambdax[NUMVERTICES];
-	IssmDouble    lambday[NUMVERTICES];
-
-	/*Get dof list: */
+	int  i;
+	int* doflist=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
+
+	/*Fetch dof list and allocate solution vectors*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values  = xNew<IssmDouble>(numdof);
+	IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
+	IssmDouble* lambday = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
+	/*Transform solution in Cartesian Space*/
+	TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
+
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<NUMVERTICES;i++){
+	for(i=0;i<numdof;i++){
 		lambdax[i]=values[i*NDOF2+0];
 		lambday[i]=values[i*NDOF2+1];
@@ -5357,33 +5367,7 @@
 
 	/*Free ressources:*/
-	xDelete<int>(doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{*/
-void  Tria::InputUpdateFromSolutionAdjointBalancethickness(IssmDouble* solution){
-
-	const int numdof=NDOF1*NUMVERTICES;
-
-	int       i;
-	int*      doflist=NULL;
-	IssmDouble    values[numdof];
-	IssmDouble    lambda[NUMVERTICES];
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numdof;i++){
-		lambda[i]=values[i];
-		if(xIsNan<IssmDouble>(lambda[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Add vx and vy as inputs to the tria element: */
-	this->inputs->AddInput(new TriaInput(AdjointEnum,lambda,P1Enum));
-
-	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(lambdax);
+	xDelete<IssmDouble>(lambday);
 	xDelete<int>(doflist);
 }
@@ -5464,11 +5448,8 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element vector and vectors*/
 	ElementMatrix* Ke=new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
 	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Initialize Element matrix*/
 
 	/*Retrieve all inputs and parameters*/
@@ -5488,7 +5469,7 @@
 		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
 
-		TripleMultiply(basis,numdof,1,0,
+		TripleMultiply(basis,numnodes,1,0,
 					&D_scalar,1,1,0,
-					basis,1,numdof,0,
+					basis,1,numnodes,0,
 					&Ke->values[0],1);
 	}
@@ -5691,9 +5672,8 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
-	IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    B      = xNew<IssmDouble>(5*numnodes);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble     D[2][2];
@@ -5719,7 +5699,7 @@
 		D[1][0]=0.;       D[1][1]=D_scalar;
 		GetBHydro(B,&xyz_list[0][0],gauss); 
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					B,2,numdof,0,
+					B,2,numnodes,0,
 					&Ke->values[0],1);
 
@@ -5729,7 +5709,7 @@
 			D_scalar=sediment_storing*gauss->weight*Jdet;
 
-			TripleMultiply(basis,numdof,1,0,
+			TripleMultiply(basis,numnodes,1,0,
 						&D_scalar,1,1,0,
-						basis,1,numdof,0,
+						basis,1,numnodes,0,
 						&Ke->values[0],1);
 		}
@@ -5758,9 +5738,8 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element matrix and vectors*/
 	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
-	IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    B      = xNew<IssmDouble>(5*numnodes);
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble     D[2][2];
@@ -5785,7 +5764,7 @@
 		D[1][0]=0.;       D[1][1]=D_scalar;
 		GetBHydro(B,&xyz_list[0][0],gauss); 
-		TripleMultiply(B,2,numdof,1,
+		TripleMultiply(B,2,numnodes,1,
 					&D[0][0],2,2,0,
-					B,2,numdof,0,
+					B,2,numnodes,0,
 					&Ke->values[0],1);
 
@@ -5795,7 +5774,7 @@
 			D_scalar=epl_storing*gauss->weight*Jdet;
 
-			TripleMultiply(basis,numdof,1,0,
+			TripleMultiply(basis,numnodes,1,0,
 						&D_scalar,1,1,0,
-						basis,1,numdof,0,
+						basis,1,numnodes,0,
 						&Ke->values[0],1);
 		}
@@ -5824,5 +5803,4 @@
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = this->NumberofNodes();
-	int numdof   = numnodes*NDOF1;
 
 	/*Initialize Element vector*/
@@ -5850,8 +5828,8 @@
 
 		if(reCast<int,IssmDouble>(dt)){
-			for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
+			for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
 		}
 		else{
-			for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
+			for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
 		}
 	}
@@ -5999,14 +5977,13 @@
 void  Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){
 
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	int         i;
 	int        *doflist = NULL;
-	IssmDouble  enum_value;
-	IssmDouble  values[numdof];
-	GaussTria  *gauss   = NULL;
-
-	/*Get dof list: */
+	IssmDouble  value;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values = xNew<IssmDouble>(numnodes);
 
 	/*Get inputs*/
@@ -6014,20 +5991,18 @@
 
 	/*Ok, we have the values, fill in the array: */
-	/*P1 element only for now*/
-	gauss=new GaussTria();
-	for(i=0;i<NUMVERTICES;i++){
-
-		gauss->GaussVertex(i);
-
-		/*Recover dof values*/
-		enum_input->GetInputValue(&enum_value,gauss);
-		values[i]=enum_value;
-	}
-
-	solution->SetValues(numdof,doflist,values,INS_VAL);
+	GaussTria* gauss=new GaussTria();
+	for(int i=0;i<numnodes;i++){
+		gauss->GaussNode(this->element_type,i);
+
+		enum_input->GetInputValue(&value,gauss);
+		values[i]=value;
+	}
+
+	solution->SetValues(numnodes,doflist,values,INS_VAL);
 
 	/*Free ressources:*/
+	xDelete<int>(doflist);
+	xDelete<IssmDouble>(values);
 	delete gauss;
-	xDelete<int>(doflist);
 }
 /*}}}*/
@@ -6035,14 +6010,16 @@
 void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){
 
-	/*Intermediaries*/
-	const int   numdof         = NDOF1 *NUMVERTICES;
-	int        *doflist        = NULL;
-	IssmDouble  values[numdof];
-
-	/*Get dof list: */
+	/*Intermediary*/
+	int* doflist = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numdof;i++){
+	for(int i=0;i<numnodes;i++){
 		values[i]=solution[doflist[i]];
 		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
@@ -6054,4 +6031,5 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
 	xDelete<int>(doflist);
 }
@@ -6061,17 +6039,19 @@
 
 	/*Intermediaries*/
-	const int   numdof   = NDOF1 *NUMVERTICES;
 	int        *doflist  = NULL;
 	bool        converged;
-	IssmDouble  values[numdof];
-	IssmDouble  residual[numdof];
 	IssmDouble  penalty_factor, dt;
 	IssmDouble  kmax, kappa, h_max;
 
-	/*Get dof list: */
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Fetch dof list and allocate solution vector*/
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values   = xNew<IssmDouble>(numnodes);
+	IssmDouble* residual = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
-	for(int i=0;i<numdof;i++){
+	for(int i=0;i<numnodes;i++){
 		values[i]=solution[doflist[i]];
 		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
@@ -6088,5 +6068,5 @@
 		kappa=kmax*pow(10.,penalty_factor);
 
-		for(int i=0;i<NUMVERTICES;i++){
+		for(int i=0;i<numnodes;i++){
 			this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
 			if(values[i]>h_max){
@@ -6104,4 +6084,6 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(residual);
 	xDelete<int>(doflist);
 }
@@ -6145,6 +6127,6 @@
 void  Tria::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
 
-	const int  numdof         = NDOF1 *NUMVERTICES;
-	int        *doflist       = NULL;
+	const int  numdof   = NDOF1 *NUMVERTICES;
+	int        *doflist = NULL;
 	bool       isefficientlayer;
 	int        transfermethod;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15730)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15731)
@@ -242,6 +242,5 @@
 		ElementVector* CreatePVectorAdjointHoriz(void);
 		ElementVector* CreatePVectorAdjointBalancethickness(void);
-		void	  InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);
-		void	  InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution);
+		void	         InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution);
 		#endif