Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5676)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5677)
@@ -2746,13 +2746,8 @@
 
 	/* gaussian points: */
-	int     num_area_gauss;
-	int     igarea,igvert;
-	double* first_gauss_area_coord  =  NULL;
-	double* second_gauss_area_coord =  NULL;
-	double* third_gauss_area_coord  =  NULL;
-	double* vert_gauss_coord = NULL;
-	double* area_gauss_weights  =  NULL;
-	double* vert_gauss_weights  =  NULL;
-	double  gaussgrids[NUMVERTICES][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+	int     ig;
+	GaussPenta *gauss=NULL;
+	GaussTria  *gauss_tria=NULL;
+
 
 	/* specific gaussian point: */
@@ -2806,10 +2801,4 @@
 		points, same deal, which yields 3 gaussian points.*/
 
-	area_order=5;
-	num_vert_gauss=5;
-
-	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-	gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
-
 	/*Retrieve all inputs we will be needing: */
 	vx_input=inputs->GetInput(VxEnum);
@@ -2818,51 +2807,37 @@
 
 	/* Start  looping on the number of gaussian points: */
-	for (igarea=0; igarea<num_area_gauss; igarea++){
-		for (igvert=0; igvert<num_vert_gauss; igvert++){
-			/*Pick up the gaussian point: */
-			area_gauss_weight=*(area_gauss_weights+igarea);
-			vert_gauss_weight=*(vert_gauss_weights+igvert);
-			gauss_weight=area_gauss_weight*vert_gauss_weight;
-			gauss_coord[0]=*(first_gauss_area_coord+igarea); 
-			gauss_coord[1]=*(second_gauss_area_coord+igarea);
-			gauss_coord[2]=*(third_gauss_area_coord+igarea);
-			gauss_coord[3]=*(vert_gauss_coord+igvert);
-
-			/*Compute strain rate: */
-			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
-
-			/*Get viscosity: */
-			matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
-
-			/*Get B and Bprime matrices: */
-			GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); 
-			GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); 
-
-			/* Get Jacobian determinant: */
-			GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]);
-
-			/* 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: */
-			D_scalar=gauss_weight*Jdet;
-			for (i=0;i<6;i++){
-				D[i][i]=D_scalar*2*viscosity;
-			}
-			for (i=6;i<8;i++){
-				D[i][i]=-D_scalar*stokesreconditioning;
-			}
-
-			/*  Do the triple product tB*D*Bprime: */
-			TripleMultiply( &B[0][0],8,27,1,
-						&D[0][0],8,8,0,
-						&B_prime[0][0],8,27,0,
-						&Ke_gaussian[0][0],0);
-
-			/*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
-			for(i=0;i<27;i++){
-				for(j=0;j<27;j++){
-					Ke_temp[i][j]+=Ke_gaussian[i][j];
-				}
-			}
-		}
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(5,5);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/*Compute strain rate: */
+		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+
+		/*Get viscosity: */
+		matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+		/*Get B and Bprime matrices: */
+		GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
+		GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
+
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&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: */
+		D_scalar=gauss->weight*Jdet;
+		for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
+		for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
+
+		/*  Do the triple product tB*D*Bprime: */
+		TripleMultiply( &B[0][0],8,27,1,
+					&D[0][0],8,8,0,
+					&B_prime[0][0],8,27,0,
+					&Ke_gaussian[0][0],0);
+
+		/*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */
+		for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
 	}
 
@@ -2878,28 +2853,23 @@
 		}
 
-		xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
-		GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
-
-		/* Start looping on the number of gauss 2d (nodes on the bedrock) */
-		for (igarea=0; igarea<num_area_gauss; igarea++){
-			gauss_weight=*(area_gauss_weights+igarea);
-			gauss_coord[0]=*(first_gauss_area_coord+igarea); 
-			gauss_coord[1]=*(second_gauss_area_coord+igarea);
-			gauss_coord[2]=*(third_gauss_area_coord+igarea);
-			gauss_coord[3]=-1;
-
-			gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); 
-			gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
-			gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
+		/* Start  looping on the number of gaussian points: */
+		delete gauss; //gauss of previous loop
+		gauss=new GaussPenta();
+		gauss->GaussFaceTria(0,1,2,2);
+		gauss_tria=new GaussTria();
+		for (ig=gauss->begin();ig<gauss->end();ig++){
+
+			gauss->GaussPoint(ig);
+			gauss->SynchronizeGaussTria(gauss_tria);
 
 			/*Get the Jacobian determinant */
-			tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
+			tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);
 
 			/*Get L matrix if viscous basal drag present: */
-			GetLStokes(&LStokes[0][0],  gauss_coord);
-			GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord);
+			GetLStokes(&LStokes[0][0], gauss);
+			GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
 
 			/*Compute strain rate: */
-			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
+			this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 
 			/*Get viscosity at last iteration: */
@@ -2910,20 +2880,20 @@
 
 			/*Calculate DL on gauss point */
-			friction->GetAlpha2(&alpha2_gauss, gauss_coord,VxEnum,VyEnum,VzEnum);
-
-			DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d;
-			DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d;
-			DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
-			DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
-			DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2];
-			DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2];
-			DLStokes[6][6]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0];
-			DLStokes[7][7]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1];
-			DLStokes[8][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[2];
-			DLStokes[9][8]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0;
-			DLStokes[10][10]=-2*viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0;
-			DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0];
-			DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1];
-			DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2];
+			friction->GetAlpha2(&alpha2_gauss, gauss_tria,VxEnum,VyEnum,VzEnum);
+
+			DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
+			DLStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
+			DLStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+			DLStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+			DLStokes[4][4]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+			DLStokes[5][5]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+			DLStokes[6][6]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
+			DLStokes[7][7]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
+			DLStokes[8][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[2];
+			DLStokes[9][8]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0]/2.0;
+			DLStokes[10][10]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1]/2.0;
+			DLStokes[11][11]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
+			DLStokes[12][12]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
+			DLStokes[13][13]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[2];
 
 			/*  Do the triple product tL*D*L: */
@@ -2933,9 +2903,5 @@
 						&Ke_drag_gaussian[0][0],0);
 
-			for(i=0;i<numdof2d;i++){
-				for(j=0;j<numdof2d;j++){
-					Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
-				}
-			}
+			for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_temp[i][j]+=Ke_drag_gaussian[i][j];
 		}
 	
@@ -2952,10 +2918,6 @@
 
 	/*Free ressources:*/
-	xfree((void**)&first_gauss_area_coord);
-	xfree((void**)&second_gauss_area_coord);
-	xfree((void**)&third_gauss_area_coord);
-	xfree((void**)&area_gauss_weights);
-	xfree((void**)&vert_gauss_coord);
-	xfree((void**)&vert_gauss_weights);
+	delete gauss;
+	delete gauss_tria;
 	xfree((void**)&doflist);
 
@@ -3955,5 +3917,5 @@
 
 			/*Get L matrix : */
-			GetNodalFunctionsP1(l1l6, gauss);
+			GetNodalFunctionsP1(l1l6, gauss_coord);
 
 			/*Get water_pressure at gaussian point */
@@ -3965,5 +3927,5 @@
 			for(i=0;i<NUMVERTICES2D;i++){
 				for(j=0;j<3;j++){
-					Pe_temp[i*NDOF4+j]+=water_pressure*gauss_weight*Jdet2d*L[i]*bed_normal[j];
+					Pe_temp[i*NDOF4+j]+=water_pressure*gauss_weight*Jdet2d*l1l6[i]*bed_normal[j];
 				}
 			}
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 5676)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 5677)
@@ -1406,6 +1406,103 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::GetLStokes {{{1*/
+void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){
+	/*
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For grid i, Li can be expressed in the actual coordinate system
+	 * by: 
+	 *       Li=[ h    0    0   0]
+	 *	 	      [ 0    h    0   0]
+	 *		      [ 0    0    h   0]
+	 *		      [ 0    0    h   0]
+	 *	 	      [ h    0    0   0]
+	 *	 	      [ 0    h    0   0]
+	 *	 	      [ h    0    0   0]
+	 *	 	      [ 0    h    0   0]
+	 *		      [ 0    0    h   0]
+	 *		      [ 0    0    h   0]
+	 *		      [ 0    0    h   0]
+	 *	 	      [ h    0    0   0]
+	 *	 	      [ 0    h    0   0]
+	 *		      [ 0    0    h   0]
+	 * where h is the interpolation function for grid i.
+	 */
+
+	int i;
+	const int numgrids2d=3;
+	int num_dof=4;
+
+	double l1l2l3[numgrids2d];
+
+
+	/*Get l1l2l3 in actual coordinate system: */
+	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LStokes: */
+	for (i=0;i<3;i++){
+		*(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
+		*(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
+		*(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
+		*(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
+		*(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
+
+	}
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/
-void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss){
+void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
 
 	/*
@@ -1437,7 +1534,7 @@
 
 	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss_tria[0];
-	l1l2l3[1]=gauss_tria[1];
-	l1l2l3[2]=gauss_tria[2];
+	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
Index: /issm/trunk/src/c/objects/Elements/PentaRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 5676)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 5677)
@@ -70,7 +70,7 @@
 		void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
 		void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss);
-		//void GetLStokes(double* LStokes, double* gauss_tria);
-		void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss_tria, GaussPenta* gauss);
-		void GetParameterValue(double* pvalue,double* plist,GaussPenta* gauss);
+		void GetLStokes(double* LStokes, GaussPenta* gauss);
+		void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
+		void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss);
 		//void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
 		void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss);