Changeset 15625


Ignore:
Timestamp:
07/25/13 15:42:34 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: some changes necessary for P2 Elements

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

Legend:

Unmodified
Added
Removed

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

    r15621 r15625  
    845845void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
    846846
    847         int  i,count=0;
    848         int  numberofdofs=0;
    849         int* doflist=NULL;
    850 
    851         /*First, figure out size of doflist: */
    852         for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    853 
    854         /*Allocate: */
    855         doflist=xNew<int>(numberofdofs);
     847        /*Fetch number of nodes and dof for this finite element*/
     848        int numnodes = this->NumberofNodes();
     849
     850        /*First, figure out size of doflist and create it: */
     851        int numberofdofs=0;
     852        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     853
     854        /*Allocate output*/
     855        int* doflist=xNew<int>(numberofdofs);
    856856
    857857        /*Populate: */
    858         count=0;
    859         for(i=0;i<6;i++){
     858        int count=0;
     859        for(int i=0;i<numnodes;i++){
    860860                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
    861861                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    … …  
    11021102                for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
    11031103        }
     1104}
     1105/*}}}*/
     1106/*FUNCTION Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/
     1107void Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){
     1108
     1109        _assert_(pvalue);
     1110
     1111        Input *input    = inputs->GetInput(enumtype);
     1112        int    numnodes = this->NumberofNodes();
     1113
     1114        /* Start looping on the number of vertices: */
     1115        if(input){
     1116                GaussPenta* gauss=new GaussPenta();
     1117                for(int iv=0;iv<this->NumberofNodes();iv++){
     1118                        gauss->GaussNode(numnodes,iv);
     1119                        input->GetInputValue(&pvalue[iv],gauss);
     1120                }
     1121                delete gauss;
     1122        }
     1123        else{
     1124                for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
     1125        }
     1126}
     1127/*}}}*/
     1128/*FUNCTION Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype) {{{*/
     1129void Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){
     1130
     1131        _assert_(pvalue);
     1132
     1133        /*Recover input*/
     1134        Input* input=inputs->GetInput(enumtype);
     1135        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
     1136        int    numnodes = this->NumberofNodes();
     1137
     1138        /* Start looping on the number of vertices: */
     1139        GaussPenta* gauss=new GaussPenta();
     1140        for (int iv=0;iv<this->NumberofNodes();iv++){
     1141                gauss->GaussNode(numnodes,iv);
     1142                input->GetInputValue(&pvalue[iv],gauss);
     1143        }
     1144        delete gauss;
    11041145}
    11051146/*}}}*/
    … …  
    94079448void  Penta::InputUpdateFromSolutionDiagnosticHO(IssmDouble* solution){
    94089449
    9409         const int    numdof=NDOF2*NUMVERTICES;
    9410 
    9411         int    i;
    9412         IssmDouble rho_ice,g;
    9413         IssmDouble values[numdof];
    9414         IssmDouble vx[NUMVERTICES];
    9415         IssmDouble vy[NUMVERTICES];
    9416         IssmDouble vz[NUMVERTICES];
    9417         IssmDouble vel[NUMVERTICES];
    9418         IssmDouble pressure[NUMVERTICES];
    9419         IssmDouble surface[NUMVERTICES];
    9420         IssmDouble xyz_list[NUMVERTICES][3];
    9421         int*   doflist = NULL;
    9422 
    9423         /*Get dof list: */
     9450        int         i;
     9451        IssmDouble  rho_ice,g;
     9452        IssmDouble  xyz_list[NUMVERTICES][3];
     9453        int        *doflist = NULL;
     9454
     9455        /*Fetch number of nodes and dof for this finite element*/
     9456        int numnodes = this->NumberofNodes();
     9457        int numdof   = numnodes*NDOF2;
     9458
     9459        /*Fetch dof list and allocate solution vectors*/
    94249460        GetDofList(&doflist,HOApproximationEnum,GsetEnum);
    9425 
    9426         /*Get node data: */
    9427         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     9461        IssmDouble* values    = xNew<IssmDouble>(numdof);
     9462        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     9463        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     9464        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     9465        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     9466        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
     9467        IssmDouble* surface   = xNew<IssmDouble>(numnodes);
    94289468
    94299469        /*Use the dof list to index into the solution vector: */
    … …  
    94329472        /*Transform solution in Cartesian Space*/
    94339473        TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
     9474        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    94349475
    94359476        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    9436         for(i=0;i<NUMVERTICES;i++){
     9477        for(i=0;i<numnodes;i++){
    94379478                vx[i]=values[i*NDOF2+0];
    94389479                vy[i]=values[i*NDOF2+1];
    … …  
    94439484        }
    94449485
    9445         /*Get Vz*/
    9446         Input* vz_input=inputs->GetInput(VzEnum);
    9447         if (vz_input){
    9448                 GetInputListOnVertices(&vz[0],VzEnum);
    9449         }
    9450         else{
    9451                 for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
    9452         }
    9453 
    9454         /*Now Compute vel*/
    9455         for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     9486        /*Get Vz and compute vel*/
     9487        GetInputListOnNodes(&vz[0],VzEnum,0.);
     9488        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    94569489
    94579490        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    … …  
    94759508
    94769509        /*Free ressources:*/
     9510        xDelete<IssmDouble>(surface);
     9511        xDelete<IssmDouble>(pressure);
     9512        xDelete<IssmDouble>(vel);
     9513        xDelete<IssmDouble>(vz);
     9514        xDelete<IssmDouble>(vy);
     9515        xDelete<IssmDouble>(vx);
     9516        xDelete<IssmDouble>(values);
    94779517        xDelete<int>(doflist);
    94789518}

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

    r15613 r15625  
    192192                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    193193                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
     194                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
     195                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    194196                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    195197                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);

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

    r15621 r15625  
    34253425        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    34263426        IssmDouble* values    = xNew<IssmDouble>(numdof);
    3427         IssmDouble* vx        = xNew<IssmDouble>(numdof);
    3428         IssmDouble* vy        = xNew<IssmDouble>(numdof);
    3429         IssmDouble* vz        = xNew<IssmDouble>(numdof);
    3430         IssmDouble* vel       = xNew<IssmDouble>(numdof);
    3431         IssmDouble* pressure  = xNew<IssmDouble>(numdof);
    3432         IssmDouble* thickness = xNew<IssmDouble>(numdof);
     3427        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     3428        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     3429        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     3430        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     3431        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
     3432        IssmDouble* thickness = xNew<IssmDouble>(numnodes);
    34333433
    34343434        /*Use the dof list to index into the solution vector: */

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp

    r15564 r15625  
    141141        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    142142        IssmDouble velocity[numnodes][DOFVELOCITY];
    143 
    144143        _assert_(this->NumberofNodes()==6); //Check Tria too
    145144
    … …  
    313312        /*Get B matrix: */
    314313        GetBHO(&B[0][0], xyz_list, gauss);
    315 
    316314        _assert_(this->NumberofNodes()==6); //Check Tria too
    317315

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r15538 r15625  
    366366        /*update static arrays*/
    367367        switch(iv){
    368                 case 0:
    369                         coord1=1; coord2=0; coord3=0; coord4= -1;
    370                         break;
    371                 case 1:
    372                         coord1=0; coord2=1; coord3=0; coord4= -1;
    373                         break;
    374                 case 2:
    375                         coord1=0; coord2=0; coord3=1; coord4= -1;
    376                         break;
    377                 case 3:
    378                         coord1=1; coord2=0; coord3=0; coord4= +1;
    379                         break;
    380                 case 4:
    381                         coord1=0; coord2=1; coord3=0; coord4= +1;
    382                         break;
    383                 case 5:
    384                         coord1=0; coord2=0; coord3=1; coord4= +1;
    385                         break;
    386                 default:
    387                         _error_("vertex index should be in [0 5]");
     368                case 0: coord1=1.; coord2=0.; coord3=0.; coord4= -1.; break;
     369                case 1: coord1=0.; coord2=1.; coord3=0.; coord4= -1.; break;
     370                case 2: coord1=0.; coord2=0.; coord3=1.; coord4= -1.; break;
     371                case 3: coord1=1.; coord2=0.; coord3=0.; coord4= +1.; break;
     372                case 4: coord1=0.; coord2=1.; coord3=0.; coord4= +1.; break;
     373                case 5: coord1=0.; coord2=0.; coord3=1.; coord4= +1.; break;
     374                default: _error_("vertex index should be in [0 5]");
    388375
    389376        }
    … …  
    405392        else{
    406393                _error_("Tria not supported yet");
     394        }
     395
     396}
     397/*}}}*/
     398/*FUNCTION GaussPenta::GaussNode{{{*/
     399void GaussPenta::GaussNode(int numnodes,int iv){
     400
     401        /*in debugging mode: check that the default constructor has been called*/
     402        _assert_(numgauss==-1);
     403
     404        /*update static arrays*/
     405        switch(numnodes){
     406                case 6: //P1 Lagrange element
     407                        switch(iv){
     408                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     409                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     410                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     411                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     412                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     413                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     414                                default: _error_("node index should be in [0 5]");
     415                        }
     416                        break;
     417                case 15: //P2 Lagrange element
     418                        switch(iv){
     419                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     420                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     421                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     422                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     423                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     424                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     425
     426                                case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
     427                                case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
     428                                case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
     429
     430                                case  9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
     431                                case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
     432                                case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
     433                                case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
     434                                case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
     435                                case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
     436                                default: _error_("node index should be in [0 5]");
     437                        }
     438                        break;
     439                default: _error_("supported number of nodes are 6 and 15");
    407440        }
    408441

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r15538 r15625  
    4444                void GaussPoint(int ig);
    4545                void GaussVertex(int iv);
     46                void GaussNode(int numnodes,int iv);
    4647                void GaussFaceTria(int index1, int index2, int index3, int order);
    4748                void GaussCenter(void);
Note: See TracChangeset for help on using the changeset viewer.