Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 10522)
+++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 10523)
@@ -445,4 +445,8 @@
 	NearestInterpEnum,
 	/*}}}*/
+	/*Coordinate Systems{{{1*/
+	XYEnum,
+	XYZPEnum,
+	/*}}}*/
 	/*Options{{{1*/
 	OptionEnum,
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10522)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 10523)
@@ -5745,5 +5745,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,NDOF2);
+	TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum);
 
 	/*Clean up and return*/
@@ -5877,5 +5877,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -5951,5 +5951,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -6043,5 +6043,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
 
 	/*Clean up and return*/
@@ -6114,5 +6114,5 @@
 
 	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
-	//TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
+	//TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
 	
 	/*Clean up and return*/
@@ -6735,5 +6735,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -6829,5 +6829,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
 
 	/*Clean up and return*/
@@ -6897,5 +6897,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
 
 	/*Clean up and return*/
@@ -7256,5 +7256,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,NDOF2); /*2D: only the first 3 nodes are taken*/
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
@@ -7515,5 +7515,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -7841,5 +7841,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF4);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYZPEnum);
 
 	/*Ok, we have vx and vy in values, fill in all arrays: */
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10522)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 10523)
@@ -2707,5 +2707,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -2775,5 +2775,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
+	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -2851,5 +2851,5 @@
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
 
 	/*Clean up and return*/
@@ -3009,5 +3009,5 @@
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
+	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 10522)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 10523)
@@ -530,5 +530,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,NDOF2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,XYEnum);
 
 	/*Clean up and return*/
@@ -647,5 +647,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF2);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum);
 
 	/*Clean up and return*/
@@ -725,5 +725,5 @@
 
 	/*Transform load vector*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF4);
+	TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum);
 
 	/*Clean up and return*/
Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 10522)
+++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp	(revision 10523)
@@ -560,5 +560,5 @@
 
 	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,NDOF4);
+	TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,XYZPEnum);
 
 	/*Clean up and return*/
Index: /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp	(revision 10522)
+++ /issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp	(revision 10523)
@@ -5,7 +5,8 @@
 #include <math.h>
 
-void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dofspernode){
+void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int* cs_array){
 
-	int     i;
+	int     i,counter;
+	int     numdofs           = 0;
 	double  x_norm,y_norm;
 	double *transform         = NULL;
@@ -13,11 +14,19 @@
 	double  coord_system[3][3];
 
-	/*Some checks*/
+	/*Some checks in debugging mode*/
 	_assert_(numnodes && nodes);
-	if(dofspernode!=2 && dofspernode!=3 && dofspernode!=4) _error_("list of dofs per node supported: 2, 3 and 4");
+
+	/*Get total number of dofs*/
+	for(i=0;i<numnodes;i++){
+		switch(cs_array[i]){
+			case XYEnum:   numdofs+=2; break;
+			case XYZPEnum: numdofs+=4; break;
+			default: _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
+		}
+	}
 
 	/*Allocate and initialize transform matrix*/
-	transform=(double*)xmalloc(numnodes*dofspernode*numnodes*dofspernode*sizeof(double));
-	for(i=0;i<numnodes*dofspernode*numnodes*dofspernode;i++) transform[i]=0.0;
+	transform=(double*)xmalloc(numdofs*numdofs*sizeof(double));
+	for(i=0;i<numdofs*numdofs;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
@@ -30,38 +39,34 @@
 	 * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
 	 * associated to this node*/
+	counter=0;
 	for(i=0;i<numnodes;i++){
 		nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
-		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[(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(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(dofspernode==4){
-			/*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
-			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;
+		switch(cs_array[i]){
+			case XYEnum:
+				/*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[(numdofs)*(counter+0) + counter+0] = coord_system[0][0]/x_norm;
+				transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1]/y_norm;
+				transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0]/x_norm;
+				transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1]/y_norm;
+				counter+=2;
+				break;
+			case XYZPEnum:
+				/*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
+				transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
+				transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
+				transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
+				transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
+				transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
+				transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
+				transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
+				transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
+				transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
+				transform[(numdofs)*(counter+3) + counter+3] = 1.0;
+				counter+=4;
+				break;
+			default:
+				_error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
 		}
 	}
Index: /issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp	(revision 10522)
+++ /issm/trunk/src/c/shared/Elements/TransformLoadVectorCoord.cpp	(revision 10523)
@@ -4,9 +4,33 @@
 #include "./elements.h"
 
-void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int dofspernode){
+void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum){
+	int* cs_array=NULL;
+
+	/*All nodes have the same Coordinate System*/
+	cs_array=(int*)xmalloc(numnodes*sizeof(int));
+	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
+
+	/*Call core*/
+	TransformLoadVectorCoord(pe,nodes,numnodes,cs_array);
+
+	/*Clean-up*/
+	xfree((void**)&cs_array);
+}
+
+void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array){
 
 	int     i,j;
+	int     numdofs   = 0;
 	double *transform = NULL;
 	double *values    = NULL;
+
+	/*Get total number of dofs*/
+	for(i=0;i<numnodes;i++){
+		switch(cs_array[i]){
+			case XYEnum:   numdofs+=2; break;
+			case XYZPEnum: numdofs+=4; break;
+			default: _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
+		}
+	}
 
 	/*Copy current load vector*/
@@ -15,8 +39,8 @@
 
 	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
 
 	/*Transform matrix: R^T*pe */
-	MatrixMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,1,
+	MatrixMultiply(transform,numdofs,numdofs,1,
 				values,pe->nrows,1,0,
 				&pe->values[0],0);
Index: /issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp	(revision 10522)
+++ /issm/trunk/src/c/shared/Elements/TransformSolutionCoord.cpp	(revision 10523)
@@ -4,20 +4,45 @@
 #include "./elements.h"
 
-void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int dofspernode){
+void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int cs_enum){
+
+	int* cs_array=NULL;
+
+	/*All nodes have the same Coordinate System*/
+	cs_array=(int*)xmalloc(numnodes*sizeof(int));
+	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
+
+	/*Call core*/
+	TransformSolutionCoord(solution,nodes,numnodes,cs_array);
+
+	/*Clean-up*/
+	xfree((void**)&cs_array);
+}
+
+void TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int* cs_array){
 
 	int     i,j;
+	int     numdofs   = 0;
 	double *transform = NULL;
 	double *values    = NULL;
 
+	/*Get total number of dofs*/
+	for(i=0;i<numnodes;i++){
+		switch(cs_array[i]){
+			case XYEnum:   numdofs+=2; break;
+			case XYZPEnum: numdofs+=4; break;
+			default: _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
+		}
+	}
+
 	/*Copy current solution vector*/
-	values=(double*)xmalloc(numnodes*dofspernode*sizeof(double));
-	for(i=0;i<numnodes*dofspernode;i++) values[i]=solution[i];
+	values=(double*)xmalloc(numdofs*sizeof(double));
+	for(i=0;i<numdofs;i++) values[i]=solution[i];
 
 	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
 
 	/*Transform matrix: R*U */
-	MatrixMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,0,
-				values,numnodes*dofspernode,1,0,
+	MatrixMultiply(transform,numdofs,numdofs,0,
+				values,numdofs,1,0,
 				&solution[0],0);
 
Index: /issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp
===================================================================
--- /issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp	(revision 10522)
+++ /issm/trunk/src/c/shared/Elements/TransformStiffnessMatrixCoord.cpp	(revision 10523)
@@ -4,9 +4,34 @@
 #include "./elements.h"
 
-void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int dofspernode){
+void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){
+
+	int* cs_array=NULL;
+
+	/*All nodes have the same Coordinate System*/
+	cs_array=(int*)xmalloc(numnodes*sizeof(int));
+	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
+
+	/*Call core*/
+	TransformStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);
+
+	/*Clean-up*/
+	xfree((void**)&cs_array);
+}
+
+void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){
 
 	int     i,j;
+	int     numdofs   = 0;
 	double *transform = NULL;
 	double *values    = NULL;
+
+	/*Get total number of dofs*/
+	for(i=0;i<numnodes;i++){
+		switch(cs_array[i]){
+			case XYEnum:   numdofs+=2; break;
+			case XYZPEnum: numdofs+=4; break;
+			default: _error_("Coordinate system %s not supported yet",EnumToStringx(cs_array[i]));
+		}
+	}
 
 	/*Copy current stiffness matrix*/
@@ -15,10 +40,10 @@
 
 	/*Get Coordinate Systems transform matrix*/
-	CoordinateSystemTransform(&transform,nodes,numnodes,dofspernode);
+	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
 
 	/*Transform matrix: R^T*Ke*R */
-	TripleMultiply(transform,numnodes*dofspernode,numnodes*dofspernode,1,
+	TripleMultiply(transform,numdofs,numdofs,1,
 				values,Ke->nrows,Ke->ncols,0,
-				transform,numnodes*dofspernode,numnodes*dofspernode,0,
+				transform,numdofs,numdofs,0,
 				&Ke->values[0],0);
 
Index: /issm/trunk/src/c/shared/Elements/elements.h
===================================================================
--- /issm/trunk/src/c/shared/Elements/elements.h	(revision 10522)
+++ /issm/trunk/src/c/shared/Elements/elements.h	(revision 10523)
@@ -17,9 +17,12 @@
 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 dofspernode);
+void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int* cs_array);
 #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);
+void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
+void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
+void   TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum);
+void   TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array);
+void   TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int cs_enum);
+void   TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int* cs_array);
 #endif