Index: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5635)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5636)
@@ -1931,11 +1931,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* parameters: */
@@ -2008,31 +2003,17 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-		/*Compute misfit at gaussian point: */
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
-
-		/*compute SurfaceRelVelMisfit*/
-		Jelem+=misfit*Jdet*gauss_weight;
-	}
-
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-
-	/*Return: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
+		Jelem+=misfit*Jdet*gauss->weight;
+	}
+
+	/*clean up and Return: */
+	delete gauss;
 	return Jelem;
 }
@@ -2062,11 +2043,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* parameters: */
@@ -2083,5 +2059,4 @@
 	/*inputs: */
 	bool onwater;
-	
 	double  meanvel, epsvel;
 
@@ -2137,31 +2112,17 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-		/*Compute misfit at gaussian point: */
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
-
-		/*compute SurfaceLogVelMisfit*/
-		Jelem+=misfit*Jdet*gauss_weight;
-	}
-
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-		
-	/*Return: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
+		Jelem+=misfit*Jdet*gauss->weight;
+	}
+
+	/*clean-up and Return: */
+	delete gauss;
 	return Jelem;
 }
@@ -2191,11 +2152,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* parameters: */
@@ -2268,31 +2224,17 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-		/*Compute misfit at gaussian point: */
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
-
-		/*compute SurfaceLogVxVyMisfit*/
-		Jelem+=misfit*Jdet*gauss_weight;
-	}
-
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-		
-	/*Return: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
+		Jelem+=misfit*Jdet*gauss->weight;
+	}
+
+	/*clean-up and Return: */
+	delete gauss;
 	return Jelem;
 }
@@ -2322,11 +2264,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* parameters: */
@@ -2398,31 +2335,17 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-
-		/*Compute misfit at gaussian point: */
-		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
-
-		/*compute SurfaceAverageVelMisfit*/
-		Jelem+=misfit*Jdet*gauss_weight;
-	}
-
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
-
-	/*Return: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+		TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
+		Jelem+=misfit*Jdet*gauss->weight;
+	}
+
+	/*clean-up and Return: */
+		delete gauss;
 	return Jelem;
 }
@@ -2434,5 +2357,4 @@
 	int row;
 	int vertices_ids[3];
-
 
 	/*recover pointer: */
@@ -2779,11 +2701,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* matrices: */
@@ -2839,6 +2756,8 @@
 	if(artdiff){
 		//Get the Jacobian determinant
-		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+		gauss=new GaussTria();
+		gauss->GaussCenter();
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
+		delete gauss;
 
 		//Build K matrix (artificial diffusivity matrix)
@@ -2850,33 +2769,28 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
 
 		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
 
 		/*Get B  and B prime matrix: */
-		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
-		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
 
 		//Get vx, vy and their derivatives at gauss point
-		vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);
-		vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);
-
-		vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
-		vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
+		vxaverage_input->GetParameterValue(&vx,gauss);
+		vyaverage_input->GetParameterValue(&vy,gauss);
+
+		vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
+		vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
 
-		DL_scalar=gauss_weight*Jdettria;
+		DL_scalar=gauss->weight*Jdettria;
 
 		//Create DL and DLprime matrix
@@ -2924,8 +2838,6 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
+	/*Clean up*/
+	delete gauss;
 	xfree((void**)&doflist);
 }
@@ -2945,11 +2857,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* matrices: */
@@ -2977,7 +2884,4 @@
 	this->parameters->FindParam(&dim,DimEnum);
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/*Retrieve all inputs we will be needed*/
 	if(dim==2){
@@ -2987,24 +2891,22 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
 
 		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
 
 		/*Get B  and B prime matrix: */
 		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
-		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
-		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
+		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
 
 		//Get vx, vy and their derivatives at gauss point
-		vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]);
-		vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]);
-
-		DL_scalar=-gauss_weight*Jdettria;
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+
+		DL_scalar=-gauss->weight*Jdettria;
 
 		DLprime[0][0]=DL_scalar*vx;
@@ -3025,8 +2927,6 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
+	/*Clean up and return*/
+	delete gauss;
 	xfree((void**)&doflist);
 }
@@ -3046,11 +2946,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* matrices: */
@@ -3129,6 +3024,8 @@
 	if(artdiff){
 		//Get the Jacobian determinant
-		gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+		gauss=new GaussTria();
+		gauss->GaussCenter();
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
+		delete gauss;
 
 		//Build K matrix (artificial diffusivity matrix)
@@ -3140,33 +3037,28 @@
 	}
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
 
 		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
+		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
 
 		/*Get B  and B prime matrix: */
-		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
-		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
 
 		//Get vx, vy and their derivatives at gauss point
-		vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]);
-		vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]);
-
-		vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
-		vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
+		vxaverage_input->GetParameterValue(&vx,gauss);
+		vyaverage_input->GetParameterValue(&vy,gauss);
+
+		vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
+		vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
 
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
 
-		DL_scalar=gauss_weight*Jdettria;
+		DL_scalar=gauss->weight*Jdettria;
 
 		DLprime[0][0]=DL_scalar*nx;
@@ -3204,8 +3096,6 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
+	/*Clean up and return*/
+	delete gauss;
 	xfree((void**)&doflist);
 }
@@ -3224,11 +3114,6 @@
 
 	/* gaussian points: */
-	int     num_gauss,ig;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* gauss_weights           =  NULL;
-	double  gauss_weight;
-	double  gauss_l1l2l3[3];
+	int     ig;
+	GaussTria *gauss=NULL;
 
 	/* material data: */
@@ -3281,7 +3166,4 @@
 	GetDofList(&doflist,MacAyealApproximationEnum);
 
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
-
 	/*Retrieve all inputs we will be needing: */
 	thickness_input=inputs->GetInput(ThicknessEnum);
@@ -3292,18 +3174,15 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for (ig=0; ig<num_gauss; ig++){
-		/*Pick up the gaussian point: */
-		gauss_weight=*(gauss_weights+ig);
-		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
-		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
-		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
-
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
 
 		/*Compute thickness at gaussian point: */
-		thickness_input->GetParameterValue(&thickness, gauss_l1l2l3);
+		thickness_input->GetParameterValue(&thickness, gauss);
 
 		/*Get strain rate from velocity: */
-		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input);
-		this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,vxold_input,vyold_input);
+		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
 
 		/*Get viscosity: */
@@ -3312,10 +3191,10 @@
 
 		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
 
 		/* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
 			onto this scalar matrix, so that we win some computational time: */
 		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2*newviscosity*thickness*gauss_weight*Jdet;
+		D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
 
 		for (i=0;i<3;i++){
@@ -3324,6 +3203,6 @@
 
 		/*Get B and Bprime matrices: */
-		GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
-		GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
 
 		/*  Do the triple product tB*D*Bprime: */
@@ -3346,8 +3225,7 @@
 	}
 
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&gauss_weights);
+
+	/*Clean up and return*/
+	delete gauss;
 	xfree((void**)&doflist);
 }
@@ -6949,5 +6827,4 @@
 /*}}}1*/
 /*FUNCTION Tria::SurfaceNormal{{{1*/
-
 void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
 
@@ -6972,5 +6849,4 @@
 	*(surface_normal+1)=normal[1]/normal_norm;
 	*(surface_normal+2)=normal[2]/normal_norm;
-
-}
-/*}}}*/
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5635)
+++ issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5636)
@@ -149,4 +149,33 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::GetBPrognostic{{{1*/
+void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss){
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For grid i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ h ]
+	 *          [ h ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
+	 */
+
+	int i;
+	const int NDOF1=1;
+	const int numgrids=3;
+
+	double l1l2l3[numgrids];
+
+
+	/*Get dh1dh2dh3 in actual coordinate system: */
+	GetNodalFunctions(&l1l2l3[0],gauss);
+
+	/*Build B_prog: */
+	for (i=0;i<numgrids;i++){
+		*(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
+		*(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
+	}
+}
+/*}}}*/
 /*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss){
@@ -184,6 +213,71 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
+void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){
+
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For grid i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ 2*dh/dx    dh/dy ]
+	 *                [   dh/dx  2*dh/dy ]
+	 *                [   dh/dy    dh/dx ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
+	 */
+
+	int i;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i]; 
+	}
+}
+/*}}}*/
 /*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
 void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss){
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For grid i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ dh/dx ]
+	 *                [ dh/dy ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
+	 */
+
+	int i;
+	const int NDOF1=1;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i]; 
+		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i]; 
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
+void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss){
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For grid i, Bi' can be expressed in the actual coordinate system
Index: issm/trunk/src/c/objects/Elements/TriaRef.h
===================================================================
--- issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 5635)
+++ issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 5636)
@@ -28,6 +28,9 @@
 		void GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss);
 		void GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss);
+		void GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss);
 		void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);
+		void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss);
 		void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);
+		void GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss);
 		void GetL(double* L, double* xyz_list,double* gauss,int numdof);
 		void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
