Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15716)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15717)
@@ -6606,21 +6606,18 @@
 ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
 
-	/*Constants*/
-	const int    numnodes=NUMVERTICES;
-
 	/*Intermediaries */
 	int        dim;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble Jdettria,dt,vx,vy;
-	IssmDouble basis[NUMVERTICES];
-	IssmDouble B[2][NUMVERTICES];
-	IssmDouble Bprime[2][NUMVERTICES];
-	IssmDouble DL[2][2]={0.0};
-	IssmDouble DLprime[2][2]={0.0};
-	IssmDouble DL_scalar;
-	GaussTria  *gauss=NULL;
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+	IssmDouble Jdettria,D_scalar,dt,vx,vy;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+	IssmDouble     D[2][2];
 
 	/*Retrieve all inputs and parameters*/
@@ -6640,5 +6637,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	GaussTria* gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -6649,29 +6646,33 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		DL_scalar=gauss->weight*Jdettria;
-
-		TripleMultiply(&basis[0],1,numnodes,1,
-					&DL_scalar,1,1,0,
-					&basis[0],1,numnodes,0,
+		GetNodalFunctions(basis,gauss);
+
+		D_scalar=gauss->weight*Jdettria;
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D_scalar,1,1,0,
+					basis,1,numnodes,0,
 					&Ke->values[0],1);
 
 		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
-		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
-
-		DL_scalar=-dt*gauss->weight*Jdettria;
-
-		DLprime[0][0]=DL_scalar*vx;
-		DLprime[1][1]=DL_scalar*vy;
-
-		TripleMultiply( &B[0][0],2,numnodes,1,
-					&DLprime[0][0],2,2,0,
-					&Bprime[0][0],2,numnodes,0,
+		GetBPrognostic(Bprime, &xyz_list[0][0], gauss);
+		GetBprimePrognostic(B, &xyz_list[0][0], gauss);
+
+		D_scalar=-dt*gauss->weight*Jdettria;
+		D[0][0]=D_scalar*vx;
+		D[0][1]=0.;
+		D[1][0]=0.;
+		D[1][1]=D_scalar*vy;
+
+		TripleMultiply(B,2,numnodes,1,
+					&D[0][0],2,2,0,
+					Bprime,2,numnodes,0,
 					&Ke->values[0],1);
 	}
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	return Ke;
@@ -6744,16 +6745,15 @@
 ElementVector* Tria::CreatePVectorPrognostic_DG(void){
 
-	/*Constants*/
-	const int    numnodes=NUMVERTICES;
-
 	/*Intermediaries */
 	IssmDouble Jdettria,dt;
 	IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble basis[NUMVERTICES];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -6765,5 +6765,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	GaussTria* gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -6771,5 +6771,5 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
+		GetNodalFunctions(basis,gauss);
 
 		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
@@ -6781,4 +6781,5 @@
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return pe;
@@ -7055,19 +7056,18 @@
 ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
 
-	/*Constants*/
-	const int  numnodes=NUMVERTICES;
-
 	/*Intermediaries*/
-	int        i,j,dim;
-	IssmDouble vx,vy,Jdettria;
+	int        dim;
+	IssmDouble vx,vy,D_scalar,Jdettria;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble B[2][NUMVERTICES];
-	IssmDouble Bprime[2][NUMVERTICES];
-	IssmDouble DL[2][2]={0.0};
-	IssmDouble DL_scalar;
-	GaussTria  *gauss=NULL;
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
+	IssmDouble     D[2][2];
 
 	/*Retrieve all inputs and parameters*/
@@ -7078,5 +7078,5 @@
 
 	/*Start looping on the number of gaussian points:*/
-	gauss=new GaussTria(2);
+	GaussTria* gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -7085,21 +7085,26 @@
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
 		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
-		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
+		GetBPrognostic(Bprime,&xyz_list[0][0],gauss);
+		GetBprimePrognostic(B,&xyz_list[0][0],gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 
-		DL_scalar=-gauss->weight*Jdettria;
-		DL[0][0]=DL_scalar*vx;
-		DL[1][1]=DL_scalar*vy;
-
-		TripleMultiply( &B[0][0],2,numnodes,1,
-					&DL[0][0],2,2,0,
-					&Bprime[0][0],2,numnodes,0,
+		D_scalar=-gauss->weight*Jdettria;
+		D[0][0]=D_scalar*vx;
+		D[0][1]=0.;
+		D[1][0]=0.;
+		D[1][1]=D_scalar*vy;
+
+		TripleMultiply(B,2,numnodes,1,
+					&D[0][0],2,2,0,
+					Bprime,2,numnodes,0,
 					&Ke->values[0],1);
 	}
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	return Ke;
@@ -7165,16 +7170,14 @@
 ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
 
-	/*Constants*/
-	const int    numnodes=NUMVERTICES;
-
 	/*Intermediaries */
-	int        i,j;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
-	IssmDouble basis[NUMVERTICES];
-	GaussTria* gauss=NULL;
-
-	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -7185,5 +7188,5 @@
 
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
+	GaussTria* gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -7195,10 +7198,11 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetNodalFunctions(&basis[0],gauss);
-
-		for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
+		GetNodalFunctions(basis,gauss);
+
+		for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
 	}
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return pe;