Changeset 10367


Ignore:
Timestamp:
10/31/11 10:33:50 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added transformation of coordinate systems in diagnostic

Location:
issm/trunk/src/c
Files:
2 added
10 edited

Legend:

Unmodified
Added
Removed

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

    r10308 r10367  
    184184                                        ./shared/Numerics/cross.cpp\
    185185                                        ./shared/Numerics/extrema.cpp\
     186                                        ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\
    186187                                        ./shared/Numerics/UnitConversion.cpp\
    187188                                        ./shared/Numerics/PetscOptionsFromAnalysis.cpp\
    … …  
    198199                                        ./shared/Elements/GetGlobalDofList.cpp\
    199200                                        ./shared/Elements/GetNumberOfDofs.cpp\
     201                                        ./shared/Elements/CoordinateSystemTransform.cpp\
    200202                                        ./shared/String/sharedstring.h\
    201203                                        ./toolkits/petsc\

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

    r10355 r10367  
    28302830        ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
    28312831        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2832        TransformStiffnessMatrixCoord(Ke,2);
    28322833       
    28332834        /*clean-up and return*/
    … …  
    30313032        }
    30323033
     3034        /*Transform coordinate system*/
     3035        TransformLoadVectorCoord(pe,2);
     3036
    30333037        /*Clean up and return*/
    30343038        delete gauss;
    … …  
    31883192        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    31893193        for(i=0;i<NUMVERTICES;i++){
    3190                 vx[i]=values[i*NDOF2+0];
    3191                 vy[i]=values[i*NDOF2+1];
     3194                /*TEMP*/
     3195                vel[0]=values[i*NDOF2+0];
     3196                vel[1]=values[i*NDOF2+1];
     3197                TransformSolutionCoord(&vel[0],2,i);
     3198                vx[i]=vel[0];
     3199                vy[i]=vel[1];
    31923200
    31933201                /*Check solution*/
    … …  
    32813289        /*Free ressources:*/
    32823290        xfree((void**)&doflist);
     3291}
     3292/*}}}*/
     3293/*FUNCTION Tria::TransformStiffnessMatrixCoord{{{1*/
     3294void Tria::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim){
     3295
     3296        int     i,j;
     3297        int     numnodes          = NUMVERTICES;
     3298        double *transform         = NULL;
     3299        double *values            = NULL;
     3300
     3301        /*Copy current stiffness matrix*/
     3302        values=(double*)xmalloc(Ke->nrows*Ke->ncols*sizeof(double));
     3303        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];
     3304
     3305        /*Get Coordinate Systems transform matrix*/
     3306        CoordinateSystemTransform(&transform,nodes,numnodes,dim);
     3307
     3308        /*Transform matrix: T*Ke*T^t */
     3309        TripleMultiply(transform,numnodes*dim,numnodes*dim,1,
     3310                                values,Ke->nrows,Ke->ncols,0,
     3311                                transform,numnodes*dim,numnodes*dim,0,
     3312                                &Ke->values[0],0);
     3313
     3314        /*Free Matrix*/
     3315        xfree((void**)&transform);
     3316        xfree((void**)&values);
     3317}
     3318/*}}}*/
     3319/*FUNCTION Tria::TransformLoadVectorCoord{{{1*/
     3320void Tria::TransformLoadVectorCoord(ElementVector* pe,int dim){
     3321
     3322        int     i,j;
     3323        int     numnodes          = NUMVERTICES;
     3324        double *transform         = NULL;
     3325        double *values            = NULL;
     3326
     3327        /*Copy current load vector*/
     3328        values=(double*)xmalloc(pe->nrows*sizeof(double));
     3329        for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];
     3330
     3331        /*Get Coordinate Systems transform matrix*/
     3332        CoordinateSystemTransform(&transform,nodes,numnodes,dim);
     3333
     3334        /*Transform matrix: T*pe */
     3335        MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
     3336                                  values,pe->nrows,1,0,
     3337                                  &pe->values[0],0);
     3338
     3339        /*Free Matrix*/
     3340        xfree((void**)&transform);
     3341        xfree((void**)&values);
     3342}
     3343/*}}}*/
     3344/*FUNCTION Tria::TransformSolutionCoord{{{1*/
     3345void Tria::TransformSolutionCoord(double* vel,int dim,int node_index){
     3346
     3347        int     i,j;
     3348        double *transform         = NULL;
     3349        double *values            = NULL;
     3350
     3351        /*Copy current solution vector*/
     3352        values=(double*)xmalloc(dim*sizeof(double));
     3353        for(i=0;i<dim;i++) values[i]=vel[i];
     3354
     3355        /*Get Coordinate Systems transform matrix*/
     3356        CoordinateSystemTransform(&transform,&nodes[node_index],1,dim);
     3357
     3358        /*Transform matrix: T*pe */
     3359        MatrixMultiply(transform,dim,dim,0,
     3360                                values,dim,1,0,
     3361                                vel,0);
     3362
     3363        /*Free Matrix*/
     3364        xfree((void**)&transform);
     3365        xfree((void**)&values);
    32833366}
    32843367/*}}}*/

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

    r10355 r10367  
    9292                void   GetVectorFromInputs(Vec vector, int name_enum);
    9393                void   GetVectorFromResults(Vec vector,int name_enum);
    94                
    9594                void   InputArtificialNoise(int enum_type,double min, double max);
    9695                bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    … …  
    103102                void   DeleteResults(void);
    104103                void   MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
    105                
    106104                void   Migration(double* sheet_ungrounding);
    107105                void   ShelfSync();
    108106                void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding);
    109107                void   MigrateGroundingline();
    110                
    111108                void   RequestedOutput(int output_enum,int step,double time);
    112109                void   ListResultsEnums(int** results_enums,int* num_results);
    … …  
    193190                void           GetInputValue(double* pvalue,Node* node,int enumtype);
    194191                void           GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
    195                 void           InputUpdateFromSolutionOneDof(double* solution,int enum_type);
    196                 void           InputUpdateFromSolutionPrognostic(double* solution);
    197                 bool           IsInput(int name);
    198                 void           SetClone(int* minranks);
    199                 void           SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
     192                void             InputUpdateFromSolutionOneDof(double* solution,int enum_type);
     193                void             InputUpdateFromSolutionPrognostic(double* solution);
     194                bool             IsInput(int name);
     195                void             SetClone(int* minranks);
     196                void             SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
    200197
    201198               
    … …  
    211208                void      InputUpdateFromSolutionDiagnosticHoriz( double* solution);
    212209                void      InputUpdateFromSolutionDiagnosticHutter( double* solution);
     210                void    TransformStiffnessMatrixCoord(ElementMatrix* Ke,int dim);
     211                void    TransformLoadVectorCoord(ElementVector* pe,int dim);
     212                void    TransformSolutionCoord(double* vel,int dim,int node_index);
    213213                #endif
    214214

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

    r10135 r10367  
    528528                }
    529529        }
     530        /*Transform load vector*/
     531        TransformLoadVectorCoord(pe,2);
    530532
    531533        /*Clean up and return*/
    532534        delete gauss;
    533535        return pe;
     536}
     537/*}}}*/
     538/*FUNCTION Icefront::TransformLoadVectorCoord{{{1*/
     539void Icefront::TransformLoadVectorCoord(ElementVector* pe,int dim){
     540
     541        int     numnodes          = 2;
     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,dim);
     551
     552        /*Transform matrix: T*pe */
     553        MatrixMultiply(transform,numnodes*dim,numnodes*dim,1,
     554                                values,pe->nrows,1,0,
     555                                &pe->values[0],0);
     556
     557        /*Free Matrix*/
     558        xfree((void**)&transform);
     559        xfree((void**)&values);
    534560}
    535561/*}}}*/

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

    r9883 r10367  
    8080                /*}}}*/
    8181                /*Load management: {{{1*/
    82                 void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
     82                void GetDofList(int** pdoflist,int approximation_enum,int setenum);
    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 dim);
    8586                #ifdef _HAVE_CONTROL_
    8687                ElementVector* CreatePVectorAdjointHoriz(void);

  • issm/trunk/src/c/objects/Node.cpp

    r10245 r10367  
    4343        this->sid=node_sid;
    4444        this->analysis_type=analysis_type;
    45         for(k=0;k<3;k++) for(l=0;l<3;l++) this->referential[k][l]=1.0;
     45
     46        /*Initialize coord_system: Identity matrix by default*/
     47        for(k=0;k<3;k++) for(l=0;l<3;l++) this->coord_system[k][l]=0.0;
     48        for(k=0;k<3;k++) this->coord_system[k][k]=1.0;
    4649
    4750        /*indexing:*/
    … …  
    5457        //intialize inputs, and add as many inputs per element as requested:
    5558        this->inputs=new Inputs();
    56         if (iomodel->Data(MeshVertexonbedEnum))      this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index]));
    57         if (iomodel->Data(MeshVertexonsurfaceEnum))  this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index]));
    58         if (iomodel->Data(MaskVertexonfloatingiceEnum)) this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index]));
    59         if (iomodel->Data(MaskVertexongroundediceEnum)) this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index]));
    60         if (iomodel->numbernodetoelementconnectivity) this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,(IssmInt)iomodel->numbernodetoelementconnectivity[io_index]));
    61         if (analysis_type==DiagnosticHorizAnalysisEnum) this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index]));
     59        if (iomodel->Data(MeshVertexonbedEnum))
     60         this->inputs->AddInput(new BoolInput(MeshVertexonbedEnum,(IssmBool)iomodel->Data(MeshVertexonbedEnum)[io_index]));
     61        if (iomodel->Data(MeshVertexonsurfaceEnum))
     62         this->inputs->AddInput(new BoolInput(MeshVertexonsurfaceEnum,(IssmBool)iomodel->Data(MeshVertexonsurfaceEnum)[io_index]));
     63        if (iomodel->Data(MaskVertexonfloatingiceEnum))
     64         this->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,(IssmBool)iomodel->Data(MaskVertexonfloatingiceEnum)[io_index]));
     65        if (iomodel->Data(MaskVertexongroundediceEnum))
     66         this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,(IssmBool)iomodel->Data(MaskVertexongroundediceEnum)[io_index]));
     67        if (iomodel->numbernodetoelementconnectivity)
     68         this->inputs->AddInput(new IntInput(NumberNodeToElementConnectivityEnum,(IssmInt)iomodel->numbernodetoelementconnectivity[io_index]));
     69        if (analysis_type==DiagnosticHorizAnalysisEnum)
     70         this->inputs->AddInput(new IntInput(ApproximationEnum,(IssmInt)iomodel->Data(FlowequationVertexEquationEnum)[io_index]));
    6271       
    6372        /*set single point constraints: */
    … …  
    7281
    7382        /*Diagnostic Horiz*/
     83        #ifdef _HAVE_DIAGNOSTIC_
    7484        if (analysis_type==DiagnosticHorizAnalysisEnum){
     85
     86                /*Coordinate system provided, convert to coord_system matrix*/
     87                _assert_(iomodel->Data(DiagnosticReferentialEnum));
     88                XZvectorsToCoordinateSystem(&this->coord_system[0][0],iomodel->Data(DiagnosticReferentialEnum)+io_index*6);
     89
    7590                if (dim==3){
    7691                        /*We have a  3d mesh, we may have collapsed elements, hence dead nodes. Freeze them out: */
    77                         if (!iomodel->Data(MeshVertexonbedEnum)) _error_("iomodel->nodeonbed is NULL");
    78                         if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL");
     92                        _assert_(iomodel->Data(MeshVertexonbedEnum));
     93                        _assert_(iomodel->Data(FlowequationVertexEquationEnum));
    7994                        if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){
    8095                                for(k=1;k<=gsize;k++) this->FreezeDof(k);
    … …  
    97112                        }
    98113                }
    99                
    100         }
     114        }
     115        #endif
    101116
    102117        /*Diagnostic Hutter*/
    … …  
    202217        memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
    203218        memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
    204         memcpy(marshalled_dataset,&referential,9*sizeof(double));marshalled_dataset+=9*sizeof(double); 
     219        memcpy(marshalled_dataset,&coord_system,9*sizeof(double));marshalled_dataset+=9*sizeof(double); 
    205220
    206221        /*marshall objects: */
    … …  
    248263        memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
    249264        memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
    250         memcpy(&referential,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double);
     265        memcpy(&coord_system,marshalled_dataset,9*sizeof(double));marshalled_dataset+=9*sizeof(double);
    251266       
    252267        /*demarshall objects: */
    … …  
    507522}
    508523/*}}}*/
     524#ifdef _HAVE_DIAGNOSTIC_
     525/*FUNCTION Node::GetCoordinateSystem{{{1*/
     526void Node::GetCoordinateSystem(double* coord_system_out){
     527
     528        /*Copy coord_system*/
     529        for(int k=0;k<3;k++) for(int l=0;l<3;l++) coord_system_out[3*k+l]=this->coord_system[k][l];
     530
     531}
     532/*}}}*/
     533#endif
    509534/*FUNCTION Node::SetVertexDof {{{1*/
    510535void   Node::SetVertexDof(int in_dof){

  • issm/trunk/src/c/objects/Node.h

    r10143 r10367  
    3030                Inputs*        inputs; //properties of this node
    3131                int            analysis_type;
    32                 double         referential[3][3];
     32                double         coord_system[3][3];
    3333
    3434                /*Node constructors, destructors {{{1*/
    … …  
    6565                /*}}}*/
    6666                /*Node numerical routines {{{1*/
    67                 void  Configure(DataSet* nodes,Vertices* vertices);
    68                 void  CreateNodalConstraints(Vec ys);
    69                 void  SetCurrentConfiguration(DataSet* nodes,Vertices* vertices);
    70                 int   Sid(void);
    71                 int   GetVertexDof(void);
    72                 int   GetVertexId(void);
    73                 void  SetVertexDof(int in_dof);
    74                 bool  InAnalysis(int analysis_type);
    75                 int   GetApproximation();
    76                 int   GetNumberOfDofs(int approximation_enum,int setenum);
    77                 int   IsClone();
    78                 void  ApplyConstraint(int dof,double value);
    79                 void  RelaxConstraint(int dof);
    80                 void  DofInSSet(int dof);
    81                 void  DofInFSet(int dof);
    82                 int   GetDof(int dofindex,int setenum);
    83                 void  CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s);
    84                 int   GetConnectivity();
    85                 void  GetDofList(int* poutdoflist,int approximation_enum,int setenum);
    86                 void  GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
    87                 int   GetDofList1(void);
    88                 int   GetSidList(void);
     67                void   Configure(DataSet* nodes,Vertices* vertices);
     68                void   CreateNodalConstraints(Vec ys);
     69                void   SetCurrentConfiguration(DataSet* nodes,Vertices* vertices);
     70                int    Sid(void);
     71                int    GetVertexDof(void);
     72                int    GetVertexId(void);
     73#ifdef _HAVE_DIAGNOSTIC_
     74                void   GetCoordinateSystem(double* coord_system_out);
     75#endif
     76                void   SetVertexDof(int in_dof);
     77                bool   InAnalysis(int analysis_type);
     78                int    GetApproximation();
     79                int    GetNumberOfDofs(int approximation_enum,int setenum);
     80                int    IsClone();
     81                void   ApplyConstraint(int dof,double value);
     82                void   RelaxConstraint(int dof);
     83                void   DofInSSet(int dof);
     84                void   DofInFSet(int dof);
     85                int    GetDof(int dofindex,int setenum);
     86                void   CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s);
     87                int    GetConnectivity();
     88                void   GetDofList(int* poutdoflist,int approximation_enum,int setenum);
     89                void   GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
     90                int    GetDofList1(void);
     91                int    GetSidList(void);
    8992                double GetX();
    9093                double GetY();
    9194                double GetZ();
    9295                double GetSigma();
    93                 int   IsOnBed();
    94                 int   IsOnSurface();
    95                 void  FreezeDof(int dof);
    96                 int   IsFloating();
    97                 int   IsGrounded();
    98                 void  UpdateSpcs(double* ys);
    99                 void  VecMerge(Vec ug, double* vector_serial,int setnum);
    100                 void  VecReduce(Vec vector, double* ug_serial,int setnum);
     96                int    IsOnBed();
     97                int    IsOnSurface();
     98                void   FreezeDof(int dof);
     99                int    IsFloating();
     100                int    IsGrounded();
     101                void   UpdateSpcs(double* ys);
     102                void   VecMerge(Vec ug, double* vector_serial,int setnum);
     103                void   VecReduce(Vec vector, double* ug_serial,int setnum);
    101104               
    102105                /*}}}*/

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

    r7848 r10367  
    1717int*   GetLocalDofList( Node** nodes,int numnodes,int setenum,int approximation_enum);
    1818int*   GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation_enum);
     19void   CoordinateSystemTransform(double** ptransform,Node** nodes,int numnodes,int dimension);
    1920
    2021inline void printarray(double* array,int lines,int cols=1){

  • issm/trunk/src/c/shared/Numerics/numerics.h

    r10144 r10367  
    2828double UnitConversion(double value, int direction_enum, int type_enum);
    2929void   PetscOptionsFromAnalysis(Parameters* parameters,int analysis_type);
     30void   XZvectorsToCoordinateSystem(double* T,double* xzvectors);
    3031
    3132#endif //ifndef _NUMERICS_H_

  • issm/trunk/src/c/solutions/issm.cpp

    r9761 r10367  
    107107                #ifdef _HAVE_CONTROL_
    108108                #ifdef _HAVE_TAO_
    109                 controltao_core(femmodel);
     109                //controltao_core(femmodel);
     110                control_core(femmodel);
    110111                #else
    111112                control_core(femmodel);
Note: See TracChangeset for help on using the changeset viewer.