source: issm/oecreview/Archive/16133-16554/ISSM-16349-16350.diff

Last change on this file was 16556, checked in by Mathieu Morlighem, 11 years ago

NEW: added Archive/16133-16554
File size: 36.7 KB

  • ../trunk-jpl/src/c/classes/Materials/Matice.cpp

     
    291291        *pviscosity=viscosity;
    292292}
    293293/*}}}*/
     294/*FUNCTION Matice::GetViscosity2dvertical {{{*/
     295void  Matice::GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* epsilon){
     296        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
     297                                                                           (1-D) B
     298          viscosity= --------------------------------------
     299                                                  2[ exx^2+eyy^2+ 2exy^2]^[(n-1)/2n]
     300
     301          where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
     302          vector, and n the flow law exponent.
     303
     304          If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
     305          return 10^14, initial viscosity.
     306          */
     307
     308        /*output: */
     309        IssmDouble viscosity;
     310
     311        /*input strain rate: */
     312        IssmDouble exx,eyy,exy;
     313
     314        /*Intermediary: */
     315        IssmDouble A,e;
     316        IssmDouble B,D,n;
     317
     318        /*Get B and n*/
     319        B=GetB();
     320        n=GetN();
     321        D=GetD();
     322
     323        if (n==1){
     324                /*Viscous behaviour! viscosity=B: */
     325                viscosity=(1-D)*B/2;
     326        }
     327        else{
     328                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
     329                        viscosity=0.5*pow(10.,14);
     330                }
     331                else{
     332                        /*Retrive strain rate components: */
     333                        exx=epsilon[0];
     334                        eyy=epsilon[1];
     335                        exy=epsilon[2];
     336
     337                        /*Build viscosity: viscosity=B/(2*A^e) */
     338                        A=exx*exx+eyy*eyy+2.*exy*exy;
     339                        if(A==0.){
     340                                /*Maxiviscositym viscosity for 0 shear areas: */
     341                                viscosity=2.5*2.e+17;
     342                        }
     343                        else{
     344                                e=(n-1.)/(2.*n);
     345                                viscosity=(1.-D)*B/(2.*pow(A,e));
     346                        }
     347                }
     348        }
     349
     350        /*Checks in debugging mode*/
     351        if(viscosity<=0) _error_("Negative viscosity");
     352        _assert_(B>0);
     353        _assert_(n>0);
     354        _assert_(D>=0 && D<1);
     355
     356        /*Return: */
     357        *pviscosity=viscosity;
     358}
     359/*}}}*/
    294360/*FUNCTION Matice::GetViscosity3d {{{*/
    295361void  Matice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){
    296362
    … …  
    792858        if(iomodel->meshtype==Mesh2DhorizontalEnum){
    793859
    794860                /*Intermediaries*/
    795                 const int num_vertices = 3; //Tria has 3 vertices
    796                 IssmDouble    nodeinputs[num_vertices];
    797                 IssmDouble    cmmininputs[num_vertices];
    798                 IssmDouble    cmmaxinputs[num_vertices];
     861                const int  num_vertices = 3; //Tria has 3 vertices
     862                IssmDouble nodeinputs[num_vertices];
     863                IssmDouble cmmininputs[num_vertices];
     864                IssmDouble cmmaxinputs[num_vertices];
    799865
    800866                /*Get B*/
    801867                if (iomodel->Data(MaterialsRheologyBEnum)) {
    … …  
    866932                /*Get D:*/
    867933                if (iomodel->Data(DamageDEnum)) {
    868934                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
    869                         this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum));
     935                        this->inputs->AddInput(new TriaInput(DamageDEnum,nodeinputs,P1Enum));
    870936                }
    871937        }
    872938        #ifdef _HAVE_3D_

  • ../trunk-jpl/src/c/classes/Materials/Material.h

     
    2626                virtual void       Configure(Elements* elements)=0;
    2727                virtual void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum)=0;
    2828                virtual void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
     29                virtual void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
    2930                virtual void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0;
    3031                virtual void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0;
    3132                virtual void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;

  • ../trunk-jpl/src/c/classes/Materials/Matice.h

     
    5353                void   GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum);
    5454                void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    5555                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
     56                void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon);
    5657                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
    5758                void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon);
    5859                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);

  • ../trunk-jpl/src/c/classes/Materials/Matpar.h

     
    7979                void   InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    8080                /*}}}*/
    8181                /*Material virtual functions resolution: {{{*/
    82                 void   InputDuplicate(int original_enum,int new_enum);
    83                 void   Configure(Elements* elements);
    84                 void   GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
     82                void       InputDuplicate(int original_enum,int new_enum);
     83                void       Configure(Elements* elements);
     84                void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
    8585                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
     86                void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
    8687                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");};
    8788                void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};
    8889                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};

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

     
    1919/*}}}*/
    2020
    2121/*Element macros*/
    22 #define NUMVERTICES 3
     22#define NUMVERTICES   3
     23#define NUMVERTICES1D 2
     24
    2325/*Constructors/destructor/copy*/
    2426/*FUNCTION Tria::Tria(){{{*/
    2527Tria::Tria(){
    … …  
    218220        switch(analysis_type){
    219221                #ifdef _HAVE_STRESSBALANCE_
    220222                case StressbalanceAnalysisEnum:
    221                         return CreateKMatrixStressbalanceSSA();
     223                        int approximation;
     224                        inputs->GetInputValue(&approximation,ApproximationEnum);
     225                        switch(approximation){
     226                                case SSAApproximationEnum:
     227                                        return CreateKMatrixStressbalanceSSA();
     228                                case FSApproximationEnum:
     229                                        return CreateKMatrixStressbalanceFS();
     230                                case NoneApproximationEnum:
     231                                        return NULL;
     232                                default:
     233                                        _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
     234                        }
    222235                        break;
    223236                case StressbalanceSIAAnalysisEnum:
    224237                        return CreateKMatrixStressbalanceSIA();
    … …  
    372385        switch(analysis_type){
    373386                #ifdef _HAVE_STRESSBALANCE_
    374387                case StressbalanceAnalysisEnum:
    375                         return CreatePVectorStressbalanceSSA();
     388                        int approximation;
     389                        inputs->GetInputValue(&approximation,ApproximationEnum);
     390                        switch(approximation){
     391                                case SSAApproximationEnum:
     392                                        return CreatePVectorStressbalanceSSA();
     393                                case FSApproximationEnum:
     394                                        return CreatePVectorStressbalanceFS();
     395                                case NoneApproximationEnum:
     396                                        return NULL;
     397                                default:
     398                                        _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
     399                        }
    376400                        break;
    377401                case StressbalanceSIAAnalysisEnum:
    378402                        return CreatePVectorStressbalanceSIA();
    … …  
    16781702        switch(analysis_type){
    16791703                #ifdef _HAVE_STRESSBALANCE_
    16801704                case StressbalanceAnalysisEnum:
    1681                         InputUpdateFromSolutionStressbalanceHoriz(solution);
     1705                        int approximation;
     1706                        inputs->GetInputValue(&approximation,ApproximationEnum);
     1707                        if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
     1708                                InputUpdateFromSolutionStressbalanceFS(solution);
     1709                        }
     1710                        else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){
     1711                                InputUpdateFromSolutionStressbalanceHoriz(solution);
     1712                        }
     1713                        else{
     1714                                _error_("approximation not supported yet");
     1715                        }
    16821716                        break;
    16831717                case StressbalanceSIAAnalysisEnum:
    16841718                        InputUpdateFromSolutionStressbalanceHoriz(solution);
    … …  
    20552089        return onbed;
    20562090}
    20572091/*}}}*/
     2092/*FUNCTION Tria::HasEdgeOnBed {{{*/
     2093bool Tria::HasEdgeOnBed(){
     2094
     2095        IssmDouble values[NUMVERTICES];
     2096        IssmDouble sum;
     2097
     2098        /*Retrieve all inputs and parameters*/
     2099        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
     2100        sum = values[0]+values[1]+values[2];
     2101
     2102        _assert_(sum==0. || sum==1. || sum==2.);
     2103
     2104        if(sum>1.){
     2105                return true;
     2106        }
     2107        else{
     2108                return false;
     2109        }
     2110}
     2111/*}}}*/
     2112/*FUNCTION Tria::EdgeOnBedIndices{{{*/
     2113void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
     2114
     2115        bool found=false;
     2116        IssmDouble values[NUMVERTICES];
     2117
     2118        /*Retrieve all inputs and parameters*/
     2119        GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
     2120
     2121        for(int i=0;i<NUMVERTICES;i++){
     2122                if(values[i]==1.){
     2123                        if(found){
     2124                                *pindex2 = i;
     2125                                return;
     2126                        }
     2127                        else{
     2128                                *pindex1 = i;
     2129                        }
     2130                }
     2131        }
     2132
     2133        _error_("Could not find 2 vertices on bed");
     2134}
     2135/*}}}*/
    20582136/*FUNCTION Tria::IsFloating {{{*/
    20592137bool   Tria::IsFloating(){
    20602138
    … …  
    30363114#endif
    30373115
    30383116#ifdef _HAVE_STRESSBALANCE_
     3117/*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/
     3118ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){
     3119
     3120        ElementMatrix* Ke1 = NULL;
     3121        ElementMatrix* Ke2 = NULL;
     3122        ElementMatrix* Ke  = NULL;
     3123
     3124        /*compute all stiffness matrices for this element*/
     3125        Ke1=CreateKMatrixStressbalanceFSViscous();
     3126        Ke2=CreateKMatrixStressbalanceFSFriction();
     3127        Ke =new ElementMatrix(Ke1,Ke2);
     3128
     3129        /*clean-up and return*/
     3130        delete Ke1;
     3131        delete Ke2;
     3132        return Ke;
     3133
     3134}
     3135/*}}}*/
     3136/*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/
     3137ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){
     3138
     3139        /*Intermediaries */
     3140        int        i,approximation;
     3141        IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
     3142        IssmDouble xyz_list[NUMVERTICES][3];
     3143        IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     3144        GaussTria *gauss=NULL;
     3145
     3146        /*Fetch number of nodes and dof for this finite element*/
     3147        int vnumnodes = this->NumberofNodesVelocity();
     3148        int pnumnodes = this->NumberofNodesPressure();
     3149        int numdof    = vnumnodes*NDOF2 + pnumnodes*NDOF1;
     3150
     3151        /*Prepare coordinate system list*/
     3152        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3153        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
     3154        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3155
     3156        /*Initialize Element matrix and vectors*/
     3157        ElementMatrix* Ke     = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
     3158        IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
     3159        IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
     3160        IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
     3161
     3162        /*Retrieve all inputs and parameters*/
     3163        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3164        parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3165        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3166        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3167
     3168        /* Start  looping on the number of gaussian points: */
     3169        gauss=new GaussTria(5);
     3170        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3171
     3172                gauss->GaussPoint(ig);
     3173
     3174                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3175                GetBFS(B,&xyz_list[0][0],gauss);
     3176                GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
     3177
     3178                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3179                material->GetViscosity2dvertical(&viscosity,&epsilon[0]);
     3180
     3181                D_scalar=gauss->weight*Jdet;
     3182                for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity;
     3183                for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning;
     3184
     3185                TripleMultiply(B,5,numdof,1,
     3186                                        D,5,5,0,
     3187                                        Bprime,5,numdof,0,
     3188                                        Ke->values,1);
     3189        }
     3190
     3191        /*Transform Coordinate System*/
     3192        TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
     3193
     3194        /*Clean up and return*/
     3195        delete gauss;
     3196        xDelete<int>(cs_list);
     3197        xDelete<IssmDouble>(B);
     3198        xDelete<IssmDouble>(Bprime);
     3199        xDelete<IssmDouble>(D);
     3200        return Ke;
     3201}
     3202/*}}}*/
     3203/*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/
     3204ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){
     3205
     3206        /*Intermediaries */
     3207        int         i,j;
     3208        int         analysis_type,approximation;
     3209        IssmDouble  alpha2,Jdet2d;
     3210        IssmDouble  FSreconditioning,viscosity;
     3211        IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     3212        IssmDouble  xyz_list[NUMVERTICES][3];
     3213        IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
     3214        Friction   *friction = NULL;
     3215        GaussTria  *gauss    = NULL;
     3216
     3217        /*Initialize Element matrix and return if necessary*/
     3218        if(IsFloating() || !HasEdgeOnBed()) return NULL;
     3219
     3220        _error_("STOP");
     3221}
     3222/*}}}*/
    30393223/*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/
    30403224ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){
    30413225
    … …  
    32253409        return Ke;
    32263410}
    32273411/*}}}*/
     3412/*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/
     3413ElementVector* Tria::CreatePVectorStressbalanceFS(void){
     3414
     3415        ElementVector* pe1;
     3416        ElementVector* pe2;
     3417        ElementVector* pe3;
     3418        ElementVector* pe;
     3419
     3420        /*compute all stiffness matrices for this element*/
     3421        pe1=CreatePVectorStressbalanceFSViscous();
     3422        pe2=CreatePVectorStressbalanceFSShelf();
     3423        pe3=CreatePVectorStressbalanceFSFront();
     3424        pe =new ElementVector(pe1,pe2,pe3);
     3425
     3426        /*clean-up and return*/
     3427        delete pe1;
     3428        delete pe2;
     3429        delete pe3;
     3430        return pe;
     3431}
     3432/*}}}*/
     3433/*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/
     3434ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){
     3435
     3436        /*Intermediaries */
     3437        int         i;
     3438        IssmDouble  ls[NUMVERTICES];
     3439        IssmDouble  xyz_list[NUMVERTICES][3];
     3440        bool        isfront;
     3441
     3442        /*Retrieve all inputs and parameters*/
     3443        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
     3444
     3445        /*If the level set is awlays <=0, there is no ice front here*/
     3446        isfront = false;
     3447        if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
     3448                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
     3449                        isfront = true;
     3450                }
     3451        }
     3452
     3453        /*If no front, return NULL*/
     3454        if(!isfront) return NULL;
     3455
     3456        _error_("STOP");
     3457}
     3458/*}}}*/
     3459/*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/
     3460ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){
     3461
     3462        /*Intermediaries*/
     3463        int        i;
     3464        int        approximation;
     3465        IssmDouble Jdet,gravity,rho_ice;
     3466        IssmDouble forcex,forcey;
     3467        IssmDouble xyz_list[NUMVERTICES][3];
     3468        GaussTria *gauss=NULL;
     3469
     3470        /*Fetch number of nodes and dof for this finite element*/
     3471        int vnumnodes = this->NumberofNodesVelocity();
     3472        int pnumnodes = this->NumberofNodesPressure();
     3473
     3474        /*Prepare coordinate system list*/
     3475        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3476        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
     3477        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3478
     3479        /*Initialize Element matrix and vectors*/
     3480        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
     3481        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3482
     3483        /*Retrieve all inputs and parameters*/
     3484        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3485        Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     3486        Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     3487        rho_ice = matpar->GetRhoIce();
     3488        gravity = matpar->GetG();
     3489
     3490        /* Start  looping on the number of gaussian points: */
     3491        gauss=new GaussTria(5);
     3492        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3493
     3494                gauss->GaussPoint(ig);
     3495
     3496                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     3497                GetNodalFunctionsVelocity(vbasis, gauss);
     3498
     3499                loadingforcex_input->GetInputValue(&forcex,gauss);
     3500                loadingforcey_input->GetInputValue(&forcey,gauss);
     3501
     3502                for(i=0;i<vnumnodes;i++){
     3503                        pe->values[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     3504                        pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     3505                        pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     3506                }
     3507        }
     3508
     3509        /*Transform coordinate system*/
     3510        TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
     3511
     3512        /*Clean up and return*/
     3513        delete gauss;
     3514        xDelete<int>(cs_list);
     3515        xDelete<IssmDouble>(vbasis);
     3516        return pe;
     3517}
     3518/*}}}*/
     3519/*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/
     3520ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){
     3521
     3522        /*Intermediaries*/
     3523        int         i,j;
     3524        int         indices[2];
     3525        IssmDouble  gravity,rho_water,bed,water_pressure;
     3526        IssmDouble  normal_vel,vx,vy,dt;
     3527        IssmDouble      xyz_list_seg[NUMVERTICES1D][3];
     3528        IssmDouble  xyz_list[NUMVERTICES][3];
     3529        IssmDouble      normal[2];
     3530        IssmDouble  Jdet;
     3531
     3532        /*Initialize Element vector and return if necessary*/
     3533        if(!HasEdgeOnBed() || !IsFloating()) return NULL;
     3534
     3535        /*Fetch number of nodes and dof for this finite element*/
     3536        int vnumnodes = this->NumberofNodesVelocity();
     3537        int pnumnodes = this->NumberofNodesPressure();
     3538
     3539        /*Prepare coordinate system list*/
     3540        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3541        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
     3542        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3543
     3544        /*Initialize Element matrix and vectors*/
     3545        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
     3546        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3547
     3548        /*Retrieve all inputs and parameters*/
     3549        rho_water=matpar->GetRhoWater();
     3550        gravity=matpar->GetG();
     3551        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3552        Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input);
     3553
     3554        /*Get vertex indices that lie on bed*/
     3555        this->EdgeOnBedIndices(&indices[0],&indices[1]);
     3556
     3557        for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
     3558
     3559        /* Start looping on the number of gauss 1d (nodes on the bedrock) */
     3560        GaussTria* gauss=new GaussTria(indices[0],indices[1],2);
     3561        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3562
     3563                gauss->GaussPoint(ig);
     3564
     3565                GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
     3566                GetNodalFunctionsVelocity(vbasis, gauss);
     3567
     3568                GetSegmentNormal(&normal[0],xyz_list_seg);
     3569                _assert_(normal[1]<0.);
     3570                bed_input->GetInputValue(&bed, gauss);
     3571                water_pressure=gravity*rho_water*bed;
     3572
     3573                for(i=0;i<vnumnodes;i++){
     3574                        pe->values[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0];
     3575                        pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1];
     3576                }
     3577        }
     3578
     3579        /*Transform coordinate system*/
     3580        TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
     3581
     3582        /*Clean up and return*/
     3583        delete gauss;
     3584        xDelete<int>(cs_list);
     3585        xDelete<IssmDouble>(vbasis);
     3586        return pe;
     3587}
     3588/*}}}*/
    32283589/*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/
    32293590ElementVector* Tria::CreatePVectorStressbalanceSSA(){
    32303591
    … …  
    36954056
    36964057}
    36974058/*}}}*/
     4059/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/
     4060void  Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
     4061
     4062        int          i;
     4063        int*         vdoflist=NULL;
     4064        int*         pdoflist=NULL;
     4065        IssmDouble   FSreconditioning;
     4066
     4067        /*Fetch number of nodes and dof for this finite element*/
     4068        int vnumnodes = this->NumberofNodesVelocity();
     4069        int pnumnodes = this->NumberofNodesPressure();
     4070        int vnumdof   = vnumnodes*NDOF2;
     4071        int pnumdof   = pnumnodes*NDOF1;
     4072
     4073        /*Initialize values*/
     4074        IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
     4075        IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
     4076        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
     4077        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
     4078        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
     4079        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
     4080
     4081        /*Get dof list: */
     4082        GetDofListVelocity(&vdoflist,GsetEnum);
     4083        GetDofListPressure(&pdoflist,GsetEnum);
     4084
     4085        /*Use the dof list to index into the solution vector: */
     4086        for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
     4087        for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
     4088
     4089        /*Transform solution in Cartesian Space*/
     4090        TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYEnum);
     4091
     4092        /*Ok, we have vx and vy in values, fill in all arrays: */
     4093        for(i=0;i<vnumnodes;i++){
     4094                vx[i] = vvalues[i*NDOF2+0];
     4095                vy[i] = vvalues[i*NDOF2+1];
     4096                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     4097                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     4098        }
     4099        for(i=0;i<pnumnodes;i++){
     4100                pressure[i] = pvalues[i];
     4101                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
     4102        }
     4103
     4104        /*Recondition pressure and compute vel: */
     4105        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     4106        for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
     4107        for(i=0;i<vnumnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
     4108
     4109        /*Now, we have to move the previous inputs  to old
     4110         * status, otherwise, we'll wipe them off: */
     4111        this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
     4112        this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
     4113        this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
     4114
     4115        /*Add vx and vy as inputs to the tria element: */
     4116        this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
     4117        this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
     4118        this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
     4119        this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
     4120
     4121        /*Free ressources:*/
     4122        xDelete<IssmDouble>(pressure);
     4123        xDelete<IssmDouble>(vel);
     4124        xDelete<IssmDouble>(vy);
     4125        xDelete<IssmDouble>(vx);
     4126        xDelete<IssmDouble>(vvalues);
     4127        xDelete<IssmDouble>(pvalues);
     4128        xDelete<int>(vdoflist);
     4129        xDelete<int>(pdoflist);
     4130}
     4131/*}}}*/
    36984132/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceSIA {{{*/
    36994133void  Tria::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
    37004134

  • ../trunk-jpl/src/c/classes/Elements/TriaRef.h

     
    2222                void SetElementType(int type,int type_counter);
    2323
    2424                /*Numerics*/
     25                void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
    2526                void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
    2627                void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
     28                void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
    2729                void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    2830                void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    2931                void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    … …  
    3537                void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
    3638                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
    3739                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
     40                void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
     41                void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
    3842                void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
    3943                void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
    4044                void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
    4145                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
     46                void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
     47                void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
    4248                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
    4349                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
    4450                void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);

  • ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp

     
    333333
    334334        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
    335335         * For node i, Bvi can be expressed in the actual coordinate system
    336          * by:     Bvi=[ dh/dx          0              0      ]
     336         * by:     Bvi=[ dh/dx          0             0      ]
    337337         *                                      [   0           dh/dy           0      ]
    338          *                                      [   0             0           dh/dy    ]
     338         *                                      [   0             0           dh/dz    ]
    339339         *                                      [ 1/2*dh/dy    1/2*dh/dx        0      ]
    340340         *                                      [ 1/2*dh/dz       0         1/2*dh/dx  ]
    341341         *                                      [   0          1/2*dh/dz    1/2*dh/dy  ]

  • ../trunk-jpl/src/c/classes/Elements/Tria.h

     
    8383                int         GetNumberOfNodes(void);
    8484                int         Sid();
    8585                bool        IsOnBed();
     86                bool        HasEdgeOnBed();
     87                void        EdgeOnBedIndices(int* pindex1,int* pindex);
    8688                bool        IsFloating();
    8789                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
    8890                bool        NoIceInElement();
    … …  
    232234                ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void);
    233235                ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void);
    234236                ElementMatrix* CreateKMatrixStressbalanceSIA(void);
     237                ElementMatrix* CreateKMatrixStressbalanceFS(void);
     238                ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
     239                ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
    235240                ElementVector* CreatePVectorStressbalanceSSA(void);
    236241                ElementVector* CreatePVectorStressbalanceSSADrivingStress(void);
    237242                ElementVector* CreatePVectorStressbalanceSSAFront(void);
    238243                ElementVector* CreatePVectorStressbalanceSIA(void);
     244                ElementVector* CreatePVectorStressbalanceFS(void);
     245                ElementVector* CreatePVectorStressbalanceFSFront(void);
     246                ElementVector* CreatePVectorStressbalanceFSViscous(void);
     247                void           PVectorGLSstabilization(ElementVector* pe);
     248                ElementVector* CreatePVectorStressbalanceFSShelf(void);
    239249                ElementMatrix* CreateJacobianStressbalanceSSA(void);
    240250                void      GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
    241251                void      GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
    242252                void      GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
    243253                void      InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
     254                void      InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
    244255                void      InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
    245256                #endif
    246257

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

     
    1065710657        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    1065810658        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
    1065910659        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    10660         IssmDouble* pressure = xNew<IssmDouble>(vnumnodes);
     10660        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
    1066110661
    1066210662        /*Get dof list: */
    1066310663        GetDofListVelocity(&vdoflist,GsetEnum);

  • ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp

     
    199199        xDelete<IssmDouble>(basis);
    200200}
    201201/*}}}*/
     202/*FUNCTION TriaRef::GetBFS {{{*/
     203void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
     204        /*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
     205         * For node i, Bvi can be expressed in the actual coordinate system
     206         * by:     Bvi=[ dphi/dx          0        ]
     207         *                                       [   0           dphi/dy     ]
     208         *                                       [ 1/2*dphi/dy    1/2*dphi/dx]
     209         *                                       [   0             0         ]
     210         *                                       [ dphi/dx         dphi/dy   ]
     211         *
     212         * by:    Bpi=[  0    ]
     213         *                                      [  0    ]
     214         *                                      [  0    ]
     215         *                                      [ phi_p ]
     216         *                                      [  0    ]
     217         *      where phi is the interpolation function for node i.
     218         *      Same thing for Bb except the last column that does not exist.
     219         */
     220
     221        /*Fetch number of nodes for this finite element*/
     222        int pnumnodes = this->NumberofNodesPressure();
     223        int vnumnodes = this->NumberofNodesVelocity();
     224
     225        /*Get nodal functions derivatives*/
     226        IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
     227        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     228        GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     229        GetNodalFunctionsPressure(pbasis,gauss);
     230
     231        /*Build B: */
     232        for(int i=0;i<vnumnodes;i++){
     233                B[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
     234                B[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
     235                B[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
     236                B[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
     237                B[(2*vnumnodes+pnumnodes)*2+2*i+0] = .5*vdbasis[1*vnumnodes+i];
     238                B[(2*vnumnodes+pnumnodes)*2+2*i+1] = .5*vdbasis[0*vnumnodes+i];
     239                B[(2*vnumnodes+pnumnodes)*3+2*i+0] = 0.;
     240                B[(2*vnumnodes+pnumnodes)*3+2*i+1] = 0.;
     241                B[(2*vnumnodes+pnumnodes)*4+2*i+0] = vdbasis[0*vnumnodes+i];
     242                B[(2*vnumnodes+pnumnodes)*4+2*i+1] = vdbasis[1*vnumnodes+i];
     243        }
     244        for(int i=0;i<pnumnodes;i++){
     245                B[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
     246                B[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
     247                B[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
     248                B[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = pbasis[i];
     249                B[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = 0.;
     250        }
     251
     252        /*Clean up*/
     253        xDelete<IssmDouble>(vdbasis);
     254        xDelete<IssmDouble>(pbasis);
     255}
     256/*}}}*/
     257/*FUNCTION TriaRef::GetBprimeFS {{{*/
     258void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss){
     259        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
     260         *      For node i, Bi' can be expressed in the actual coordinate system
     261         *      by:
     262         *                      Bvi' = [  dphi/dx     0     ]
     263         *                                       [     0      dphi/dy ]
     264         *                                       [  dphi/dy   dphi/dx ]
     265         *                                       [  dphi/dx   dphi/dy ]
     266         *                                       [     0      0       ]
     267         *
     268         * by:    Bpi=[  0  ]
     269         *                                      [  0  ]
     270         *                                      [  0  ]
     271         *                                      [  0  ]
     272         *                                      [ phi ]
     273         *      where phi is the interpolation function for node i.
     274         */
     275
     276        /*Fetch number of nodes for this finite element*/
     277        int pnumnodes = this->NumberofNodesPressure();
     278        int vnumnodes = this->NumberofNodesVelocity();
     279
     280        /*Get nodal functions derivatives*/
     281        IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
     282        IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
     283        GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     284        GetNodalFunctionsPressure(pbasis,gauss);
     285
     286        /*Build B_prime: */
     287        for(int i=0;i<vnumnodes;i++){
     288                B_prime[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
     289                B_prime[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
     290                B_prime[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
     291                B_prime[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
     292                B_prime[(2*vnumnodes+pnumnodes)*2+2*i+0] = vdbasis[1*vnumnodes+i];
     293                B_prime[(2*vnumnodes+pnumnodes)*2+2*i+1] = vdbasis[0*vnumnodes+i];
     294                B_prime[(2*vnumnodes+pnumnodes)*3+2*i+0] = vdbasis[0*vnumnodes+i];
     295                B_prime[(2*vnumnodes+pnumnodes)*3+2*i+1] = vdbasis[1*vnumnodes+i];
     296                B_prime[(2*vnumnodes+pnumnodes)*4+2*i+0] = 0.;
     297                B_prime[(2*vnumnodes+pnumnodes)*4+2*i+1] = 0.;
     298        }
     299        for(int i=0;i<pnumnodes;i++){
     300                B_prime[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
     301                B_prime[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
     302                B_prime[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
     303                B_prime[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = 0.;
     304                B_prime[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = pbasis[i];
     305        }
     306
     307        /*Clean up*/
     308        xDelete<IssmDouble>(vdbasis);
     309        xDelete<IssmDouble>(pbasis);
     310}
     311/*}}}*/
    202312/*FUNCTION TriaRef::GetBMasstransport{{{*/
    203313void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
    204314        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    … …  
    457567        }
    458568}
    459569/*}}}*/
     570/*FUNCTION TriaRef::GetNodalFunctionsVelocity{{{*/
     571void TriaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussTria* gauss){
     572        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     573
     574        switch(this->element_type){
     575                case P1P1Enum:
     576                        this->element_type = P1Enum;
     577                        this->GetNodalFunctions(basis,gauss);
     578                        this->element_type = P1P1Enum;
     579                        return;
     580                case P1P1GLSEnum:
     581                        this->element_type = P1Enum;
     582                        this->GetNodalFunctions(basis,gauss);
     583                        this->element_type = P1P1GLSEnum;
     584                        return;
     585                case MINIcondensedEnum:
     586                        this->element_type = P1bubbleEnum;
     587                        this->GetNodalFunctions(basis,gauss);
     588                        this->element_type = MINIcondensedEnum;
     589                        return;
     590                case MINIEnum:
     591                        this->element_type = P1bubbleEnum;
     592                        this->GetNodalFunctions(basis,gauss);
     593                        this->element_type = MINIEnum;
     594                        return;
     595                case TaylorHoodEnum:
     596                        this->element_type = P2Enum;
     597                        this->GetNodalFunctions(basis,gauss);
     598                        this->element_type = TaylorHoodEnum;
     599                        return;
     600                default:
     601                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     602        }
     603}
     604/*}}}*/
     605/*FUNCTION TriaRef::GetNodalFunctionsPressure{{{*/
     606void TriaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussTria* gauss){
     607        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     608
     609        switch(this->element_type){
     610                case P1P1Enum:
     611                        this->element_type = P1Enum;
     612                        this->GetNodalFunctions(basis,gauss);
     613                        this->element_type = P1P1Enum;
     614                        return;
     615                case P1P1GLSEnum:
     616                        this->element_type = P1Enum;
     617                        this->GetNodalFunctions(basis,gauss);
     618                        this->element_type = P1P1GLSEnum;
     619                        return;
     620                case MINIcondensedEnum:
     621                        this->element_type = P1Enum;
     622                        this->GetNodalFunctions(basis,gauss);
     623                        this->element_type = MINIcondensedEnum;
     624                        return;
     625                case MINIEnum:
     626                        this->element_type = P1Enum;
     627                        this->GetNodalFunctions(basis,gauss);
     628                        this->element_type = MINIEnum;
     629                        return;
     630                case TaylorHoodEnum:
     631                        this->element_type = P1Enum;
     632                        this->GetNodalFunctions(basis,gauss);
     633                        this->element_type = TaylorHoodEnum;
     634                        return;
     635                default:
     636                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     637        }
     638}
     639/*}}}*/
    460640/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
    461641void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
    462642        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    528708
    529709}
    530710/*}}}*/
     711/*FUNCTION TriaRef::GetNodalFunctionsDerivativesPressure{{{*/
     712void TriaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
     713        switch(this->element_type){
     714                case P1P1Enum:
     715                        this->element_type = P1Enum;
     716                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     717                        this->element_type = P1P1Enum;
     718                        return;
     719                case P1P1GLSEnum:
     720                        this->element_type = P1Enum;
     721                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     722                        this->element_type = P1P1GLSEnum;
     723                        return;
     724                case MINIcondensedEnum:
     725                        this->element_type = P1Enum;
     726                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     727                        this->element_type = MINIcondensedEnum;
     728                        return;
     729                case MINIEnum:
     730                        this->element_type = P1Enum;
     731                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     732                        this->element_type = MINIEnum;
     733                        return;
     734                case TaylorHoodEnum:
     735                        this->element_type = P1Enum;
     736                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     737                        this->element_type = TaylorHoodEnum;
     738                        return;
     739                default:
     740                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     741        }
     742}
     743/*}}}*/
     744/*FUNCTION TriaRef::GetNodalFunctionsDerivativesVelocity{{{*/
     745void TriaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
     746        switch(this->element_type){
     747                case P1P1Enum:
     748                        this->element_type = P1Enum;
     749                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     750                        this->element_type = P1P1Enum;
     751                        return;
     752                case P1P1GLSEnum:
     753                        this->element_type = P1Enum;
     754                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     755                        this->element_type = P1P1GLSEnum;
     756                        return;
     757                case MINIcondensedEnum:
     758                        this->element_type = P1bubbleEnum;
     759                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     760                        this->element_type = MINIcondensedEnum;
     761                        return;
     762                case MINIEnum:
     763                        this->element_type = P1bubbleEnum;
     764                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     765                        this->element_type = MINIEnum;
     766                        return;
     767                case TaylorHoodEnum:
     768                        this->element_type = P2Enum;
     769                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     770                        this->element_type = TaylorHoodEnum;
     771                        return;
     772                default:
     773                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     774        }
     775}
     776/*}}}*/
    531777/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
    532778void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){
    533779        /*This routine returns the values of the nodal functions derivatives  (with respect to the
Note: See TracBrowser for help on using the repository browser.