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/Elements
Files:
3 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: */
Note: See TracChangeset for help on using the changeset viewer.