Changeset 25514


Ignore:
Timestamp:
09/02/20 12:00:13 (5 years ago)
Author:
Cheng Gong
Message:

NEW: Galerkin Least Squares stabilization for P1P1 elements in FS is added together with a Nitsches method for weakly imposing Dirichlet boundary condition on the bottom surface of the ice. Current issue: the nonlinear iteration may fail to converge if the number of vertical layers is too small(~5).

Location:
issm/trunk-jpl/src
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r25451 r25514  
    929929        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isFS",FlowequationIsFSEnum));
    930930        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.fe_FS",FlowequationFeFSEnum));
     931        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.isNitscheBC",FlowequationIsNitscheEnum));
     932        parameters->AddObject(iomodel->CopyConstantObject("md.flowequation.FSNitscheGamma",FeFSNitscheGammaEnum));
    931933        parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.restol",StressbalanceRestolEnum));
    932934        parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.reltol",StressbalanceReltolEnum));
    … …  
    34333435        /*Get type of algorithm*/
    34343436        int fe_FS;
     3437        bool isNitsche;
     3438
    34353439        element->FindParam(&fe_FS,FlowequationFeFSEnum);
     3440        element->FindParam(&isNitsche,FlowequationIsNitscheEnum);
    34363441
    34373442        /*compute all stiffness matrices for this element*/
    … …  
    34413446        else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
    34423447         Ke1=CreateKMatrixFSViscousLA(element);
     3448        else if(fe_FS==P1P1GLSEnum)
     3449         Ke1=CreateKMatrixFSViscousGLS(element);
    34433450        else
    34443451         Ke1=CreateKMatrixFSViscous(element);
    34453452
    3446         ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
     3453        ElementMatrix* Ke2;
     3454        if (isNitsche) {
     3455                Ke2 = CreateKMatrixFSFrictionNitsche(element);
     3456        }
     3457        else {
     3458                Ke2 = CreateKMatrixFSFriction(element);
     3459        }
    34473460        ElementMatrix* Ke3=CreateKMatrixFSShelf(element);
    34483461        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    … …  
    35643577                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    35653578
    3566                 if(dim==3){
     3579                if(dim==2 || dim==3){
    35673580                        /*Stress balance*/
    35683581                        for(int i=0;i<vnumnodes;i++){
    35693582                                for(int j=0;j<vnumnodes;j++){
    3570                                         Ke->values[(3*i+0)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
    3571                                         Ke->values[(3*i+0)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
    3572                                         Ke->values[(3*i+0)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
    3573                                         Ke->values[(3*i+1)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
    3574                                         Ke->values[(3*i+1)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
    3575                                         Ke->values[(3*i+1)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
    3576                                         Ke->values[(3*i+2)*numdof+3*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
    3577                                         Ke->values[(3*i+2)*numdof+3*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[2*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
    3578                                         Ke->values[(3*i+2)*numdof+3*j+2] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i] + 2.*vdbasis[2*vnumnodes+j]*vdbasis[2*vnumnodes+i]);
     3583                                        for (int p=0;p<dim;p++){
     3584                                                for (int q=0;q<dim;q++){
     3585                                                        /* diagonal only */
     3586                                                        Ke->values[(dim*i+p)*numdof+dim*j+p] += gauss->weight*Jdet*viscosity*(vdbasis[q*vnumnodes+j]*vdbasis[q*vnumnodes+i]);
     3587                                                        /* All the entries */
     3588                                                        Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet*viscosity*(vdbasis[p*vnumnodes+j]*vdbasis[q*vnumnodes+i]);
     3589                                                }
     3590                                        }
    35793591                                }
    35803592                                for(int k=0;k<pnumnodes;k++){
    3581                                         Ke->values[(3*i+0)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
    3582                                         Ke->values[(3*i+1)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
    3583                                         Ke->values[(3*i+2)*numdof+3*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[2*vnumnodes+i]);
     3593                                        for (int p=0;p<dim;p++){
     3594                                                Ke->values[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[p*vnumnodes+i]);
     3595                                        }
    35843596                                }
    35853597                        }
    … …  
    35873599                        for(int k=0;k<pnumnodes;k++){
    35883600                                for(int j=0;j<vnumnodes;j++){
    3589                                         Ke->values[(3*vnumnodes+k)*numdof+3*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
    3590                                         Ke->values[(3*vnumnodes+k)*numdof+3*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
    3591                                         Ke->values[(3*vnumnodes+k)*numdof+3*j+2] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[2*vnumnodes+j]*pbasis[k]);
    3592                                 }
    3593                         }
    3594                 }
    3595                 else{
    3596                         /*Stress balance*/
    3597                         for(int i=0;i<vnumnodes;i++){
    3598                                 for(int j=0;j<vnumnodes;j++){
    3599                                         Ke->values[(2*i+0)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(2.*vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
    3600                                         Ke->values[(2*i+0)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
    3601                                         Ke->values[(2*i+1)*numdof+2*j+0] += gauss->weight*Jdet*viscosity*(vdbasis[1*vnumnodes+j]*vdbasis[0*vnumnodes+i]);
    3602                                         Ke->values[(2*i+1)*numdof+2*j+1] += gauss->weight*Jdet*viscosity*(vdbasis[0*vnumnodes+j]*vdbasis[0*vnumnodes+i] + 2.*vdbasis[1*vnumnodes+j]*vdbasis[1*vnumnodes+i]);
    3603                                 }
    3604                                 for(int k=0;k<pnumnodes;k++){
    3605                                         Ke->values[(2*i+0)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[0*vnumnodes+i]);
    3606                                         Ke->values[(2*i+1)*numdof+2*vnumnodes+k] += gauss->weight*Jdet*FSreconditioning*(-pbasis[k]*vdbasis[1*vnumnodes+i]);
    3607                                 }
    3608                         }
    3609                         /*Incompressibility*/
    3610                         for(int k=0;k<pnumnodes;k++){
    3611                                 for(int j=0;j<vnumnodes;j++){
    3612                                         Ke->values[(2*vnumnodes+k)*numdof+2*j+0] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[0*vnumnodes+j]*pbasis[k]);
    3613                                         Ke->values[(2*vnumnodes+k)*numdof+2*j+1] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[1*vnumnodes+j]*pbasis[k]);
     3601                                        for (int p=0;p<dim;p++){
     3602                                                Ke->values[(dim*vnumnodes+k)*numdof+dim*j+p] += gauss->weight*Jdet*(-FSreconditioning*vdbasis[p*vnumnodes+j]*pbasis[k]);
     3603                                        }
    36143604                                }
    36153605                        }
    … …  
    37303720        xDelete<IssmDouble>(pbasis);
    37313721        return Ke;
     3722}/*}}}*/
     3723
     3724ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousGLS(Element* element){/*{{{*/
     3725
     3726        /*Intermediaries*/
     3727        int         i,dim;
     3728        IssmDouble  viscosity,FSreconditioning,Jdet;
     3729        IssmDouble *xyz_list = NULL;
     3730
     3731        /*Get problem dimension*/
     3732        element->FindParam(&dim,DomainDimensionEnum);
     3733
     3734        /*Fetch number of nodes and dof for this finite element*/
     3735        int vnumnodes = element->NumberofNodesVelocity();
     3736        int pnumnodes = element->NumberofNodesPressure();
     3737        int numdof    = vnumnodes*dim + pnumnodes;
     3738
     3739        /* only work for P1 elements */
     3740        if (dim == 2) {_assert_(vnumnodes==3);}
     3741        else if (dim == 3) {_assert_(vnumnodes==6);}
     3742        else {_error_("GLS is not implemented except for 2D and 3D problems.");}
     3743
     3744        /*Prepare coordinate system list*/
     3745        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3746        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3747        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3748        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3749
     3750        /*Initialize Element matrix and vectors*/
     3751        ElementMatrix* Ke   = element->NewElementMatrix(FSvelocityEnum);
     3752        IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
     3753        IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
     3754
     3755        /*Retrieve all inputs and parameters*/
     3756        element->GetVerticesCoordinates(&xyz_list);
     3757        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3758        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     3759        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     3760        Input* vz_input = NULL;
     3761        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     3762       
     3763        /* prepare viscosity gradient for GLS */
     3764        IssmDouble NodalViscosity[6];
     3765        IssmDouble gradViscos[3];
     3766        IssmDouble etapq, s, Tau, mk, hk;
     3767        IssmDouble hx, hy, hz;
     3768        IssmDouble SU[3*(3+1)*6];
     3769    Gauss* vert = element->NewGauss();
     3770        /* Compute the nodal values of the viscosity */
     3771        for(int i=0;i<vnumnodes;i++){
     3772        vert->GaussNode(element->element_type, i);
     3773                element->material->ViscosityFS(&NodalViscosity[i],dim,xyz_list,vert,vx_input,vy_input,vz_input);
     3774        }
     3775        delete vert;
     3776
     3777        /* Start  looping on the number of gaussian points: */
     3778        Gauss* gauss = element->NewGauss(5);
     3779        while(gauss->next()){
     3780                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3781                element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     3782                element->NodalFunctionsPressure(pbasis,gauss);
     3783                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3784
     3785                /* compute viscosity gradient at the gaussian point */
     3786                element->ValueP1DerivativesOnGauss(&gradViscos[0],NodalViscosity,xyz_list,gauss);
     3787
     3788                /*  weight*detJ */
     3789                s = gauss->weight*Jdet;
     3790
     3791                /* GLS Stabilization */
     3792                mk = 1.0 /3.0;
     3793                element->ElementSizes(&hx,&hy,&hz);
     3794                hk = max(max(hx, hy), hz);
     3795                // if(!element->IsOnDirichletBoundary()) {
     3796                        Tau = - mk * hk * hk * 0.125 / viscosity;
     3797
     3798
     3799                for (int q=0;q<dim;q++){
     3800                        /* column 1-3 for the velocities */
     3801                        for (int p=0;p<dim;p++){
     3802                                for(int i=0;i<vnumnodes;i++){
     3803                                        SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
     3804                                }
     3805                        }
     3806                        /* add the diagnal components */
     3807                        for(int i=0;i<vnumnodes;i++){
     3808                                for (int p=0;p<dim;p++){
     3809                                        SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
     3810                                }
     3811                        }
     3812                        /* column 4 for the pressure */
     3813                        for(int i=0;i<vnumnodes;i++){
     3814                                        SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
     3815                        }
     3816                }
     3817
     3818                if(dim==2 || dim==3){
     3819                        /*Stress balance*/
     3820                        for(int i=0;i<vnumnodes;i++){
     3821                                for(int j=0;j<vnumnodes;j++){
     3822                                        for (int p=0;p<dim;p++){
     3823                                                for (int q=0;q<dim;q++){
     3824                                                        /* diagonal only */
     3825                                                        Ke->values[(dim*i+p)*numdof+dim*j+p] += s*viscosity*(vdbasis[q*vnumnodes+j]*vdbasis[q*vnumnodes+i]);
     3826                                                        /* All the entries */
     3827                                                        Ke->values[(dim*i+p)*numdof+dim*j+q] += s*viscosity*(vdbasis[p*vnumnodes+j]*vdbasis[q*vnumnodes+i]);
     3828                                                }
     3829                                        }
     3830                                }
     3831                                for(int k=0;k<pnumnodes;k++){
     3832                                        for (int p=0;p<dim;p++){
     3833                                                Ke->values[(dim*i+p)*numdof+dim*vnumnodes+k] += s*FSreconditioning*(-pbasis[k]*vdbasis[p*vnumnodes+i]);
     3834                                        }
     3835                                }
     3836                        }
     3837                        /*Incompressibility*/
     3838                        for(int k=0;k<pnumnodes;k++){
     3839                                for(int j=0;j<vnumnodes;j++){
     3840                                        for (int p=0;p<dim;p++){
     3841                                                Ke->values[(dim*vnumnodes+k)*numdof+dim*j+p] += s*(-FSreconditioning*vdbasis[p*vnumnodes+j]*pbasis[k]);
     3842                                        }
     3843                                }
     3844                        }
     3845
     3846                        /* GLS */
     3847                        for(int i=0;i<numdof;i++){
     3848                                for(int j=0;j<numdof;j++){
     3849                                        for (int p=0;p<dim;p++){
     3850                                                        Ke->values[i*numdof+j] += Tau * s * SU[p*numdof+i]*SU[p*numdof+j];
     3851                                        }
     3852                                }
     3853                        }
     3854                }
     3855        }
     3856
     3857        /*Transform Coordinate System*/
     3858        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     3859
     3860        /*Clean up and return*/
     3861        delete gauss;
     3862        xDelete<IssmDouble>(xyz_list);
     3863        xDelete<IssmDouble>(pbasis);
     3864        xDelete<IssmDouble>(vdbasis);
     3865        xDelete<int>(cs_list);
     3866        return Ke;
     3867}/*}}}*/
     3868ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousGLS(Element* element){/*{{{*/
     3869
     3870        int         i,dim;
     3871        IssmDouble  Jdet,forcex,forcey,forcez;
     3872        IssmDouble *xyz_list = NULL;
     3873
     3874        /*Get problem dimension*/
     3875        element->FindParam(&dim,DomainDimensionEnum);
     3876
     3877        /*Fetch number of nodes and dof for this finite element*/
     3878        int vnumnodes = element->NumberofNodesVelocity();
     3879        int pnumnodes = element->NumberofNodesPressure();
     3880
     3881        /*Prepare coordinate system list*/
     3882        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3883        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3884        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3885        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3886
     3887        /*Initialize vectors*/
     3888        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     3889        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3890        IssmDouble*    vdbasis = xNew<IssmDouble>(dim*vnumnodes);
     3891
     3892        /*Retrieve all inputs and parameters*/
     3893        element->GetVerticesCoordinates(&xyz_list);
     3894        IssmDouble  rho_ice =element->FindParam(MaterialsRhoIceEnum);
     3895        IssmDouble  gravity =element->FindParam(ConstantsGEnum);
     3896        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     3897        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     3898        Input*      loadingforcez_input=NULL;
     3899        if(dim==3){
     3900                loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     3901        }
     3902
     3903        /* prepare viscosity gradient for GLS */
     3904        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     3905        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     3906        Input* vz_input = NULL;
     3907        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     3908
     3909        IssmDouble viscosity,FSreconditioning;
     3910        IssmDouble NodalViscosity[6];
     3911        IssmDouble gradViscos[3];
     3912        IssmDouble etapq, s, Tau, mk, hk;
     3913        IssmDouble hx, hy, hz;
     3914        IssmDouble SU[3*(3+1)*6];
     3915    Gauss* vert = element->NewGauss();
     3916       
     3917        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     3918        /* Compute the nodal values of the viscosity */
     3919        for(int i=0;i<vnumnodes;i++){
     3920        vert->GaussNode(element->element_type, i);
     3921                element->material->ViscosityFS(&NodalViscosity[i],dim,xyz_list,vert,vx_input,vy_input,vz_input);
     3922        }
     3923        delete vert;
     3924       
     3925        /* Start  looping on the number of gaussian points: */
     3926        Gauss* gauss=element->NewGauss(5);
     3927        while(gauss->next()){
     3928
     3929                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3930                element->NodalFunctionsVelocity(vbasis,gauss);
     3931                element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     3932
     3933                loadingforcex_input->GetInputValue(&forcex,gauss);
     3934                loadingforcey_input->GetInputValue(&forcey,gauss);
     3935                if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
     3936
     3937                /* compute viscosity gradient at the gaussian point */
     3938                element->ValueP1DerivativesOnGauss(&gradViscos[0],NodalViscosity,xyz_list,gauss);
     3939                /*  weight*detJ */
     3940                s = gauss->weight*Jdet;
     3941                /* GLS Stabilization */
     3942                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3943                mk = 1.0 /3.0;
     3944                element->ElementSizes(&hx,&hy,&hz);
     3945                hk = max(max(hx, hy), hz);
     3946                // if(!element->IsOnDirichletBoundary()) {
     3947                        Tau = - mk * hk * hk * 0.125 / viscosity;
     3948                // }
     3949                // else {
     3950                //      Tau = 0.0;
     3951                // }
     3952
     3953                for (int q=0;q<dim;q++){
     3954                        /* column 1-3 for the velocities */
     3955                        for (int p=0;p<dim;p++){
     3956                                for(int i=0;i<vnumnodes;i++){
     3957                                        SU[q*vnumnodes*(dim+1)+i*dim+p] = (-gradViscos[p])*vdbasis[q*vnumnodes+i];
     3958                                }
     3959                        }
     3960                        /* add the diagnal components */
     3961                        for(int i=0;i<vnumnodes;i++){
     3962                                for (int p=0;p<dim;p++){
     3963                                        SU[q*vnumnodes*(dim+1)+i*dim+q] += (-gradViscos[p])*vdbasis[p*vnumnodes+i];
     3964                                }
     3965                        }
     3966                        /* column 4 for the pressure */
     3967                        for(int i=0;i<vnumnodes;i++){
     3968                                        SU[q*vnumnodes*(dim+1)+dim*vnumnodes+i] = FSreconditioning*vdbasis[q*vnumnodes+i];
     3969                        }
     3970                }
     3971
     3972                for(i=0;i<vnumnodes;i++){
     3973                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     3974                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     3975                        if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
     3976
     3977                        pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     3978                }
     3979
     3980                for(int i=0;i<vnumnodes*(dim+1);i++){
     3981                        pe->values[i] += Tau*Jdet*gauss->weight*rho_ice*(-gravity)*SU[(dim-1)*(dim+1)*vnumnodes+i];
     3982                }
     3983        }
     3984
     3985        /*Transform coordinate system*/
     3986        element->TransformLoadVectorCoord(pe,cs_list);
     3987
     3988        /*Clean up and return*/
     3989        delete gauss;
     3990        xDelete<int>(cs_list);
     3991        xDelete<IssmDouble>(vbasis);
     3992        xDelete<IssmDouble>(vdbasis);
     3993        xDelete<IssmDouble>(xyz_list);
     3994        return pe;
    37323995}/*}}}*/
    37333996
    … …  
    43354598                delete pe4;
    43364599        }
     4600        else if(fe_FS==P1P1GLSEnum) {
     4601                ElementVector* pe1=CreatePVectorFSViscousGLS(element);
     4602                ElementVector* pe2=CreatePVectorFSShelf(element);
     4603                ElementVector* pe3=CreatePVectorFSFront(element);
     4604                pe =new ElementVector(pe1,pe2,pe3);
     4605                delete pe1;
     4606                delete pe2;
     4607                delete pe3;
     4608        }
    43374609        else{
    43384610                ElementVector* pe1=CreatePVectorFSViscous(element);
    … …  
    43964668                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
    43974669                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
    4398                         if(dim==3){
    4399                                 pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
    4400                                 pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    4401                         }
    4402                         else{
    4403                                 pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    4404                         }
     4670                        if(dim==3) pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
     4671
     4672                        pe->values[i*dim+dim-1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    44054673                }
    44064674        }
    … …  
    44154683        xDelete<IssmDouble>(xyz_list);
    44164684        return pe;
     4685}/*}}}*/
     4686ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFrictionNitsche(Element* element){/*{{{*/
     4687
     4688        if(element->IsFloating() || !element->IsOnBase()) return NULL;
     4689
     4690        /*If on water or not FS, skip stiffness: */
     4691        int approximation;
     4692        element->GetInputValue(&approximation,ApproximationEnum);
     4693        if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     4694
     4695        /*Intermediaries*/
     4696        bool        mainlyfloating;
     4697        int         dim,domaintype;
     4698        int         friction_style,point1;
     4699        IssmDouble  alpha2,Jdet,fraction1,fraction2;
     4700        IssmDouble  viscosity, FSreconditioning;
     4701        IssmDouble  gllevelset,phi=1.;
     4702        IssmDouble *xyz_list_base = NULL;
     4703        IssmDouble *xyz_list = NULL;
     4704        Gauss*      gauss         = NULL;
     4705        /*coefficient of Nitsche's method*/
     4706        IssmDouble gamma, hx, hy, hz;
     4707        IssmDouble pnu, qnv;
     4708
     4709        /*Get problem dimension*/
     4710        element->FindParam(&dim,DomainDimensionEnum);
     4711        element->FindParam(&domaintype,DomainTypeEnum);
     4712        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     4713
     4714        /*Fetch number of nodes and dof for this finite element*/
     4715        int vnumnodes = element->NumberofNodesVelocity();
     4716        int pnumnodes = element->NumberofNodesPressure();
     4717        int numdof    = vnumnodes*dim + pnumnodes;
     4718
     4719        /*Initialize Element matrix and vectors*/
     4720        ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
     4721        IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
     4722        IssmDouble* vdbasis = xNew<IssmDouble>(dim*vnumnodes);
     4723        IssmDouble* pbasis  = xNew<IssmDouble>(pnumnodes);
     4724        IssmDouble* vdbasisn = xNew<IssmDouble>(vnumnodes);
     4725        IssmDouble* nsigma = xNew<IssmDouble>(numdof);
     4726        IssmDouble  bed_normal[3], bed_tangent[6];
     4727
     4728        /*Prepare coordinate system list*/
     4729        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     4730        if(dim==2) for(int i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     4731        else       for(int i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     4732        for(int i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     4733
     4734        /*Retrieve all inputs and parameters*/
     4735        element->GetVerticesCoordinatesBase(&xyz_list_base);
     4736        element->GetVerticesCoordinates(&xyz_list);
     4737        element->FindParam(&friction_style,GroundinglineFrictionInterpolationEnum);
     4738        Input* gllevelset_input = NULL;
     4739        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     4740        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     4741        Input* vz_input = NULL;
     4742        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     4743
     4744        /*build friction object, used later on: */
     4745        Friction* friction=new Friction(element,dim==3?3:1);
     4746
     4747        /*Recover portion of element that is grounded*/
     4748        if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base);
     4749        if(friction_style==SubelementFriction2Enum){
     4750                if(domaintype==Domain2DverticalEnum) _error_("Subelement Friction 2 not implemented yet for Flowline");
     4751                gllevelset_input=element->GetInput(MaskOceanLevelsetEnum); _assert_(gllevelset_input);
     4752                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     4753                gauss=element->NewGaussBase(3);
     4754        }
     4755        else{
     4756                gauss=element->NewGaussBase(3);
     4757        }
     4758
     4759        element->NormalBase(&bed_normal[0],xyz_list_base);
     4760        element->TangentBase(&bed_tangent[0], &bed_normal[0]);
     4761        element->ElementSizes(&hx,&hy,&hz);
     4762        element->FindParam(&gamma,FeFSNitscheGammaEnum);
     4763
     4764        gamma = gamma / sqrt(hx*hx+hy*hy+hz*hz);
     4765
     4766        /* Start  looping on the number of gaussian points: */
     4767        while(gauss->next()){
     4768                friction->GetAlpha2(&alpha2,gauss);
     4769                if(friction_style==SubelementFriction1Enum) alpha2=phi*alpha2;
     4770                else if(friction_style==SubelementFriction2Enum){
     4771                        gllevelset_input->GetInputValue(&gllevelset, gauss);
     4772                        if(gllevelset<0.) alpha2=0.;
     4773                }
     4774                else if(friction_style==NoFrictionOnPartiallyFloatingEnum){
     4775                        if (phi<0.99999999) alpha2=0.;
     4776                }
     4777                else  _error_("friction interpolation "<<EnumToStringx(friction_style)<<" not implemented yet");
     4778
     4779                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     4780                element->NodalFunctionsVelocity(vbasis,gauss);
     4781                /* The full element is needed for calculating derivative of the basis functions */
     4782                element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     4783                element->NodalFunctionsPressure(pbasis,gauss);
     4784                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     4785
     4786                /*Frictions*/
     4787                for(int i=0;i<vnumnodes;i++){
     4788                        for(int j=0;j<vnumnodes;j++){
     4789                                for (int d=0; d<2; d++) {
     4790                                        for (int p=0; p<dim; p++) {
     4791                                                for (int q=0; q<dim; q++) {
     4792                                                        Ke->values[(dim*i+p)*numdof+dim*j+q] += alpha2*gauss->weight*Jdet*vbasis[i]*vbasis[j]*bed_tangent[d*3+p]*bed_tangent[d*3+q];
     4793                                                }
     4794                                        }
     4795                                }
     4796                        }
     4797                }       
     4798                /* -------- Nitsche terms -------- */
     4799                /* n\cdot\nabla\phi*/
     4800                for(int i=0;i<vnumnodes;i++){
     4801                        vdbasisn[i] = 0.0;
     4802                        for (int d=0; d<dim; d++) {
     4803                                vdbasisn[i] += bed_normal[d]*vdbasis[d*vnumnodes+i];
     4804                        }
     4805                }
     4806                /* Boundary terms for integration by parts.
     4807                        The coefficient matrix of n*sigma*n-gamma*n*u is stored in the following order:
     4808                        rows--dimensions, columns--number of nodes.
     4809                        If we consider nsigma as a 1d vector, it has exactly the same order as the unknown vector.
     4810                */
     4811                for(int i=0;i<vnumnodes;i++){
     4812                        for(int d=0;d<dim;d++){
     4813                                nsigma[d*vnumnodes+i] = 2.0 * viscosity * bed_normal[d]*vdbasisn[i];
     4814                        }
     4815                }
     4816                for(int i=0;i<pnumnodes;i++){   
     4817                        nsigma[dim*vnumnodes+i] = -FSreconditioning*pbasis[i];
     4818                }
     4819
     4820                /*velocity components A11 */
     4821                for(int i=0;i<vnumnodes;i++){
     4822                        for (int p=0; p<dim; p++) {
     4823                                for(int j=0;j<vnumnodes;j++){
     4824                                        for (int q=0; q<dim; q++) {
     4825                                                /* gamma*(n*u)*(n*v) */
     4826                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += gauss->weight*Jdet * (-gamma) * bed_normal[q] * vbasis[j]*bed_normal[p] * vbasis[i];
     4827                                                /* sigma(u)*(n*v) */
     4828                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* nsigma[p*vnumnodes+i] * bed_normal[q] * vbasis[j];
     4829                                                /* sigma(v)*(n*u) */
     4830                                                Ke->values[(dim*i+p)*numdof+dim*j+q] += (- gauss->weight*Jdet)* bed_normal[p] * vbasis[i] * nsigma[q*vnumnodes+j];
     4831                                        }
     4832                                }
     4833                        }
     4834                }       
     4835
     4836                /* pressure x velocity  component A12 */
     4837                for(int k=0;k<pnumnodes;k++){
     4838                        pnu = nsigma[dim*vnumnodes+k];
     4839                        for(int i=0;i<vnumnodes;i++){
     4840                                for (int p=0;p<dim;p++) {
     4841                                        Ke->values[(dim*i+p)*numdof+dim*vnumnodes+k] += gauss->weight*Jdet * pnu * (- bed_normal[p] * vbasis[i]);
     4842                                }
     4843                        }
     4844                }
     4845                /* velocity x pressure component A21 */
     4846                for(int k=0;k<pnumnodes;k++){
     4847                        qnv = nsigma[dim*vnumnodes+k];
     4848                        for(int j=0;j<vnumnodes;j++){
     4849                                for (int q=0;q<dim;q++) {
     4850                                        Ke->values[(dim*vnumnodes+k)*numdof+dim*j+q] += gauss->weight*Jdet * ( - bed_normal[q] * vbasis[j]) * qnv;
     4851                                }
     4852                        }
     4853                }
     4854        }
     4855
     4856        /*Transform Coordinate System*/
     4857        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     4858
     4859        /*Clean up and return*/
     4860        delete gauss;
     4861        delete friction;
     4862        xDelete<IssmDouble>(xyz_list_base);
     4863        xDelete<IssmDouble>(xyz_list);
     4864        xDelete<IssmDouble>(vbasis);
     4865        xDelete<IssmDouble>(vdbasis);
     4866        xDelete<IssmDouble>(pbasis);
     4867        xDelete<IssmDouble>(vdbasisn);
     4868        xDelete<IssmDouble>(nsigma);
     4869        xDelete<int>(cs_list);
     4870        return Ke;
    44174871}/*}}}*/
    44184872#endif

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r25451 r25514  
    7373                ElementMatrix* CreateKMatrixFS(Element* element);
    7474                ElementMatrix* CreateKMatrixFSFriction(Element* element);
     75                ElementMatrix* CreateKMatrixFSFrictionNitsche(Element* element);
    7576                ElementMatrix* CreateKMatrixFSShelf(Element* element);
    7677                ElementMatrix* CreateKMatrixFSViscous(Element* element);
     78                ElementMatrix* CreateKMatrixFSViscousGLS(Element* element);
    7779                ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
    7880                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
    … …  
    8587                ElementVector* CreatePVectorFSStress(Element* element);
    8688                ElementVector* CreatePVectorFSViscous(Element* element);
     89                ElementVector* CreatePVectorFSViscousGLS(Element* element);
    8790                ElementVector* CreatePVectorFSViscousLA(Element* element);
    8891                ElementVector* CreatePVectorFSViscousXTH(Element* element);

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

    r25508 r25514  
    343343                virtual void        StressIntensityFactor(void)=0;
    344344                virtual IssmDouble SurfaceArea(void)=0;
     345                virtual void       TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal){_error_("not implemented yet");};
    345346                virtual int        TensorInterpolation()=0;
    346347                virtual IssmDouble TimeAdapt()=0;

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

    r25508 r25514  
    36793679                return S;
    36803680        }
     3681}
     3682/*}}}*/
     3683void       Penta::TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal){/*{{{*/
     3684        /*
     3685         To compute the two tangent bed_tangent[0:2] and bed_tangent[3:5] from the given normal vecotr.
     3686        */
     3687        IssmDouble n1, n2, n3;
     3688        IssmDouble tangent_norm;
     3689
     3690        n1 = fabs(bed_normal[0]);
     3691        n2 = fabs(bed_normal[1]);
     3692        n3 = fabs(bed_normal[2]);
     3693
     3694        /* For the first tangent, on x-y plane in most cases*/
     3695        if ((n1<=n3) && (n2<=n3)) {
     3696                bed_tangent[0] = 0.0;
     3697                bed_tangent[1] = -bed_normal[2];
     3698                bed_tangent[2] = bed_normal[1];
     3699        }
     3700        else {
     3701                bed_tangent[0] = -bed_normal[1];
     3702                bed_tangent[1] = bed_normal[0];
     3703                bed_tangent[2] = 0.0;
     3704        }
     3705        tangent_norm = sqrt(bed_tangent[0]*bed_tangent[0]+bed_tangent[1]*bed_tangent[1]+bed_tangent[2]*bed_tangent[2]);
     3706        for(int i=0;i<3;i++) bed_tangent[i] = bed_tangent[i]/tangent_norm;
     3707        /* The second tangent*/
     3708        bed_tangent[3] = bed_normal[1]*bed_tangent[2] - bed_normal[2]*bed_tangent[1];
     3709        bed_tangent[4] = bed_normal[2]*bed_tangent[0] - bed_normal[0]*bed_tangent[2];
     3710        bed_tangent[5] = bed_normal[0]*bed_tangent[1] - bed_normal[1]*bed_tangent[0];
     3711        tangent_norm = sqrt(bed_tangent[3]*bed_tangent[3]+bed_tangent[4]*bed_tangent[4]+bed_tangent[5]*bed_tangent[5]);
     3712        for(int i=3;i<6;i++) bed_tangent[i] = bed_tangent[i]/tangent_norm;
     3713
     3714        IssmDouble checksum = 0.0;
     3715        for (int i=0;i<3;i++) {
     3716                checksum += bed_normal[i]*bed_tangent[i];
     3717                checksum += bed_normal[i]*bed_tangent[i+3];
     3718                checksum += bed_tangent[i]*bed_tangent[i+3];
     3719        }
     3720        _assert_(fabs(checksum)<1e-10);
    36813721}
    36823722/*}}}*/

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

    r25508 r25514  
    178178                void           StrainRateperpendicular();
    179179                IssmDouble     SurfaceArea(void);
     180                void               TangentBase(IssmDouble* bed_tangent,IssmDouble* bed_normal);
    180181                int            TensorInterpolation(){_error_("not implemented yet");};
    181182                IssmDouble     TimeAdapt();

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

    r25446 r25514  
    571571        switch(finiteelement){
    572572                case P1Enum: case P1DGEnum:
     573                case P1P1GLSEnum: case P1P1Enum: /* added to allow P1-P1 GLS */
    573574                        switch(iv){
    574575                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;

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

    r25451 r25514  
    567567                        break;
    568568                case P1Enum: case P1DGEnum:
     569                case P1P1GLSEnum: case P1P1Enum:/* added to allow P1-P1 GLS */
    569570                        switch(iv){
    570571                                case 0: coord1=1.; coord2=0.; coord3=0.; break;

  • issm/trunk-jpl/src/c/cores/stressbalance_core.cpp

    r24240 r25514  
    1919        bool       dakota_analysis,control_analysis;
    2020        int        domaintype;
    21         bool       isSIA,isSSA,isL1L2,isHO,isFS;
     21        bool       isSIA,isSSA,isL1L2,isHO,isFS,isNitsche;
    2222        bool       save_results;
    2323        int        solution_type;
    … …  
    3333        femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
    3434        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
     35        femmodel->parameters->FindParam(&isNitsche,FlowequationIsNitscheEnum);
    3536        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    3637        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    4647                bedslope_core(femmodel);
    4748                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
    48                 ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     49                if (!isNitsche) {
     50                        ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     51                }
    4952
    5053                /*We need basal melt rates for the shelf dampening*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r24098 r25514  
    6666                                numnodes = iomodel->numberofvertices+iomodel->numberofedges;
    6767                                break;
     68                        case P1P1Enum: case P1P1GLSEnum:
     69                                /*P1 velocity + P1 pressure element with GLS stabilization*/
     70                                numnodes = (iomodel->numberofvertices) + iomodel->numberofvertices;
     71                                break;
    6872                        case MINIEnum: case MINIcondensedEnum:
    6973                                /*P1+ velocity (bubble statically condensed), P1 pressure*/
    … …  
    127131                                FacesPartitioning(iomodel);
    128132                                numnodes = iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces;
     133                                break;
     134                        case P1P1Enum: case P1P1GLSEnum:
     135                                /*P1 velocity + P1 pressure element with GLS stabilization*/
     136                                numnodes = (iomodel->numberofvertices) + iomodel->numberofvertices;
    129137                                break;
    130138                        case MINIEnum: case MINIcondensedEnum:
    … …  
    207215                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
    208216                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     217                                        break;
     218                                case P1P1Enum: case P1P1GLSEnum:
     219                                        element_numnodes = 3+3;
     220                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     221                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     222                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     223                                        for(int n=0;n<3;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     224                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elements[3*i+0]-1;
     225                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elements[3*i+1]-1;
     226                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elements[3*i+2]-1;
     227                                        for(int n=3;n<6;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
    209228                                        break;
    210229                                case MINIEnum: case MINIcondensedEnum:
    … …  
    408427                                        element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/
    409428                                        element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/
     429                                        break;
     430                                case P1P1Enum: case P1P1GLSEnum:
     431                                        /*P1-P1 very similar to MINI, but no DOF on the element
     432                                        TODO: add if(!approximations) or not?*/
     433                                        element_numnodes = 6+6;
     434                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     435                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     436                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     437                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     438                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     439                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     440                                        if(!approximations) for(int n=0;n<6;n ++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     441                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elements[6*i+0]-1;
     442                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elements[6*i+1]-1;
     443                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elements[6*i+2]-1;
     444                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elements[6*i+3]-1;
     445                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elements[6*i+4]-1;
     446                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elements[6*i+5]-1;
     447                                        if(!approximations) for(int n=6;n<12;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
    410448                                        break;
    411449                                case MINIEnum: case MINIcondensedEnum:

  • issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp

    r25379 r25514  
    2424                /*include switch for elements with multiple different sets of nodes*/
    2525                switch(element->GetElementType()){
     26                        case P1P1GLSEnum: case P1P1Enum:/* added to allow P1-P1 GLS */
    2627                        case MINIEnum:case MINIcondensedEnum:
    2728                        case TaylorHoodEnum:case XTaylorHoodEnum:case LATaylorHoodEnum:

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r25467 r25514  
    148148syn keyword cConstant FlowequationIsSIAEnum
    149149syn keyword cConstant FlowequationIsSSAEnum
     150syn keyword cConstant FlowequationIsNitscheEnum
     151syn keyword cConstant FeFSNitscheGammaEnum
    150152syn keyword cConstant FrictionCouplingEnum
    151153syn keyword cConstant FrictionDeltaEnum

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25467 r25514  
    142142        FlowequationIsSIAEnum,
    143143        FlowequationIsSSAEnum,
     144        FlowequationIsNitscheEnum,
     145        FeFSNitscheGammaEnum,
    144146        FrictionCouplingEnum,
    145147        FrictionDeltaEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r25467 r25514  
    150150                case FlowequationIsSIAEnum : return "FlowequationIsSIA";
    151151                case FlowequationIsSSAEnum : return "FlowequationIsSSA";
     152                case FlowequationIsNitscheEnum : return "FlowequationIsNitsche";
     153                case FeFSNitscheGammaEnum : return "FeFSNitscheGamma";
    152154                case FrictionCouplingEnum : return "FrictionCoupling";
    153155                case FrictionDeltaEnum : return "FrictionDelta";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r25467 r25514  
    153153              else if (strcmp(name,"FlowequationIsSIA")==0) return FlowequationIsSIAEnum;
    154154              else if (strcmp(name,"FlowequationIsSSA")==0) return FlowequationIsSSAEnum;
     155              else if (strcmp(name,"FlowequationIsNitsche")==0) return FlowequationIsNitscheEnum;
     156              else if (strcmp(name,"FeFSNitscheGamma")==0) return FeFSNitscheGammaEnum;
    155157              else if (strcmp(name,"FrictionCoupling")==0) return FrictionCouplingEnum;
    156158              else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
    … …  
    258260              else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
    259261              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
    260               else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    261               else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
     265              if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
     266              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
     267              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    266268              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    267269              else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
    … …  
    381383              else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum;
    382384              else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
    383               else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
    384               else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
     388              if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
     389              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
     390              else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
    389391              else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
    390392              else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
    … …  
    504506              else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
    505507              else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
    506               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    507               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
     511              if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
     512              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     513              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    512514              else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
    513515              else if (strcmp(name,"BasalStressy")==0) return BasalStressyEnum;
    … …  
    627629              else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum;
    628630              else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
    629               else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
    630               else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
     634              if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
     635              else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
     636              else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
    635637              else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum;
    636638              else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
    … …  
    750752              else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
    751753              else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
    752               else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
    753               else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
     757              if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
     758              else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
     759              else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
    758760              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    759761              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    … …  
    873875              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    874876              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    875               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    876               else if (strcmp(name,"Vy")==0) return VyEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
     880              if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     881              else if (strcmp(name,"Vy")==0) return VyEnum;
     882              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    881883              else if (strcmp(name,"VyObs")==0) return VyObsEnum;
    882884              else if (strcmp(name,"Vz")==0) return VzEnum;
    … …  
    996998              else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum;
    997999              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    998               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    999               else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
     1003              if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     1004              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
     1005              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
    10041006              else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
    10051007              else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
    … …  
    11191121              else if (strcmp(name,"Gset")==0) return GsetEnum;
    11201122              else if (strcmp(name,"Gsl")==0) return GslEnum;
    1121               else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
    1122               else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Hook")==0) return HookEnum;
     1126              if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
     1127              else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
     1128              else if (strcmp(name,"Hook")==0) return HookEnum;
    11271129              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    11281130              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    … …  
    12421244              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    12431245              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    1244               else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    1245               else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
     1249              if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
     1250              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
     1251              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    12501252              else if (strcmp(name,"P1xP4")==0) return P1xP4Enum;
    12511253              else if (strcmp(name,"P2")==0) return P2Enum;
    … …  
    13651367              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    13661368              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    1367               else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    1368               else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
     1372              if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
     1373              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
     1374              else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
    13731375              else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
    13741376              else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r24861 r25514  
    1111                isHO                           = 0;
    1212                isFS                           = 0;
     13                isNitscheBC                    = 0;
     14      FSNitscheGamma                 = 1e6;
    1315                fe_SSA                         = '';
    1416                fe_HO                          = '';
    … …  
    98100                        md = checkfield(md,'fieldname','flowequation.isHO','numel',[1],'values',[0 1]);
    99101                        md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
     102                        md = checkfield(md,'fieldname','flowequation.isNitscheBC','numel',[1],'values',[0 1]);
     103                        md = checkfield(md,'fieldname','flowequation.FSNitscheGamma','numel',[1], '>=', 0.);
    100104                        md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'});
    101105                        md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P1xP4','P2xP4'});
    … …  
    144148                        fielddisplay(self,'isHO','is the Higher-Order (HO) approximation used ?');
    145149                        fielddisplay(self,'isFS','are the Full-FS (FS) equations used ?');
     150                        fielddisplay(self,'isNitscheBC','is weakly imposed condition used?');
     151                        fielddisplay(self,'FSNitscheGamma','Gamma value for the Nitsche term, by default gamma=1e6?');
    146152                        fielddisplay(self,'fe_SSA','Finite Element for SSA  ''P1'', ''P1bubble'' ''P1bubblecondensed'' ''P2''');
    147153                        fielddisplay(self,'fe_HO', 'Finite Element for HO   ''P1'' ''P1bubble'' ''P1bubblecondensed'' ''P1xP2'' ''P2xP1'' ''P2''');
    … …  
    160166                        WriteData(fid,prefix,'object',self,'fieldname','isHO','format','Boolean');
    161167                        WriteData(fid,prefix,'object',self,'fieldname','isFS','format','Boolean');
     168                        WriteData(fid,prefix,'object',self,'fieldname','isNitscheBC','format','Boolean');
     169                        WriteData(fid,prefix,'object',self,'fieldname','FSNitscheGamma','format','Double');
    162170                        WriteData(fid,prefix,'object',self,'fieldname','fe_SSA','data',self.fe_SSA,'format','String');
    163171                        WriteData(fid,prefix,'object',self,'fieldname','fe_HO' ,'data',self.fe_HO,'format','String');
    … …  
    182190                        writejsdouble(fid,[modelname '.flowequation.isHO'],self.isHO);
    183191                        writejsdouble(fid,[modelname '.flowequation.isFS'],self.isFS);
    184                         writejsstring(fid,[modelname '.flowequation.fe_SSA'],self.fe_SSA);
     192         writejsstring(fid,[modelname '.flowequation.isNitscheBC'],self.isNitscheBC);
     193         writejsstring(fid,[modelname '.flowequation.FSNitscheGamma'],self.FSNitscheGamma);
     194         writejsstring(fid,[modelname '.flowequation.fe_SSA'],self.fe_SSA);
    185195                        writejsstring(fid,[modelname '.flowequation.fe_HO'],self.fe_HO);
    186196                        writejsstring(fid,[modelname '.flowequation.fe_FS'],self.fe_FS);

  • issm/trunk-jpl/src/m/classes/frictionschoof.m

    r24666 r25514  
    4545                        disp('Schoof sliding law parameters:');
    4646                        disp('   Schoof''s sliding law reads:');
    47                         disp('                         C |u_b|^(m-1)                ');
     47                        disp('                         C^2 |u_b|^(m-1)                ');
    4848                        disp('      tau_b = - _____________________________   u_b   ');
    49                         disp('               (1+(C/(Cmax N))^1/m |u_b| )^m          ');
     49                        disp('               (1+(C^2/(Cmax N))^1/m |u_b| )^m          ');
    5050                        disp(' ');
    5151                        fielddisplay(self,'C','friction coefficient [SI]');
Note: See TracChangeset for help on using the changeset viewer.