Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 5624)
+++ /issm/trunk/src/c/Makefile.am	(revision 5625)
@@ -77,4 +77,6 @@
 					./objects/Bamg/Mesh.cpp\
 					./objects/Bamg/Mesh.h\
+					./objects/Gauss/GaussTria.h\
+					./objects/Gauss/GaussTria.cpp\
 					./objects/Update.h\
 					./objects/Element.h\
@@ -630,4 +632,6 @@
 					./objects/Bamg/Mesh.h\
 					./objects/Bamg/Mesh.cpp\
+					./objects/Gauss/GaussTria.h\
+					./objects/Gauss/GaussTria.cpp\
 					./objects/Update.h\
 					./objects/Element.h\
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5624)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5625)
@@ -4580,74 +4580,46 @@
 void  Tria::CreateKMatrixSlope(Mat Kgg){
 
-	/* local declarations */
-	int             i,j;
-
-	/* node data: */
-	const int    numgrids=3;
+	/*constants: */
+	const int    numnodes=3;
 	const int    NDOF1=1;
-	const int    numdof=NDOF1*numgrids;
-	double       xyz_list[numgrids][3];
-	int*         doflist=NULL;
-	
-	/* 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];
-
-	/* matrices: */
-	double L[1][3];
-	double DL_scalar;
-
-	/* local element matrices: */
-	double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
-	double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-	
-	double Jdet;
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
+	const int    numdof=NDOF1*numnodes;
+
+	/* Intermediaries */
+	int        i,j,ig;
+	double     DL_scalar,Jdet;
+	double     xyz_list[numnodes][3];
+	double     L[1][3];
+	double     Ke[numdof][numdof] = {0.0};
+	double     Ke_g[numdof][numdof];
+	GaussTria *gauss = NULL;
+	int       *doflist = NULL;
+
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes);
 	GetDofList(&doflist);
 
-	/* 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);
-
+	/* Start looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
 		
-		/*Get L matrix: */
-		GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
-		
-		DL_scalar=gauss_weight*Jdet;
-
-		/*  Do the triple producte tL*D*L: */
-		TripleMultiply( &L[0][0],1,3,1,
+		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+		DL_scalar=gauss->weight*Jdet;
+
+		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1);
+
+		TripleMultiply(&L[0][0],1,3,1,
 					&DL_scalar,1,1,0,
 					&L[0][0],1,3,0,
-					&Ke_gg_gaussian[0][0],0);
-
-		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
-		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-	} //for (ig=0; ig<num_gauss; ig++
-
-	/*Add Ke_gg to global matrix Kgg: */
-	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+					&Ke_g[0][0],0);
+
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j];
+	}
+
+	/*Add Ke to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,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);
 }
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5624)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 5625)
@@ -180,5 +180,5 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetL {{{1*/
+/*FUNCTION TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof) {{{1*/
 void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){
 	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
@@ -221,5 +221,44 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetJacobian {{{1*/
+/*FUNCTION TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof) {{{1*/
+void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){
+	/*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: 
+	 *       numdof=1: 
+	 *                 Li=h;
+	 *       numdof=2:
+	 *                 Li=[ h   0 ]
+	 *                    [ 0   h ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
+	 */
+
+	int i;
+	const int NDOF2=2;
+	const int numgrids=3;
+	double l1l2l3[3];
+
+	/*Get l1l2l3 in actual coordinate system: */
+	GetNodalFunctions(l1l2l3,gauss);
+
+	/*Build L: */
+	if(numdof==1){
+		for (i=0;i<numgrids;i++){
+			L[i]=l1l2l3[i]; 
+		}
+	}
+	else{
+		for (i=0;i<numgrids;i++){
+			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
+			*(L+numdof*numgrids*0+numdof*i+1)=0;
+			*(L+numdof*numgrids*1+numdof*i)=0;
+			*(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
+		}
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss) {{{1*/
 void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){
 	/*The Jacobian is constant over the element, discard the gaussian points. 
@@ -244,7 +283,45 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetJacobianDeterminant2d {{{1*/
+/*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/
+void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){
+	/*The Jacobian is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	const int NDOF2=2;
+	const int numgrids=3;
+	double x1,y1,x2,y2,x3,y3;
+
+	x1=*(xyz_list+numgrids*0+0);
+	y1=*(xyz_list+numgrids*0+1);
+	x2=*(xyz_list+numgrids*1+0);
+	y2=*(xyz_list+numgrids*1+1);
+	x3=*(xyz_list+numgrids*2+0);
+	y3=*(xyz_list+numgrids*2+1);
+
+
+	*(J+NDOF2*0+0)=0.5*(x2-x1);
+	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
+	*(J+NDOF2*0+1)=0.5*(y2-y1);
+	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss) {{{1*/
 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){
 
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+	double J[2][2];
+
+	/*Get Jacobian*/
+	GetJacobian(&J[0][0],xyz_list,gauss);
+
+	/*Get Determinant*/
+	Matrix2x2Determinant(Jdet,&J[0][0]);
+	if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/
+void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){
 	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
 	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
@@ -296,5 +373,5 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetNodalFunctions {{{1*/
+/*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss) {{{1*/
 void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
@@ -303,4 +380,14 @@
 	l1l2l3[1]=gauss[1];
 	l1l2l3[2]=gauss[2];
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss) {{{1*/
+void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	l1l2l3[0]=gauss->coord1;
+	l1l2l3[1]=gauss->coord2;
+	l1l2l3[2]=gauss->coord3;
 
 }
Index: /issm/trunk/src/c/objects/Elements/TriaRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 5624)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 5625)
@@ -7,4 +7,6 @@
 #ifndef _TRIAREF_H_
 #define _TRIAREF_H_
+
+class GaussTria;
 
 class TriaRef{
@@ -27,10 +29,14 @@
 		void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);
 		void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);
-		void GetL(double* L, double* xyz_list, double* gauss,int numdof);
+		void GetL(double* L, double* xyz_list,double* gauss,int numdof);
+		void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
 		void GetJacobian(double* J, double* xyz_list,double* gauss);
+		void GetJacobian(double* J, double* xyz_list,GaussTria* gauss);
 		void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);
+		void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);
 		void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
 		void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
 		void GetNodalFunctions(double* l1l2l3, double* gauss);
+		void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
 		void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
 		void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);
Index: /issm/trunk/src/c/objects/Gauss/GaussTria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Gauss/GaussTria.cpp	(revision 5625)
+++ /issm/trunk/src/c/objects/Gauss/GaussTria.cpp	(revision 5625)
@@ -0,0 +1,145 @@
+/*!\file GaussTria.c
+ * \brief: implementation of the GaussTria object
+ */
+
+/*Include files: {{{1*/
+#include "./../objects.h"
+/*}}}*/
+
+/*GaussTria constructors and destructors:*/
+/*FUNCTION GaussTria::GaussTria() {{{1*/
+GaussTria::GaussTria(){
+
+	numgauss=-1;
+
+	weights=NULL;
+	coords1=NULL;
+	coords2=NULL;
+	coords3=NULL;
+
+	weight=UNDEF;
+	coord1=UNDEF;
+	coord2=UNDEF;
+	coord3=UNDEF;
+}
+/*}}}*/
+/*FUNCTION GaussTria::GaussTria(int order) {{{1*/
+GaussTria::GaussTria(int order){
+
+	/*Get gauss points*/
+	GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
+
+	/*Initialize static fields as undefinite*/
+	weight=UNDEF;
+	coord1=UNDEF;
+	coord2=UNDEF;
+	coord3=UNDEF;
+
+}
+/*}}}*/
+/*FUNCTION GaussTria::GaussTria(int order,int node1,int node2) {{{1*/
+GaussTria::GaussTria(int order,int node1,int node2){
+
+	ISSMERROR("not implemented yet");
+
+	/*Get gauss points*/
+	GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
+
+	/*Initialize static fields as undefinite*/
+	weight=UNDEF;
+	coord1=UNDEF;
+	coord2=UNDEF;
+	coord3=UNDEF;
+
+}
+/*}}}*/
+/*FUNCTION GaussTria::~GaussTria(){{{1*/
+GaussTria::~GaussTria(){
+	xfree((void**)&weights);
+	xfree((void**)&coords1);
+	xfree((void**)&coords2);
+	xfree((void**)&coords3);
+}
+/*}}}*/
+
+/*Methods*/
+/*FUNCTION GaussTria::Echo{{{1*/
+void GaussTria::Echo(void){
+
+	printf("GaussTria:\n");
+	printf("   numgauss: %i\n",numgauss);
+
+	if (weights){
+	 printf("   weights = ["); 
+	 for(int i=0;i<numgauss;i++) printf(" %g\n",weights[i]);
+	 printf("]\n");
+	}
+	else printf("weights = NULL\n");
+	if (coords1){
+	 printf("   coords1 = ["); 
+	 for(int i=0;i<numgauss;i++) printf(" %g\n",coords1[i]);
+	 printf("]\n");
+	}
+	else printf("coords1 = NULL\n");
+	if (coords2){
+	 printf("   coords2 = ["); 
+	 for(int i=0;i<numgauss;i++) printf(" %g\n",coords2[i]);
+	 printf("]\n");
+	}
+	else printf("coords2 = NULL\n");
+	if (coords3){
+	 printf("   coords3 = ["); 
+	 for(int i=0;i<numgauss;i++) printf(" %g\n",coords3[i]);
+	 printf("]\n");
+	}
+	else printf("coords3 = NULL\n");
+
+	printf("   weight = %g\n",weight);
+	printf("   coord1 = %g\n",coord1);
+	printf("   coord2 = %g\n",coord2);
+	printf("   coord3 = %g\n",coord3);
+
+}
+/*}}}*/
+/*FUNCTION GaussTria::GaussPoint{{{1*/
+void GaussTria::GaussPoint(int ig){
+
+	/*Check input in debugging mode*/
+	 ISSMASSERT(ig>=0 && ig< numgauss);
+
+	 /*update static arrays*/
+	 weight=weights[ig];
+	 coord1=coords1[ig];
+	 coord2=coords2[ig];
+	 coord3=coords3[ig];
+
+}
+/*}}}*/
+/*FUNCTION GaussTria::begin{{{1*/
+int GaussTria::begin(void){
+
+	/*Check that this has been initialized*/
+	ISSMASSERT(numgauss>0);
+	ISSMASSERT(weights);
+	ISSMASSERT(coords1);
+	ISSMASSERT(coords2);
+	ISSMASSERT(coords3);
+
+	/*return first gauss index*/
+	return 0;
+}
+/*}}}*/
+/*FUNCTION GaussTria::end{{{1*/
+int GaussTria::end(void){
+
+	/*Check that this has been initialized*/
+	ISSMASSERT(numgauss>0);
+	ISSMASSERT(weights);
+	ISSMASSERT(coords1);
+	ISSMASSERT(coords2);
+	ISSMASSERT(coords3);
+
+	/*return last gauss index +1*/
+	return numgauss;
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Gauss/GaussTria.h
===================================================================
--- /issm/trunk/src/c/objects/Gauss/GaussTria.h	(revision 5625)
+++ /issm/trunk/src/c/objects/Gauss/GaussTria.h	(revision 5625)
@@ -0,0 +1,42 @@
+/*!\file GaussTria.h
+ * \brief: header file for node object
+ */
+
+#ifndef _GAUSSTRIA_H_
+#define _GAUSSTRIA_H_
+
+/*Headers:*/
+/*{{{1*/
+#include "./../../shared/shared.h"
+/*}}}*/
+
+class GaussTria{
+
+	private:
+		int numgauss;
+		double* weights;
+		double* coords1;
+		double* coords2;
+		double* coords3;
+
+	public:
+		double weight;
+		double coord1;
+		double coord2;
+		double coord3;
+		
+	public:
+
+		/*GaussTria constructors, destructors*/
+		GaussTria();
+		GaussTria(int order);
+		GaussTria(int order,int node1, int node2);
+		~GaussTria();
+
+		/*Methods*/
+		int  begin(void);
+		int  end(void);
+		void Echo(void);
+		void GaussPoint(int ig);
+};
+#endif  /* _GAUSSTRIA_H_ */
Index: /issm/trunk/src/c/objects/objects.h
===================================================================
--- /issm/trunk/src/c/objects/objects.h	(revision 5624)
+++ /issm/trunk/src/c/objects/objects.h	(revision 5625)
@@ -22,4 +22,7 @@
 /*Constraints: */
 #include "./Constraints/Spc.h"
+
+/*Gauss*/
+#include "./Gauss/GaussTria.h"
 
 /*Loads: */