Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16018)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16019)
@@ -5587,5 +5587,5 @@
 	/*Intermediaries */
 	IssmDouble diffusivity;
-	IssmDouble Jdettria,D_scalar,dt,h;
+	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vx,vy,vel,dvxdx,dvydy;
 	IssmDouble dvx[2],dvy[2];
@@ -5622,5 +5622,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -5630,5 +5630,5 @@
 		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		TripleMultiply(basis,1,numnodes,1,
@@ -5642,5 +5642,5 @@
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
-		D_scalar=dt*gauss->weight*Jdettria;
+		D_scalar=dt*gauss->weight*Jdet;
 
 		D[0][0]=D_scalar*dvxdx;
@@ -5812,5 +5812,5 @@
 	/*Intermediaries */
 	int        i;
-	IssmDouble Jdettria,dt;
+	IssmDouble Jdet,dt;
 	IssmDouble basal_melting_g;
 	IssmDouble old_watercolumn_g;
@@ -5840,5 +5840,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis, gauss);
 
@@ -5847,8 +5847,8 @@
 
 		if(reCast<int,IssmDouble>(dt)){
-			for(i=0;i<numnodes;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]+=Jdet*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
 		}
 		else{
-			for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
+			for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*basal_melting_g*basis[i];
 		}
 	}
@@ -6311,5 +6311,5 @@
 	int        stabilization;
 	int        dim;
-	IssmDouble Jdettria,D_scalar,dt,h;
+	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
 	IssmDouble dvx[2],dvy[2];
@@ -6349,5 +6349,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6357,5 +6357,5 @@
 		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		TripleMultiply(basis,1,numnodes,1,
@@ -6369,5 +6369,5 @@
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
-		D_scalar=dt*gauss->weight*Jdettria;
+		D_scalar=dt*gauss->weight*Jdet;
 
 		D[0][0]=D_scalar*dvxdx;
@@ -6430,5 +6430,5 @@
 	int        dim;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdettria,D_scalar,dt,vx,vy;
+	IssmDouble Jdet,D_scalar,dt,vx,vy;
 
 	/*Fetch number of nodes for this finite element*/
@@ -6466,8 +6466,8 @@
 		vyaverage_input->GetInputValue(&vy,gauss);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		TripleMultiply(basis,1,numnodes,1,
@@ -6480,5 +6480,5 @@
 		GetBprimeMasstransport(B, &xyz_list[0][0], gauss);
 
-		D_scalar=-dt*gauss->weight*Jdettria;
+		D_scalar=-dt*gauss->weight*Jdet;
 		D[0][0]=D_scalar*vx;
 		D[0][1]=0.;
@@ -6506,5 +6506,5 @@
 	int        stabilization;
 	int        dim;
-	IssmDouble Jdettria,D_scalar,dt,h;
+	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6534,5 +6534,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6540,5 +6540,5 @@
 		vy_input->GetInputValue(&vy,gauss);
 
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		TripleMultiply(basis,1,numnodes,1,
@@ -6550,5 +6550,5 @@
 		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
 
-		D_scalar=dt*gauss->weight*Jdettria;
+		D_scalar=dt*gauss->weight*Jdet;
 
 		D[0][0]=D_scalar*vx;
@@ -6604,5 +6604,5 @@
 	int        stabilization;
 	int        dim;
-	IssmDouble Jdettria,D_scalar,dt,h;
+	IssmDouble Jdet,D_scalar,dt,h;
 	IssmDouble vel,vx,vy;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6632,5 +6632,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6638,5 +6638,5 @@
 		vy_input->GetInputValue(&vy,gauss);
 
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		TripleMultiply(basis,1,numnodes,1,
@@ -6648,5 +6648,5 @@
 		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
 
-		D_scalar=dt*gauss->weight*Jdettria;
+		D_scalar=dt*gauss->weight*Jdet;
 
 		D[0][0]=D_scalar*vx;
@@ -6713,5 +6713,5 @@
 
 	/*Intermediaries */
-	IssmDouble Jdettria,dt;
+	IssmDouble Jdet,dt;
 	IssmDouble ms,mb,mb_correction,thickness;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6739,5 +6739,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6750,5 +6750,5 @@
 		 mb_correction=0.;
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
 	}
 
@@ -6763,5 +6763,5 @@
 
 	/*Intermediaries */
-	IssmDouble Jdettria,dt;
+	IssmDouble Jdet,dt;
 	IssmDouble ms,mb,thickness;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6787,5 +6787,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6794,5 +6794,5 @@
 		thickness_input->GetInputValue(&thickness,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
 	}
 
@@ -6807,5 +6807,5 @@
 
 	/*Intermediaries */
-	IssmDouble Jdettria,dt;
+	IssmDouble Jdet,dt;
 	IssmDouble ms,surface,vz;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6832,5 +6832,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6839,5 +6839,5 @@
 		surface_input->GetInputValue(&surface,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
 	}
 
@@ -6852,5 +6852,5 @@
 
 	/*Intermediaries */
-	IssmDouble Jdettria,dt;
+	IssmDouble Jdet,dt;
 	IssmDouble mb,mb_correction,bed,vz;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6878,5 +6878,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
@@ -6889,5 +6889,5 @@
 		 mb_correction=0.;
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
 	}
 
@@ -7062,5 +7062,5 @@
 	int        stabilization;
 	int        i,j,dim;
-	IssmDouble Jdettria,vx,vy,dvxdx,dvydy,vel,h;
+	IssmDouble Jdet,vx,vy,dvxdx,dvydy,vel,h;
 	IssmDouble D_scalar;
 	IssmDouble dvx[2],dvy[2];
@@ -7099,5 +7099,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetBMasstransport(B,&xyz_list[0][0],gauss);
 		GetBprimeMasstransport(Bprime,&xyz_list[0][0],gauss);
@@ -7110,5 +7110,5 @@
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
-		D_scalar=gauss->weight*Jdettria;
+		D_scalar=gauss->weight*Jdet;
 
 		D[0][0]=D_scalar*dvxdx;
@@ -7170,5 +7170,5 @@
 	/*Intermediaries*/
 	int        dim;
-	IssmDouble vx,vy,D_scalar,Jdettria;
+	IssmDouble vx,vy,D_scalar,Jdet;
 	IssmDouble xyz_list[NUMVERTICES][3];
 
@@ -7195,5 +7195,5 @@
 		gauss->GaussPoint(ig);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		/*WARNING: B and Bprime are inverted compared to usual masstransport!!!!*/
 		GetBMasstransport(Bprime,&xyz_list[0][0],gauss);
@@ -7203,5 +7203,5 @@
 		vy_input->GetInputValue(&vy,gauss);
 
-		D_scalar=-gauss->weight*Jdettria;
+		D_scalar=-gauss->weight*Jdet;
 		D[0][0]=D_scalar*vx;
 		D[0][1]=0.;
@@ -7228,10 +7228,10 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
+	IssmDouble dhdt_g,mb_g,ms_g,Jdet;
 	IssmDouble h,gamma,thickness;
 	IssmDouble hnx,hny,dhnx[2],dhny[2];
 	IssmDouble vel,vx,vy,dvxdx,dvydy;
 	IssmDouble D[2][2];
-	int        stabilization=1;
+	int        stabilization=0;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7286,5 +7286,5 @@
 		gamma=h/(2.*thickness+1.e-10);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 		GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
@@ -7292,5 +7292,5 @@
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j] += gauss->weight*Jdettria*(
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
 							(basis[i]+gamma*(basis[i]*(dhnx[0]+dhny[1]) + dbasis[0*numnodes+i]*hnx + dbasis[1*numnodes+i]*hny))*
 							(basis[j]*(dhnx[0]+dhny[1])  + dbasis[0*numnodes+j]*hnx + dbasis[1*numnodes+j]*hny)
@@ -7318,8 +7318,8 @@
 		if(stabilization==1 || stabilization==2){
 			GetBprimeMasstransport(B,&xyz_list[0][0],gauss);
-			D[0][0]=gauss->weight*Jdettria*D[0][0];
-			D[1][0]=gauss->weight*Jdettria*D[1][0];
-			D[0][1]=gauss->weight*Jdettria*D[0][1];
-			D[1][1]=gauss->weight*Jdettria*D[1][1];
+			D[0][0]=gauss->weight*Jdet*D[0][0];
+			D[1][0]=gauss->weight*Jdet*D[1][0];
+			D[0][1]=gauss->weight*Jdet*D[0][1];
+			D[1][1]=gauss->weight*Jdet*D[1][1];
 			TripleMultiply(B,2,numnodes,1,
 						&D[0][0],2,2,0,
@@ -7349,5 +7349,5 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble D[2][2];
-	IssmDouble l=15.;
+	IssmDouble l=12.;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7355,6 +7355,7 @@
 
 	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke    = CreateMassMatrix();
-	IssmDouble*    B     = xNew<IssmDouble>(2*numnodes);
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters);
+	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -7368,24 +7369,25 @@
 		gauss->GaussPoint(ig);
 
-		GetBMasstransport(B,&xyz_list[0][0], gauss);
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		thickness_input->GetInputValue(&thickness,gauss);
 		if(thickness<50.) thickness=50.;
 
-		D_scalar=gauss->weight*Jdet*(l*thickness)*(l*thickness);
-		D[0][0]=D_scalar;
-	   D[1][0]=0.;
-		D[0][1]=0.;
-		D[1][1]=D_scalar;
-
-		TripleMultiply(B,2,numnodes,1,
-					&D[0][0],2,2,0,
-					B,2,numnodes,0,
-					&Ke->values[0],1);
+		GetNodalFunctions(basis,gauss);
+		GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
+
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
+							basis[i]*basis[j]
+							+(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+							);
+			}
+		}
 	}
 
 	/*Clean up and return*/
 	delete gauss;
-	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(basis);
 	return Ke;
 }
@@ -7410,5 +7412,5 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
+	IssmDouble dhdt_g,mb_g,ms_g,Jdet;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7435,8 +7437,8 @@
 		dhdt_input->GetInputValue(&dhdt_g,gauss);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
 	}
 
@@ -7452,5 +7454,5 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble mb_g,ms_g,dhdt_g,Jdettria;
+	IssmDouble mb_g,ms_g,dhdt_g,Jdet;
 
 	/*Fetch number of nodes and dof for this finite element*/
@@ -7477,8 +7479,8 @@
 		dhdt_input->GetInputValue(&dhdt_g,gauss);
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
 	}
 
@@ -7494,5 +7496,5 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble dhdt_g,mb_g,ms_g,Jdettria;
+	IssmDouble dhdt_g,mb_g,ms_g,Jdet;
 	IssmDouble h,gamma,thickness;
 	IssmDouble hnx,hny,dhnx[2],dhny[2];
@@ -7548,5 +7550,5 @@
 		if(thickness<50.) thickness=50.;
 
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 		GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
@@ -7560,5 +7562,5 @@
 
 		for(int i=0;i<numnodes;i++){
-			pe->values[i]+=Jdettria*gauss->weight*(ms_g-mb_g-dhdt_g)*(
+			pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*(
 						basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i])
 						);
@@ -7583,5 +7585,5 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdettria;
+	IssmDouble Jdet;
 	IssmDouble thickness,slope[2];
 	IssmDouble taud_x;
@@ -7608,9 +7610,10 @@
 		surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
 		taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0];
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		//taud_x = slope[0];
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*taud_x*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_x*basis[i];
 	}
 
@@ -7626,5 +7629,5 @@
 	/*Intermediaries */
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdettria;
+	IssmDouble Jdet;
 	IssmDouble thickness,slope[2];
 	IssmDouble taud_y;
@@ -7651,9 +7654,10 @@
 		surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
 		taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1];
-
-		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
+		//taud_y = slope[1];
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		GetNodalFunctions(basis,gauss);
 
-		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*taud_y*basis[i];
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_y*basis[i];
 	}