Index: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10452)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10480)
@@ -5986,4 +5986,5 @@
 	ElementMatrix* Ke1=CreateKMatrixDiagnosticStokesViscous();
 	ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction();
+	//ElementMatrix* Ke2=CreateKMatrixDiagnosticStokesFriction2();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
@@ -6136,9 +6137,111 @@
 
 	/*Transform Coordinate System*/
+	if(id==1) printarray(&bed_normal[0],1,3);
+	if(id==1) printarray(&nodes[0]->coord_system[0][0],3,3);
+	if(id==1) printf("Before rotation\n");
+	if(id==1) printarray(Ke->values,4,Ke->ncols);
 	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
+	if(id==1) printf("After rotation\n");
+	if(id==1) printarray(Ke->values,4,Ke->ncols);
+	if(id==1) printf("New:\n");
+	if(id==1) Ke=CreateKMatrixDiagnosticStokesFriction2();
+
+	//if(id==56) printarray(Ke->values,1,9);
+	//if(id==56) printarray(&nodes[0]->coord_system[0][0],3,3);
+	//if(id==56) _error_("STOP");
 
 	/*Clean up and return*/
 	delete gauss;
 	delete friction;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction2{{{1*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction2(void){
+
+	/*Constants*/
+	const int numdof=NUMVERTICES*NDOF4;
+	const int numdof2d=NUMVERTICES2D*NDOF4;
+
+	/*Intermediaries */
+	int        i,j,ig;
+	int        analysis_type,approximation;
+	double     alpha2,Jdet2d;
+	double     stokesreconditioning,viscosity;
+	double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+	double     xyz_list[NUMVERTICES][3];
+	double	  xyz_list_tria[NUMVERTICES2D][3];
+	double     LStokes2[2][numdof2d];
+	//double     LprimeStokes2[4][numdof2d];
+	double     DLStokes[2][2]={0.0};
+	double     Ke_drag_gaussian[numdof2d][numdof2d];
+	Friction*  friction=NULL;
+	GaussPenta *gauss=NULL;
+
+	/*If on water or not Stokes, skip stiffness: */
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+	if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
+	for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+
+	/*build friction object, used later on: */
+	friction=new Friction("3d",inputs,matpar,analysis_type);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(0,1,2,2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
+		GetLStokes3(&LStokes2[0][0], gauss);
+		//GetLprimeStokes2(&LprimeStokes2[0][0], &xyz_list[0][0], gauss);
+
+		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+		matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
+
+		DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
+		DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
+		//DLStokes[2][2] = +2*viscosity*gauss->weight*Jdet2d;
+		//DLStokes[3][3] = -stokesreconditioning*gauss->weight*Jdet2d;
+
+		//TripleMultiply( &LStokes2[0][0],4,numdof2d,1,
+		//			&DLStokes[0][0],4,4,0,
+		//			&LprimeStokes2[0][0],4,numdof2d,0,
+		//			&Ke_drag_gaussian[0][0],0);
+		TripleMultiply( &LStokes2[0][0],2,numdof2d,1,
+					&DLStokes[0][0],2,2,0,
+					&LStokes2[0][0],2,numdof2d,0,
+					&Ke_drag_gaussian[0][0],0);
+
+		for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof+j]+=Ke_drag_gaussian[i][j];
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+	//TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
+	
+	//printf("New:\n");
+	//printarray(&Ke_drag_gaussian[0][0],numdof2d,numdof2d);
+	//printf("Old:\n");
+	//Ke=CreateKMatrixDiagnosticStokesFriction();
+	//_error_("STOP");
+	//if(id==1) printarray(Ke->values,4,Ke->ncols);
+	//if(id==1) printarray(Ke->values,1,9);
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+
+	//return CreateKMatrixDiagnosticStokesFriction();
 	return Ke;
 }
Index: issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.h	(revision 10452)
+++ issm/trunk/src/c/objects/Elements/Penta.h	(revision 10480)
@@ -223,4 +223,5 @@
 		ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
 		ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
+		ElementMatrix* CreateKMatrixDiagnosticStokesFriction2(void);
 		ElementMatrix* CreateKMatrixDiagnosticVert(void);
 		ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
Index: issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 10452)
+++ issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 10480)
@@ -651,4 +651,141 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::GetLStokes2{{{1*/
+void PentaRef::GetLStokes2(double* LStokes, GaussPenta* gauss){
+	/*
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
+	 * by: 
+	 *       Li=[ h 0 0 0]
+	 *	 	      [ 0 h 0 0]
+	 *		      [ 0 0 h h]
+	 *		      [ 0 0 0 0]
+	 * where h is the interpolation function for node i.
+	 */
+
+	int i;
+	int num_dof=4;
+
+	double l1l2l3[NUMNODESP1_2d];
+
+
+	/*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*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0.;
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLStokes3{{{1*/
+void PentaRef::GetLStokes3(double* LStokes, GaussPenta* gauss){
+	/*
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
+	 * by: 
+	 *       Li=[ h 0 ]
+	 *	 	      [ 0 h ]
+	 *		      [ 0 0 ]
+	 *		      [ 0 0 ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	int i;
+	int num_dof=4;
+
+	double l1l2l3[NUMNODESP1_2d];
+
+
+	/*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*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLprimeStokes2{{{1*/
+void PentaRef::GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss){
+	/*
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
+	 * by: 
+	 *       Li=[ h 0 0 0]
+	 *	 	      [ 0 h 0 0]
+	 *		      [ 0 0 dh\dz 0]
+	 *		      [ 0 0 0 h]
+	 * where h is the interpolation function for node i.
+	 */
+
+	int i;
+	int num_dof=4;
+
+	double l1l2l3[NUMNODESP1_2d];
+	double dh1dh6[3][NUMNODESP1];
+
+
+	/*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;
+
+	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+
+	/*Build LStokes: */
+	for (i=0;i<3;i++){
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=dh1dh6[2][i];
+		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0.;
+
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+0)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0.;
+		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=l1l2l3[i];
+	}
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/
 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
@@ -658,18 +795,33 @@
 	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Lpi=[ 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  dh/dz 0]
-	 *		       [ 0    0  dh/dz 0]
-	 *		       [ 0    0  dh/dz 0]
-	 *		       [dh/dz 0  dh/dx 0]
-	 *		       [ 0 dh/dz dh/dy 0]
-	 *           [ 0    0    0   h]
-	 *           [ 0    0    0   h]
-	 *           [ 0    0    0   h]
+	 *       Lpi=[ h    0    0   0]1
+	 *		       [ 0    h    0   0]2
+	 *		       [ h    0    0   0]3
+	 *		       [ 0    h    0   0]4
+	 *		       [ 0    0    h   0]5
+	 *		       [ 0    0    h   0]6
+	 *		       [ 0    0  dh/dz 0]7
+	 *		       [ 0    0  dh/dz 0]8
+	 *		       [ 0    0  dh/dz 0]9
+	 *		       [dh/dz 0  dh/dx 0]0
+	 *		       [ 0 dh/dz dh/dy 0]1
+	 *           [ 0    0    0   h]2
+	 *           [ 0    0    0   h]3
+	 *           [ 0    0    0   h]4
+	 *
+	 *       Li=[ h    0    0   0]1
+	 *	 	      [ 0    h    0   0]2
+	 *		      [ 0    0    h   0]3
+	 *		      [ 0    0    h   0]4
+	 *	 	      [ h    0    0   0]5
+	 *	 	      [ 0    h    0   0]6
+	 *	 	      [ h    0    0   0]7
+	 *	 	      [ 0    h    0   0]8
+	 *		      [ 0    0    h   0]9
+	 *		      [ 0    0    h   0]0
+	 *		      [ 0    0    h   0]1
+	 *	 	      [ h    0    0   0]2
+	 *	 	      [ 0    h    0   0]3
+	 *		      [ 0    0    h   0]4
 	 * where h is the interpolation function for node i.
 	 */
Index: issm/trunk/src/c/objects/Elements/PentaRef.h
===================================================================
--- issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 10452)
+++ issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 10480)
@@ -51,4 +51,7 @@
 		void GetL(double* L, GaussPenta* gauss,int numdof);
 		void GetLStokes(double* LStokes, GaussPenta* gauss);
+		void GetLStokes2(double* LStokes, GaussPenta* gauss);
+		void GetLStokes3(double* LStokes, GaussPenta* gauss);
+		void GetLprimeStokes2(double* LStokes,double* xyz_list,GaussPenta* gauss);
 		void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
 		void GetLMacAyealStokes(double* LMacAyealStokes, GaussPenta* gauss);
Index: issm/trunk/src/c/shared/Elements/elements.h
===================================================================
--- issm/trunk/src/c/shared/Elements/elements.h	(revision 10452)
+++ issm/trunk/src/c/shared/Elements/elements.h	(revision 10480)
@@ -28,5 +28,5 @@
 	for(int i=0;i<lines;i++){  
 		printf("   [ ");
-		for(int j=0;j<cols;j++) printf(" %7.5g ",array[i*cols+j]);
+		for(int j=0;j<cols;j++) printf(" %12.7g ",array[i*cols+j]);
 		printf(" ]\n");
 	}  
