Changeset 10407


Ignore:
Timestamp:
11/01/11 11:40:58 (13 years ago)
Author:
Mathieu Morlighem
Message:

Simplified Coordinate system transformation

Location:
issm/trunk/src/c
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r10370 r10407  
    1818#}}}
    1919
    20 #Our sources
     20#sources
    2121#Core sources{{{1
    2222core_sources   = ./include/macros.h\
    … …  
    482482                                              ./modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp \
    483483                                              ./modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
    484                                               ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp
     484                                              ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp \
     485                                                        ./shared/Elements/TransformLoadVectorCoord.cpp \
     486                                                        ./shared/Elements/TransformStiffnessMatrixCoord.cpp \
     487                                                        ./shared/Elements/TransformSolutionCoord.cpp
    485488diagnostic_psources =./solutions/diagnostic_core.cpp\
    486489                                              ./solvers/solver_stokescoupling_nonlinear.cpp

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r10400 r10407  
    56515651
    56525652        /*Transform Coordinate System*/
    5653         tria->TransformStiffnessMatrixCoord(Ke,2);
     5653        TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,NDOF2);
    56545654
    56555655        /*Clean up and return*/
    … …  
    57835783
    57845784        /*Transform Coordinate System*/
    5785         TransformStiffnessMatrixCoord(Ke,2);
     5785        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
    57865786
    57875787        /*Clean up and return*/
    … …  
    58575857
    58585858        /*Transform Coordinate System*/
    5859         TransformStiffnessMatrixCoord(Ke,2);
     5859        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
    58605860
    58615861        /*Clean up and return*/
    … …  
    59525952
    59535953        /*Transform Coordinate System*/
    5954         TransformStiffnessMatrixCoord(Ke,4);
     5954        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
    59555955
    59565956        /*Clean up and return*/
    … …  
    60406040
    60416041        /*Transform Coordinate System*/
    6042         TransformStiffnessMatrixCoord(Ke,4);
     6042        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF4);
    60436043
    60446044        /*Clean up and return*/
    … …  
    66616661
    66626662        /*Transform coordinate system*/
    6663         TransformLoadVectorCoord(pe,2);
     6663        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
    66646664
    66656665        /*Clean up and return*/
    … …  
    67576757
    67586758        /*Transform coordinate system*/
    6759         TransformLoadVectorCoord(pe,4);
     6759        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
    67606760
    67616761        /*Clean up and return*/
    … …  
    68256825
    68266826        /*Transform coordinate system*/
    6827         TransformLoadVectorCoord(pe,4);
     6827        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF4);
    68286828
    68296829        /*Clean up and return*/
    … …  
    71857185
    71867186        /*Transform solution in Cartesian Space*/
    7187         TransformSolutionCoord(&values[0],&values0[0],2,true);
     7187        TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,NDOF2); /*2D: only the first 3 nodes are taken*/
    71887188
    71897189        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    … …  
    74257425        double rho_ice,g;
    74267426        double values[numdof];
    7427         double values0[numdof];
    74287427        double vx[NUMVERTICES];
    74297428        double vy[NUMVERTICES];
    … …  
    74427441
    74437442        /*Use the dof list to index into the solution vector: */
    7444         for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
     7443        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    74457444
    74467445        /*Transform solution in Cartesian Space*/
    7447         TransformSolutionCoord(&values[0],&values0[0],2);
     7446        TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
    74487447
    74497448        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    77567755        int     i;
    77577756        double  values[numdof];
    7758         double  values0[numdof];
    77597757        double  vx[NUMVERTICES];
    77607758        double  vy[NUMVERTICES];
    … …  
    77697767
    77707768        /*Use the dof list to index into the solution vector: */
    7771         for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
     7769        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    77727770
    77737771        /*Transform solution in Cartesian Space*/
    7774         TransformSolutionCoord(&values[0],&values0[0],NDOF4);
     7772        TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF4);
    77757773
    77767774        /*Ok, we have vx and vy in values, fill in all arrays: */
    … …  
    78117809}
    78127810/*}}}*/
    7813 /*FUNCTION Penta::TransformStiffnessMatrixCoord{{{1*/
    7814 void Penta::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
    7815 
    7816         int     i,j;
    7817         int     numnodes          = NUMVERTICES;
    7818         double *transform         = NULL;
    7819         double *values            = NULL;
    7820 
    7821         /*Copy current stiffness matrix*/
    7822         values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
    7823         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];
    7824 
    7825         /*Get Coordinate Systems transform matrix*/
    7826         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    7827 
    7828         /*Transform matrix: T*Ke*T^t */
    7829         TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
    7830                                 values,Ke->nrows,Ke->ncols,0,
    7831                                 transform,numnodes*dim,numnodes*dim,0,
    7832                                 &Ke->values[0],0);
    7833 
    7834         /*Free Matrix*/
    7835         xfree((void**)&transform);
    7836         xfree((void**)&values);
    7837 }
    7838 /*}}}*/
    7839 /*FUNCTION Penta::TransformLoadVectorCoord{{{1*/
    7840 void Penta::TransformLoadVectorCoord(ElementVector* pe,int dim){
    7841 
    7842         int     i,j;
    7843         int     numnodes          = NUMVERTICES;
    7844         double *transform         = NULL;
    7845         double *values            = NULL;
    7846 
    7847         /*Copy current load vector*/
    7848         values=(double*)xmalloc(pe->nrows*sizeof(double));
    7849         for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
    7850 
    7851         /*Get Coordinate Systems transform matrix*/
    7852         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    7853 
    7854         /*Transform matrix: T^t*pe */
    7855         MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
    7856                                 values,pe->nrows,1,0,
    7857                                 &pe->values[0],0);
    7858 
    7859         /*Free Matrix*/
    7860         xfree((void**)&transform);
    7861         xfree((void**)&values);
    7862 }
    7863 /*}}}*/
    7864 /*FUNCTION Penta::TransformSolutionCoord{{{1*/
    7865 void Penta::TransformSolutionCoord(double* solution,double* solution0,int dim,bool is2d){
    7866 
    7867         int     i,j;
    7868         int     numnodes;
    7869         double *transform = NULL;
    7870 
    7871         /*Get Coordinate Systems transform matrix*/
    7872         if(is2d) numnodes=NUMVERTICES2D;
    7873         else     numnodes=NUMVERTICES;
    7874         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    7875 
    7876         /*Transform matrix: T*U */
    7877         MatrixMultiply(transform,numnodes*dim,numnodes*dim,0,
    7878                                 solution0,numnodes*dim,1,0,
    7879                                 &solution[0],0);
    7880 
    7881         /*Free Matrix*/
    7882         xfree((void**)&transform);
    7883 }
    7884 /*}}}*/
    78857811#endif
    78867812

  • issm/trunk/src/c/objects/Elements/Penta.h

    r10377 r10407  
    257257                ElementVector* CreatePVectorDiagnosticVertVolume(void);
    258258                ElementVector* CreatePVectorDiagnosticVertBase(void);
    259                 void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
    260                 void           TransformLoadVectorCoord(ElementVector* pe,int dim);
    261                 void           TransformSolutionCoord(double* solution,double* solution0,int dim,bool is2d=false);
    262259                #endif
    263260

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r10400 r10407  
    27072707
    27082708        /*Transform Coordinate System*/
    2709         TransformStiffnessMatrixCoord(Ke,2);
     2709        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
    27102710
    27112711        /*Clean up and return*/
    … …  
    27752775
    27762776        /*Transform Coordinate System*/
    2777         TransformStiffnessMatrixCoord(Ke,2);
     2777        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,NDOF2);
    27782778
    27792779        /*Clean up and return*/
    … …  
    28512851
    28522852        /*Transform coordinate system*/
    2853         TransformLoadVectorCoord(pe,2);
     2853        TransformLoadVectorCoord(pe,nodes,NUMVERTICES,NDOF2);
    28542854
    28552855        /*Clean up and return*/
    … …  
    29942994        int*      doflist=NULL;
    29952995        double    rho_ice,g;
    2996         double    values0[numdof];
    29972996        double    values[numdof];
    29982997        double    vx[NUMVERTICES];
    … …  
    30073006
    30083007        /*Use the dof list to index into the solution vector: */
    3009         for(i=0;i<numdof;i++) values0[i]=solution[doflist[i]];
     3008        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    30103009
    30113010        /*Transform solution in Cartesian Space*/
    3012         TransformSolutionCoord(&values[0],&values0[0],2);
     3011        TransformSolutionCoord(&values[0],nodes,NUMVERTICES,NDOF2);
    30133012
    30143013        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    31073106        /*Free ressources:*/
    31083107        xfree((void**)&doflist);
    3109 }
    3110 /*}}}*/
    3111 /*FUNCTION Tria::TransformStiffnessMatrixCoord{{{1*/
    3112 void Tria::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
    3113 
    3114         int     i,j;
    3115         int     numnodes          = NUMVERTICES;
    3116         double *transform         = NULL;
    3117         double *values            = NULL;
    3118 
    3119         /*Copy current stiffness matrix*/
    3120         values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
    3121         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];
    3122 
    3123         /*Get Coordinate Systems transform matrix*/
    3124         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    3125 
    3126         /*Transform matrix: T*Ke*T^t */
    3127         TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
    3128                                 values,Ke->nrows,Ke->ncols,0,
    3129                                 transform,numnodes*dim,numnodes*dim,0,
    3130                                 &Ke->values[0],0);
    3131 
    3132         /*Free Matrix*/
    3133         xfree((void**)&transform);
    3134         xfree((void**)&values);
    3135 }
    3136 /*}}}*/
    3137 /*FUNCTION Tria::TransformLoadVectorCoord{{{1*/
    3138 void Tria::TransformLoadVectorCoord(ElementVector* pe,int dim){
    3139 
    3140         int     i,j;
    3141         int     numnodes          = NUMVERTICES;
    3142         double *transform         = NULL;
    3143         double *values            = NULL;
    3144 
    3145         /*Copy current load vector*/
    3146         values=(double*)xmalloc(pe->nrows*sizeof(double));
    3147         for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
    3148 
    3149         /*Get Coordinate Systems transform matrix*/
    3150         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    3151 
    3152         /*Transform matrix: T^t*pe */
    3153         MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
    3154                                   values,pe->nrows,1,0,
    3155                                   &pe->values[0],0);
    3156 
    3157         /*Free Matrix*/
    3158         xfree((void**)&transform);
    3159         xfree((void**)&values);
    3160 }
    3161 /*}}}*/
    3162 /*FUNCTION Tria::TransformSolutionCoord{{{1*/
    3163 void Tria::TransformSolutionCoord(double* solution,double* solution0,int dim){
    3164 
    3165         int     i,j;
    3166         int     numnodes          = NUMVERTICES;
    3167         double *transform         = NULL;
    3168 
    3169         /*Get Coordinate Systems transform matrix*/
    3170         CoordinateSystemTransform(&transform,nodes,numnodes,dim);
    3171 
    3172         /*Transform matrix: T*U */
    3173         MatrixMultiply(transform,numnodes*dim,numnodes*dim,0,
    3174                                 solution0,numnodes*dim,1,0,
    3175                                 &solution[0],0);
    3176 
    3177         /*Free Matrix*/
    3178         xfree((void**)&transform);
    31793108}
    31803109/*}}}*/

  • issm/trunk/src/c/objects/Elements/Tria.h

    r10377 r10407  
    203203                void      InputUpdateFromSolutionDiagnosticHoriz( double* solution);
    204204                void      InputUpdateFromSolutionDiagnosticHutter( double* solution);
    205                 void    TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
    206                 void    TransformLoadVectorCoord(ElementVector* pe,int dim);
    207                 void    TransformSolutionCoord(double* solution,double* solution0,int dim);
    208205                #endif
    209206

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r10391 r10407  
    530530
    531531        /*Transform load vector*/
    532         TransformLoadVectorCoord(pe,NUMVERTICESSEG,2);
     532        TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,NDOF2);
    533533
    534534        /*Clean up and return*/
    … …  
    537537}
    538538/*}}}*/
    539 /*FUNCTION Icefront::TransformLoadVectorCoord{{{1*/
    540 void Icefront::TransformLoadVectorCoord(ElementVector* pe,int numnodes,int dimension){
    541 
    542         double *transform         = NULL;
    543         double *values            = NULL;
    544 
    545         /*Copy current load matrix*/
    546         values=(double*)xmalloc(pe->nrows*sizeof(double));
    547         for(int i=0;i<pe->nrows;i++) values[i]=pe->values[i];
    548 
    549         /*Get Coordinate Systems transform matrix*/
    550         CoordinateSystemTransform(&transform,nodes,numnodes,dimension);
    551 
    552         /*Transform matrix: T*pe */
    553         MatrixMultiply(transform,numnodes*dimension,numnodes*dimension,1,
    554                                 values,pe->nrows,1,0,
    555                                 &pe->values[0],0);
    556 
    557         /*Free Matrix*/
    558         xfree((void**)&transform);
    559         xfree((void**)&values);
    560 }
    561 /*}}}*/
    562539#endif
     540
    563541#ifdef _HAVE_CONTROL_
    564542/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{1*/
    … …  
    669647
    670648        /*Transform load vector*/
    671         TransformLoadVectorCoord(pe,NUMVERTICESQUA,2);
     649        TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF2);
    672650
    673651        /*Clean up and return*/
    … …  
    747725
    748726        /*Transform load vector*/
    749         TransformLoadVectorCoord(pe,NUMVERTICESQUA,4);
     727        TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,NDOF4);
    750728
    751729        /*Clean up and return*/

  • issm/trunk/src/c/objects/Loads/Icefront.h

    r10381 r10407  
    8383                void GetSegmentNormal(double* normal,double xyz_list[2][3]);
    8484                void GetQuadNormal(double* normal,double xyz_list[4][3]);
    85                 void TransformLoadVectorCoord(ElementVector* pe,int numnodes,int dimension);
    8685                #ifdef _HAVE_CONTROL_
    8786                ElementVector* CreatePVectorAdjointHoriz(void);

  • issm/trunk/src/c/objects/Loads/Pengrid.cpp

    r10395 r10407  
    560560
    561561        /*Transform Coordinate System*/
    562         TransformStiffnessMatrixCoord(Ke,NDOF4);
     562        TransformStiffnessMatrixCoord(Ke,&node,NUMVERTICES,NDOF4);
    563563
    564564        /*Clean up and return*/
    … …  
    691691}
    692692/*}}}1*/
    693 /*FUNCTION Pengrid::TransformStiffnessMatrixCoord{{{1*/
    694 void Pengrid::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
    695 
    696         int     i,j;
    697         int     numnodes          = 1;
    698         double *transform         = NULL;
    699         double *values            = NULL;
    700 
    701         /*Copy current stiffness matrix*/
    702         values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
    703         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];
    704 
    705         /*Get Coordinate Systems transform matrix*/
    706         CoordinateSystemTransform(&transform,&node,numnodes,dim);
    707 
    708         /*Transform matrix: T*Ke*T^t */
    709         TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
    710                                 values,Ke->nrows,Ke->ncols,0,
    711                                 transform,numnodes*dim,numnodes*dim,0,
    712                                 &Ke->values[0],0);
    713 
    714         /*Free Matrix*/
    715         xfree((void**)&transform);
    716         xfree((void**)&values);
    717 }
    718 /*}}}*/
    719693/*FUNCTION Pengrid::UpdateInputs {{{1*/
    720694void  Pengrid::UpdateInputs(double* solution){

  • issm/trunk/src/c/objects/Loads/Pengrid.h

    r10395 r10407  
    9090                void  UpdateInputs(double* solution);
    9191                void  ResetConstraint(void);
    92                 void  TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
    9392                /*}}}*/
    9493

  • issm/trunk/src/c/shared/Elements/CoordinateSystemTransform.cpp

    r10393 r10407  
    55#include <math.h>
    66
    7 void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension){
     7void CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dofspernode){
    88
    99        int     i;
    … …  
    1515        /*Some checks*/
    1616        _assert_(numnodes && nodes);
    17         if(dimension!=2 && dimension!=3 && dimension!=4) _error_("only 2d and 3d coordinate systems are supported");
     17        if(dofspernode!=2 && dofspernode!=3 && dofspernode!=4) _error_("list of dofs per node supported: 2, 3 and 4");
    1818
    1919        /*Allocate and initialize transform matrix*/
    20         transform=(double*)xmalloc(numnodes*dimension*numnodes*dimension*sizeof(double));
    21         for(i=0;i<numnodes*dimension*numnodes*dimension;i++) transform[i]=0.0;
     20        transform=(double*)xmalloc(numnodes*dofspernode*numnodes*dofspernode*sizeof(double));
     21        for(i=0;i<numnodes*dofspernode*numnodes*dofspernode;i++) transform[i]=0.0;
    2222
    2323        /*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
    … …  
    3232        for(i=0;i<numnodes;i++){
    3333                nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
    34                 if(dimension==2){
     34                if(dofspernode==2){
    3535                        /*We remove the z component, we need to renormalize x and y*/
    3636                        x_norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]);
    3737                        y_norm = sqrt( coord_system[0][1]*coord_system[0][1] + coord_system[1][1]*coord_system[1][1]);
    38                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0]/x_norm;
    39                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1]/y_norm;
    40                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0]/x_norm;
    41                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1]/y_norm;
     38                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0]/x_norm;
     39                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1]/y_norm;
     40                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0]/x_norm;
     41                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1]/y_norm;
    4242                }
    43                 else if(dimension==3){
    44                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0];
    45                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1];
    46                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+2] = coord_system[0][2];
    47                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0];
    48                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1];
    49                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+2] = coord_system[1][2];
    50                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+0] = coord_system[2][0];
    51                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+1] = coord_system[2][1];
    52                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+2] = coord_system[2][2];
     43                else if(dofspernode==3){
     44                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0];
     45                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1];
     46                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+2] = coord_system[0][2];
     47                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0];
     48                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1];
     49                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+2] = coord_system[1][2];
     50                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+0] = coord_system[2][0];
     51                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+1] = coord_system[2][1];
     52                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+2] = coord_system[2][2];
    5353                }
    54                 else if(dimension==4){
     54                else if(dofspernode==4){
    5555                        /*Only the first 3 coordinates are changed (x,y,z), leave the others (P) unchanged*/
    56                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+0] = coord_system[0][0];
    57                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+1] = coord_system[0][1];
    58                         transform[(dimension*numnodes)*(i*dimension+0) + i*dimension+2] = coord_system[0][2];
    59                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+0] = coord_system[1][0];
    60                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+1] = coord_system[1][1];
    61                         transform[(dimension*numnodes)*(i*dimension+1) + i*dimension+2] = coord_system[1][2];
    62                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+0] = coord_system[2][0];
    63                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+1] = coord_system[2][1];
    64                         transform[(dimension*numnodes)*(i*dimension+2) + i*dimension+2] = coord_system[2][2];
    65                         transform[(dimension*numnodes)*(i*dimension+3) + i*dimension+3] = 1.0;
     56                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+0] = coord_system[0][0];
     57                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+1] = coord_system[0][1];
     58                        transform[(dofspernode*numnodes)*(i*dofspernode+0) + i*dofspernode+2] = coord_system[0][2];
     59                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+0] = coord_system[1][0];
     60                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+1] = coord_system[1][1];
     61                        transform[(dofspernode*numnodes)*(i*dofspernode+1) + i*dofspernode+2] = coord_system[1][2];
     62                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+0] = coord_system[2][0];
     63                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+1] = coord_system[2][1];
     64                        transform[(dofspernode*numnodes)*(i*dofspernode+2) + i*dofspernode+2] = coord_system[2][2];
     65                        transform[(dofspernode*numnodes)*(i*dofspernode+3) + i*dofspernode+3] = 1.0;
    6666                }
    6767        }

  • issm/trunk/src/c/shared/Elements/elements.h

    r10367 r10407  
    1717int*   GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum);
    1818int*   GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation_enum);
    19 void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension);
     19void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dofspernode);
     20#ifdef _HAVE_DIAGNOSTIC_
     21void   TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int dofspernode);
     22void   TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int dofspernode);
     23void   TransformSolutionCoord(double* solution,Node** nodes,int numnodes,int dofspernode);
     24#endif
    2025
    2126inline void printarray(double* array,int lines,int cols=1){
    … …  
    2328        for(int i=0;i<lines;i++){ 
    2429                printf("   [ ");
    25                 for(int j=0;j<cols;j++) printf(" %12.6g ",array[i*cols+j]);
     30                for(int j=0;j<cols;j++) printf(" %7.5g ",array[i*cols+j]);
    2631                printf(" ]\n");
    2732        } 
Note: See TracChangeset for help on using the changeset viewer.