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;inrows*Ke->ncols*sizeof(double)); - for(i=0;inrows;i++) for(j=0;jncols;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;inrows;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;inrows*Ke->ncols*sizeof(double)); - for(i=0;inrows;i++) for(j=0;jncols;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;inrows;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;inrows;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;inrows;i++) for(j=0;jncols;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 -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;iGetCoordinateSystem(&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;inrows;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;inrows*Ke->ncols*sizeof(double)); + for(i=0;inrows;i++) for(j=0;jncols;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