Changeset 17255


Ignore:
Timestamp:
02/11/14 13:41:23 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moved matrix rotation from Node.cpp to Element.cpp

Location:
issm/trunk-jpl/src/c/classes
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r17248 r17255  
    5353        return false;
    5454}/*}}}*/
     55void Element::CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
     56
     57        int         i,counter;
     58        int         numdofs   = 0;
     59        IssmDouble  norm;
     60        IssmDouble *transform = NULL;
     61        IssmDouble  coord_system[3][3];
     62
     63        /*Some checks in debugging mode*/
     64        _assert_(numnodes && nodes_list);
     65
     66        /*Get total number of dofs*/
     67        for(i=0;i<numnodes;i++){
     68                switch(cs_array[i]){
     69                        case PressureEnum: numdofs+=1; break;
     70                        case XYEnum:       numdofs+=2; break;
     71                        case XYZEnum:      numdofs+=3; break;
     72                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     73                }
     74        }
     75
     76        /*Allocate and initialize transform matrix*/
     77        transform=xNew<IssmDouble>(numdofs*numdofs);
     78        for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;
     79
     80        /*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
     81         *for 3 nodes:
     82
     83         *     | T1 0  0 |
     84         * Q = | 0  T2 0 |
     85         *     | 0  0  T3|
     86         *
     87         * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
     88         * associated to this node*/
     89        counter=0;
     90        for(i=0;i<numnodes;i++){
     91                nodes_list[i]->GetCoordinateSystem(&coord_system[0][0]);
     92                switch(cs_array[i]){
     93                        case PressureEnum:
     94                                /*DO NOT change anything*/
     95                                transform[(numdofs)*(counter) + counter] = 1.;
     96                                counter+=1;
     97                                break;
     98                        case XYEnum:
     99                                /*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
     100                                norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
     101                                transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
     102                                transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
     103                                transform[(numdofs)*(counter+1) + counter+0] =   coord_system[1][0]/norm;
     104                                transform[(numdofs)*(counter+1) + counter+1] =   coord_system[0][0]/norm;
     105                                counter+=2;
     106                                break;
     107                        case XYZEnum:
     108                                /*The 3 coordinates are changed (x,y,z)*/
     109                                transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
     110                                transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
     111                                transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
     112                                transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
     113                                transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
     114                                transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
     115                                transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
     116                                transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
     117                                transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
     118                                counter+=3;
     119                                break;
     120                        default:
     121                                _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     122                }
     123        }
     124
     125        /*Assign output pointer*/
     126        *ptransform=transform;
     127}
     128/*}}}*/
    55129void Element::DeleteMaterials(void){/*{{{*/
    56130        delete this->material;
    … …  
    537611}/*}}}*/
    538612void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
    539         ::TransformInvStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum);
    540 }/*}}}*/
    541 void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
    542         ::TransformInvStiffnessMatrixCoord(Ke,nodes_list,numnodes,transformenum_list);
     613
     614        /*All nodes have the same Coordinate System*/
     615        int numnodes  = this->GetNumberOfNodes();
     616        int* cs_array = xNew<int>(numnodes);
     617        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
     618
     619        /*Call core*/
     620        TransformInvStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
     621
     622        /*Clean-up*/
     623        xDelete<int>(cs_array);
     624}/*}}}*/
     625void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
     626
     627        int         i,j;
     628        int         numdofs   = 0;
     629        IssmDouble *transform = NULL;
     630        IssmDouble *values    = NULL;
     631
     632        /*Get total number of dofs*/
     633        for(i=0;i<numnodes;i++){
     634                switch(cs_array[i]){
     635                        case PressureEnum: numdofs+=1; break;
     636                        case XYEnum:       numdofs+=2; break;
     637                        case XYZEnum:      numdofs+=3; break;
     638                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     639                }
     640        }
     641
     642        /*Copy current stiffness matrix*/
     643        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
     644        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];
     645
     646        /*Get Coordinate Systems transform matrix*/
     647        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
     648
     649        /*Transform matrix: R*Ke*R^T */
     650        TripleMultiply(transform,numdofs,numdofs,0,
     651                                values,Ke->nrows,Ke->ncols,0,
     652                                transform,numdofs,numdofs,1,
     653                                &Ke->values[0],0);
     654
     655        /*Free Matrix*/
     656        xDelete<IssmDouble>(transform);
     657        xDelete<IssmDouble>(values);
    543658}/*}}}*/
    544659void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/
    545         ::TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),transformenum);
    546 }/*}}}*/
    547 void Element::TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list){/*{{{*/
    548         ::TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),transformenum_list);
     660
     661        /*All nodes have the same Coordinate System*/
     662        int  numnodes = this->GetNumberOfNodes();
     663        int* cs_array = xNew<int>(numnodes);
     664        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
     665
     666        /*Call core*/
     667        this->TransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array);
     668
     669        /*Clean-up*/
     670        xDelete<int>(cs_array);
     671}/*}}}*/
     672void Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/
     673
     674        this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array);
     675
     676}/*}}}*/
     677void Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
     678
     679        int         i;
     680        int         numdofs   = 0;
     681        IssmDouble *transform = NULL;
     682        IssmDouble *values    = NULL;
     683
     684        /*Get total number of dofs*/
     685        for(i=0;i<numnodes;i++){
     686                switch(cs_array[i]){
     687                        case PressureEnum: numdofs+=1; break;
     688                        case XYEnum:       numdofs+=2; break;
     689                        case XYZEnum:      numdofs+=3; break;
     690                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     691                }
     692        }
     693
     694        /*Copy current load vector*/
     695        values=xNew<IssmDouble>(pe->nrows);
     696        for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
     697
     698        /*Get Coordinate Systems transform matrix*/
     699        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
     700
     701        /*Transform matrix: R^T*pe */
     702        MatrixMultiply(transform,numdofs,numdofs,1,
     703                                values,pe->nrows,1,0,
     704                                &pe->values[0],0);
     705
     706        /*Free Matrices*/
     707        xDelete<IssmDouble>(transform);
     708        xDelete<IssmDouble>(values);
    549709}/*}}}*/
    550710void Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/
    551         ::TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum);
     711
     712        /*All nodes have the same Coordinate System*/
     713        int  numnodes = this->GetNumberOfNodes();
     714        int* cs_array = xNew<int>(numnodes);
     715        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
     716
     717        /*Call core*/
     718        this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);
     719
     720        /*Clean-up*/
     721        xDelete<int>(cs_array);
    552722}/*}}}*/
    553723void Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/
    554         ::TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
     724        this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
    555725}/*}}}*/
    556726void Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/
    557         ::TransformSolutionCoord(values,this->nodes,numnodes,transformenum);
    558 }/*}}}*/
    559 void Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list){/*{{{*/
    560         ::TransformSolutionCoord(values,this->nodes,numnodes,transformenum_list);
    561 }/*}}}*/
    562 void Element::TransformLoadVectorCoord(ElementVector* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
    563         ::TransformLoadVectorCoord(Ke,nodes_list,numnodes,transformenum_list);
     727        this->TransformSolutionCoord(values,this->nodes,numnodes,transformenum);
     728}/*}}}*/
     729void Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/
     730        this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array);
     731}/*}}}*/
     732void Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/
     733        this->TransformSolutionCoord(values,nodes_list,numnodes,transformenum);
     734}/*}}}*/
     735void Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
     736
     737        int         i;
     738        int         numdofs   = 0;
     739        IssmDouble *transform = NULL;
     740        IssmDouble *values    = NULL;
     741
     742        /*Get total number of dofs*/
     743        for(i=0;i<numnodes;i++){
     744                switch(cs_array[i]){
     745                        case PressureEnum: numdofs+=1; break;
     746                        case XYEnum:       numdofs+=2; break;
     747                        case XYZEnum:      numdofs+=3; break;
     748                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     749                }
     750        }
     751
     752        /*Copy current solution vector*/
     753        values=xNew<IssmDouble>(numdofs);
     754        for(i=0;i<numdofs;i++) values[i]=solution[i];
     755
     756        /*Get Coordinate Systems transform matrix*/
     757        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
     758
     759        /*Transform matrix: R*U */
     760        MatrixMultiply(transform,numdofs,numdofs,0,
     761                                values,numdofs,1,0,
     762                                &solution[0],0);
     763
     764        /*Free Matrices*/
     765        xDelete<IssmDouble>(transform);
     766        xDelete<IssmDouble>(values);
    564767}/*}}}*/
    565768void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
    566         ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum);
     769
     770        /*All nodes have the same Coordinate System*/
     771        int  numnodes = this->GetNumberOfNodes();
     772        int* cs_array = xNew<int>(numnodes);
     773        for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;
     774
     775        /*Call core*/
     776        this->TransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);
     777
     778        /*Clean-up*/
     779        xDelete<int>(cs_array);
    567780}/*}}}*/
    568781void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/
    569         ::TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
    570 }/*}}}*/
    571 void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list){/*{{{*/
    572         ::TransformStiffnessMatrixCoord(Ke,nodes_list,numnodes,transformenum_list);
    573 }/*}}}*/
     782        this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
     783}/*}}}*/
     784void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/
     785
     786        int         numdofs = 0;
     787        IssmDouble *transform = NULL;
     788        IssmDouble *values    = NULL;
     789
     790        /*Get total number of dofs*/
     791        for(int i=0;i<numnodes;i++){
     792                switch(cs_array[i]){
     793                        case PressureEnum: numdofs+=1; break;
     794                        case XYEnum:       numdofs+=2; break;
     795                        case XYZEnum:      numdofs+=3; break;
     796                        default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
     797                }
     798        }
     799
     800        /*Copy current stiffness matrix*/
     801        values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
     802        for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
     803
     804        /*Get Coordinate Systems transform matrix*/
     805        CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);
     806
     807        /*Transform matrix: R^T*Ke*R */
     808        TripleMultiply(transform,numdofs,numdofs,1,
     809                                values,Ke->nrows,Ke->ncols,0,
     810                                transform,numdofs,numdofs,0,
     811                                &Ke->values[0],0);
     812
     813        /*Free Matrix*/
     814        xDelete<IssmDouble>(transform);
     815        xDelete<IssmDouble>(values);
     816}/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r17242 r17255  
    5454                bool       AllActive(void);
    5555                bool       AnyActive(void);
     56                void       CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
    5657                void       DeleteMaterials(void);
    5758                void       FindParam(bool* pvalue,int paramenum);
    … …  
    7374                void       StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    7475                IssmDouble TMeltingPoint(IssmDouble pressure);
     76                void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum);
     77                void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
     78                void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
     79                void       TransformLoadVectorCoord(ElementVector* pe,int cs_enum);
     80                void       TransformLoadVectorCoord(ElementVector* pe,int* cs_array);
     81                void       TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum);
     82                void       TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array);
     83                void       TransformSolutionCoord(IssmDouble* solution,int cs_enum);
     84                void       TransformSolutionCoord(IssmDouble* solution,int* cs_array);
     85                void       TransformSolutionCoord(IssmDouble* solution,int numnodes,int cs_enum);
     86                void       TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array);
     87                void       TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum);
     88                void       TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array);
     89                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum);
     90                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* cs_array);
     91                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
     92                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
    7593                void       ViscousHeatingCreateInput(void);
    7694                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    … …  
    8199                void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    82100                void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    83                 void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
    84                 void       TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list);
    85                 void       TransformLoadVectorCoord(ElementVector* pe,int transformenum);
    86                 void       TransformLoadVectorCoord(ElementVector* pe,int* transformenum_list);
    87101                void       TransformLoadVectorCoord(ElementVector* pe,int numnodes,int transformenum){_error_("not implemented yet");};/*Tiling only*/
    88102                void       TransformLoadVectorCoord(ElementVector* pe,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/
    89                 void       TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* transformenum_list);
    90                 void       TransformSolutionCoord(IssmDouble* values,int transformenum);
    91                 void       TransformSolutionCoord(IssmDouble* values,int* transformenum_list);
    92                 void       TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum);
    93                 void       TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list);
    94                 void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
    95                 void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list);
    96                 void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum);
    97103                void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/
    98                 void       TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* transformenum_list);
    99104
    100105                /*Virtual functions*/

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r17236 r17255  
    896896}
    897897/*}}}*/
    898 void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){/*{{{*/
    899 
    900         int* cs_array=NULL;
    901 
    902         /*All nodes have the same Coordinate System*/
    903         cs_array=xNew<int>(numnodes);
    904         for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
    905 
    906         /*Call core*/
    907         TransformInvStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);
    908 
    909         /*Clean-up*/
    910         xDelete<int>(cs_array);
    911 }
    912 /*}}}*/
    913 void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){/*{{{*/
    914 
    915         int     i,j;
    916         int     numdofs   = 0;
    917         IssmDouble *transform = NULL;
    918         IssmDouble *values    = NULL;
    919 
    920         /*Get total number of dofs*/
    921         for(i=0;i<numnodes;i++){
    922                 switch(cs_array[i]){
    923                         case PressureEnum: numdofs+=1; break;
    924                         case XYEnum:       numdofs+=2; break;
    925                         case XYZEnum:      numdofs+=3; break;
    926                         default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    927                 }
    928         }
    929 
    930         /*Copy current stiffness matrix*/
    931         values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
    932         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];
    933 
    934         /*Get Coordinate Systems transform matrix*/
    935         CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
    936 
    937         /*Transform matrix: R*Ke*R^T */
    938         TripleMultiply(transform,numdofs,numdofs,0,
    939                                 values,Ke->nrows,Ke->ncols,0,
    940                                 transform,numdofs,numdofs,1,
    941                                 &Ke->values[0],0);
    942 
    943         /*Free Matrix*/
    944         xDelete<IssmDouble>(transform);
    945         xDelete<IssmDouble>(values);
    946 }
    947 /*}}}*/
    948 void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum){/*{{{*/
    949 
    950         int* cs_array=NULL;
    951 
    952         /*All nodes have the same Coordinate System*/
    953         cs_array=xNew<int>(numnodes);
    954         for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
    955 
    956         /*Call core*/
    957         TransformLoadVectorCoord(pe,nodes,numnodes,cs_array);
    958 
    959         /*Clean-up*/
    960         xDelete<int>(cs_array);
    961 }
    962 /*}}}*/
    963 void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array){/*{{{*/
    964 
    965         int     i;
    966         int     numdofs   = 0;
    967         IssmDouble *transform = NULL;
    968         IssmDouble *values    = NULL;
    969 
    970         /*Get total number of dofs*/
    971         for(i=0;i<numnodes;i++){
    972                 switch(cs_array[i]){
    973                         case PressureEnum: numdofs+=1; break;
    974                         case XYEnum:       numdofs+=2; break;
    975                         case XYZEnum:      numdofs+=3; break;
    976                         default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    977                 }
    978         }
    979 
    980         /*Copy current load vector*/
    981         values=xNew<IssmDouble>(pe->nrows);
    982         for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
    983 
    984         /*Get Coordinate Systems transform matrix*/
    985         CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
    986 
    987         /*Transform matrix: R^T*pe */
    988         MatrixMultiply(transform,numdofs,numdofs,1,
    989                                 values,pe->nrows,1,0,
    990                                 &pe->values[0],0);
    991 
    992         /*Free Matrices*/
    993         xDelete<IssmDouble>(transform);
    994         xDelete<IssmDouble>(values);
    995 }
    996 /*}}}*/
    997 void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum){/*{{{*/
    998 
    999         int* cs_array=NULL;
    1000 
    1001         /*All nodes have the same Coordinate System*/
    1002         cs_array=xNew<int>(numnodes);
    1003         for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
    1004 
    1005         /*Call core*/
    1006         TransformSolutionCoord(solution,nodes,numnodes,cs_array);
    1007 
    1008         /*Clean-up*/
    1009         xDelete<int>(cs_array);
    1010 }
    1011 /*}}}*/
    1012 void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array){/*{{{*/
    1013 
    1014         int     i;
    1015         int     numdofs   = 0;
    1016         IssmDouble *transform = NULL;
    1017         IssmDouble *values    = NULL;
    1018 
    1019         /*Get total number of dofs*/
    1020         for(i=0;i<numnodes;i++){
    1021                 switch(cs_array[i]){
    1022                         case PressureEnum: numdofs+=1; break;
    1023                         case XYEnum:       numdofs+=2; break;
    1024                         case XYZEnum:      numdofs+=3; break;
    1025                         default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    1026                 }
    1027         }
    1028 
    1029         /*Copy current solution vector*/
    1030         values=xNew<IssmDouble>(numdofs);
    1031         for(i=0;i<numdofs;i++) values[i]=solution[i];
    1032 
    1033         /*Get Coordinate Systems transform matrix*/
    1034         CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
    1035 
    1036         /*Transform matrix: R*U */
    1037         MatrixMultiply(transform,numdofs,numdofs,0,
    1038                                 values,numdofs,1,0,
    1039                                 &solution[0],0);
    1040 
    1041         /*Free Matrices*/
    1042         xDelete<IssmDouble>(transform);
    1043         xDelete<IssmDouble>(values);
    1044 }
    1045 /*}}}*/
    1046 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){/*{{{*/
    1047 
    1048         int* cs_array=NULL;
    1049 
    1050         /*All nodes have the same Coordinate System*/
    1051         cs_array=xNew<int>(numnodes);
    1052         for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;
    1053 
    1054         /*Call core*/
    1055         TransformStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);
    1056 
    1057         /*Clean-up*/
    1058         xDelete<int>(cs_array);
    1059 }
    1060 /*}}}*/
    1061 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){/*{{{*/
    1062 
    1063         int         numdofs = 0;
    1064         IssmDouble *transform = NULL;
    1065         IssmDouble *values    = NULL;
    1066 
    1067         /*Get total number of dofs*/
    1068         for(int i=0;i<numnodes;i++){
    1069                 switch(cs_array[i]){
    1070                         case PressureEnum: numdofs+=1; break;
    1071                         case XYEnum:       numdofs+=2; break;
    1072                         case XYZEnum:      numdofs+=3; break;
    1073                         default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    1074                 }
    1075         }
    1076 
    1077         /*Copy current stiffness matrix*/
    1078         values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
    1079         for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
    1080 
    1081         /*Get Coordinate Systems transform matrix*/
    1082         CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);
    1083 
    1084         /*Transform matrix: R^T*Ke*R */
    1085         TripleMultiply(transform,numdofs,numdofs,1,
    1086                                 values,Ke->nrows,Ke->ncols,0,
    1087                                 transform,numdofs,numdofs,0,
    1088                                 &Ke->values[0],0);
    1089 
    1090         /*Free Matrix*/
    1091         xDelete<IssmDouble>(transform);
    1092         xDelete<IssmDouble>(values);
    1093 }
    1094 /*}}}*/
    1095 void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array){/*{{{*/
    1096 
    1097         int         i,counter;
    1098         int         numdofs   = 0;
    1099         IssmDouble  norm;
    1100         IssmDouble *transform = NULL;
    1101         IssmDouble  coord_system[3][3];
    1102 
    1103         /*Some checks in debugging mode*/
    1104         _assert_(numnodes && nodes);
    1105 
    1106         /*Get total number of dofs*/
    1107         for(i=0;i<numnodes;i++){
    1108                 switch(cs_array[i]){
    1109                         case PressureEnum: numdofs+=1; break;
    1110                         case XYEnum:       numdofs+=2; break;
    1111                         case XYZEnum:      numdofs+=3; break;
    1112                         default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    1113                 }
    1114         }
    1115 
    1116         /*Allocate and initialize transform matrix*/
    1117         transform=xNew<IssmDouble>(numdofs*numdofs);
    1118         for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;
    1119 
    1120         /*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
    1121          *for 3 nodes:
    1122 
    1123          *     | T1 0  0 |
    1124          * Q = | 0  T2 0 |
    1125          *     | 0  0  T3|
    1126          *
    1127          * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
    1128          * associated to this node*/
    1129         counter=0;
    1130         for(i=0;i<numnodes;i++){
    1131                 nodes[i]->GetCoordinateSystem(&coord_system[0][0]);
    1132                 switch(cs_array[i]){
    1133                         case PressureEnum:
    1134                                 /*DO NOT change anything*/
    1135                                 transform[(numdofs)*(counter) + counter] = 1.;
    1136                                 counter+=1;
    1137                                 break;
    1138                         case XYEnum:
    1139                                 /*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
    1140                                 norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
    1141                                 transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
    1142                                 transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
    1143                                 transform[(numdofs)*(counter+1) + counter+0] =   coord_system[1][0]/norm;
    1144                                 transform[(numdofs)*(counter+1) + counter+1] =   coord_system[0][0]/norm;
    1145                                 counter+=2;
    1146                                 break;
    1147                         case XYZEnum:
    1148                                 /*The 3 coordinates are changed (x,y,z)*/
    1149                                 transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
    1150                                 transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
    1151                                 transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
    1152                                 transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
    1153                                 transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
    1154                                 transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
    1155                                 transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
    1156                                 transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
    1157                                 transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
    1158                                 counter+=3;
    1159                                 break;
    1160                         default:
    1161                                 _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
    1162                 }
    1163         }
    1164 
    1165         /*Assign output pointer*/
    1166         *ptransform=transform;
    1167 }
    1168 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Node.h

    r17236 r17255  
    8989int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation);
    9090int  GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation);
    91 void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
    92 void TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
    93 void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int cs_enum);
    94 void TransformLoadVectorCoord(ElementVector* pe,Node** nodes,int numnodes,int* cs_array);
    95 void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int cs_enum);
    96 void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array);
    97 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
    98 void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
    99 void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array);
    10091
    10192#endif  /* _NODE_H_ */
Note: See TracChangeset for help on using the changeset viewer.