Index: /issm/trunk/src/c/objects/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Tria.cpp	(revision 601)
+++ /issm/trunk/src/c/objects/Tria.cpp	(revision 602)
@@ -551,4 +551,6 @@
 	/* matrices: */
 	double L[numgrids];
+	double B[2][numgrids];
+	double Bprime[2][numgrids];
 	double DL[2][2]={0.0};
 	double DLprime[2][2]={0.0};
@@ -556,15 +558,17 @@
 	double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 
 	double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
-	double Ke_gg_thickness[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
+	double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
 	
 	double Jdettria;
 	
 	/*input parameters for structural analysis (diagnostic): */
-	double  vxvy_list[numgrids][2]={0,0};
-	double  vx_list[numgrids]={0,0};
-	double  vy_list[numgrids]={0,0};
+	double  vxvy_list[numgrids][2]={0.0};
+	double  vx_list[numgrids]={0.0};
+	double  vy_list[numgrids]={0.0};
 	double  dvx[2];
 	double  dvy[2];
 	double  vx,vy;
+	double  dvxdx,dvydy;
 	double  v_gauss[2]={0.0};
 	double  K[2][2]={0.0};
@@ -636,6 +640,6 @@
 		
 		/*Get B  and B prime matrix: */
-		//GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
-		//GetBprime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
+		GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
 		
 		//Get vx, vy and their derivatives at gauss point
@@ -646,12 +650,12 @@
 		GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
 
-		//dvxdx=dvx[0];
-		//dvydy=dvy[1];
+		dvxdx=dvx[0];
+		dvydy=dvy[1];
 
 		DL_scalar=dt*gauss_weight*Jdettria;
 
 		//Create DL and DLprime matrix
-		//DL[0][0]=DL_scalar*dvxdx;
-		//DL[1][1]=DL_scalar*dvydy;
+		DL[0][0]=DL_scalar*dvxdx;
+		DL[1][1]=DL_scalar*dvydy;
 
 		DLprime[0][0]=DL_scalar*vx;
@@ -661,15 +665,18 @@
 		//Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
 
-		/*TripleMultiply( &B[0][0],numdof,numdof,1,
-					  &DL_scalar,1,1,0,
-					  &L[0],1,numdof,0,
-					  &Ke_gg_gaussian[0][0],0);*/
-
-
-
-
+		TripleMultiply( &B[0][0],2,numdof,1,
+					  &DL[0][0],2,2,0,
+					  &B[0][0],2,numdof,0,
+					  &Ke_gg_thickness1[0][0],0);
+
+		TripleMultiply( &B[0][0],2,numdof,1,
+					  &DLprime[0][0],2,2,0,
+					  &Bprime[0][0],2,numdof,0,
+					  &Ke_gg_thickness2[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( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j];
+		for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j];
 		
 		#ifdef _DEBUGELEMENTS_
@@ -696,11 +703,4 @@
 	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
 
-
-	/*Do not forget to include friction: */
-	if(!shelf){
-		//CreateKMatrixPrognosticFriction(Kgg,inputs,analysis_type,sub_analysis_type);
-	}
-
-
 	#ifdef _DEBUGELEMENTS_
 	if(my_rank==RANK && id==ELID){ 
@@ -766,5 +766,5 @@
 	double  thickness_g;
 	double  dt;
-	int     dofs[1]={1};
+	int     dofs[1]={0};
 	int     found=0;
 
@@ -1706,4 +1706,76 @@
 	}
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetB_prog"
+
+void Tria::GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3){
+
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For grid i, Bi can be expressed in the basic coordinate system
+	 * by: 
+	 *       Bi_basic=[ 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 basic coordinate system: */
+	GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
+
+	#ifdef _DEBUG_ 
+	for (i=0;i<3;i++){
+		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
+	}
+	#endif
+
+	/*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];
+	}
+}
+
+#undef __FUNCT__ 
+#define __FUNCT__ "Tria::GetBPrime_prog"
+
+void Tria::GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3){
+
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For grid i, Bi' can be expressed in the basic coordinate system
+	 * by: 
+	 *       Bi_prime__basic=[ 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 basic coordinate system: */
+	double dh1dh2dh3_basic[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in basic coordinates system : */
+	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh2dh3_basic[0][i]; 
+		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh2dh3_basic[1][i]; 
+	}
+}
+
 
 #undef __FUNCT__ 
Index: /issm/trunk/src/c/objects/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Tria.h	(revision 601)
+++ /issm/trunk/src/c/objects/Tria.h	(revision 602)
@@ -84,4 +84,6 @@
 		void  GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3);
 		void  GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof);
+		void  GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3);
+		void  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
 		void  GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);
 		void  GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3);