Index: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am	(revision 10406)
+++ issm/trunk/src/c/Makefile.am	(revision 10407)
@@ -18,5 +18,5 @@
 #}}}
 
-#Our sources
+#sources
 #Core sources{{{1
 core_sources   = ./include/macros.h\
@@ -482,5 +482,8 @@
 					      ./modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp \
 					      ./modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
-					      ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp
+					      ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp \
+							./shared/Elements/TransformLoadVectorCoord.cpp \
+							./shared/Elements/TransformStiffnessMatrixCoord.cpp \
+							./shared/Elements/TransformSolutionCoord.cpp
 diagnostic_psources =./solutions/diagnostic_core.cpp\
 					      ./solvers/solver_stokescoupling_nonlinear.cpp
Index: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10406)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10407)
@@ -5651,5 +5651,5 @@
 
 	/*Transform Coordinate System*/
-	tria->TransformStiffnessMatrixCoord(Ke,2);
+	TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,NDOF2);
 
 	/*Clean up and return*/
@@ -5783,5 +5783,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -5857,5 +5857,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -5952,5 +5952,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,4);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
 
 	/*Clean up and return*/
@@ -6040,5 +6040,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,4);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
 
 	/*Clean up and return*/
@@ -6661,5 +6661,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -6757,5 +6757,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
 
 	/*Clean up and return*/
@@ -6825,5 +6825,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
 
 	/*Clean up and return*/
@@ -7185,5 +7185,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],&values0[0],2,true);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,NDOF2); /*2D: only the first 3 nodes are taken*/
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
@@ -7425,5 +7425,4 @@
 	double rho_ice,g;
 	double values[numdof];
-	double values0[numdof];
 	double vx[NUMVERTICES];
 	double vy[NUMVERTICES];
@@ -7442,8 +7441,8 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],&values0[0],2);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -7756,5 +7755,4 @@
 	int     i;
 	double  values[numdof];
-	double  values0[numdof];
 	double  vx[NUMVERTICES];
 	double  vy[NUMVERTICES];
@@ -7769,8 +7767,8 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],&values0[0],NDOF4);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF4);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
@@ -7811,76 +7809,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::TransformStiffnessMatrixCoord{{{1*/
-void Penta::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
-
-	int     i,j;
-	int     numnodes          = NUMVERTICES;
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current stiffness matrix*/
-	values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
-	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T*Ke*T^t */
-	TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
-				values,Ke->nrows,Ke->ncols,0,
-				transform,numnodes*dim,numnodes*dim,0,
-				&Ke->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
-/*FUNCTION Penta::TransformLoadVectorCoord{{{1*/
-void Penta::TransformLoadVectorCoord(ElementVector* pe,int dim){
-
-	int     i,j;
-	int     numnodes          = NUMVERTICES;
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current load vector*/
-	values=(double*)xmalloc(pe->nrows*sizeof(double));
-	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T^t*pe */
-	MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
-				values,pe->nrows,1,0,
-				&pe->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
-/*FUNCTION Penta::TransformSolutionCoord{{{1*/
-void Penta::TransformSolutionCoord(double* solution,double* solution0,int dim,bool is2d){
-
-	int     i,j;
-	int     numnodes;
-	double *transform = NULL;
-
-	/*Get Coordinate Systems transform matrix*/
-	if(is2d) numnodes=NUMVERTICES2D;
-	else     numnodes=NUMVERTICES;
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T*U */
-	MatrixMultiply(transform,numnodes*dim,numnodes*dim,0,
-				solution0,numnodes*dim,1,0,
-				&solution[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-}
-/*}}}*/
 #endif
 
Index: issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.h	(revision 10406)
+++ issm/trunk/src/c/objects/Elements/Penta.h	(revision 10407)
@@ -257,7 +257,4 @@
 		ElementVector* CreatePVectorDiagnosticVertVolume(void);
 		ElementVector* CreatePVectorDiagnosticVertBase(void);
-		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
-		void           TransformLoadVectorCoord(ElementVector* pe,int dim);
-		void           TransformSolutionCoord(double* solution,double* solution0,int dim,bool is2d=false);
 		#endif
 
Index: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10406)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10407)
@@ -2707,5 +2707,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -2775,5 +2775,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -2851,5 +2851,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
 
 	/*Clean up and return*/
@@ -2994,5 +2994,4 @@
 	int*      doflist=NULL;
 	double    rho_ice,g;
-	double    values0[numdof];
 	double    values[numdof];
 	double    vx[NUMVERTICES];
@@ -3007,8 +3006,8 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],&values0[0],2);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -3107,74 +3106,4 @@
 	/*Free ressources:*/
 	xfree((void**)&doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::TransformStiffnessMatrixCoord{{{1*/
-void Tria::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
-
-	int     i,j;
-	int     numnodes          = NUMVERTICES;
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current stiffness matrix*/
-	values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
-	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T*Ke*T^t */
-	TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
-				values,Ke->nrows,Ke->ncols,0,
-				transform,numnodes*dim,numnodes*dim,0,
-				&Ke->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
-/*FUNCTION Tria::TransformLoadVectorCoord{{{1*/
-void Tria::TransformLoadVectorCoord(ElementVector* pe,int dim){
-
-	int     i,j;
-	int     numnodes          = NUMVERTICES;
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current load vector*/
-	values=(double*)xmalloc(pe->nrows*sizeof(double));
-	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T^t*pe */
-	MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
-				  values,pe->nrows,1,0,
-				  &pe->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
-/*FUNCTION Tria::TransformSolutionCoord{{{1*/
-void Tria::TransformSolutionCoord(double* solution,double* solution0,int dim){
-
-	int     i,j;
-	int     numnodes          = NUMVERTICES;
-	double *transform         = NULL;
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dim);
-
-	/*Transform matrix: T*U */
-	MatrixMultiply(transform,numnodes*dim,numnodes*dim,0,
-				solution0,numnodes*dim,1,0,
-				&solution[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
 }
 /*}}}*/
Index: issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.h	(revision 10406)
+++ issm/trunk/src/c/objects/Elements/Tria.h	(revision 10407)
@@ -203,7 +203,4 @@
 		void	  InputUpdateFromSolutionDiagnosticHoriz( double* solution);
 		void	  InputUpdateFromSolutionDiagnosticHutter( double* solution);
-		void    TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
-		void    TransformLoadVectorCoord(ElementVector* pe,int dim);
-		void    TransformSolutionCoord(double* solution,double* solution0,int dim);
 		#endif
 
Index: issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 10406)
+++ issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 10407)
@@ -530,5 +530,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,NUMVERTICESSEG,2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,NDOF2);
 
 	/*Clean up and return*/
@@ -537,28 +537,6 @@
 }
 /*}}}*/
-/*FUNCTION Icefront::TransformLoadVectorCoord{{{1*/
-void Icefront::TransformLoadVectorCoord(ElementVector* pe,int numnodes,int dimension){
-
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current load matrix*/
-	values=(double*)xmalloc(pe->nrows*sizeof(double));
-	for(int i=0;i<pe->nrows;i++) values[i]=pe->values[i];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dimension);
-
-	/*Transform matrix: T*pe */
-	MatrixMultiply(transform,numnodes*dimension,numnodes*dimension,1,
-				values,pe->nrows,1,0,
-				&pe->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
 #endif
+
 #ifdef _HAVE_CONTROL_
 /*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{1*/
@@ -669,5 +647,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,NUMVERTICESQUA,2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF2);
 
 	/*Clean up and return*/
@@ -747,5 +725,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,NUMVERTICESQUA,4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF4);
 
 	/*Clean up and return*/
Index: issm/trunk/src/c/objects/Loads/Icefront.h
===================================================================
--- issm/trunk/src/c/objects/Loads/Icefront.h	(revision 10406)
+++ issm/trunk/src/c/objects/Loads/Icefront.h	(revision 10407)
@@ -83,5 +83,4 @@
 		void GetSegmentNormal(double* normal,double xyz_list[2][3]);
 		void GetQuadNormal(double* normal,double xyz_list[4][3]);
-		void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int dimension);
 		#ifdef _HAVE_CONTROL_
 		ElementVector* CreatePVectorAdjointHoriz(void);
Index: issm/trunk/src/c/objects/Loads/Pengrid.cpp
===================================================================
--- issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 10406)
+++ issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 10407)
@@ -560,5 +560,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,NDOF4);
+	TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,NDOF4);
 
 	/*Clean up and return*/
@@ -691,30 +691,4 @@
 }
 /*}}}1*/
-/*FUNCTION Pengrid::TransformStiffnessMatrixCoord{{{1*/
-void Pengrid::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
-
-	int     i,j;
-	int     numnodes          = 1;
-	double *transform         = NULL;
-	double *values            = NULL;
-
-	/*Copy current stiffness matrix*/
-	values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
-	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
-
-	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,&node,numnodes,dim);
-
-	/*Transform matrix: T*Ke*T^t */
-	TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
-				values,Ke->nrows,Ke->ncols,0,
-				transform,numnodes*dim,numnodes*dim,0,
-				&Ke->values[0],0);
-
-	/*Free Matrix*/
-	xfree((void**)&transform);
-	xfree((void**)&values);
-}
-/*}}}*/
 /*FUNCTION Pengrid::UpdateInputs {{{1*/
 void  Pengrid::UpdateInputs(double* solution){
Index: issm/trunk/src/c/objects/Loads/Pengrid.h
===================================================================
--- issm/trunk/src/c/objects/Loads/Pengrid.h	(revision 10406)
+++ issm/trunk/src/c/objects/Loads/Pengrid.h	(revision 10407)
@@ -90,5 +90,4 @@
 		void  UpdateInputs(double* solution);
 		void  ResetConstraint(void);
-		void  TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
 		/*}}}*/
 
Index: issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp
===================================================================
--- issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp	(revision 10406)
+++ issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp	(revision 10407)
@@ -5,5 +5,5 @@
 #include <math.h>
 
-void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension){
+void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dofspernode){
 
 	int     i;
@@ -15,9 +15,9 @@
 	/*Some checks*/
 	_assert_(numnodes && nodes);
-	if(dimension!=2 && dimension!=3 && dimension!=4) _error_("only 2d and 3d coordinate systems are supported");
+	if(dofspernode!=2 && dofspernode!=3 && dofspernode!=4) _error_("list of dofs per node supported: 2, 3 and 4");
 
 	/*Allocate and initialize transform matrix*/
-	transform=(double*)xmalloc(numnodes*dimension*numnodes*dimension*sizeof(double));
-	for(i=0;i<numnodes*dimension*numnodes*dimension;i++) transform[i]=0.0;
+	transform=(double*)xmalloc(numnodes*dofspernode*numnodes*dofspernode*sizeof(double));
+	for(i=0;i<numnodes*dofspernode*numnodes*dofspernode;i++) transform[i]=0.0;
 
 	/*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
@@ -32,36 +32,36 @@
 	for(i=0;i<numnodes;i++){
 		nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
-		if(dimension==2){
+		if(dofspernode==2){
 			/*We remove the z component, we need to renormalize x and y*/
 			x_norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]);
 			y_norm = sqrt( coord_system[0][1]*coord_system[0][1] + coord_system[1][1]*coord_system[1][1]);
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0]/x_norm;
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1]/y_norm;
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0]/x_norm;
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1]/y_norm;
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0]/x_norm;
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1]/y_norm;
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0]/x_norm;
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1]/y_norm;
 		}
-		else if(dimension==3){
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0];
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1];
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+2] = coord_system[0][2];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+2] = coord_system[1][2];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+0] = coord_system[2][0];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+1] = coord_system[2][1];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+2] = coord_system[2][2];
+		else if(dofspernode==3){
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+2] = coord_system[0][2];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+2] = coord_system[1][2];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+0] = coord_system[2][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+1] = coord_system[2][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+2] = coord_system[2][2];
 		}
-		else if(dimension==4){
+		else if(dofspernode==4){
 			/*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0];
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1];
-			transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+2] = coord_system[0][2];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1];
-			transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+2] = coord_system[1][2];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+0] = coord_system[2][0];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+1] = coord_system[2][1];
-			transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+2] = coord_system[2][2];
-			transform[(dimension*numnodes)*(i*dimension+3) + i*dimension+3] = 1.0;
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+2] = coord_system[0][2];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+2] = coord_system[1][2];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+0] = coord_system[2][0];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+1] = coord_system[2][1];
+			transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+2] = coord_system[2][2];
+			transform[(dofspernode*numnodes)*(i*dofspernode+3) + i*dofspernode+3] = 1.0;
 		}
 	}
Index: issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp
===================================================================
--- issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp	(revision 10407)
+++ issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp	(revision 10407)
@@ -0,0 +1,27 @@
+/*!\file:  TransformLoadVectorCoord.cpp
+ * \brief transform load vector coordinate system
+ */ 
+#include "./elements.h"
+
+void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int dofspernode){
+
+	int     i,j;
+	double *transform = NULL;
+	double *values    = NULL;
+
+	/*Copy current load vector*/
+	values=(double*)xmalloc(pe->nrows*sizeof(double));
+	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
+
+	/*Get Coordinate Systems transform matrix*/
+	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+
+	/*Transform matrix: R^T*pe */
+	MatrixMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,1,
+				values,pe->nrows,1,0,
+				&pe->values[0],0);
+
+	/*Free Matrix*/
+	xfree((void**)&transform);
+	xfree((void**)&values);
+}
Index: issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp
===================================================================
--- issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp	(revision 10407)
+++ issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp	(revision 10407)
@@ -0,0 +1,27 @@
+/*!\file:  TransformSolutionCoord.cpp
+ * \brief transform solution vector coordinate system
+ */ 
+#include "./elements.h"
+
+void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int dofspernode){
+
+	int     i,j;
+	double *transform = NULL;
+	double *values    = NULL;
+
+	/*Copy current solution vector*/
+	values=(double*)xmalloc(numnodes*dofspernode*sizeof(double));
+	for(i=0;i<numnodes*dofspernode;i++) values[i]=solution[i];
+
+	/*Get Coordinate Systems transform matrix*/
+	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+
+	/*Transform matrix: R*U */
+	MatrixMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,0,
+				values,numnodes*dofspernode,1,0,
+				&solution[0],0);
+
+	/*Free Matrix*/
+	xfree((void**)&transform);
+	xfree((void**)&values);
+}
Index: issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp
===================================================================
--- issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp	(revision 10407)
+++ issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp	(revision 10407)
@@ -0,0 +1,28 @@
+/*!\file:  TransformStiffnessMatrixCoord.cpp
+ * \brief transform stiffness matrix coordinate system
+ */ 
+#include "./elements.h"
+
+void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int dofspernode){
+
+	int     i,j;
+	double *transform = NULL;
+	double *values    = NULL;
+
+	/*Copy current stiffness matrix*/
+	values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
+	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
+
+	/*Get Coordinate Systems transform matrix*/
+	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+
+	/*Transform matrix: R^T*Ke*R */
+	TripleMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,1,
+				values,Ke->nrows,Ke->ncols,0,
+				transform,numnodes*dofspernode,numnodes*dofspernode,0,
+				&Ke->values[0],0);
+
+	/*Free Matrix*/
+	xfree((void**)&transform);
+	xfree((void**)&values);
+}
Index: issm/trunk/src/c/shared/Elements/elements.h
===================================================================
--- issm/trunk/src/c/shared/Elements/elements.h	(revision 10406)
+++ issm/trunk/src/c/shared/Elements/elements.h	(revision 10407)
@@ -17,5 +17,10 @@
 int*   GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum);
 int*   GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation_enum);
-void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension);
+void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dofspernode);
+#ifdef _HAVE_DIAGNOSTIC_
+void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int dofspernode);
+void   TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int dofspernode);
+void   TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int dofspernode);
+#endif
 
 inline void printarray(double* array,int lines,int cols=1){
@@ -23,5 +28,5 @@
 	for(int i=0;i<lines;i++){  
 		printf("   [ ");
-		for(int j=0;j<cols;j++) printf(" %12.6g ",array[i*cols+j]);
+		for(int j=0;j<cols;j++) printf(" %7.5g ",array[i*cols+j]);
 		printf(" ]\n");
 	}  
