Ignore:
Timestamp:
07/23/13 14:23:07 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: hutter-> SIA macayeal->SSA pattyn->HO stokes->FS

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

Legend:

Unmodified
Added
Removed

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

    r15556 r15564  
    217217        IssmDouble  xyz_list[NUMVERTICES][3];
    218218        IssmDouble  xyz_list_tria[3][3];
    219         IssmDouble  rho_ice,gravity,stokesreconditioning;
     219        IssmDouble  rho_ice,gravity,FSreconditioning;
    220220        IssmDouble  pressure,viscosity,Jdet2d;
    221221        IssmDouble  bed_normal[3];
    … …  
    234234        /*Check analysis_types*/
    235235        if (analysis_type!=DiagnosticHorizAnalysisEnum) _error_("Not supported yet!");
    236         if (approximation!=StokesApproximationEnum) _error_("Not supported yet!");
     236        if (approximation!=FSApproximationEnum) _error_("Not supported yet!");
    237237
    238238        /*retrieve some parameters: */
    239         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     239        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    240240
    241241        if(!IsOnBed()){
    … …  
    267267                /*Compute strain rate viscosity and pressure: */
    268268                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    269                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     269                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    270270                pressure_input->GetInputValue(&pressure,gauss);
    271271
    272272                /*Compute Stress*/
    273                 sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
    274                 sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning;
    275                 sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning;
     273                sigma_xx=2*viscosity*epsilon[0]-pressure*FSreconditioning; // sigma = nu eps - pressure
     274                sigma_yy=2*viscosity*epsilon[1]-pressure*FSreconditioning;
     275                sigma_zz=2*viscosity*epsilon[2]-pressure*FSreconditioning;
    276276                sigma_xy=2*viscosity*epsilon[3];
    277277                sigma_xz=2*viscosity*epsilon[4];
    … …  
    415415                        Ke=CreateKMatrixAdjointHoriz();
    416416                        break;
    417                 case DiagnosticHutterAnalysisEnum:
    418                         Ke=CreateKMatrixDiagnosticHutter();
     417                case DiagnosticSIAAnalysisEnum:
     418                        Ke=CreateKMatrixDiagnosticSIA();
    419419                        break;
    420420                case DiagnosticVertAnalysisEnum:
    … …  
    546546                        pe=CreatePVectorDiagnosticHoriz();
    547547                        break;
    548                 case DiagnosticHutterAnalysisEnum:
    549                         pe=CreatePVectorDiagnosticHutter();
     548                case DiagnosticSIAAnalysisEnum:
     549                        pe=CreatePVectorDiagnosticSIA();
    550550                        break;
    551551                case DiagnosticVertAnalysisEnum:
    … …  
    11901190                int approximation;
    11911191                inputs->GetInputValue(&approximation,ApproximationEnum);
    1192                 if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    1193                         GetSolutionFromInputsDiagnosticStokes(solution);
    1194                 }
    1195                 else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
     1192                if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
     1193                        GetSolutionFromInputsDiagnosticFS(solution);
     1194                }
     1195                else if (approximation==MacAyealApproximationEnum || approximation==HOApproximationEnum || approximation==SIAApproximationEnum){
    11961196                        GetSolutionFromInputsDiagnosticHoriz(solution);
    11971197                }
    1198                 else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
     1198                else if (approximation==MacAyealHOApproximationEnum || approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
    11991199                        return; //the elements around will create the solution
    12001200                }
    12011201                break;
    1202         case DiagnosticHutterAnalysisEnum:
    1203                 GetSolutionFromInputsDiagnosticHutter(solution);
     1202        case DiagnosticSIAAnalysisEnum:
     1203                GetSolutionFromInputsDiagnosticSIA(solution);
    12041204                break;
    12051205        case DiagnosticVertAnalysisEnum:
    … …  
    12491249}
    12501250/*}}}*/
    1251 /*FUNCTION Penta::GetStrainRate3dPattyn{{{*/
    1252 void Penta::GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
    1253         /*Compute the 3d Blatter/PattynStrain Rate (5 components):
     1251/*FUNCTION Penta::GetStrainRate3dHO{{{*/
     1252void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
     1253        /*Compute the 3d Blatter/HOStrain Rate (5 components):
    12541254         *
    12551255         * epsilon=[exx eyy exy exz eyz]
    … …  
    12711271
    12721272        /*Get strain rate assuming that epsilon has been allocated*/
    1273         vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
    1274         vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     1273        vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
     1274        vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
    12751275
    12761276        /*Sum all contributions*/
    … …  
    19731973                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum));
    19741974                }
    1975                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynApproximationEnum){
    1976                         this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum));
    1977                 }
    1978                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealPattynApproximationEnum){
    1979                         this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum));
    1980                 }
    1981                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HutterApproximationEnum){
    1982                         this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum));
     1975                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOApproximationEnum){
     1976                        this->inputs->AddInput(new IntInput(ApproximationEnum,HOApproximationEnum));
     1977                }
     1978                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealHOApproximationEnum){
     1979                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealHOApproximationEnum));
     1980                }
     1981                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SIAApproximationEnum){
     1982                        this->inputs->AddInput(new IntInput(ApproximationEnum,SIAApproximationEnum));
    19831983                }
    19841984                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){
    19851985                        this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum));
    19861986                }
    1987                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==StokesApproximationEnum){
    1988                         this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));
    1989                 }
    1990                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealStokesApproximationEnum){
    1991                         this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum));
    1992                 }
    1993                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynStokesApproximationEnum){
    1994                         this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum));
     1987                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==FSApproximationEnum){
     1988                        this->inputs->AddInput(new IntInput(ApproximationEnum,FSApproximationEnum));
     1989                }
     1990                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealFSApproximationEnum){
     1991                        this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealFSApproximationEnum));
     1992                }
     1993                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOFSApproximationEnum){
     1994                        this->inputs->AddInput(new IntInput(ApproximationEnum,HOFSApproximationEnum));
    19951995                }
    19961996                else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){
    … …  
    20312031                InputUpdateFromSolutionDiagnosticHoriz( solution);
    20322032                break;
    2033         case DiagnosticHutterAnalysisEnum:
    2034                 InputUpdateFromSolutionDiagnosticHutter( solution);
     2033        case DiagnosticSIAAnalysisEnum:
     2034                InputUpdateFromSolutionDiagnosticSIA( solution);
    20352035                break;
    20362036        case DiagnosticVertAnalysisEnum:
    … …  
    20422042                int approximation;
    20432043                inputs->GetInputValue(&approximation,ApproximationEnum);
    2044                 if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    2045                         InputUpdateFromSolutionAdjointStokes( solution);
     2044                if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
     2045                        InputUpdateFromSolutionAdjointFS( solution);
    20462046                }
    20472047                else{
    … …  
    26512651}
    26522652/*}}}*/
    2653 /*FUNCTION Penta::ReduceMatrixStokes {{{*/
    2654 void Penta::ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
     2653/*FUNCTION Penta::ReduceMatrixFS {{{*/
     2654void Penta::ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
    26552655
    26562656        int    i,j;
    … …  
    26932693}
    26942694/*}}}*/
    2695 /*FUNCTION Penta::ReduceVectorStokes {{{*/
    2696 void Penta::ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
     2695/*FUNCTION Penta::ReduceVectorFS {{{*/
     2696void Penta::ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
    26972697
    26982698        int    i,j;
    … …  
    27862786        IssmDouble xz_plane[6];
    27872787
    2788         /*For Stokes only: we want the CS to be tangential to the bedrock*/
     2788        /*For FS only: we want the CS to be tangential to the bedrock*/
    27892789        inputs->GetInputValue(&approximation,ApproximationEnum);
    2790         if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return;
     2790        if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
    27912791
    27922792        /*Get slope on each node*/
    … …  
    28832883        if(analysis_type==DiagnosticHorizAnalysisEnum){
    28842884                inputs->GetInputValue(&approximation,ApproximationEnum);
    2885                 if(approximation==MacAyealPattynApproximationEnum || approximation==MacAyealStokesApproximationEnum){
     2885                if(approximation==MacAyealHOApproximationEnum || approximation==MacAyealFSApproximationEnum){
    28862886                        parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
    28872887                        o_nz += numlayers*3;
    … …  
    31023102        int        stabilization;
    31033103        bool       dakota_analysis;
    3104         bool       isstokes;
     3104        bool       isFS;
    31053105        IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
    31063106
    … …  
    31093109        iomodel->Constant(&stabilization,PrognosticStabilizationEnum);
    31103110        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    3111         iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
     3111        iomodel->Constant(&isFS,FlowequationIsFSEnum);
    31123112        iomodel->Constant(&beta,MaterialsBetaEnum);
    31133113        iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
    … …  
    31753175                                        this->inputs->AddInput(new PentaInput(QmuPressureEnum,nodeinputs,P1Enum));
    31763176                                }
    3177                                 if(isstokes){
     3177                                if(isFS){
    31783178                                        this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
    31793179                                        this->inputs->AddInput(new PentaInput(PressurePicardEnum,nodeinputs,P1Enum));
    31803180                                }
    31813181                        }
    3182                         if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynStokesApproximationEnum){
    3183                                 /*Create VzPattyn and VzStokes Enums*/
    3184                                 if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
    3185                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
    3186                                         this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
    3187                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]);
    3188                                         this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
     3182                        if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==HOFSApproximationEnum){
     3183                                /*Create VzHO and VzFS Enums*/
     3184                                if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
     3185                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
     3186                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
     3187                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
     3188                                        this->inputs->AddInput(new PentaInput(VzHOEnum,nodeinputs,P1Enum));
    31893189                                }
    31903190                                else{
    31913191                                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3192                                         this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
    3193                                         this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
     3192                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
     3193                                        this->inputs->AddInput(new PentaInput(VzHOEnum,nodeinputs,P1Enum));
    31943194                                }
    31953195                        }
    3196                         if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealStokesApproximationEnum){
    3197                                 /*Create VzMacAyeal and VzStokes Enums*/
    3198                                 if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
    3199                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
    3200                                         this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
    3201                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]);
     3196                        if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealFSApproximationEnum){
     3197                                /*Create VzMacAyeal and VzFS Enums*/
     3198                                if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
     3199                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
     3200                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
     3201                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
    32023202                                        this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
    32033203                                }
    32043204                                else{
    32053205                                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3206                                         this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
     3206                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
    32073207                                        this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
    32083208                                }
    … …  
    32893289
    32903290                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    3291                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3291                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    32923292                GetPhi(&phi, &epsilon[0], viscosity);
    32933293
    … …  
    40674067
    40684068                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    4069                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     4069                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    40704070                GetPhi(&phi, &epsilon[0], viscosity);
    40714071
    … …  
    43234323
    43244324                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    4325                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     4325                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    43264326                GetPhi(&phi, &epsilon[0], viscosity);
    43274327
    … …  
    47764776                case MacAyealApproximationEnum:
    47774777                        return CreateKMatrixAdjointMacAyeal2d();
    4778                 case PattynApproximationEnum:
    4779                         return CreateKMatrixAdjointPattyn();
    4780                 case StokesApproximationEnum:
    4781                         return CreateKMatrixAdjointStokes();
     4778                case HOApproximationEnum:
     4779                        return CreateKMatrixAdjointHO();
     4780                case FSApproximationEnum:
     4781                        return CreateKMatrixAdjointFS();
    47824782                case NoneApproximationEnum:
    47834783                        return NULL;
    … …  
    48214821}
    48224822/*}}}*/
    4823 /*FUNCTION Penta::CreateKMatrixAdjointPattyn{{{*/
    4824 ElementMatrix* Penta::CreateKMatrixAdjointPattyn(void){
     4823/*FUNCTION Penta::CreateKMatrixAdjointHO{{{*/
     4824ElementMatrix* Penta::CreateKMatrixAdjointHO(void){
    48254825
    48264826        /*Intermediaries */
    … …  
    48384838        GaussPenta *gauss=NULL;
    48394839
    4840         /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
     4840        /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
    48414841        parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    4842         ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
     4842        ElementMatrix* Ke=CreateKMatrixDiagnosticHO();
    48434843        if(incomplete_adjoint) return Ke;
    48444844
    … …  
    48574857                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    48584858
    4859                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     4859                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    48604860                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    48614861                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    … …  
    48864886}
    48874887/*}}}*/
    4888 /*FUNCTION Penta::CreateKMatrixAdjointStokes{{{*/
    4889 ElementMatrix* Penta::CreateKMatrixAdjointStokes(void){
     4888/*FUNCTION Penta::CreateKMatrixAdjointFS{{{*/
     4889ElementMatrix* Penta::CreateKMatrixAdjointFS(void){
    48904890
    48914891        /*Constants*/
    … …  
    49064906        GaussPenta *gauss=NULL;
    49074907
    4908         /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
     4908        /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
    49094909        parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
    4910         ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
     4910        ElementMatrix* Ke=CreateKMatrixDiagnosticFS();
    49114911        if(incomplete_adjoint) return Ke;
    49124912
    … …  
    49264926                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    49274927
    4928                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     4928                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    49294929                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    49304930                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    … …  
    49734973                case MacAyealApproximationEnum:
    49744974                        return CreatePVectorAdjointMacAyeal();
    4975                 case PattynApproximationEnum:
    4976                         return CreatePVectorAdjointPattyn();
     4975                case HOApproximationEnum:
     4976                        return CreatePVectorAdjointHO();
    49774977                case NoneApproximationEnum:
    49784978                        return NULL;
    4979                 case StokesApproximationEnum:
    4980                         return CreatePVectorAdjointStokes();
     4979                case FSApproximationEnum:
     4980                        return CreatePVectorAdjointFS();
    49814981                default:
    49824982                        _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    49984998}
    49994999/*}}}*/
    5000 /*FUNCTION Penta::CreatePVectorAdjointPattyn{{{*/
    5001 ElementVector* Penta::CreatePVectorAdjointPattyn(void){
     5000/*FUNCTION Penta::CreatePVectorAdjointHO{{{*/
     5001ElementVector* Penta::CreatePVectorAdjointHO(void){
    50025002
    50035003        if (!IsOnSurface()) return NULL;
    … …  
    50125012}
    50135013/*}}}*/
    5014 /*FUNCTION Penta::CreatePVectorAdjointStokes{{{*/
    5015 ElementVector* Penta::CreatePVectorAdjointStokes(void){
     5014/*FUNCTION Penta::CreatePVectorAdjointFS{{{*/
     5015ElementVector* Penta::CreatePVectorAdjointFS(void){
    50165016
    50175017        if (!IsOnSurface()) return NULL;
    … …  
    50195019        /*Call Tria function*/
    50205020        Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
    5021         ElementVector* pe=tria->CreatePVectorAdjointStokes();
     5021        ElementVector* pe=tria->CreatePVectorAdjointFS();
    50225022        delete tria->material; delete tria;
    50235023
    … …  
    50595059                                        GradjDragMacAyeal(gradient,control_index);
    50605060                                        break;
    5061                                 case PattynApproximationEnum:
    5062                                         GradjDragPattyn(gradient,control_index);
     5061                                case HOApproximationEnum:
     5062                                        GradjDragHO(gradient,control_index);
    50635063                                        break;
    5064                                 case StokesApproximationEnum:
    5065                                         GradjDragStokes(gradient,control_index);
     5064                                case FSApproximationEnum:
     5065                                        GradjDragFS(gradient,control_index);
    50665066                                        break;
    50675067                                case NoneApproximationEnum:
    … …  
    50795079                                        GradjBbarMacAyeal(gradient,control_index);
    50805080                                        break;
    5081                                 case PattynApproximationEnum:
    5082                                         GradjBbarPattyn(gradient,control_index);
     5081                                case HOApproximationEnum:
     5082                                        GradjBbarHO(gradient,control_index);
    50835083                                        break;
    5084                                 case StokesApproximationEnum:
    5085                                         GradjBbarStokes(gradient,control_index);
     5084                                case FSApproximationEnum:
     5085                                        GradjBbarFS(gradient,control_index);
    50865086                                        break;
    50875087                                case NoneApproximationEnum:
    … …  
    51435143
    51445144} /*}}}*/
    5145 /*FUNCTION Penta::GradjDragPattyn {{{*/
    5146 void  Penta::GradjDragPattyn(Vector<IssmDouble>* gradient,int control_index){
     5145/*FUNCTION Penta::GradjDragHO {{{*/
     5146void  Penta::GradjDragHO(Vector<IssmDouble>* gradient,int control_index){
    51475147
    51485148        int        i,j;
    … …  
    52145214}
    52155215/*}}}*/
    5216 /*FUNCTION Penta::GradjDragStokes {{{*/
    5217 void  Penta::GradjDragStokes(Vector<IssmDouble>* gradient,int control_index){
     5216/*FUNCTION Penta::GradjDragFS {{{*/
     5217void  Penta::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){
    52185218
    52195219        int        i,j;
    … …  
    53245324
    53255325} /*}}}*/
    5326 /*FUNCTION Penta::GradjBbarPattyn {{{*/
    5327 void  Penta::GradjBbarPattyn(Vector<IssmDouble>* gradient,int control_index){
     5326/*FUNCTION Penta::GradjBbarHO {{{*/
     5327void  Penta::GradjBbarHO(Vector<IssmDouble>* gradient,int control_index){
    53285328
    53295329        /*Gradient is computed on bed only (Bbar)*/
    … …  
    53415341        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    53425342} /*}}}*/
    5343 /*FUNCTION Penta::GradjBbarStokes {{{*/
    5344 void  Penta::GradjBbarStokes(Vector<IssmDouble>* gradient,int control_index){
     5343/*FUNCTION Penta::GradjBbarFS {{{*/
     5344void  Penta::GradjBbarFS(Vector<IssmDouble>* gradient,int control_index){
    53455345
    53465346        /*Gradient is computed on bed only (Bbar)*/
    … …  
    54035403}
    54045404/*}}}*/
    5405 /*FUNCTION Penta::InputUpdateFromSolutionAdjointStokes {{{*/
    5406 void  Penta::InputUpdateFromSolutionAdjointStokes(IssmDouble* solution){
     5405/*FUNCTION Penta::InputUpdateFromSolutionAdjointFS {{{*/
     5406void  Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){
    54075407
    54085408        const int    numdof=NDOF4*NUMVERTICES;
    … …  
    59515951
    59525952        switch(approximation){
    5953                 case StokesApproximationEnum:
    5954                         return CreateDVectorDiagnosticStokes();
     5953                case FSApproximationEnum:
     5954                        return CreateDVectorDiagnosticFS();
    59555955                default:
    5956                         return NULL; //no need for doftypes outside of stokes approximation
    5957         }
    5958 }
    5959 /*}}}*/
    5960 /*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{*/
    5961 ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
     5956                        return NULL; //no need for doftypes outside of FS approximation
     5957        }
     5958}
     5959/*}}}*/
     5960/*FUNCTION Penta::CreateDVectorDiagnosticFS{{{*/
     5961ElementVector* Penta::CreateDVectorDiagnosticFS(void){
    59625962
    59635963        /*output: */
    … …  
    59695969        /*Initialize Element vector and return if necessary*/
    59705970        inputs->GetInputValue(&approximation,ApproximationEnum);
    5971         if(approximation!=StokesApproximationEnum) return NULL;
    5972 
    5973         De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     5971        if(approximation!=FSApproximationEnum) return NULL;
     5972
     5973        De=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    59745974
    59755975        for (i=0;i<NUMVERTICES;i++){
    … …  
    59835983}
    59845984/*}}}*/
    5985 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{*/
    5986 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattyn(void){
     5985/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHO{{{*/
     5986ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHO(void){
    59875987
    59885988        /*compute all stiffness matrices for this element*/
    5989         ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous();
    5990         ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction();
     5989        ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealHOViscous();
     5990        ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealHOFriction();
    59915991        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    59925992
    … …  
    59975997}
    59985998/*}}}*/
    5999 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynViscous{{{*/
    6000 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){
     5999/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHOViscous{{{*/
     6000ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHOViscous(void){
    60016001
    60026002        /*Constants*/
    … …  
    60236023        int         cs_list[numnodes];
    60246024
    6025         /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */
     6025        /*Find penta on bed as HO must be coupled to the dofs on the bed: */
    60266026        Penta* pentabase=GetBasalElement();
    60276027        Tria*  tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    … …  
    60376037        /*Initialize Element matrix*/
    60386038        ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    6039         ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,PattynApproximationEnum);
     6039        ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,HOApproximationEnum);
    60406040        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    60416041        delete Ke1; delete Ke2;
    … …  
    60586058
    60596059                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6060                 GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
     6060                GetBMacAyealHO(&B[0][0], &xyz_list[0][0], gauss);
    60616061                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    60626062
    6063                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6064                 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     6063                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     6064                this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    60656065                material->GetViscosity3d(&viscosity, &epsilon[0]);
    60666066                material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    60906090}
    60916091/*}}}*/
    6092 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynFriction{{{*/
    6093 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){
     6092/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHOFriction{{{*/
     6093ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHOFriction(void){
    60946094
    60956095        /*Constants*/
    … …  
    61166116        if(IsFloating() || !IsOnBed()) return NULL;
    61176117        ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    6118         ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     6118        ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    61196119        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    61206120        delete Ke1; delete Ke2;
    … …  
    61496149
    61506150                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    6151                 GetBPattynFriction(&L[0][0],gauss);
     6151                GetBHOFriction(&L[0][0],gauss);
    61526152
    61536153                DL_scalar=alpha2*gauss->weight*Jdet2d;
    … …  
    61756175}
    61766176/*}}}*/
    6177 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokes{{{*/
    6178 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokes(void){
     6177/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFS{{{*/
     6178ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFS(void){
    61796179
    61806180        /*compute all stiffness matrices for this element*/
    6181         ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealStokesViscous();
    6182         ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealStokesFriction();
     6181        ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealFSViscous();
     6182        ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealFSFriction();
    61836183        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    61846184
    … …  
    61896189}
    61906190/*}}}*/
    6191 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesViscous{{{*/
    6192 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesViscous(void){
     6191/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFSViscous{{{*/
     6192ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFSViscous(void){
    61936193
    61946194        /*Constants*/
    … …  
    62016201        int         i,j;
    62026202        IssmDouble Jdet;
    6203         IssmDouble viscosity,stokesreconditioning; //viscosity
     6203        IssmDouble viscosity,FSreconditioning; //viscosity
    62046204        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    62056205        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    62206220        int         cs_list[numnodes];
    62216221
    6222         /*Find penta on bed as stokes must be coupled to the dofs on the bed: */
     6222        /*Find penta on bed as FS must be coupled to the dofs on the bed: */
    62236223        Penta* pentabase=GetBasalElement();
    62246224        Tria* tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    … …  
    62346234        /*Initialize Element matrix and return if necessary*/
    62356235        ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    6236         ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6236        ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,FSApproximationEnum);
    62376237        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    62386238        delete Ke1; delete Ke2;
    … …  
    62406240        /* Get node coordinates and dof list: */
    62416241        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6242         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     6242        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    62436243        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    62446244        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    … …  
    62546254
    62556255                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6256                 GetBMacAyealStokes(&B[0][0], &xyz_list[0][0], gauss);
    6257                 tria->GetBprimeMacAyealStokes(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    6258                 tria->GetBMacAyealStokes(&B2[0][0], &xyz_list[0][0], gauss_tria);
    6259                 GetBprimeMacAyealStokes(&Bprime2[0][0], &xyz_list[0][0], gauss);
     6256                GetBMacAyealFS(&B[0][0], &xyz_list[0][0], gauss);
     6257                tria->GetBprimeMacAyealFS(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
     6258                tria->GetBMacAyealFS(&B2[0][0], &xyz_list[0][0], gauss_tria);
     6259                GetBprimeMacAyealFS(&Bprime2[0][0], &xyz_list[0][0], gauss);
    62606260
    62616261                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6262                 material->GetViscosity3dStokes(&viscosity, &epsilon[0]);
     6262                material->GetViscosity3dFS(&viscosity, &epsilon[0]);
    62636263
    62646264                D_scalar=2*viscosity*gauss->weight*Jdet;
    62656265                for (i=0;i<3;i++) D[i][i]=D_scalar;
    6266                 D[3][3]=-gauss->weight*Jdet*stokesreconditioning;
     6266                D[3][3]=-gauss->weight*Jdet*FSreconditioning;
    62676267                for (i=0;i<3;i++) D2[i][i]=D_scalar;
    62686268
    … …  
    62936293}
    62946294/*}}}*/
    6295 /*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction {{{*/
    6296 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesFriction(void){
     6295/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFSFriction {{{*/
     6296ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFSFriction(void){
    62976297
    62986298        /*Constants*/
    … …  
    63076307        int        i,j;
    63086308        int        analysis_type,approximation;
    6309         IssmDouble stokesreconditioning;
     6309        IssmDouble FSreconditioning;
    63106310        IssmDouble viscosity,alpha2_gauss,Jdet2d;
    63116311        IssmDouble bed_normal[3];
    … …  
    63136313        IssmDouble xyz_list[NUMVERTICES][3];
    63146314        IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    6315         IssmDouble LMacAyealStokes[8][numdof2dm];
    6316         IssmDouble LprimeMacAyealStokes[8][numdof2d];
    6317         IssmDouble DLMacAyealStokes[8][8]={0.0};
    6318         IssmDouble LStokesMacAyeal[4][numdof2d];
    6319         IssmDouble LprimeStokesMacAyeal[4][numdof2dm];
    6320         IssmDouble DLStokesMacAyeal[4][4]={0.0};
     6315        IssmDouble LMacAyealFS[8][numdof2dm];
     6316        IssmDouble LprimeMacAyealFS[8][numdof2d];
     6317        IssmDouble DLMacAyealFS[8][8]={0.0};
     6318        IssmDouble LFSMacAyeal[4][numdof2d];
     6319        IssmDouble LprimeFSMacAyeal[4][numdof2dm];
     6320        IssmDouble DLFSMacAyeal[4][4]={0.0};
    63216321        IssmDouble Ke_drag_gaussian[numdof2dm][numdof2d];
    63226322        IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm];
    … …  
    63266326        int         cs_list[numnodes];
    63276327
    6328         /*If on water or not Stokes, skip stiffness: */
     6328        /*If on water or not FS, skip stiffness: */
    63296329        inputs->GetInputValue(&approximation,ApproximationEnum);
    63306330        if(IsFloating() || !IsOnBed()) return NULL;
    63316331        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
    6332         ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6332        ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    63336333        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    63346334        delete Ke1; delete Ke2;
    … …  
    63456345        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    63466346        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    6347         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     6347        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    63486348        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    63496349        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    63616361
    63626362                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    6363                 GetLMacAyealStokes(&LMacAyealStokes[0][0], gauss);
    6364                 GetLprimeMacAyealStokes(&LprimeMacAyealStokes[0][0], &xyz_list[0][0], gauss);
    6365                 GetLStokesMacAyeal(&LStokesMacAyeal[0][0], gauss);
    6366                 GetLprimeStokesMacAyeal(&LprimeStokesMacAyeal[0][0], &xyz_list[0][0], gauss);
     6363                GetLMacAyealFS(&LMacAyealFS[0][0], gauss);
     6364                GetLprimeMacAyealFS(&LprimeMacAyealFS[0][0], &xyz_list[0][0], gauss);
     6365                GetLFSMacAyeal(&LFSMacAyeal[0][0], gauss);
     6366                GetLprimeFSMacAyeal(&LprimeFSMacAyeal[0][0], &xyz_list[0][0], gauss);
    63676367
    63686368                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6369                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     6369                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    63706370
    63716371                BedNormal(&bed_normal[0],xyz_list_tria);
    63726372                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    63736373
    6374                 DLMacAyealStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    6375                 DLMacAyealStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
    6376                 DLMacAyealStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    6377                 DLMacAyealStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    6378                 DLMacAyealStokes[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
    6379                 DLMacAyealStokes[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
    6380                 DLMacAyealStokes[6][6]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
    6381                 DLMacAyealStokes[7][7]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
    6382 
    6383                 DLStokesMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
    6384                 DLStokesMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
    6385                 DLStokesMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
    6386                 DLStokesMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
    6387 
    6388                 TripleMultiply( &LMacAyealStokes[0][0],8,numdof2dm,1,
    6389                                         &DLMacAyealStokes[0][0],8,8,0,
    6390                                         &LprimeMacAyealStokes[0][0],8,numdof2d,0,
     6374                DLMacAyealFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     6375                DLMacAyealFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     6376                DLMacAyealFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     6377                DLMacAyealFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     6378                DLMacAyealFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
     6379                DLMacAyealFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
     6380                DLMacAyealFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0];
     6381                DLMacAyealFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1];
     6382
     6383                DLFSMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
     6384                DLFSMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
     6385                DLFSMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
     6386                DLFSMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
     6387
     6388                TripleMultiply( &LMacAyealFS[0][0],8,numdof2dm,1,
     6389                                        &DLMacAyealFS[0][0],8,8,0,
     6390                                        &LprimeMacAyealFS[0][0],8,numdof2d,0,
    63916391                                        &Ke_drag_gaussian[0][0],0);
    63926392
    6393                 TripleMultiply( &LStokesMacAyeal[0][0],4,numdof2d,1,
    6394                                         &DLStokesMacAyeal[0][0],4,4,0,
    6395                                         &LprimeStokesMacAyeal[0][0],4,numdof2dm,0,
     6393                TripleMultiply( &LFSMacAyeal[0][0],4,numdof2d,1,
     6394                                        &DLFSMacAyeal[0][0],4,4,0,
     6395                                        &LprimeFSMacAyeal[0][0],4,numdof2dm,0,
    63966396                                        &Ke_drag_gaussian2[0][0],0);
    63976397
    … …  
    64096409}
    64106410/*}}}*/
    6411 /*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{*/
    6412 ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){
     6411/*FUNCTION Penta::CreateKMatrixCouplingHOFS{{{*/
     6412ElementMatrix* Penta::CreateKMatrixCouplingHOFS(void){
    64136413
    64146414        /*Constants*/
    … …  
    64326432
    64336433        /*compute all stiffness matrices for this element*/
    6434         ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
    6435         ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     6434        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
     6435        ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    64366436        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    64376437        delete Ke1;
    64386438        delete Ke2;
    6439         Ke1=CreateKMatrixDiagnosticPattyn(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
    6440         Ke2=CreateKMatrixDiagnosticStokes(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum);
     6439        Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
     6440        Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum);
    64416441
    64426442        for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
    … …  
    64696469                case L1L2ApproximationEnum:
    64706470                        return CreateKMatrixDiagnosticL1L2();
    6471                 case PattynApproximationEnum:
    6472                         return CreateKMatrixDiagnosticPattyn();
    6473                 case StokesApproximationEnum:
    6474                         return CreateKMatrixDiagnosticStokes();
    6475                 case HutterApproximationEnum:
     6471                case HOApproximationEnum:
     6472                        return CreateKMatrixDiagnosticHO();
     6473                case FSApproximationEnum:
     6474                        return CreateKMatrixDiagnosticFS();
     6475                case SIAApproximationEnum:
    64766476                        return NULL;
    64776477                case NoneApproximationEnum:
    64786478                        return NULL;
    6479                 case MacAyealPattynApproximationEnum:
    6480                         return CreateKMatrixDiagnosticMacAyealPattyn();
    6481                 case MacAyealStokesApproximationEnum:
    6482                         return CreateKMatrixDiagnosticMacAyealStokes();
    6483                 case PattynStokesApproximationEnum:
    6484                         return CreateKMatrixDiagnosticPattynStokes();
     6479                case MacAyealHOApproximationEnum:
     6480                        return CreateKMatrixDiagnosticMacAyealHO();
     6481                case MacAyealFSApproximationEnum:
     6482                        return CreateKMatrixDiagnosticMacAyealFS();
     6483                case HOFSApproximationEnum:
     6484                        return CreateKMatrixDiagnosticHOFS();
    64856485                default:
    64866486                        _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    64886488}
    64896489/*}}}*/
    6490 /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{*/
    6491 ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
     6490/*FUNCTION Penta::CreateKMatrixDiagnosticSIA{{{*/
     6491ElementMatrix* Penta::CreateKMatrixDiagnosticSIA(void){
    64926492
    64936493        /*Constants*/
    … …  
    66106610        IssmDouble  viscosity , oldviscosity, newviscosity, viscosity_overshoot;
    66116611        IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
    6612         IssmDouble  epsilons[6];                    //6 for stokes
     6612        IssmDouble  epsilons[6];                    //6 for FS
    66136613        IssmDouble  xyz_list[NUMVERTICES][3];
    66146614        IssmDouble  B[3][numdof2d];
    … …  
    66226622        GaussTria  *gauss_tria = NULL;
    66236623
    6624         /*Find penta on bed as this is a macayeal elements: */
     6624        /*Find penta on bed as this is a SSA elements: */
    66256625        pentabase=GetBasalElement();
    66266626        tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    … …  
    66516651                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    66526652
    6653                 if(approximation==MacAyealPattynApproximationEnum){
    6654                         this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6655                         this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     6653                if(approximation==MacAyealHOApproximationEnum){
     6654                        this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     6655                        this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    66566656                        material->GetViscosity3d(&viscosity, &epsilon[0]);
    66576657                        material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    66596659                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    66606660                }
    6661                 else if (approximation==MacAyealStokesApproximationEnum){
     6661                else if (approximation==MacAyealFSApproximationEnum){
    66626662                        this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    6663                         material->GetViscosity3dStokes(&newviscosity,&epsilons[0]);
     6663                        material->GetViscosity3dFS(&newviscosity,&epsilons[0]);
    66646664                }
    66656665                else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
    … …  
    67046704}
    67056705/*}}}*/
    6706 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealPattyn{{{*/
    6707 ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealPattyn(void){
     6706/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealHO{{{*/
     6707ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealHO(void){
    67086708
    67096709        /*compute all stiffness matrices for this element*/
    67106710        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3d();
    6711         ElementMatrix* Ke2=CreateKMatrixDiagnosticPattyn();
    6712         ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealPattyn();
     6711        ElementMatrix* Ke2=CreateKMatrixDiagnosticHO();
     6712        ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealHO();
    67136713        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    67146714
    … …  
    67206720}
    67216721/*}}}*/
    6722 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealStokes{{{*/
    6723 ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealStokes(void){
     6722/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealFS{{{*/
     6723ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealFS(void){
    67246724
    67256725        /*compute all stiffness matrices for this element*/
    67266726        ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3d();
    6727         ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes();
    6728         ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealStokes();
     6727        ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
     6728        ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealFS();
    67296729        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    67306730
    … …  
    67706770        GaussTria  *gauss_tria = NULL;
    67716771
    6772         /*Find penta on bed as this is a macayeal elements: */
     6772        /*Find penta on bed as this is a SSA elements: */
    67736773        pentabase=GetBasalElement();
    67746774        tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    … …  
    68356835}
    68366836/*}}}*/
    6837 /*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{*/
    6838 ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){
     6837/*FUNCTION Penta::CreateKMatrixDiagnosticHO{{{*/
     6838ElementMatrix* Penta::CreateKMatrixDiagnosticHO(void){
    68396839
    68406840        /*compute all stiffness matrices for this element*/
    6841         ElementMatrix* Ke1=CreateKMatrixDiagnosticPattynViscous();
    6842         ElementMatrix* Ke2=CreateKMatrixDiagnosticPattynFriction();
     6841        ElementMatrix* Ke1=CreateKMatrixDiagnosticHOViscous();
     6842        ElementMatrix* Ke2=CreateKMatrixDiagnosticHOFriction();
    68436843        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
    68446844
    … …  
    68506850}
    68516851/*}}}*/
    6852 /*FUNCTION Penta::CreateKMatrixDiagnosticPattynViscous{{{*/
    6853 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){
     6852/*FUNCTION Penta::CreateKMatrixDiagnosticHOViscous{{{*/
     6853ElementMatrix* Penta::CreateKMatrixDiagnosticHOViscous(void){
    68546854
    68556855        /*Constants*/
    … …  
    68716871
    68726872        /*Initialize Element matrix*/
    6873         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     6873        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    68746874
    68756875        /*Retrieve all inputs and parameters*/
    … …  
    68896889
    68906890                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6891                 GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
    6892                 GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
    6893 
    6894                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6895                 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
     6891                GetBHO(&B[0][0], &xyz_list[0][0], gauss);
     6892                GetBprimeHO(&Bprime[0][0], &xyz_list[0][0], gauss);
     6893
     6894                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     6895                this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    68966896                material->GetViscosity3d(&viscosity, &epsilon[0]);
    68976897                material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    … …  
    69156915}
    69166916/*}}}*/
    6917 /*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{*/
    6918 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
     6917/*FUNCTION Penta::CreateKMatrixDiagnosticHOFriction{{{*/
     6918ElementMatrix* Penta::CreateKMatrixDiagnosticHOFriction(void){
    69196919
    69206920        /*Constants*/
    … …  
    69376937        if(IsFloating() || !IsOnBed()) return NULL;
    69386938
    6939         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     6939        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    69406940
    69416941        /*Retrieve all inputs and parameters*/
    … …  
    69626962
    69636963                GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
    6964                 GetBPattynFriction(&L[0][0],gauss);
     6964                GetBHOFriction(&L[0][0],gauss);
    69656965
    69666966                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    … …  
    69856985}
    69866986/*}}}*/
    6987 /*FUNCTION Penta::CreateKMatrixDiagnosticPattynStokes{{{*/
    6988 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynStokes(void){
     6987/*FUNCTION Penta::CreateKMatrixDiagnosticHOFS{{{*/
     6988ElementMatrix* Penta::CreateKMatrixDiagnosticHOFS(void){
    69896989
    69906990        /*compute all stiffness matrices for this element*/
    6991         ElementMatrix* Ke1=CreateKMatrixDiagnosticPattyn();
    6992         ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes();
    6993         ElementMatrix* Ke3=CreateKMatrixCouplingPattynStokes();
     6991        ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
     6992        ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
     6993        ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
    69946994        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    69956995
    … …  
    70017001}
    70027002/*}}}*/
    7003 /*FUNCTION Penta::CreateKMatrixDiagnosticStokes{{{*/
    7004 ElementMatrix* Penta::CreateKMatrixDiagnosticStokes(void){
    7005 
    7006         int fe_stokes;
     7003/*FUNCTION Penta::CreateKMatrixDiagnosticFS{{{*/
     7004ElementMatrix* Penta::CreateKMatrixDiagnosticFS(void){
     7005
     7006        int fe_FS;
    70077007        ElementMatrix* Ke1;
    70087008        ElementMatrix* Ke2;
    70097009        ElementMatrix* Ke;
    7010         parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum);
    7011 
    7012         switch(fe_stokes){
     7010        parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
     7011
     7012        switch(fe_FS){
    70137013                case 0:
    70147014                        /*compute all stiffness matrices for this element*/
    7015                         Ke1=CreateKMatrixDiagnosticStokesViscous();
    7016                         Ke2=CreateKMatrixDiagnosticStokesFriction();
     7015                        Ke1=CreateKMatrixDiagnosticFSViscous();
     7016                        Ke2=CreateKMatrixDiagnosticFSFriction();
    70177017                        Ke =new ElementMatrix(Ke1,Ke2);
    70187018                        break;
    70197019                case 1:
    70207020                        /*compute all stiffness matrices for this element*/
    7021                         Ke1=CreateKMatrixDiagnosticStokesGLSViscous();
    7022                         Ke2=CreateKMatrixDiagnosticStokesFriction();
     7021                        Ke1=CreateKMatrixDiagnosticFSGLSViscous();
     7022                        Ke2=CreateKMatrixDiagnosticFSFriction();
    70237023                        Ke =new ElementMatrix(Ke1,Ke2);
    70247024                        break;
    70257025                default:
    7026                         _error_("Finite element" << fe_stokes << " not supported yet");
     7026                        _error_("Finite element" << fe_FS << " not supported yet");
    70277027        }
    70287028
    … …  
    70347034}
    70357035/*}}}*/
    7036 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesViscous {{{*/
    7037 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){
     7036/*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/
     7037ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){
    70387038
    70397039        /*Intermediaries */
    70407040        int        i,approximation;
    7041         IssmDouble Jdet,viscosity,stokesreconditioning;
     7041        IssmDouble Jdet,viscosity,FSreconditioning;
    70427042        IssmDouble xyz_list[NUMVERTICES][3];
    70437043        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    70497049        GaussPenta *gauss=NULL;
    70507050
    7051         /*If on water or not Stokes, skip stiffness: */
     7051        /*If on water or not FS, skip stiffness: */
    70527052        inputs->GetInputValue(&approximation,ApproximationEnum);
    7053         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    7054         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7053        if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     7054        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    70557055
    70567056        /*Retrieve all inputs and parameters*/
    70577057        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7058         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7058        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    70597059        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    70607060        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    70687068
    70697069                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7070                 GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
    7071                 GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss);
     7070                GetBFS(&B[0][0],&xyz_list[0][0],gauss);
     7071                GetBprimeFS(&B_prime[0][0],&xyz_list[0][0],gauss);
    70727072
    70737073                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7074                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7074                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    70757075
    70767076                D_scalar=gauss->weight*Jdet;
    70777077                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    7078                 for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     7078                for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    70797079
    70807080                TripleMultiply( &B[0][0],8,27,1,
    … …  
    70857085
    70867086        /*Condensation*/
    7087         ReduceMatrixStokes(Ke->values, &Ke_temp[0][0]);
     7087        ReduceMatrixFS(Ke->values, &Ke_temp[0][0]);
    70887088        //Ke->Echo();
    70897089        //_error_("stop");
    … …  
    70977097}
    70987098/*}}}*/
    7099 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSViscous {{{*/
    7100 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){
     7099/*FUNCTION Penta::CreateKMatrixDiagnosticFSGLSViscous {{{*/
     7100ElementMatrix* Penta::CreateKMatrixDiagnosticFSGLSViscous(void){
    71017101
    71027102        int        numdof  = NUMVERTICES*NDOF4;
    … …  
    71047104        /*Intermediaries */
    71057105        int        i,j,approximation;
    7106         IssmDouble Jdet,viscosity,stokesreconditioning,diameter,rigidity;
     7106        IssmDouble Jdet,viscosity,FSreconditioning,diameter,rigidity;
    71077107        IssmDouble xyz_list[NUMVERTICES][3];
    71087108        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    71287128        int c=3; //index of pressure
    71297129
    7130         /*If on water or not Stokes, skip stiffness: */
     7130        /*If on water or not FS, skip stiffness: */
    71317131        inputs->GetInputValue(&approximation,ApproximationEnum);
    7132         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    7133         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7132        if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     7133        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    71347134
    71357135        /*Retrieve all inputs and parameters*/
    71367136        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7137         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7137        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    71387138        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    71397139        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    71657165
    71667166                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7167                 GetBStokesGLS(&B[0][0],&xyz_list[0][0],gauss);
    7168                 GetBprimeStokesGLS(&B_prime[0][0],&xyz_list[0][0],gauss);
     7167                GetBFSGLS(&B[0][0],&xyz_list[0][0],gauss);
     7168                GetBprimeFSGLS(&B_prime[0][0],&xyz_list[0][0],gauss);
    71697169
    71707170                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7171                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7171                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    71727172
    71737173                D_scalar=gauss->weight*Jdet;
    71747174                for (i=0;i<6;i++) D[i][i]=D_scalar*2.*2.*viscosity;
    7175                 for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     7175                for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    71767176
    71777177                TripleMultiply( &B[0][0],8,24,1,
    … …  
    72357235}
    72367236/*}}}*/
    7237 /*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/
    7238 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
     7237/*FUNCTION Penta::CreateKMatrixDiagnosticFSFriction{{{*/
     7238ElementMatrix* Penta::CreateKMatrixDiagnosticFSFriction(void){
    72397239
    72407240        /*Constants*/
    … …  
    72467246        int        analysis_type,approximation;
    72477247        IssmDouble alpha2,Jdet2d;
    7248         IssmDouble stokesreconditioning,viscosity;
     7248        IssmDouble FSreconditioning,viscosity;
    72497249        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    72507250        IssmDouble xyz_list[NUMVERTICES][3];
    72517251        IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    7252         IssmDouble LStokes[2][numdof2d];
    7253         IssmDouble DLStokes[2][2]={0.0};
     7252        IssmDouble LFS[2][numdof2d];
     7253        IssmDouble DLFS[2][2]={0.0};
    72547254        IssmDouble Ke_drag_gaussian[numdof2d][numdof2d];
    72557255        Friction*  friction=NULL;
    72567256        GaussPenta *gauss=NULL;
    72577257
    7258         /*If on water or not Stokes, skip stiffness: */
     7258        /*If on water or not FS, skip stiffness: */
    72597259        inputs->GetInputValue(&approximation,ApproximationEnum);
    7260         if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
    7261         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7260        if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return NULL;
     7261        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    72627262
    72637263        /*Retrieve all inputs and parameters*/
    72647264        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    72657265        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    7266         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7266        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    72677267        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    72687268        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    72807280
    72817281                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    7282                 GetLStokes(&LStokes[0][0], gauss);
     7282                GetLFS(&LFS[0][0], gauss);
    72837283
    72847284                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7285                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7285                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    72867286
    72877287                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
    72887288
    7289                 DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
    7290                 DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
    7291 
    7292                 TripleMultiply( &LStokes[0][0],2,numdof2d,1,
    7293                                         &DLStokes[0][0],2,2,0,
    7294                                         &LStokes[0][0],2,numdof2d,0,
     7289                DLFS[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
     7290                DLFS[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
     7291
     7292                TripleMultiply( &LFS[0][0],2,numdof2d,1,
     7293                                        &DLFS[0][0],2,2,0,
     7294                                        &LFS[0][0],2,numdof2d,0,
    72957295                                        &Ke_drag_gaussian[0][0],0);
    72967296
    … …  
    74127412}
    74137413/*}}}*/
    7414 /*FUNCTION Penta::CreatePVectorCouplingMacAyealStokes {{{*/
    7415 ElementVector* Penta::CreatePVectorCouplingMacAyealStokes(void){
     7414/*FUNCTION Penta::CreatePVectorCouplingMacAyealFS {{{*/
     7415ElementVector* Penta::CreatePVectorCouplingMacAyealFS(void){
    74167416
    74177417        /*compute all load vectors for this element*/
    7418         ElementVector* pe1=CreatePVectorCouplingMacAyealStokesViscous();
    7419         ElementVector* pe2=CreatePVectorCouplingMacAyealStokesFriction();
     7418        ElementVector* pe1=CreatePVectorCouplingMacAyealFSViscous();
     7419        ElementVector* pe2=CreatePVectorCouplingMacAyealFSFriction();
    74207420        ElementVector* pe =new ElementVector(pe1,pe2);
    74217421
    … …  
    74267426}
    74277427/*}}}*/
    7428 /*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesViscous {{{*/
    7429 ElementVector* Penta::CreatePVectorCouplingMacAyealStokesViscous(void){
     7428/*FUNCTION Penta::CreatePVectorCouplingMacAyealFSViscous {{{*/
     7429ElementVector* Penta::CreatePVectorCouplingMacAyealFSViscous(void){
    74307430
    74317431        /*Constants*/
    … …  
    74367436        int         approximation;
    74377437        IssmDouble  viscosity,Jdet;
    7438         IssmDouble  stokesreconditioning;
     7438        IssmDouble  FSreconditioning;
    74397439        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    74407440        IssmDouble  dw[3];
    … …  
    74467446        /*Initialize Element vector and return if necessary*/
    74477447        inputs->GetInputValue(&approximation,ApproximationEnum);
    7448         if(approximation!=MacAyealStokesApproximationEnum) return NULL;
    7449         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7448        if(approximation!=MacAyealFSApproximationEnum) return NULL;
     7449        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    74507450
    74517451        /*Retrieve all inputs and parameters*/
    74527452        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7453         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7453        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    74547454        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    74557455        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    74567456        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    7457         Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzmacayeal_input);
     7457        Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzSSA_input);
    74587458
    74597459        /* Start  looping on the number of gaussian points: */
    … …  
    74677467                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    74687468
    7469                 vzmacayeal_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     7469                vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    74707470
    74717471                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7472                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7472                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    74737473
    74747474                for(i=0;i<NUMVERTICES;i++){
    … …  
    74767476                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    74777477                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    7478                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
     7478                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    74797479                }
    74807480        }
    … …  
    74887488}
    74897489/*}}}*/
    7490 /*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesFriction{{{*/
    7491 ElementVector* Penta::CreatePVectorCouplingMacAyealStokesFriction(void){
     7490/*FUNCTION Penta::CreatePVectorCouplingMacAyealFSFriction{{{*/
     7491ElementVector* Penta::CreatePVectorCouplingMacAyealFSFriction(void){
    74927492
    74937493        /*Constants*/
    … …  
    74987498        int         approximation,analysis_type;
    74997499        IssmDouble  Jdet,Jdet2d;
    7500         IssmDouble  stokesreconditioning;
     7500        IssmDouble  FSreconditioning;
    75017501        IssmDouble      bed_normal[3];
    75027502        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    75137513        if(!IsOnBed() || IsFloating()) return NULL;
    75147514        inputs->GetInputValue(&approximation,ApproximationEnum);
    7515         if(approximation!=MacAyealStokesApproximationEnum) return NULL;
    7516         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7515        if(approximation!=MacAyealFSApproximationEnum) return NULL;
     7516        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    75177517
    75187518        /*Retrieve all inputs and parameters*/
    75197519        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    75207520        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    7521         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7521        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    75227522        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    75237523        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    75247524        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    7525         Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzmacayeal_input);
     7525        Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzSSA_input);
    75267526
    75277527        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    … …  
    75397539                GetNodalFunctionsP1(basis, gauss);
    75407540
    7541                 vzmacayeal_input->GetInputValue(&w, gauss);
    7542                 vzmacayeal_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     7541                vzSSA_input->GetInputValue(&w, gauss);
     7542                vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    75437543
    75447544                BedNormal(&bed_normal[0],xyz_list_tria);
    75457545                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7546                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7546                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    75477547                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    75487548
    … …  
    75637563}
    75647564/*}}}*/
    7565 /*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{*/
    7566 ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
     7565/*FUNCTION Penta::CreatePVectorCouplingHOFS {{{*/
     7566ElementVector* Penta::CreatePVectorCouplingHOFS(void){
    75677567
    75687568        /*compute all load vectors for this element*/
    7569         ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
    7570         ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
     7569        ElementVector* pe1=CreatePVectorCouplingHOFSViscous();
     7570        ElementVector* pe2=CreatePVectorCouplingHOFSFriction();
    75717571        ElementVector* pe =new ElementVector(pe1,pe2);
    75727572
    … …  
    75777577}
    75787578/*}}}*/
    7579 /*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{*/
    7580 ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
     7579/*FUNCTION Penta::CreatePVectorCouplingHOFSViscous {{{*/
     7580ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){
    75817581
    75827582        /*Constants*/
    … …  
    75877587        int         approximation;
    75887588        IssmDouble  viscosity,Jdet;
    7589         IssmDouble  stokesreconditioning;
     7589        IssmDouble  FSreconditioning;
    75907590        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    75917591        IssmDouble  dw[3];
    … …  
    75977597        /*Initialize Element vector and return if necessary*/
    75987598        inputs->GetInputValue(&approximation,ApproximationEnum);
    7599         if(approximation!=PattynStokesApproximationEnum) return NULL;
    7600         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7599        if(approximation!=HOFSApproximationEnum) return NULL;
     7600        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    76017601
    76027602        /*Retrieve all inputs and parameters*/
    76037603        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7604         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7604        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    76057605        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    76067606        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    76077607        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    7608         Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
     7608        Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
    76097609
    76107610        /* Start  looping on the number of gaussian points: */
    … …  
    76187618                GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    76197619
    7620                 vzpattyn_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     7620                vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    76217621
    76227622                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7623                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7623                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    76247624
    76257625                for(i=0;i<NUMVERTICES;i++){
    … …  
    76277627                        pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    76287628                        pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    7629                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
     7629                        pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    76307630                }
    76317631        }
    … …  
    76397639}
    76407640/*}}}*/
    7641 /*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{*/
    7642 ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
     7641/*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/
     7642ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){
    76437643
    76447644        /*Constants*/
    … …  
    76497649        int         approximation,analysis_type;
    76507650        IssmDouble  Jdet,Jdet2d;
    7651         IssmDouble  stokesreconditioning;
     7651        IssmDouble  FSreconditioning;
    76527652        IssmDouble      bed_normal[3];
    76537653        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    76647664        if(!IsOnBed() || IsFloating()) return NULL;
    76657665        inputs->GetInputValue(&approximation,ApproximationEnum);
    7666         if(approximation!=PattynStokesApproximationEnum) return NULL;
    7667         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     7666        if(approximation!=HOFSApproximationEnum) return NULL;
     7667        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    76687668
    76697669        /*Retrieve all inputs and parameters*/
    76707670        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    76717671        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    7672         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     7672        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    76737673        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    76747674        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    76757675        Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    7676         Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
     7676        Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
    76777677
    76787678        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    … …  
    76907690                GetNodalFunctionsP1(basis, gauss);
    76917691
    7692                 vzpattyn_input->GetInputValue(&w, gauss);
    7693                 vzpattyn_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     7692                vzHO_input->GetInputValue(&w, gauss);
     7693                vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    76947694
    76957695                BedNormal(&bed_normal[0],xyz_list_tria);
    76967696                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7697                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     7697                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    76987698                friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
    76997699
    … …  
    77237723                case MacAyealApproximationEnum:
    77247724                        return CreatePVectorDiagnosticMacAyeal();
    7725                 case PattynApproximationEnum:
    7726                         return CreatePVectorDiagnosticPattyn();
     7725                case HOApproximationEnum:
     7726                        return CreatePVectorDiagnosticHO();
    77277727                case L1L2ApproximationEnum:
    77287728                        return CreatePVectorDiagnosticL1L2();
    7729                 case HutterApproximationEnum:
     7729                case SIAApproximationEnum:
    77307730                        return NULL;
    77317731                case NoneApproximationEnum:
    77327732                        return NULL;
    7733                 case StokesApproximationEnum:
    7734                         return CreatePVectorDiagnosticStokes();
    7735                 case MacAyealPattynApproximationEnum:
    7736                         return CreatePVectorDiagnosticMacAyealPattyn();
    7737                 case MacAyealStokesApproximationEnum:
    7738                         return CreatePVectorDiagnosticMacAyealStokes();
    7739                 case PattynStokesApproximationEnum:
    7740                         return CreatePVectorDiagnosticPattynStokes();
     7733                case FSApproximationEnum:
     7734                        return CreatePVectorDiagnosticFS();
     7735                case MacAyealHOApproximationEnum:
     7736                        return CreatePVectorDiagnosticMacAyealHO();
     7737                case MacAyealFSApproximationEnum:
     7738                        return CreatePVectorDiagnosticMacAyealFS();
     7739                case HOFSApproximationEnum:
     7740                        return CreatePVectorDiagnosticHOFS();
    77417741                default:
    77427742                        _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    … …  
    77447744}
    77457745/*}}}*/
    7746 /*FUNCTION Penta::CreatePVectorDiagnosticMacAyealPattyn{{{*/
    7747 ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){
     7746/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealHO{{{*/
     7747ElementVector* Penta::CreatePVectorDiagnosticMacAyealHO(void){
    77487748
    77497749        /*compute all load vectors for this element*/
    77507750        ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
    7751         ElementVector* pe2=CreatePVectorDiagnosticPattyn();
     7751        ElementVector* pe2=CreatePVectorDiagnosticHO();
    77527752        ElementVector* pe =new ElementVector(pe1,pe2);
    77537753
    … …  
    77587758}
    77597759/*}}}*/
    7760 /*FUNCTION Penta::CreatePVectorDiagnosticMacAyealStokes{{{*/
    7761 ElementVector* Penta::CreatePVectorDiagnosticMacAyealStokes(void){
     7760/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealFS{{{*/
     7761ElementVector* Penta::CreatePVectorDiagnosticMacAyealFS(void){
    77627762
    77637763        /*compute all load vectors for this element*/
    77647764        ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
    7765         ElementVector* pe2=CreatePVectorDiagnosticStokes();
    7766         ElementVector* pe3=CreatePVectorCouplingMacAyealStokes();
     7765        ElementVector* pe2=CreatePVectorDiagnosticFS();
     7766        ElementVector* pe3=CreatePVectorCouplingMacAyealFS();
    77677767        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    77687768
    … …  
    77747774}
    77757775/*}}}*/
    7776 /*FUNCTION Penta::CreatePVectorDiagnosticPattynStokes{{{*/
    7777 ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){
     7776/*FUNCTION Penta::CreatePVectorDiagnosticHOFS{{{*/
     7777ElementVector* Penta::CreatePVectorDiagnosticHOFS(void){
    77787778
    77797779        /*compute all load vectors for this element*/
    7780         ElementVector* pe1=CreatePVectorDiagnosticPattyn();
    7781         ElementVector* pe2=CreatePVectorDiagnosticStokes();
    7782         ElementVector* pe3=CreatePVectorCouplingPattynStokes();
     7780        ElementVector* pe1=CreatePVectorDiagnosticHO();
     7781        ElementVector* pe2=CreatePVectorDiagnosticFS();
     7782        ElementVector* pe3=CreatePVectorCouplingHOFS();
    77837783        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    77847784
    … …  
    77907790}
    77917791/*}}}*/
    7792 /*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{*/
    7793 ElementVector* Penta::CreatePVectorDiagnosticHutter(void){
     7792/*FUNCTION Penta::CreatePVectorDiagnosticSIA{{{*/
     7793ElementVector* Penta::CreatePVectorDiagnosticSIA(void){
    77947794
    77957795        /*Intermediaries*/
    … …  
    79037903}
    79047904/*}}}*/
    7905 /*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{*/
    7906 ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){
     7905/*FUNCTION Penta::CreatePVectorDiagnosticHO{{{*/
     7906ElementVector* Penta::CreatePVectorDiagnosticHO(void){
    79077907
    79087908        /*compute all load vectors for this element*/
    7909         ElementVector* pe1=CreatePVectorDiagnosticPattynDrivingStress();
    7910         ElementVector* pe2=CreatePVectorDiagnosticPattynFront();
     7909        ElementVector* pe1=CreatePVectorDiagnosticHODrivingStress();
     7910        ElementVector* pe2=CreatePVectorDiagnosticHOFront();
    79117911        ElementVector* pe =new ElementVector(pe1,pe2);
    79127912
    … …  
    79177917}
    79187918/*}}}*/
    7919 /*FUNCTION Penta::CreatePVectorDiagnosticPattynDrivingStress{{{*/
    7920 ElementVector* Penta::CreatePVectorDiagnosticPattynDrivingStress(void){
     7919/*FUNCTION Penta::CreatePVectorDiagnosticHODrivingStress{{{*/
     7920ElementVector* Penta::CreatePVectorDiagnosticHODrivingStress(void){
    79217921
    79227922        /*Constants*/
    … …  
    79337933
    79347934        /*Initialize Element vector*/
    7935         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
     7935        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    79367936
    79377937        /*Retrieve all inputs and parameters*/
    … …  
    79657965}
    79667966/*}}}*/
    7967 /*FUNCTION Penta::CreatePVectorDiagnosticPattynFront{{{*/
    7968 ElementVector* Penta::CreatePVectorDiagnosticPattynFront(void){
     7967/*FUNCTION Penta::CreatePVectorDiagnosticHOFront{{{*/
     7968ElementVector* Penta::CreatePVectorDiagnosticHOFront(void){
    79697969
    79707970        /*Intermediaries */
    … …  
    80008000
    80018001        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8002         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,PattynApproximationEnum);
     8002        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
    80038003        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    80048004        rho_water=matpar->GetRhoWater();
    … …  
    80438043}
    80448044/*}}}*/
    8045 /*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{*/
    8046 ElementVector* Penta::CreatePVectorDiagnosticStokes(void){
    8047 
    8048         int fe_stokes;
     8045/*FUNCTION Penta::CreatePVectorDiagnosticFS {{{*/
     8046ElementVector* Penta::CreatePVectorDiagnosticFS(void){
     8047
     8048        int fe_FS;
    80498049        ElementVector* pe1;
    80508050        ElementVector* pe2;
    80518051        ElementVector* pe3;
    80528052        ElementVector* pe;
    8053         parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum);
    8054 
    8055         switch(fe_stokes){
     8053        parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
     8054
     8055        switch(fe_FS){
    80568056                case 0:
    80578057                        /*compute all stiffness matrices for this element*/
    8058                         pe1=CreatePVectorDiagnosticStokesViscous();
    8059                         pe2=CreatePVectorDiagnosticStokesShelf();
    8060                         pe3=CreatePVectorDiagnosticStokesFront();
     8058                        pe1=CreatePVectorDiagnosticFSViscous();
     8059                        pe2=CreatePVectorDiagnosticFSShelf();
     8060                        pe3=CreatePVectorDiagnosticFSFront();
    80618061                        pe =new ElementVector(pe1,pe2,pe3);
    80628062                        break;
    80638063                case 1:
    80648064                        /*compute all stiffness matrices for this element*/
    8065                         pe1=CreatePVectorDiagnosticStokesGLSViscous();
    8066                         pe2=CreatePVectorDiagnosticStokesShelf();
    8067                         pe3=CreatePVectorDiagnosticStokesFront();
     8065                        pe1=CreatePVectorDiagnosticFSGLSViscous();
     8066                        pe2=CreatePVectorDiagnosticFSShelf();
     8067                        pe3=CreatePVectorDiagnosticFSFront();
    80688068                        pe =new ElementVector(pe1,pe2,pe3);
    80698069                        break;
    80708070                default:
    8071                         _error_("Finite element" << fe_stokes << " not supported yet");
     8071                        _error_("Finite element" << fe_FS << " not supported yet");
    80728072        }
    80738073
    … …  
    80808080}
    80818081/*}}}*/
    8082 /*FUNCTION Penta::CreatePVectorDiagnosticStokesViscous {{{*/
    8083 ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){
     8082/*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/
     8083ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){
    80848084
    80858085        /*Constants*/
    … …  
    80908090        int        approximation;
    80918091        IssmDouble Jdet,viscosity;
    8092         IssmDouble gravity,rho_ice,stokesreconditioning;
     8092        IssmDouble gravity,rho_ice,FSreconditioning;
    80938093        IssmDouble forcex,forcey,forcez;
    80948094        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    81068106        /*Initialize Element vector and return if necessary*/
    81078107        inputs->GetInputValue(&approximation,ApproximationEnum);
    8108         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    8109         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     8108        if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     8109        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    81108110
    81118111        /*Retrieve all inputs and parameters*/
    8112         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     8112        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    81138113        rho_ice=matpar->GetRhoIce();
    81148114        gravity=matpar->GetG();
    … …  
    81288128
    81298129                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8130                 GetBStokes(&B[0][0],&xyz_list[0][0],gauss);
    8131                 GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss);
     8130                GetBFS(&B[0][0],&xyz_list[0][0],gauss);
     8131                GetBprimeFS(&B_prime[0][0],&xyz_list[0][0], gauss);
    81328132                GetNodalFunctionsMINI(&l1l7[0], gauss);
    81338133
    81348134                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8135                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     8135                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    81368136
    81378137                loadingforcex_input->GetInputValue(&forcex, gauss);
    … …  
    81518151                D_scalar=gauss->weight*Jdet;
    81528152                for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
    8153                 for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
     8153                for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
    81548154
    81558155                TripleMultiply(&B[0][0],8,numdofbubble,1,
    … …  
    81608160
    81618161        /*Condensation*/
    8162         ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
     8162        ReduceVectorFS(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
    81638163
    81648164        /*Transform coordinate system*/
    … …  
    81708170}
    81718171/*}}}*/
    8172 /*FUNCTION Penta::CreatePVectorDiagnosticStokesFront{{{*/
    8173 ElementVector* Penta::CreatePVectorDiagnosticStokesFront(void){
     8172/*FUNCTION Penta::CreatePVectorDiagnosticFSFront{{{*/
     8173ElementVector* Penta::CreatePVectorDiagnosticFSFront(void){
    81748174
    81758175        /*Intermediaries */
    … …  
    82058205
    82068206        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8207         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,StokesApproximationEnum);
     8207        ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSApproximationEnum);
    82088208        rho_water=matpar->GetRhoWater();
    82098209        rho_ice  =matpar->GetRhoIce();
    … …  
    82478247}
    82488248/*}}}*/
    8249 /*FUNCTION Penta::CreatePVectorDiagnosticStokesGLSViscous {{{*/
    8250 ElementVector* Penta::CreatePVectorDiagnosticStokesGLSViscous(void){
     8249/*FUNCTION Penta::CreatePVectorDiagnosticFSGLSViscous {{{*/
     8250ElementVector* Penta::CreatePVectorDiagnosticFSGLSViscous(void){
    82518251
    82528252        /*Constants*/
    … …  
    82578257        int        approximation;
    82588258        IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity;
    8259         IssmDouble forcex,forcey,forcez,diameter,stokesreconditioning;
     8259        IssmDouble forcex,forcey,forcez,diameter,FSreconditioning;
    82608260        IssmDouble xyz_list[NUMVERTICES][3];
    82618261        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    82768276        /*Initialize Element vector and return if necessary*/
    82778277        inputs->GetInputValue(&approximation,ApproximationEnum);
    8278         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    8279         parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
    8280         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     8278        if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     8279        parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
     8280        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    82818281
    82828282        /*Retrieve all inputs and parameters*/
    … …  
    83188318                GetNodalFunctionsP1(&l1l6[0], gauss);
    83198319                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8320                 material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     8320                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    83218321
    83228322                loadingforcex_input->GetInputValue(&forcex, gauss);
    … …  
    83728372}
    83738373/*}}}*/
    8374 /*FUNCTION Penta::CreatePVectorDiagnosticStokesShelf{{{*/
    8375 ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){
     8374/*FUNCTION Penta::CreatePVectorDiagnosticFSShelf{{{*/
     8375ElementVector* Penta::CreatePVectorDiagnosticFSShelf(void){
    83768376
    83778377        /*Intermediaries*/
    … …  
    83928392        inputs->GetInputValue(&approximation,ApproximationEnum);
    83938393        this->parameters->FindParam(&shelf_dampening,DiagnosticShelfDampeningEnum);
    8394         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
    8395         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     8394        if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
     8395        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    83968396
    83978397        /*Retrieve all inputs and parameters*/
    … …  
    84778477        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    84788478        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    8479         Input* vzstokes_input=NULL;
    8480         if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
    8481                 vzstokes_input=inputs->GetInput(VzStokesEnum); _assert_(vzstokes_input);
     8479        Input* vzFS_input=NULL;
     8480        if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
     8481                vzFS_input=inputs->GetInput(VzFSEnum); _assert_(vzFS_input);
    84828482        }
    84838483
    … …  
    84938493                vx_input->GetInputDerivativeValue(&du[0],&xyz_list[0][0],gauss);
    84948494                vy_input->GetInputDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
    8495                 if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
    8496                         vzstokes_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
     8495                if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
     8496                        vzFS_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    84978497                        dwdz=dw[2];
    84988498                }
    … …  
    85398539        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    85408540        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    8541         Input* vzstokes_input=NULL;
    8542         if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
    8543                 vzstokes_input=inputs->GetInput(VzStokesEnum);       _assert_(vzstokes_input);
     8541        Input* vzFS_input=NULL;
     8542        if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
     8543                vzFS_input=inputs->GetInput(VzFSEnum);       _assert_(vzFS_input);
    85448544        }
    85458545
    … …  
    85548554                vx_input->GetInputValue(&vx, gauss);
    85558555                vy_input->GetInputValue(&vy, gauss);
    8556                 if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
    8557                         vzstokes_input->GetInputValue(&vz, gauss);
     8556                if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
     8557                        vzFS_input->GetInputValue(&vz, gauss);
    85588558                }
    85598559                else vz=0;
    … …  
    85818581        switch(approximation){
    85828582                case MacAyealApproximationEnum:
    8583                         return CreateJacobianDiagnosticMacayeal2d();
    8584                 case PattynApproximationEnum:
    8585                         return CreateJacobianDiagnosticPattyn();
    8586                 case StokesApproximationEnum:
    8587                         return CreateJacobianDiagnosticStokes();
     8583                        return CreateJacobianDiagnosticSSA2d();
     8584                case HOApproximationEnum:
     8585                        return CreateJacobianDiagnosticHO();
     8586                case FSApproximationEnum:
     8587                        return CreateJacobianDiagnosticFS();
    85888588                case NoneApproximationEnum:
    85898589                        return NULL;
    … …  
    85938593}
    85948594/*}}}*/
    8595 /*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{*/
    8596 ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){
     8595/*FUNCTION Penta::CreateJacobianDiagnosticSSA2d{{{*/
     8596ElementMatrix* Penta::CreateJacobianDiagnosticSSA2d(void){
    85978597
    85988598        /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
    … …  
    86168616        /*Call Tria function*/
    86178617        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    8618         ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal();
     8618        ElementMatrix* Ke=tria->CreateJacobianDiagnosticSSA();
    86198619        delete tria->material; delete tria;
    86208620
    … …  
    86278627}
    86288628/*}}}*/
    8629 /*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{*/
    8630 ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
     8629/*FUNCTION Penta::CreateJacobianDiagnosticHO{{{*/
     8630ElementMatrix* Penta::CreateJacobianDiagnosticHO(void){
    86318631
    86328632        /*Constants*/
    … …  
    86458645        GaussPenta *gauss=NULL;
    86468646
    8647         /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
    8648         ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
     8647        /*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
     8648        ElementMatrix* Ke=CreateKMatrixDiagnosticHO();
    86498649
    86508650        /*Retrieve all inputs and parameters*/
    … …  
    86628662                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    86638663
    8664                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     8664                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    86658665                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    86668666                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    … …  
    86918691}
    86928692/*}}}*/
    8693 /*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{*/
    8694 ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
     8693/*FUNCTION Penta::CreateJacobianDiagnosticFS{{{*/
     8694ElementMatrix* Penta::CreateJacobianDiagnosticFS(void){
    86958695
    86968696        /*Constants*/
    … …  
    87108710        GaussPenta *gauss=NULL;
    87118711
    8712         /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
    8713         ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
     8712        /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
     8713        ElementMatrix* Ke=CreateKMatrixDiagnosticFS();
    87148714
    87158715        /*Retrieve all inputs and parameters*/
    … …  
    87288728                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    87298729
    8730                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     8730                this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    87318731                material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    87328732                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    … …  
    87848784
    87858785        /*If the element is a coupling, do nothing: every node is also on an other elements
    8786          * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
     8786         * (as coupling is between MacAyeal and HO) so the other element will take care of it*/
    87878787        GetDofList(&doflist,approximation,GsetEnum);
    87888788
    … …  
    88088808}
    88098809/*}}}*/
    8810 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{*/
    8811 void  Penta::GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution){
     8810/*FUNCTION Penta::GetSolutionFromInputsDiagnosticSIA{{{*/
     8811void  Penta::GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution){
    88128812
    88138813        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    88778877}
    88788878/*}}}*/
    8879 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{*/
    8880 void  Penta::GetSolutionFromInputsDiagnosticStokes(Vector<IssmDouble>* solution){
     8879/*FUNCTION Penta::GetSolutionFromInputsDiagnosticFS{{{*/
     8880void  Penta::GetSolutionFromInputsDiagnosticFS(Vector<IssmDouble>* solution){
    88818881
    88828882        const int    numdof=NDOF4*NUMVERTICES;
    … …  
    88858885        int*         doflist=NULL;
    88868886        IssmDouble       vx,vy,vz,p;
    8887         IssmDouble       stokesreconditioning;
     8887        IssmDouble       FSreconditioning;
    88888888        IssmDouble       values[numdof];
    88898889        GaussPenta   *gauss;
    88908890
    88918891        /*Get dof list: */
    8892         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
     8892        GetDofList(&doflist,FSApproximationEnum,GsetEnum);
    88938893        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    88948894        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    … …  
    88978897
    88988898        /*Recondition pressure: */
    8899         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     8899        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    89008900
    89018901        /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
    … …  
    89118911                values[i*NDOF4+1]=vy;
    89128912                values[i*NDOF4+2]=vz;
    8913                 values[i*NDOF4+3]=p/stokesreconditioning;
     8913                values[i*NDOF4+3]=p/FSreconditioning;
    89148914        }
    89158915
    … …  
    89618961
    89628962        /* Get eps_b*/
    8963         vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
    8964         vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     8963        vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
     8964        vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
    89658965        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    89668966        eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
    … …  
    90159015                return;
    90169016        }
    9017         else if (approximation==PattynApproximationEnum){
    9018                 InputUpdateFromSolutionDiagnosticPattyn(solution);
    9019         }
    9020         else if (approximation==PattynStokesApproximationEnum){
    9021                 InputUpdateFromSolutionDiagnosticPattynStokes(solution);
    9022         }
    9023         else if (approximation==MacAyealStokesApproximationEnum){
    9024                 InputUpdateFromSolutionDiagnosticMacAyealStokes(solution);
    9025         }
    9026         else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    9027                 InputUpdateFromSolutionDiagnosticStokes(solution);
    9028         }
    9029         else if (approximation==MacAyealPattynApproximationEnum){
    9030                 InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution);
     9017        else if (approximation==HOApproximationEnum){
     9018                InputUpdateFromSolutionDiagnosticHO(solution);
     9019        }
     9020        else if (approximation==HOFSApproximationEnum){
     9021                InputUpdateFromSolutionDiagnosticHOFS(solution);
     9022        }
     9023        else if (approximation==MacAyealFSApproximationEnum){
     9024                InputUpdateFromSolutionDiagnosticMacAyealFS(solution);
     9025        }
     9026        else if (approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
     9027                InputUpdateFromSolutionDiagnosticFS(solution);
     9028        }
     9029        else if (approximation==MacAyealHOApproximationEnum){
     9030                InputUpdateFromSolutionDiagnosticMacAyealHO(solution);
    90319031        }
    90329032}
    … …  
    91139113}
    91149114/*}}}*/
    9115 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{*/
    9116 void  Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(IssmDouble* solution){
     9115/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealHO {{{*/
     9116void  Penta::InputUpdateFromSolutionDiagnosticMacAyealHO(IssmDouble* solution){
    91179117
    91189118        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    91219121        int     i;
    91229122        IssmDouble  rho_ice,g;
    9123         IssmDouble  macayeal_values[numdof];
    9124         IssmDouble  pattyn_values[numdof];
     9123        IssmDouble  SSA_values[numdof];
     9124        IssmDouble  HO_values[numdof];
    91259125        IssmDouble  vx[NUMVERTICES];
    91269126        IssmDouble  vy[NUMVERTICES];
    … …  
    91349134        Penta   *penta   = NULL;
    91359135
    9136         /*OK, we have to add results of this element for pattyn
    9137          * and results from the penta at base for macayeal. Now recover results*/
     9136        /*OK, we have to add results of this element for HO
     9137         * and results from the penta at base for SSA. Now recover results*/
    91389138        penta=GetBasalElement();
    91399139
    9140         /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
    9141         GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
     9140        /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
     9141        GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
    91429142        penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
    91439143
    … …  
    91479147        /*Use the dof list to index into the solution vector: */
    91489148        for(i=0;i<numdof2d;i++){
    9149                 pattyn_values[i]=solution[doflistp[i]];
    9150                 macayeal_values[i]=solution[doflistm[i]];
     9149                HO_values[i]=solution[doflistp[i]];
     9150                SSA_values[i]=solution[doflistm[i]];
    91519151        }
    91529152        for(i=numdof2d;i<numdof;i++){
    9153                 pattyn_values[i]=solution[doflistp[i]];
    9154                 macayeal_values[i]=macayeal_values[i-numdof2d];
     9153                HO_values[i]=solution[doflistp[i]];
     9154                SSA_values[i]=SSA_values[i-numdof2d];
    91559155        }
    91569156
    91579157        /*Transform solution in Cartesian Space*/
    9158         TransformSolutionCoord(&macayeal_values[0],penta->nodes,NUMVERTICES,XYEnum);
    9159         TransformSolutionCoord(&pattyn_values[0],   this->nodes,NUMVERTICES,XYEnum);
     9158        TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
     9159        TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
    91609160
    91619161        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    91629162        for(i=0;i<NUMVERTICES;i++){
    9163                 vx[i]=macayeal_values[i*NDOF2+0]+pattyn_values[i*NDOF2+0];
    9164                 vy[i]=macayeal_values[i*NDOF2+1]+pattyn_values[i*NDOF2+1];
     9163                vx[i]=SSA_values[i*NDOF2+0]+HO_values[i*NDOF2+0];
     9164                vy[i]=SSA_values[i*NDOF2+1]+HO_values[i*NDOF2+1];
    91659165
    91669166                /*Check solution*/
    … …  
    91979197}
    91989198/*}}}*/
    9199 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{*/
    9200 void  Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(IssmDouble* solution){
     9199/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealFS {{{*/
     9200void  Penta::InputUpdateFromSolutionDiagnosticMacAyealFS(IssmDouble* solution){
    92019201
    92029202        const int    numdofm=NDOF2*NUMVERTICES;
    … …  
    92059205
    92069206        int     i;
    9207         IssmDouble  stokesreconditioning;
    9208         IssmDouble  macayeal_values[numdofm];
    9209         IssmDouble  stokes_values[numdofs];
     9207        IssmDouble  FSreconditioning;
     9208        IssmDouble  SSA_values[numdofm];
     9209        IssmDouble  FS_values[numdofs];
    92109210        IssmDouble  vx[NUMVERTICES];
    92119211        IssmDouble  vy[NUMVERTICES];
    92129212        IssmDouble  vz[NUMVERTICES];
    9213         IssmDouble  vzmacayeal[NUMVERTICES];
    9214         IssmDouble  vzstokes[NUMVERTICES];
     9213        IssmDouble  vzSSA[NUMVERTICES];
     9214        IssmDouble  vzFS[NUMVERTICES];
    92159215        IssmDouble  vel[NUMVERTICES];
    92169216        IssmDouble  pressure[NUMVERTICES];
    … …  
    92209220        Penta   *penta          = NULL;
    92219221
    9222         /*OK, we have to add results of this element for macayeal
    9223          * and results from the penta at base for macayeal. Now recover results*/
     9222        /*OK, we have to add results of this element for SSA
     9223         * and results from the penta at base for SSA. Now recover results*/
    92249224        penta=GetBasalElement();
    92259225
    9226         /*Get dof listof this element (macayeal dofs) and of the penta at base (macayeal dofs): */
     9226        /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
    92279227        penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
    9228         GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
    9229         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     9228        GetDofList(&doflists,FSApproximationEnum,GsetEnum);
     9229        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    92309230
    92319231        /*Get node data: */
    … …  
    92349234        /*Use the dof list to index into the solution vector: */
    92359235        for(i=0;i<numdof2d;i++){
    9236                 macayeal_values[i]=solution[doflistm[i]];
    9237                 macayeal_values[i+numdof2d]=solution[doflistm[i]];
     9236                SSA_values[i]=solution[doflistm[i]];
     9237                SSA_values[i+numdof2d]=solution[doflistm[i]];
    92389238        }
    92399239        for(i=0;i<numdofs;i++){
    9240                 stokes_values[i]=solution[doflists[i]];
     9240                FS_values[i]=solution[doflists[i]];
    92419241        }
    92429242
    92439243        /*Transform solution in Cartesian Space*/
    9244         TransformSolutionCoord(&macayeal_values[0],this->nodes,NUMVERTICES,XYEnum);
    9245         TransformSolutionCoord(&stokes_values[0],this->nodes,NUMVERTICES,XYZPEnum);
     9244        TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
     9245        TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
    92469246
    92479247        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    92489248        for(i=0;i<NUMVERTICES;i++){
    9249                 vx[i]=stokes_values[i*NDOF4+0]+macayeal_values[i*NDOF2+0];
    9250                 vy[i]=stokes_values[i*NDOF4+1]+macayeal_values[i*NDOF2+1];
    9251                 vzstokes[i]=stokes_values[i*NDOF4+2];
    9252                 pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
     9249                vx[i]=FS_values[i*NDOF4+0]+SSA_values[i*NDOF2+0];
     9250                vy[i]=FS_values[i*NDOF4+1]+SSA_values[i*NDOF2+1];
     9251                vzFS[i]=FS_values[i*NDOF4+2];
     9252                pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
    92539253
    92549254                /*Check solution*/
    92559255                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    92569256                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    9257                 if(xIsNan<IssmDouble>(vzstokes[i])) _error_("NaN found in solution vector");
     9257                if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
    92589258                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    92599259        }
    92609260
    92619261        /*Get Vz*/
    9262         Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);
    9263         if (vzmacayeal_input){
    9264                 if (vzmacayeal_input->ObjectEnum()!=PentaInputEnum){
    9265                         _error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzmacayeal_input->ObjectEnum()));
    9266                 }
    9267                 GetInputListOnVertices(&vzmacayeal[0],VzMacAyealEnum);
     9262        Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);
     9263        if (vzSSA_input){
     9264                if (vzSSA_input->ObjectEnum()!=PentaInputEnum){
     9265                        _error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzSSA_input->ObjectEnum()));
     9266                }
     9267                GetInputListOnVertices(&vzSSA[0],VzMacAyealEnum);
    92689268        }
    92699269        else{
    … …  
    92739273        /*Now Compute vel*/
    92749274        for(i=0;i<NUMVERTICES;i++) {
    9275                 vz[i]=vzmacayeal[i]+vzstokes[i];
     9275                vz[i]=vzSSA[i]+vzFS[i];
    92769276                vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    92779277        }
    … …  
    92889288        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    92899289        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    9290         this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
     9290        this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
    92919291        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    92929292        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    … …  
    93789378}
    93799379/*}}}*/
    9380 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{*/
    9381 void  Penta::InputUpdateFromSolutionDiagnosticPattyn(IssmDouble* solution){
     9380/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHO {{{*/
     9381void  Penta::InputUpdateFromSolutionDiagnosticHO(IssmDouble* solution){
    93829382
    93839383        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    93969396
    93979397        /*Get dof list: */
    9398         GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
     9398        GetDofList(&doflist,HOApproximationEnum,GsetEnum);
    93999399
    94009400        /*Get node data: */
    … …  
    94529452}
    94539453/*}}}*/
    9454 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{*/
    9455 void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(IssmDouble* solution){
     9454/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHOFS {{{*/
     9455void  Penta::InputUpdateFromSolutionDiagnosticHOFS(IssmDouble* solution){
    94569456
    94579457        const int    numdofp=NDOF2*NUMVERTICES;
    … …  
    94599459
    94609460        int    i;
    9461         IssmDouble pattyn_values[numdofp];
    9462         IssmDouble stokes_values[numdofs];
     9461        IssmDouble HO_values[numdofp];
     9462        IssmDouble FS_values[numdofs];
    94639463        IssmDouble vx[NUMVERTICES];
    94649464        IssmDouble vy[NUMVERTICES];
    94659465        IssmDouble vz[NUMVERTICES];
    9466         IssmDouble vzpattyn[NUMVERTICES];
    9467         IssmDouble vzstokes[NUMVERTICES];
     9466        IssmDouble vzHO[NUMVERTICES];
     9467        IssmDouble vzFS[NUMVERTICES];
    94689468        IssmDouble vel[NUMVERTICES];
    94699469        IssmDouble pressure[NUMVERTICES];
    94709470        IssmDouble xyz_list[NUMVERTICES][3];
    9471         IssmDouble stokesreconditioning;
     9471        IssmDouble FSreconditioning;
    94729472        int*   doflistp      = NULL;
    94739473        int*   doflists      = NULL;
    94749474        Penta  *penta        = NULL;
    94759475
    9476         /*OK, we have to add results of this element for pattyn
    9477          * and results from the penta at base for macayeal. Now recover results*/
     9476        /*OK, we have to add results of this element for HO
     9477         * and results from the penta at base for SSA. Now recover results*/
    94789478        penta=GetBasalElement();
    94799479
    9480         /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
    9481         GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
    9482         GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
    9483         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
     9480        /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
     9481        GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
     9482        GetDofList(&doflists,FSApproximationEnum,GsetEnum);
     9483        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    94849484
    94859485        /*Get node data: */
    … …  
    94879487
    94889488        /*Use the dof list to index into the solution vector: */
    9489         for(i=0;i<numdofp;i++) pattyn_values[i]=solution[doflistp[i]];
    9490         for(i=0;i<numdofs;i++) stokes_values[i]=solution[doflists[i]];
     9489        for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
     9490        for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
    94919491
    94929492        /*Transform solution in Cartesian Space*/
    9493         TransformSolutionCoord(&pattyn_values[0],this->nodes,NUMVERTICES,XYEnum);
    9494         TransformSolutionCoord(&stokes_values[0],this->nodes,NUMVERTICES,XYZPEnum);
     9493        TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
     9494        TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
    94959495
    94969496        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    94979497        for(i=0;i<NUMVERTICES;i++){
    9498                 vx[i]=stokes_values[i*NDOF4+0]+pattyn_values[i*NDOF2+0];
    9499                 vy[i]=stokes_values[i*NDOF4+1]+pattyn_values[i*NDOF2+1];
    9500                 vzstokes[i]=stokes_values[i*NDOF4+2];
    9501                 pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
     9498                vx[i]=FS_values[i*NDOF4+0]+HO_values[i*NDOF2+0];
     9499                vy[i]=FS_values[i*NDOF4+1]+HO_values[i*NDOF2+1];
     9500                vzFS[i]=FS_values[i*NDOF4+2];
     9501                pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
    95029502
    95039503                /*Check solution*/
    95049504                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    95059505                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    9506                 if(xIsNan<IssmDouble>(vzstokes[i])) _error_("NaN found in solution vector");
     9506                if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
    95079507                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    95089508        }
    95099509
    95109510        /*Get Vz*/
    9511         Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);
    9512         if (vzpattyn_input){
    9513                 if (vzpattyn_input->ObjectEnum()!=PentaInputEnum){
    9514                         _error_("Cannot compute Vel as VzPattyn is of type " << EnumToStringx(vzpattyn_input->ObjectEnum()));
    9515                 }
    9516                 GetInputListOnVertices(&vzpattyn[0],VzPattynEnum);
     9511        Input* vzHO_input=inputs->GetInput(VzHOEnum);
     9512        if (vzHO_input){
     9513                if (vzHO_input->ObjectEnum()!=PentaInputEnum){
     9514                        _error_("Cannot compute Vel as VzHO is of type " << EnumToStringx(vzHO_input->ObjectEnum()));
     9515                }
     9516                GetInputListOnVertices(&vzHO[0],VzHOEnum);
    95179517        }
    95189518        else{
    9519                 _error_("Cannot update solution as VzPattyn is not present");
     9519                _error_("Cannot update solution as VzHO is not present");
    95209520        }
    95219521
    95229522        /*Now Compute vel*/
    95239523        for(i=0;i<NUMVERTICES;i++) {
    9524                 vz[i]=vzpattyn[i]+vzstokes[i];
     9524                vz[i]=vzHO[i]+vzFS[i];
    95259525                vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    95269526        }
    … …  
    95379537        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    95389538        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    9539         this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
     9539        this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
    95409540        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    95419541        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    … …  
    95469546}
    95479547/*}}}*/
    9548 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{*/
    9549 void  Penta::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){
     9548/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticSIA {{{*/
     9549void  Penta::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
    95509550
    95519551        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    96219621        IssmDouble   vy[NUMVERTICES];
    96229622        IssmDouble   vz[NUMVERTICES];
    9623         IssmDouble   vzmacayeal[NUMVERTICES];
    9624         IssmDouble   vzpattyn[NUMVERTICES];
    9625         IssmDouble   vzstokes[NUMVERTICES];
     9623        IssmDouble   vzSSA[NUMVERTICES];
     9624        IssmDouble   vzHO[NUMVERTICES];
     9625        IssmDouble   vzFS[NUMVERTICES];
    96269626        IssmDouble   vel[NUMVERTICES];
    96279627        IssmDouble   pressure[NUMVERTICES];
    … …  
    96309630        int*     doflist      = NULL;
    96319631
    9632         /*Get the approximation and do nothing if the element in Stokes or None*/
     9632        /*Get the approximation and do nothing if the element in FS or None*/
    96339633        inputs->GetInputValue(&approximation,ApproximationEnum);
    9634         if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
     9634        if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
    96359635                return;
    96369636        }
    … …  
    96539653        GetInputListOnVertices(&vy[0],VyEnum,0.0); //default is 0
    96549654
    9655         /*Do some modifications if we actually have a PattynStokes or MacAyealStokes element*/
    9656         if(approximation==PattynStokesApproximationEnum){
    9657                 Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    9658                 if (vzstokes_input){
    9659                         if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
    9660                         GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    9661                 }
    9662                 else _error_("Cannot compute Vz as VzStokes in not present in PattynStokes element");
     9655        /*Do some modifications if we actually have a HOFS or MacAyealFS element*/
     9656        if(approximation==HOFSApproximationEnum){
     9657                Input* vzFS_input=inputs->GetInput(VzFSEnum);
     9658                if (vzFS_input){
     9659                        if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
     9660                        GetInputListOnVertices(&vzFS[0],VzFSEnum);
     9661                }
     9662                else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
    96639663                for(i=0;i<NUMVERTICES;i++){
    9664                         vzpattyn[i]=vz[i];
    9665                         vz[i]=vzpattyn[i]+vzstokes[i];
    9666                 }
    9667         }
    9668         else if(approximation==MacAyealStokesApproximationEnum){
    9669                 Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    9670                 if (vzstokes_input){
    9671                         if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
    9672                         GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    9673                 }
    9674                 else _error_("Cannot compute Vz as VzStokes in not present in MacAyealStokes element");
     9664                        vzHO[i]=vz[i];
     9665                        vz[i]=vzHO[i]+vzFS[i];
     9666                }
     9667        }
     9668        else if(approximation==MacAyealFSApproximationEnum){
     9669                Input* vzFS_input=inputs->GetInput(VzFSEnum);
     9670                if (vzFS_input){
     9671                        if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
     9672                        GetInputListOnVertices(&vzFS[0],VzFSEnum);
     9673                }
     9674                else _error_("Cannot compute Vz as VzFS in not present in MacAyealFS element");
    96759675                for(i=0;i<NUMVERTICES;i++){
    9676                         vzmacayeal[i]=vz[i];
    9677                         vz[i]=vzmacayeal[i]+vzstokes[i];
     9676                        vzSSA[i]=vz[i];
     9677                        vz[i]=vzSSA[i]+vzFS[i];
    96789678                }
    96799679        }
    … …  
    96839683
    96849684        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    9685          *so the pressure is just the pressure at the z elevation: except it this is a PattynStokes element */
    9686         if(approximation!=PattynStokesApproximationEnum &&  approximation!=MacAyealStokesApproximationEnum){
     9685         *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
     9686        if(approximation!=HOFSApproximationEnum &&  approximation!=MacAyealFSApproximationEnum){
    96879687                rho_ice=matpar->GetRhoIce();
    96889688                g=matpar->GetG();
    … …  
    96959695        this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
    96969696
    9697         if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){
     9697        if(approximation!=HOFSApproximationEnum && approximation!=MacAyealFSApproximationEnum){
    96989698                this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    96999699                this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    97009700        }
    9701         else if(approximation==PattynStokesApproximationEnum){
    9702                 this->inputs->AddInput(new PentaInput(VzPattynEnum,vzpattyn,P1Enum));
    9703         }
    9704         else if(approximation==MacAyealStokesApproximationEnum){
    9705                 this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzmacayeal,P1Enum));
     9701        else if(approximation==HOFSApproximationEnum){
     9702                this->inputs->AddInput(new PentaInput(VzHOEnum,vzHO,P1Enum));
     9703        }
     9704        else if(approximation==MacAyealFSApproximationEnum){
     9705                this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzSSA,P1Enum));
    97069706        }
    97079707        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    … …  
    97129712}
    97139713/*}}}*/
    9714 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticStokes {{{*/
    9715 void  Penta::InputUpdateFromSolutionDiagnosticStokes(IssmDouble* solution){
     9714/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticFS {{{*/
     9715void  Penta::InputUpdateFromSolutionDiagnosticFS(IssmDouble* solution){
    97169716
    97179717        const int numdof=NDOF4*NUMVERTICES;
    … …  
    97249724        IssmDouble  vel[NUMVERTICES];
    97259725        IssmDouble  pressure[NUMVERTICES];
    9726         IssmDouble  stokesreconditioning;
     9726        IssmDouble  FSreconditioning;
    97279727        int*    doflist=NULL;
    97289728
    97299729        /*Get dof list: */
    9730         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
     9730        GetDofList(&doflist,FSApproximationEnum,GsetEnum);
    97319731
    97329732        /*Use the dof list to index into the solution vector: */
    … …  
    97519751
    97529752        /*Recondition pressure and compute vel: */
    9753         this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
    9754         for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*stokesreconditioning;
     9753        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
     9754        for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*FSreconditioning;
    97559755        for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
    97569756

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

    r15538 r15564  
    145145                void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
    146146                void   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
    147                 void   GradjDragPattyn(Vector<IssmDouble>* gradient,int control_index);
    148                 void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
     147                void   GradjDragHO(Vector<IssmDouble>* gradient,int control_index);
     148                void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
    149149                void   GradjBbarMacAyeal(Vector<IssmDouble>* gradient,int control_index);
    150                 void   GradjBbarPattyn(Vector<IssmDouble>* gradient,int control_index);
    151                 void   GradjBbarStokes(Vector<IssmDouble>* gradient,int control_index);
     150                void   GradjBbarHO(Vector<IssmDouble>* gradient,int control_index);
     151                void   GradjBbarFS(Vector<IssmDouble>* gradient,int control_index);
    152152                void   GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
    153153                void   SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
    … …  
    197197                void             GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
    198198                IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
    199                 void    GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
     199                void    GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
    200200                void    GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    201201                Penta*  GetUpperElement(void);
    … …  
    215215                bool    IsOnWater(void);
    216216                IssmDouble  MinEdgeLength(IssmDouble xyz_list[6][3]);
    217                 void      ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
    218                 void      ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
     217                void      ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
     218                void      ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
    219219                void      SetClone(int* minranks);
    220220                Tria*     SpawnTria(int g0, int g1, int g2);
    … …  
    223223
    224224                #ifdef _HAVE_DIAGNOSTIC_
    225                 ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void);
    226                 ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void);
    227                 ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
    228                 ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void);
    229                 ElementMatrix* CreateKMatrixCouplingMacAyealStokesViscous(void);
    230                 ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void);
    231                 ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
     225                ElementMatrix* CreateKMatrixCouplingMacAyealHO(void);
     226                ElementMatrix* CreateKMatrixCouplingMacAyealHOViscous(void);
     227                ElementMatrix* CreateKMatrixCouplingMacAyealHOFriction(void);
     228                ElementMatrix* CreateKMatrixCouplingMacAyealFS(void);
     229                ElementMatrix* CreateKMatrixCouplingMacAyealFSViscous(void);
     230                ElementMatrix* CreateKMatrixCouplingMacAyealFSFriction(void);
     231                ElementMatrix* CreateKMatrixCouplingHOFS(void);
    232232                ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
    233233                ElementMatrix* CreateKMatrixAdjointHoriz(void);
    234234                ElementVector* CreateDVectorDiagnosticHoriz(void);
    235                 ElementVector* CreateDVectorDiagnosticStokes(void);
    236                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
     235                ElementVector* CreateDVectorDiagnosticFS(void);
     236                ElementMatrix* CreateKMatrixDiagnosticSIA(void);
    237237                ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
    238238                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);
    239239                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dViscous(void);
    240240                ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
    241                 ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void);
    242                 ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void);
     241                ElementMatrix* CreateKMatrixDiagnosticMacAyealHO(void);
     242                ElementMatrix* CreateKMatrixDiagnosticMacAyealFS(void);
    243243                ElementMatrix* CreateKMatrixDiagnosticL1L2(void);
    244244                ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void);
    245245                ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void);
    246                 ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
    247                 ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
    248                 ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
    249                 ElementMatrix* CreateKMatrixDiagnosticPattynStokes(void);
    250                 ElementMatrix* CreateKMatrixDiagnosticStokes(void);
    251                 ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
    252                 ElementMatrix* CreateKMatrixDiagnosticStokesGLSViscous(void);
    253                 ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
    254                 ElementMatrix* CreateKMatrixDiagnosticStokesGLSFriction(void);
     246                ElementMatrix* CreateKMatrixDiagnosticHO(void);
     247                ElementMatrix* CreateKMatrixDiagnosticHOViscous(void);
     248                ElementMatrix* CreateKMatrixDiagnosticHOFriction(void);
     249                ElementMatrix* CreateKMatrixDiagnosticHOFS(void);
     250                ElementMatrix* CreateKMatrixDiagnosticFS(void);
     251                ElementMatrix* CreateKMatrixDiagnosticFSViscous(void);
     252                ElementMatrix* CreateKMatrixDiagnosticFSGLSViscous(void);
     253                ElementMatrix* CreateKMatrixDiagnosticFSFriction(void);
     254                ElementMatrix* CreateKMatrixDiagnosticFSGLSFriction(void);
    255255                ElementMatrix* CreateKMatrixDiagnosticVert(void);
    256256                ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
    257257                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
    258258                ElementMatrix* CreateJacobianDiagnosticHoriz(void);
    259                 ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void);
    260                 ElementMatrix* CreateJacobianDiagnosticPattyn(void);
    261                 ElementMatrix* CreateJacobianDiagnosticStokes(void);
     259                ElementMatrix* CreateJacobianDiagnosticSSA2d(void);
     260                ElementMatrix* CreateJacobianDiagnosticHO(void);
     261                ElementMatrix* CreateJacobianDiagnosticFS(void);
    262262                void           InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solutiong);
    263263                void           InputUpdateFromSolutionDiagnosticMacAyeal( IssmDouble* solutiong);
    264                 void           InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong);
    265                 void           InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong);
     264                void           InputUpdateFromSolutionDiagnosticMacAyealHO( IssmDouble* solutiong);
     265                void           InputUpdateFromSolutionDiagnosticMacAyealFS( IssmDouble* solutiong);
    266266                void           InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong);
    267                 void           InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong);
    268                 void           InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong);
    269                 void           InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solutiong);
     267                void           InputUpdateFromSolutionDiagnosticHO( IssmDouble* solutiong);
     268                void           InputUpdateFromSolutionDiagnosticHOFS( IssmDouble* solutiong);
     269                void           InputUpdateFromSolutionDiagnosticSIA( IssmDouble* solutiong);
    270270                void           InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong);
    271                 void           InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong);
     271                void           InputUpdateFromSolutionDiagnosticFS( IssmDouble* solutiong);
    272272                void             GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong);
    273                 void             GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong);
    274                 void             GetSolutionFromInputsDiagnosticStokes(Vector<IssmDouble>* solutiong);
     273                void             GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solutiong);
     274                void             GetSolutionFromInputsDiagnosticFS(Vector<IssmDouble>* solutiong);
    275275                void             GetSolutionFromInputsDiagnosticVert(Vector<IssmDouble>* solutiong);
    276                 ElementVector* CreatePVectorCouplingMacAyealStokes(void);
    277                 ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void);
    278                 ElementVector* CreatePVectorCouplingMacAyealStokesFriction(void);
    279                 ElementVector* CreatePVectorCouplingPattynStokes(void);
    280                 ElementVector* CreatePVectorCouplingPattynStokesViscous(void);
    281                 ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
     276                ElementVector* CreatePVectorCouplingMacAyealFS(void);
     277                ElementVector* CreatePVectorCouplingMacAyealFSViscous(void);
     278                ElementVector* CreatePVectorCouplingMacAyealFSFriction(void);
     279                ElementVector* CreatePVectorCouplingHOFS(void);
     280                ElementVector* CreatePVectorCouplingHOFSViscous(void);
     281                ElementVector* CreatePVectorCouplingHOFSFriction(void);
    282282                ElementVector* CreatePVectorDiagnosticHoriz(void);
    283                 ElementVector* CreatePVectorDiagnosticHutter(void);
     283                ElementVector* CreatePVectorDiagnosticSIA(void);
    284284                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    285                 ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
    286                 ElementVector* CreatePVectorDiagnosticMacAyealStokes(void);
     285                ElementVector* CreatePVectorDiagnosticMacAyealHO(void);
     286                ElementVector* CreatePVectorDiagnosticMacAyealFS(void);
    287287                ElementVector* CreatePVectorDiagnosticL1L2(void);
    288                 ElementVector* CreatePVectorDiagnosticPattyn(void);
    289                 ElementVector* CreatePVectorDiagnosticPattynDrivingStress(void);
    290                 ElementVector* CreatePVectorDiagnosticPattynFront(void);
    291                 ElementVector* CreatePVectorDiagnosticPattynStokes(void);
    292                 ElementVector* CreatePVectorDiagnosticStokes(void);
    293                 ElementVector* CreatePVectorDiagnosticStokesFront(void);
    294                 ElementVector* CreatePVectorDiagnosticStokesViscous(void);
    295                 ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void);
    296                 ElementVector* CreatePVectorDiagnosticStokesShelf(void);
     288                ElementVector* CreatePVectorDiagnosticHO(void);
     289                ElementVector* CreatePVectorDiagnosticHODrivingStress(void);
     290                ElementVector* CreatePVectorDiagnosticHOFront(void);
     291                ElementVector* CreatePVectorDiagnosticHOFS(void);
     292                ElementVector* CreatePVectorDiagnosticFS(void);
     293                ElementVector* CreatePVectorDiagnosticFSFront(void);
     294                ElementVector* CreatePVectorDiagnosticFSViscous(void);
     295                ElementVector* CreatePVectorDiagnosticFSGLSViscous(void);
     296                ElementVector* CreatePVectorDiagnosticFSShelf(void);
    297297                ElementVector* CreatePVectorDiagnosticVert(void);
    298298                ElementVector* CreatePVectorDiagnosticVertVolume(void);
    … …  
    304304                ElementVector* CreatePVectorAdjointHoriz(void);
    305305                ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void);
    306                 ElementMatrix* CreateKMatrixAdjointPattyn(void);
    307                 ElementMatrix* CreateKMatrixAdjointStokes(void);
     306                ElementMatrix* CreateKMatrixAdjointHO(void);
     307                ElementMatrix* CreateKMatrixAdjointFS(void);
    308308                ElementVector* CreatePVectorAdjointMacAyeal(void);
    309                 ElementVector* CreatePVectorAdjointPattyn(void);
    310                 ElementVector* CreatePVectorAdjointStokes(void);
     309                ElementVector* CreatePVectorAdjointHO(void);
     310                ElementVector* CreatePVectorAdjointFS(void);
    311311                void    InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
    312                 void    InputUpdateFromSolutionAdjointStokes( IssmDouble* solutiong);
     312                void    InputUpdateFromSolutionAdjointFS( IssmDouble* solutiong);
    313313                #endif
    314314

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

    r15511 r15564  
    5454
    5555/*Reference Element numerics*/
    56 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{*/
    57 void PentaRef::GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     56/*FUNCTION PentaRef::GetBMacAyealHO {{{*/
     57void PentaRef::GetBMacAyealHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    5858        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    5959         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    8585}
    8686/*}}}*/
    87 /*FUNCTION PentaRef::GetBMacAyealStokes{{{*/
    88 void PentaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     87/*FUNCTION PentaRef::GetBMacAyealFS{{{*/
     88void PentaRef::GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    8989        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    9090         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    131131}
    132132/*}}}*/
    133 /*FUNCTION PentaRef::GetBPattyn {{{*/
    134 void PentaRef::GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     133/*FUNCTION PentaRef::GetBHO {{{*/
     134void PentaRef::GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    135135        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    136136         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    171171}
    172172/*}}}*/
    173 /*FUNCTION PentaRef::GetBprimePattyn {{{*/
    174 void PentaRef::GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
     173/*FUNCTION PentaRef::GetBprimeHO {{{*/
     174void PentaRef::GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
    175175        /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    176176         * For node i, Bi can be expressed in the actual coordinate system
    … …  
    209209}
    210210/*}}}*/
    211 /*FUNCTION PentaRef::GetBprimeMacAyealStokes{{{*/
    212 void PentaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss){
     211/*FUNCTION PentaRef::GetBprimeMacAyealFS{{{*/
     212void PentaRef::GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss){
    213213        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
    214214         * For node i, Bprimei can be expressed in the actual coordinate system
    … …  
    249249}
    250250/*}}}*/
    251 /*FUNCTION PentaRef::GetBStokes {{{*/
    252 void PentaRef::GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     251/*FUNCTION PentaRef::GetBFS {{{*/
     252void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    253253
    254254        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
    … …  
    316316}
    317317/*}}}*/
    318 /*FUNCTION PentaRef::GetBStokesGLS {{{*/
    319 void PentaRef::GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
     318/*FUNCTION PentaRef::GetBFSGLS {{{*/
     319void PentaRef::GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
    320320
    321321        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4.
    … …  
    382382}
    383383/*}}}*/
    384 /*FUNCTION PentaRef::GetBprimeStokes {{{*/
    385 void PentaRef::GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
     384/*FUNCTION PentaRef::GetBprimeFS {{{*/
     385void PentaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
    386386        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    387387         *      For node i, Bi' can be expressed in the actual coordinate system
    … …  
    449449}
    450450/*}}}*/
    451 /*FUNCTION PentaRef::GetBprimeStokesGLS {{{*/
    452 void PentaRef::GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
     451/*FUNCTION PentaRef::GetBprimeFSGLS {{{*/
     452void PentaRef::GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
    453453        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    454454         *      For node i, Bi' can be expressed in the actual coordinate system
    … …  
    618618}
    619619/*}}}*/
    620 /*FUNCTION PentaRef::GetBPattynFriction{{{*/
    621 void PentaRef::GetBPattynFriction(IssmDouble* B, GaussPenta* gauss){
     620/*FUNCTION PentaRef::GetBHOFriction{{{*/
     621void PentaRef::GetBHOFriction(IssmDouble* B, GaussPenta* gauss){
    622622        /*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2x2.
    623623         ** For node i, Bi can be expressed in the actual coordinate system
    … …  
    642642}
    643643/*}}}*/
    644 /*FUNCTION PentaRef::GetLStokes{{{*/
    645 void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){
     644/*FUNCTION PentaRef::GetLFS{{{*/
     645void PentaRef::GetLFS(IssmDouble* LFS, GaussPenta* gauss){
    646646        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    647647         * For node i, Li can be expressed in the actual coordinate system
    … …  
    655655
    656656        const int num_dof=4;
    657         IssmDouble l1l2l3[NUMNODESP1_2d];
    658 
    659         /*Get l1l2l3 in actual coordinate system: */
    660         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    661         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    662         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    663 
    664         /*Build LStokes: */
     657        IssmDouble L1L2l3[NUMNODESP1_2d];
     658
     659        /*Get L1L2l3 in actual coordinate system: */
     660        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     661        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     662        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     663
     664        /*Build LFS: */
    665665        for(int i=0;i<3;i++){
    666                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
    667                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    668                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    669                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    670 
    671                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    672                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
    673                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    674                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    675         }
    676 }
    677 /*}}}*/
    678 /*FUNCTION PentaRef::GetLprimeStokes {{{*/
    679 void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
     666                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
     667                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     668                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     669                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     670
     671                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     672                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
     673                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     674                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
     675        }
     676}
     677/*}}}*/
     678/*FUNCTION PentaRef::GetLprimeFS {{{*/
     679void PentaRef::GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
    680680        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    681681         * For node i, Lpi can be expressed in the actual coordinate system
    … …  
    715715        int num_dof=4;
    716716
    717         IssmDouble l1l2l3[NUMNODESP1_2d];
     717        IssmDouble L1L2l3[NUMNODESP1_2d];
    718718        IssmDouble dbasis[3][NUMNODESP1];
    719719
    720         /*Get l1l2l3 in actual coordinate system: */
    721         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    722         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    723         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     720        /*Get L1L2l3 in actual coordinate system: */
     721        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     722        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     723        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    724724
    725725        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    726726
    727         /*Build LprimeStokes: */
     727        /*Build LprimeFS: */
    728728        for(int i=0;i<3;i++){
    729                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = l1l2l3[i];
    730                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
    731                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
    732                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
    733                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
    734                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = l1l2l3[i];
    735                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
    736                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
    737                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = l1l2l3[i];
    738                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
    739                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
    740                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
    741                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
    742                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = l1l2l3[i];
    743                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
    744                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
    745                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
    746                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
    747                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = l1l2l3[i];
    748                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
    749                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
    750                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
    751                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = l1l2l3[i];
    752                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
    753                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
    754                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
    755                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
    756                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
    757                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
    758                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
    759                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
    760                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
    761                 LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
    762                 LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
    763                 LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
    764                 LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
    765                 LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
    766                 LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
    767                 LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
    768                 LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
    769                 LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
    770                 LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
    771                 LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
    772                 LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
    773                 LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
    774                 LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
    775                 LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
    776                 LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = l1l2l3[i];
    777                 LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
    778                 LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
    779                 LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
    780                 LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = l1l2l3[i];
    781                 LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
    782                 LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
    783                 LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
    784                 LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = l1l2l3[i];
    785         }
    786 }
    787 /*}}}*/
    788 /*FUNCTION PentaRef::GetLMacAyealStokes {{{*/
    789 void PentaRef::GetLMacAyealStokes(IssmDouble* LStokes, GaussPenta* gauss){
     729                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = L1L2l3[i];
     730                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
     731                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
     732                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
     733                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
     734                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = L1L2l3[i];
     735                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
     736                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
     737                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = L1L2l3[i];
     738                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
     739                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
     740                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
     741                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
     742                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = L1L2l3[i];
     743                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
     744                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
     745                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
     746                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
     747                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = L1L2l3[i];
     748                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
     749                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
     750                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
     751                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = L1L2l3[i];
     752                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
     753                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
     754                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
     755                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
     756                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
     757                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
     758                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
     759                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
     760                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
     761                LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
     762                LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
     763                LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
     764                LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
     765                LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
     766                LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
     767                LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
     768                LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
     769                LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
     770                LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
     771                LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
     772                LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
     773                LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
     774                LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
     775                LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
     776                LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = L1L2l3[i];
     777                LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
     778                LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
     779                LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
     780                LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = L1L2l3[i];
     781                LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
     782                LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
     783                LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
     784                LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = L1L2l3[i];
     785        }
     786}
     787/*}}}*/
     788/*FUNCTION PentaRef::GetLMacAyealFS {{{*/
     789void PentaRef::GetLMacAyealFS(IssmDouble* LFS, GaussPenta* gauss){
    790790        /*
    791791         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    … …  
    804804
    805805        int num_dof=2;
    806         IssmDouble l1l2l3[NUMNODESP1_2d];
    807 
    808         /*Get l1l2l3 in actual coordinate system: */
    809         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    810         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    811         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    812 
    813         /*Build LStokes: */
     806        IssmDouble L1L2l3[NUMNODESP1_2d];
     807
     808        /*Get L1L2l3 in actual coordinate system: */
     809        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     810        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     811        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     812
     813        /*Build LFS: */
    814814        for(int i=0;i<3;i++){
    815                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
    816                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
    817                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
    818                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
    819                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
    820                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
    821                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
    822                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
    823                 LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = l1l2l3[i];
    824                 LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
    825                 LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
    826                 LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = l1l2l3[i];
    827                 LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = l1l2l3[i];
    828                 LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
    829                 LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
    830                 LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = l1l2l3[i];
    831         }
    832 }
    833 /*}}}*/
    834 /*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
    835 void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
     815                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
     816                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
     817                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
     818                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
     819                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
     820                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
     821                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
     822                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
     823                LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = L1L2l3[i];
     824                LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
     825                LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
     826                LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = L1L2l3[i];
     827                LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = L1L2l3[i];
     828                LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
     829                LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
     830                LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = L1L2l3[i];
     831        }
     832}
     833/*}}}*/
     834/*FUNCTION PentaRef::GetLprimeMacAyealFS {{{*/
     835void PentaRef::GetLprimeMacAyealFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
    836836        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    837837         * For node i, Lpi can be expressed in the actual coordinate system
    … …  
    848848         */
    849849        int num_dof=4;
    850         IssmDouble l1l2l3[NUMNODESP1_2d];
     850        IssmDouble L1L2l3[NUMNODESP1_2d];
    851851        IssmDouble dbasis[3][NUMNODESP1];
    852852
    853         /*Get l1l2l3 in actual coordinate system: */
    854         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    855         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    856         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     853        /*Get L1L2l3 in actual coordinate system: */
     854        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     855        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     856        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    857857
    858858        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    859859
    860         /*Build LprimeStokes: */
     860        /*Build LprimeFS: */
    861861        for(int i=0;i<3;i++){
    862                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
    863                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    864                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    865                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    866                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    867                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
    868                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    869                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    870                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
    871                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
    872                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
    873                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
    874                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
    875                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
    876                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
    877                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
    878                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
    879                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
    880                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
    881                 LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
    882                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
    883                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
    884                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
    885                 LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
    886                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
    887                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
    888                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
    889                 LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = l1l2l3[i];
    890                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
    891                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
    892                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
    893                 LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = l1l2l3[i];
    894         }
    895 }
    896 /*}}}*/
    897 /*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
    898 void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){
     862                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
     863                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     864                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     865                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     866                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     867                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
     868                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     869                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
     870                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
     871                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     872                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
     873                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
     874                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     875                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
     876                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
     877                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
     878                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
     879                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
     880                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
     881                LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
     882                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
     883                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
     884                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
     885                LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
     886                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
     887                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
     888                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
     889                LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = L1L2l3[i];
     890                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
     891                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
     892                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
     893                LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = L1L2l3[i];
     894        }
     895}
     896/*}}}*/
     897/*FUNCTION PentaRef::GetLFSMacAyeal {{{*/
     898void PentaRef::GetLFSMacAyeal(IssmDouble* LFS, GaussPenta* gauss){
    899899        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    900900         * For node i, Li can be expressed in the actual coordinate system
    … …  
    908908
    909909        int num_dof=4;
    910         IssmDouble l1l2l3[NUMNODESP1_2d];
    911 
    912         /*Get l1l2l3 in actual coordinate system: */
    913         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    914         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    915         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    916 
    917         /*Build LStokes: */
     910        IssmDouble L1L2l3[NUMNODESP1_2d];
     911
     912        /*Get L1L2l3 in actual coordinate system: */
     913        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     914        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     915        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     916
     917        /*Build LFS: */
    918918        for(int i=0;i<3;i++){
    919                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
    920                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    921                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
    922                 LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
    923                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    924                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
    925                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
    926                 LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
    927                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
    928                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
    929                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
    930                 LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
    931                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
    932                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
    933                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
    934                 LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
    935         }
    936 }
    937 /*}}}*/
    938 /*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
    939 void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
     919                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
     920                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     921                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
     922                LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
     923                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     924                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
     925                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
     926                LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
     927                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
     928                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     929                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
     930                LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
     931                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     932                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
     933                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
     934                LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
     935        }
     936}
     937/*}}}*/
     938/*FUNCTION PentaRef::GetLprimeFSMacAyeal {{{*/
     939void PentaRef::GetLprimeFSMacAyeal(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
    940940        /* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    941941         * For node i, Lpi can be expressed in the actual coordinate system
    … …  
    948948         */
    949949        int num_dof=2;
    950         IssmDouble l1l2l3[NUMNODESP1_2d];
     950        IssmDouble L1L2l3[NUMNODESP1_2d];
    951951        IssmDouble dbasis[3][NUMNODESP1];
    952952
    953         /*Get l1l2l3 in actual coordinate system: */
    954         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    955         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    956         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
     953        /*Get L1L2l3 in actual coordinate system: */
     954        L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     955        L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     956        L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    957957        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    958958
    959         /*Build LprimeStokes: */
     959        /*Build LprimeFS: */
    960960        for(int i=0;i<3;i++){
    961                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
    962                 LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
    963                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
    964                 LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
    965                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
    966                 LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
    967                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
    968                 LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
     961                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
     962                LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
     963                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
     964                LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
     965                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
     966                LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
     967                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
     968                LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
    969969        }
    970970}

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

    r15425 r15564  
    3838                void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussPenta* gauss);
    3939                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussPenta* gauss);
    40                 void GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    41                 void GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    42                 void GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    43                 void GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    44                 void GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    45                 void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss);
    46                 void GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    47                 void GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
    48                 void GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
     40                void GetBMacAyealHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     41                void GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     42                void GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     43                void GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     44                void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     45                void GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss);
     46                void GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
     47                void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
     48                void GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
    4949                void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    5050                void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss);
    … …  
    5252                void GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
    5353                void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss);
    54                 void GetBPattynFriction(IssmDouble* L, GaussPenta* gauss);
    55                 void GetLStokes(IssmDouble* LStokes, GaussPenta* gauss);
    56                 void GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss);
    57                 void GetLMacAyealStokes(IssmDouble* LMacAyealStokes, GaussPenta* gauss);
    58                 void GetLprimeMacAyealStokes(IssmDouble* LprimeMacAyealStokes, IssmDouble* xyz_list, GaussPenta* gauss);
    59                 void GetLStokesMacAyeal(IssmDouble* LStokesMacAyeal, GaussPenta* gauss);
    60                 void GetLprimeStokesMacAyeal(IssmDouble* LprimeStokesMacAyeal, IssmDouble* xyz_list, GaussPenta* gauss);
     54                void GetBHOFriction(IssmDouble* L, GaussPenta* gauss);
     55                void GetLFS(IssmDouble* LFS, GaussPenta* gauss);
     56                void GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss);
     57                void GetLMacAyealFS(IssmDouble* LMacAyealFS, GaussPenta* gauss);
     58                void GetLprimeMacAyealFS(IssmDouble* LprimeMacAyealFS, IssmDouble* xyz_list, GaussPenta* gauss);
     59                void GetLFSMacAyeal(IssmDouble* LFSMacAyeal, GaussPenta* gauss);
     60                void GetLprimeFSMacAyeal(IssmDouble* LprimeFSMacAyeal, IssmDouble* xyz_list, GaussPenta* gauss);
    6161                void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, GaussPenta* gauss);
    6262                void GetInputValue(IssmDouble* pvalue,IssmDouble* plist,GaussTria* gauss){_error_("only PentaGauss are supported");};

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

    r15544 r15564  
    193193                        Ke=CreateKMatrixDiagnosticMacAyeal();
    194194                        break;
    195                 case DiagnosticHutterAnalysisEnum:
    196                         Ke=CreateKMatrixDiagnosticHutter();
     195                case DiagnosticSIAAnalysisEnum:
     196                        Ke=CreateKMatrixDiagnosticSIA();
    197197                        break;
    198198                 #endif
    … …  
    306306                        pe=CreatePVectorDiagnosticMacAyeal();
    307307                        break;
    308                 case DiagnosticHutterAnalysisEnum:
    309                         pe=CreatePVectorDiagnosticHutter();
     308                case DiagnosticSIAAnalysisEnum:
     309                        pe=CreatePVectorDiagnosticSIA();
    310310                        break;
    311311                #endif
    … …  
    423423#ifdef _HAVE_DIAGNOSTIC_
    424424                case DiagnosticHorizAnalysisEnum:
    425                         Ke=CreateJacobianDiagnosticMacayeal();
     425                        Ke=CreateJacobianDiagnosticSSA();
    426426                        break;
    427427#endif
    … …  
    12011201                GetSolutionFromInputsDiagnosticHoriz(solution);
    12021202                break;
    1203         case DiagnosticHutterAnalysisEnum:
    1204                 GetSolutionFromInputsDiagnosticHutter(solution);
     1203        case DiagnosticSIAAnalysisEnum:
     1204                GetSolutionFromInputsDiagnosticSIA(solution);
    12051205                break;
    12061206        #endif
    … …  
    15221522                        InputUpdateFromSolutionDiagnosticHoriz(solution);
    15231523                        break;
    1524                 case DiagnosticHutterAnalysisEnum:
     1524                case DiagnosticSIAAnalysisEnum:
    15251525                        InputUpdateFromSolutionDiagnosticHoriz(solution);
    15261526                        break;
    … …  
    23432343        iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
    23442344        iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
    2345         iomodel->Constant(&fe_ssa,FlowequationFeSsaEnum);
     2345        iomodel->Constant(&fe_ssa,FlowequationFeSSAEnum);
    23462346        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    23472347
    … …  
    30553055}
    30563056/*}}}*/
    3057 /*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{*/
    3058 ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
     3057/*FUNCTION Tria::CreateKMatrixDiagnosticSIA{{{*/
     3058ElementMatrix* Tria::CreateKMatrixDiagnosticSIA(void){
    30593059
    30603060        /*Intermediaries*/
    … …  
    32243224}
    32253225/*}}}*/
    3226 /*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{*/
    3227 ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
     3226/*FUNCTION Tria::CreatePVectorDiagnosticSIA{{{*/
     3227ElementVector* Tria::CreatePVectorDiagnosticSIA(void){
    32283228
    32293229        /*Intermediaries */
    … …  
    32743274}
    32753275/*}}}*/
    3276 /*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{*/
    3277 ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
     3276/*FUNCTION Tria::CreateJacobianDiagnosticSSA{{{*/
     3277ElementMatrix* Tria::CreateJacobianDiagnosticSSA(void){
    32783278
    32793279        /*Intermediaries */
    … …  
    33823382}
    33833383/*}}}*/
    3384 /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{*/
    3385 void  Tria::GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution){
     3384/*FUNCTION Tria::GetSolutionFromInputsDiagnosticSIA{{{*/
     3385void  Tria::GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution){
    33863386
    33873387        const int    numdof=NDOF2*NUMVERTICES;
    … …  
    34943494}
    34953495/*}}}*/
    3496 /*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{*/
    3497 void  Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){
     3496/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticSIA {{{*/
     3497void  Tria::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
    34983498
    34993499        int        i;
    … …  
    51375137}
    51385138/*}}}*/
    5139 /*FUNCTION Tria::CreatePVectorAdjointStokes{{{*/
    5140 ElementVector* Tria::CreatePVectorAdjointStokes(void){
     5139/*FUNCTION Tria::CreatePVectorAdjointFS{{{*/
     5140ElementVector* Tria::CreatePVectorAdjointFS(void){
    51415141
    51425142        /*Intermediaries */
    … …  
    51565156
    51575157        /*Initialize Element vector*/
    5158         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
     5158        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
    51595159
    51605160        /*Retrieve all inputs and parameters*/

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

    r15517 r15564  
    147147                void   GradjZMacAyeal(Vector<IssmDouble>* gradient,int control_index);
    148148                void   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
    149                 void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
     149                void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
    150150                void   GradjDragGradient(Vector<IssmDouble>* gradient,int weight_index,int control_index);
    151151                void   GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
    … …  
    223223                ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
    224224                ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
    225                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
     225                ElementMatrix* CreateKMatrixDiagnosticSIA(void);
    226226                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    227227                ElementVector* CreatePVectorDiagnosticMacAyealDrivingStress(void);
    228228                ElementVector* CreatePVectorDiagnosticMacAyealFront(void);
    229                 ElementVector* CreatePVectorDiagnosticHutter(void);
    230                 ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
     229                ElementVector* CreatePVectorDiagnosticSIA(void);
     230                ElementMatrix* CreateJacobianDiagnosticSSA(void);
    231231                void      GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solution);
    232                 void      GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution);
     232                void      GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution);
    233233                void      InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solution);
    234                 void      InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solution);
     234                void      InputUpdateFromSolutionDiagnosticSIA( IssmDouble* solution);
    235235                #endif
    236236
    … …  
    239239                ElementMatrix* CreateKMatrixAdjointMacAyeal(void);
    240240                ElementVector* CreatePVectorAdjointHoriz(void);
    241                 ElementVector* CreatePVectorAdjointStokes(void);
     241                ElementVector* CreatePVectorAdjointFS(void);
    242242                ElementVector* CreatePVectorAdjointBalancethickness(void);
    243243                void      InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);

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

    r15496 r15564  
    114114}
    115115/*}}}*/
    116 /*FUNCTION TriaRef::GetBMacAyealStokes {{{*/
    117 void TriaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
     116/*FUNCTION TriaRef::GetBMacAyealFS {{{*/
     117void TriaRef::GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
    118118
    119119        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    … …  
    265265}
    266266/*}}}*/
    267 /*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/
    268 void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
     267/*FUNCTION TriaRef::GetBprimeMacAyealFS {{{*/
     268void TriaRef::GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
    269269        /*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
    270270         * For node i, Bprimei can be expressed in the actual coordinate system

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

    r15326 r15564  
    2424                /*Numerics*/
    2525                void GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
    26                 void GetBMacAyealStokes(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
     26                void GetBMacAyealFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
    2727                void GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    28                 void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
     28                void GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    2929                void GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
    3030                void GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
Note: See TracChangeset for help on using the changeset viewer.