Changeset 15417


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

CHG: removed PentaP1Input and created PentaInput for quadratic elements

Location:
issm/trunk-jpl/src/c
Files:
14 edited
2 moved

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r15372 r15417  
    533533                                     ./classes/ElementResults/PentaP1ElementResult.h\
    534534                                     ./classes/ElementResults/PentaP1ElementResult.cpp\
    535                                      ./classes/Inputs/PentaP1Input.h\
    536                                      ./classes/Inputs/PentaP1Input.cpp\
     535                                     ./classes/Inputs/PentaInput.h\
     536                                     ./classes/Inputs/PentaInput.cpp\
    537537                                     ./classes/Elements/Penta.h\
    538538                                     ./classes/Elements/Penta.cpp\

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

    r15407 r15417  
    182182        if (!IsOnBed() || IsFloating()){
    183183                //empty friction:
    184                 this->inputs->AddInput(new PentaP1Input(BasalFrictionEnum,&basalfriction[0]));
     184                this->inputs->AddInput(new PentaInput(BasalFrictionEnum,&basalfriction[0],P1Enum));
    185185                return;
    186186        }
    … …  
    209209
    210210        /*Create PentaVertex input, which will hold the basal friction:*/
    211         this->inputs->AddInput(new PentaP1Input(BasalFrictionEnum,&basalfriction[0]));
     211        this->inputs->AddInput(new PentaInput(BasalFrictionEnum,&basalfriction[0],P1Enum));
    212212
    213213        /*Clean up and return*/
    … …  
    354354
    355355        /*Add Stress tensor components into inputs*/
    356         this->inputs->AddInput(new PentaP1Input(StressTensorxxEnum,&sigma_xx[0]));
    357         this->inputs->AddInput(new PentaP1Input(StressTensorxyEnum,&sigma_xy[0]));
    358         this->inputs->AddInput(new PentaP1Input(StressTensorxzEnum,&sigma_xz[0]));
    359         this->inputs->AddInput(new PentaP1Input(StressTensoryyEnum,&sigma_yy[0]));
    360         this->inputs->AddInput(new PentaP1Input(StressTensoryzEnum,&sigma_yz[0]));
    361         this->inputs->AddInput(new PentaP1Input(StressTensorzzEnum,&sigma_zz[0]));
     356        this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&sigma_xx[0],P1Enum));
     357        this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&sigma_xy[0],P1Enum));
     358        this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&sigma_xz[0],P1Enum));
     359        this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&sigma_yy[0],P1Enum));
     360        this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
     361        this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum));
    362362
    363363        /*Clean up and return*/
    … …  
    763763        for (int imonth=0;imonth<12;imonth++) {
    764764                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    765                 PentaP1Input* newmonthinput1 = new PentaP1Input(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0]);
     765                PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    766766                NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    767767
    768768                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    769                 PentaP1Input* newmonthinput2 = new PentaP1Input(SurfaceforcingsPrecipitationEnum,&tmp[0]);
     769                PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    770770                NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    771771        }
    … …  
    14201420
    14211421                        /*create static input: */
    1422                         this->inputs->AddInput(new PentaP1Input(vector_enum,nodeinputs));
     1422                        this->inputs->AddInput(new PentaInput(vector_enum,nodeinputs,P1Enum));
    14231423                }
    14241424                else if(M==numberofvertices+1){
    … …  
    14361436
    14371437                                if(t==0)transientinput=new TransientInput(vector_enum);
    1438                                 transientinput->AddTimeInput(new PentaP1Input(vector_enum,nodeinputs),time);
     1438                                transientinput->AddTimeInput(new PentaInput(vector_enum,nodeinputs,P1Enum),time);
    14391439                        }
    14401440                        this->inputs->AddInput(transientinput);
    … …  
    14901490
    14911491        /*OK, we are on bed. Initialize global inputs as 0*/
    1492         total_thickness_input =new PentaP1Input(ThicknessEnum,zeros_list);
     1492        total_thickness_input =new PentaInput(ThicknessEnum,zeros_list,P1Enum);
    14931493
    14941494        /*Now follow all the upper element from the base to the surface to integrate the input*/
    … …  
    15081508                /*If first time, initialize total_integrated_input*/
    15091509                if (step==0){
    1510                         if (original_input->ObjectEnum()==PentaP1InputEnum)
    1511                          total_integrated_input=new PentaP1Input(average_enum_type,zeros_list);
     1510                        if (original_input->ObjectEnum()==PentaInputEnum)
     1511                         total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum);
    15121512                        else if (original_input->ObjectEnum()==ControlInputEnum)
    1513                          total_integrated_input=new PentaP1Input(average_enum_type,zeros_list);
     1513                         total_integrated_input=new PentaInput(average_enum_type,zeros_list,P1Enum);
    15141514                        else if (original_input->ObjectEnum()==DoubleInputEnum)
    15151515                         total_integrated_input=new DoubleInput(average_enum_type,0.0);
    … …  
    15241524                        Helem_list[i+3]=Helem_list[i];
    15251525                }
    1526                 element_thickness_input=new PentaP1Input(ThicknessEnum,Helem_list);
     1526                element_thickness_input=new PentaInput(ThicknessEnum,Helem_list,P1Enum);
    15271527
    15281528                /*Step3: Vertically integrate A COPY of the original*/
    … …  
    17521752                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    17531753                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1754                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1754                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17551755                                        }
    17561756                                        break;
    … …  
    17601760                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    17611761                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1762                                                 this->inputs->AddInput(new ControlInput(VxEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1762                                                this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17631763                                        }
    17641764                                        break;
    … …  
    17681768                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    17691769                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1770                                                 this->inputs->AddInput(new ControlInput(VyEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1770                                                this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17711771                                        }
    17721772                                        break;
    … …  
    17761776                                                for(j=0;j<6;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    17771777                                                for(j=0;j<6;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    1778                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1778                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17791779                                        }
    17801780                                        break;
    … …  
    18301830                for(i=0;i<num_cm_responses;i++){
    18311831                        for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_cm_responses+i];
    1832                         datasetinput->inputs->AddObject(new PentaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs));
     1832                        datasetinput->inputs->AddObject(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum));
    18331833                }
    18341834
    … …  
    19831983        for(;;){
    19841984                /*Add input to the element: */
    1985                 penta->inputs->AddInput(new PentaP1Input(ThicknessEnum,newthickness));
    1986                 penta->inputs->AddInput(new PentaP1Input(SurfaceEnum,newsurface));
    1987                 penta->inputs->AddInput(new PentaP1Input(BedEnum,newbed));
     1985                penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
     1986                penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
     1987                penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
    19881988
    19891989                /*Stop if we have reached the surface*/
    … …  
    20162016
    20172017        /*Add input to the element: */
    2018         this->inputs->AddInput(new PentaP1Input(enum_type,values));
     2018        this->inputs->AddInput(new PentaInput(enum_type,values,P1Enum));
    20192019
    20202020        /*Free ressources:*/
    … …  
    20492049        for(;;){
    20502050                /*Add input to the element: */
    2051                 penta->inputs->AddInput(new PentaP1Input(enum_type,values));
     2051                penta->inputs->AddInput(new PentaInput(enum_type,values,P1Enum));
    20522052
    20532053                /*Stop if we have reached the surface*/
    … …  
    20792079                        /*update input*/
    20802080                        if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    2081                                 material->inputs->AddInput(new PentaP1Input(name,values));
     2081                                material->inputs->AddInput(new PentaInput(name,values,P1Enum));
    20822082                        }
    20832083                        else{
    2084                                 this->inputs->AddInput(new PentaP1Input(name,values));
     2084                                this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    20852085                        }
    20862086                        return;
    … …  
    20922092                        /*update input*/
    20932093                        if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    2094                                 material->inputs->AddInput(new PentaP1Input(name,values));
     2094                                material->inputs->AddInput(new PentaInput(name,values,P1Enum));
    20952095                        }
    20962096                        else{
    2097                                 this->inputs->AddInput(new PentaP1Input(name,values));
     2097                                this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    20982098                        }
    20992099                        return;
    … …  
    21092109                        }
    21102110                        /*Add input to the element: */
    2111                         this->inputs->AddInput(new PentaP1Input(name,values));
     2111                        this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    21122112
    21132113                        /*Free ressources:*/
    … …  
    21212121                        }
    21222122                        /*Add input to the element: */
    2123                         this->inputs->AddInput(new PentaP1Input(name,values));
     2123                        this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    21242124                       
    21252125                        /*Free ressources:*/
    … …  
    24642464
    24652465   /*Update inputs*/   
    2466    this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
     2466   this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    24672467   //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    24682468   this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
    … …  
    28002800        }  //end of the loop over the vertices
    28012801          /*Update inputs*/
    2802           this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0]));
     2802          this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&smb[0],P1Enum));
    28032803}
    28042804/*}}}*/
    … …  
    29742974                        if(!iomodel->Data(VxEnum)){
    29752975                                for(i=0;i<6;i++)nodeinputs[i]=0;
    2976                                 this->inputs->AddInput(new PentaP1Input(VxEnum,nodeinputs));
    2977                                 if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVxEnum,nodeinputs));
     2976                                this->inputs->AddInput(new PentaInput(VxEnum,nodeinputs,P1Enum));
     2977                                if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVxEnum,nodeinputs,P1Enum));
    29782978                        }
    29792979                        if(!iomodel->Data(VyEnum)){
    29802980                                for(i=0;i<6;i++)nodeinputs[i]=0;
    2981                                 this->inputs->AddInput(new PentaP1Input(VyEnum,nodeinputs));
    2982                                 if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVyEnum,nodeinputs));
     2981                                this->inputs->AddInput(new PentaInput(VyEnum,nodeinputs,P1Enum));
     2982                                if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVyEnum,nodeinputs,P1Enum));
    29832983                        }
    29842984                        if(!iomodel->Data(VzEnum)){
    29852985                                for(i=0;i<6;i++)nodeinputs[i]=0;
    2986                                 this->inputs->AddInput(new PentaP1Input(VzEnum,nodeinputs));
    2987                                 if(dakota_analysis) this->inputs->AddInput(new PentaP1Input(QmuVzEnum,nodeinputs));
     2986                                this->inputs->AddInput(new PentaInput(VzEnum,nodeinputs,P1Enum));
     2987                                if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVzEnum,nodeinputs,P1Enum));
    29882988                        }
    29892989                        if(!iomodel->Data(PressureEnum)){
    29902990                                for(i=0;i<6;i++)nodeinputs[i]=0;
    29912991                                if(dakota_analysis){
    2992                                         this->inputs->AddInput(new PentaP1Input(PressureEnum,nodeinputs));
    2993                                         this->inputs->AddInput(new PentaP1Input(QmuPressureEnum,nodeinputs));
     2992                                        this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
     2993                                        this->inputs->AddInput(new PentaInput(QmuPressureEnum,nodeinputs,P1Enum));
    29942994                                }
    29952995                                if(isstokes){
    2996                                         this->inputs->AddInput(new PentaP1Input(PressureEnum,nodeinputs));
    2997                                         this->inputs->AddInput(new PentaP1Input(PressurePicardEnum,nodeinputs));
     2996                                        this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
     2997                                        this->inputs->AddInput(new PentaInput(PressurePicardEnum,nodeinputs,P1Enum));
    29982998                                }
    29992999                        }
    … …  
    30023002                                if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
    30033003                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
    3004                                         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs));
     3004                                        this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
    30053005                                        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]);
    3006                                         this->inputs->AddInput(new PentaP1Input(VzPattynEnum,nodeinputs));
     3006                                        this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
    30073007                                }
    30083008                                else{
    30093009                                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3010                                         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs));
    3011                                         this->inputs->AddInput(new PentaP1Input(VzPattynEnum,nodeinputs));
     3010                                        this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
     3011                                        this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
    30123012                                }
    30133013                        }
    … …  
    30163016                                if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
    30173017                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
    3018                                         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs));
     3018                                        this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
    30193019                                        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]);
    3020                                         this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,nodeinputs));
     3020                                        this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
    30213021                                }
    30223022                                else{
    30233023                                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3024                                         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,nodeinputs));
    3025                                         this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,nodeinputs));
     3024                                        this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
     3025                                        this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
    30263026                                }
    30273027                        }
    … …  
    30313031                        /*Initialize mesh velocity*/
    30323032                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3033                         this->inputs->AddInput(new PentaP1Input(VxMeshEnum,nodeinputs));
    3034                         this->inputs->AddInput(new PentaP1Input(VyMeshEnum,nodeinputs));
    3035                         this->inputs->AddInput(new PentaP1Input(VzMeshEnum,nodeinputs));
     3033                        this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum));
     3034                        this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum));
     3035                        this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum));
    30363036                        if(dakota_analysis){
    3037                                 this->inputs->AddInput(new PentaP1Input(QmuVxMeshEnum,nodeinputs));
    3038                                 this->inputs->AddInput(new PentaP1Input(QmuVyMeshEnum,nodeinputs));
    3039                                 this->inputs->AddInput(new PentaP1Input(QmuVzMeshEnum,nodeinputs));
     3037                                this->inputs->AddInput(new PentaInput(QmuVxMeshEnum,nodeinputs,P1Enum));
     3038                                this->inputs->AddInput(new PentaInput(QmuVyMeshEnum,nodeinputs,P1Enum));
     3039                                this->inputs->AddInput(new PentaInput(QmuVzMeshEnum,nodeinputs,P1Enum));
    30403040                        }
    30413041                        break;
    … …  
    30443044                        /*Initialize mesh velocity*/
    30453045                        for(i=0;i<6;i++)nodeinputs[i]=0;
    3046                         this->inputs->AddInput(new PentaP1Input(VxMeshEnum,nodeinputs));
    3047                         this->inputs->AddInput(new PentaP1Input(VyMeshEnum,nodeinputs));
    3048                         this->inputs->AddInput(new PentaP1Input(VzMeshEnum,nodeinputs));
     3046                        this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum));
     3047                        this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum));
     3048                        this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum));
    30493049                        if (iomodel->Data(TemperatureEnum) && iomodel->Data(WaterfractionEnum) && iomodel->Data(PressureEnum)) {
    30503050                                for(i=0;i<6;i++){
    … …  
    30563056                                                +latentheat*iomodel->Data(WaterfractionEnum)[penta_vertex_ids[i]-1];
    30573057                                }
    3058                                 this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,nodeinputs));
     3058                                this->inputs->AddInput(new PentaInput(EnthalpyEnum,nodeinputs,P1Enum));
    30593059                        }
    30603060                        else _error_("temperature and waterfraction required for the enthalpy solution");
    … …  
    31143114
    31153115        /*Create PentaVertex input, which will hold the basal friction:*/
    3116         this->inputs->AddInput(new PentaP1Input(ViscousHeatingEnum,&viscousheating[0]));
     3116        this->inputs->AddInput(new PentaInput(ViscousHeatingEnum,&viscousheating[0],P1Enum));
    31173117
    31183118        /*Clean up and return*/
    … …  
    43994399        this->inputs->GetInputValue(&converged,ConvergedEnum);
    44004400        if(converged){
    4401                 this->inputs->AddInput(new PentaP1Input(TemperatureEnum,values));
     4401                this->inputs->AddInput(new PentaInput(TemperatureEnum,values,P1Enum));
    44024402
    44034403                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
    … …  
    44114411                                B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
    44124412                                for(i=0;i<numdof;i++) B[i]=B_average;
    4413                                 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));
     4413                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    44144414                                break;
    44154415                        case ArrheniusEnum:
    … …  
    44194419                                                        material->GetN());
    44204420                                for(i=0;i<numdof;i++) B[i]=B_average;
    4421                                 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));
     4421                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    44224422                                break;
    44234423                        default:
    … …  
    44274427        }
    44284428        else{
    4429                 this->inputs->AddInput(new PentaP1Input(TemperaturePicardEnum,values));
     4429                this->inputs->AddInput(new PentaInput(TemperaturePicardEnum,values,P1Enum));
    44304430        }
    44314431
    … …  
    44754475                }
    44764476
    4477                 this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values));
    4478                 this->inputs->AddInput(new PentaP1Input(WaterfractionEnum,waterfraction));
    4479                 this->inputs->AddInput(new PentaP1Input(TemperatureEnum,temperatures));
     4477                this->inputs->AddInput(new PentaInput(EnthalpyEnum,values,P1Enum));
     4478                this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
     4479                this->inputs->AddInput(new PentaInput(TemperatureEnum,temperatures,P1Enum));
    44804480
    44814481                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
    … …  
    44894489                                B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0);
    44904490                                for(i=0;i<numdof;i++) B[i]=B_average;
    4491                                 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));
     4491                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    44924492                                break;
    44934493                        case ArrheniusEnum:
    … …  
    44974497                                                        material->GetN());
    44984498                                for(i=0;i<numdof;i++) B[i]=B_average;
    4499                                 this->material->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,B));
     4499                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    45004500                                break;
    45014501                        default:
    … …  
    45054505        }
    45064506        else{
    4507                 this->inputs->AddInput(new PentaP1Input(EnthalpyPicardEnum,values));
     4507                this->inputs->AddInput(new PentaInput(EnthalpyPicardEnum,values,P1Enum));
    45084508        }
    45094509
    … …  
    45814581        GradientIndexing(&vertexpidlist[0],control_index);
    45824582        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
    4583         grad_input=new PentaP1Input(GradientEnum,grad_list);
     4583        grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
    45844584        ((ControlInput*)input)->SetGradient(grad_input);
    45854585
    … …  
    52555255
    52565256        /*Add vx and vy as inputs to the tria element: */
    5257         this->inputs->AddInput(new PentaP1Input(AdjointxEnum,lambdax));
    5258         this->inputs->AddInput(new PentaP1Input(AdjointyEnum,lambday));
    5259         this->inputs->AddInput(new PentaP1Input(AdjointzEnum,lambdaz));
    5260         this->inputs->AddInput(new PentaP1Input(AdjointpEnum,lambdap));
     5257        this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
     5258        this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
     5259        this->inputs->AddInput(new PentaInput(AdjointzEnum,lambdaz,P1Enum));
     5260        this->inputs->AddInput(new PentaInput(AdjointpEnum,lambdap,P1Enum));
    52615261
    52625262        /*Free ressources:*/
    … …  
    52925292
    52935293        /*Add vx and vy as inputs to the tria element: */
    5294         this->inputs->AddInput(new PentaP1Input(AdjointxEnum,lambdax));
    5295         this->inputs->AddInput(new PentaP1Input(AdjointyEnum,lambday));
     5294        this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
     5295        this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
    52965296
    52975297        /*Free ressources:*/
    … …  
    55835583                values[i]=vector[vertexpidlist[i]];
    55845584        }
    5585         new_input = new PentaP1Input(control_enum,values);
     5585        new_input = new PentaInput(control_enum,values,P1Enum);
    55865586
    55875587        if(control_enum==MaterialsRheologyBbarEnum){
    … …  
    56145614                case VertexEnum:
    56155615
    5616                         /*New PentaP1Input*/
     5616                        /*New PentaInput*/
    56175617                        IssmDouble values[6];
    56185618
    … …  
    56875687
    56885688                                        /*Add new inputs: */
    5689                                         this->inputs->AddInput(new PentaP1Input(ThicknessEnum,thickness));
    5690                                         this->inputs->AddInput(new PentaP1Input(BedEnum,bed));
    5691                                         this->inputs->AddInput(new PentaP1Input(SurfaceEnum,surface));
     5689                                        this->inputs->AddInput(new PentaInput(ThicknessEnum,thickness,P1Enum));
     5690                                        this->inputs->AddInput(new PentaInput(BedEnum,bed,P1Enum));
     5691                                        this->inputs->AddInput(new PentaInput(SurfaceEnum,surface,P1Enum));
    56925692
    56935693                                        /*}}}*/
    56945694                                        break;
    56955695                                default:
    5696                                         this->inputs->AddInput(new PentaP1Input(name,values));
     5696                                        this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    56975697                        }
    56985698                        break;
    … …  
    57475747
    57485748                                if(t==0) transientinput=new TransientInput(name);
    5749                                 transientinput->AddTimeInput(new PentaP1Input(name,values),time);
     5749                                transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time);
    57505750                                transientinput->Configure(parameters);
    57515751                        }
    … …  
    84398439
    84408440                /*Add vx and vy as inputs to the tria element: */
    8441                 penta->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8442                 penta->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8443                 penta->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8444                 penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8441                penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8442                penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8443                penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8444                penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    84458445
    84468446                /*Stop if we have reached the surface*/
    … …  
    85298529
    85308530        /*Add vx and vy as inputs to the tria element: */
    8531         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8532         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8533         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8534         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8531        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8532        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8533        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8534        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    85358535
    85368536        /*Free ressources:*/
    … …  
    86048604        Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);
    86058605        if (vzmacayeal_input){
    8606                 if (vzmacayeal_input->ObjectEnum()!=PentaP1InputEnum){
     8606                if (vzmacayeal_input->ObjectEnum()!=PentaInputEnum){
    86078607                        _error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzmacayeal_input->ObjectEnum()));
    86088608                }
    … …  
    86278627
    86288628        /*Add vx and vy as inputs to the tria element: */
    8629         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8630         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8631         this->inputs->AddInput(new PentaP1Input(VzEnum,vz));
    8632         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,vzstokes));
    8633         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8634         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8629        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8630        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8631        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
     8632        this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
     8633        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8634        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    86358635
    86368636        /*Free ressources:*/
    … …  
    87048704
    87058705                /*Add vx and vy as inputs to the tria element: */
    8706                 penta->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8707                 penta->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8708                 penta->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8709                 penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8706                penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8707                penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8708                penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8709                penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    87108710
    87118711                /*Stop if we have reached the surface*/
    … …  
    87858785
    87868786        /*Add vx and vy as inputs to the tria element: */
    8787         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8788         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8789         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8790         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8787        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8788        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8789        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8790        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    87918791
    87928792        /*Free ressources:*/
    … …  
    88538853        Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);
    88548854        if (vzpattyn_input){
    8855                 if (vzpattyn_input->ObjectEnum()!=PentaP1InputEnum){
     8855                if (vzpattyn_input->ObjectEnum()!=PentaInputEnum){
    88568856                        _error_("Cannot compute Vel as VzPattyn is of type " << EnumToStringx(vzpattyn_input->ObjectEnum()));
    88578857                }
    … …  
    88768876
    88778877        /*Add vx and vy as inputs to the tria element: */
    8878         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8879         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8880         this->inputs->AddInput(new PentaP1Input(VzEnum,vz));
    8881         this->inputs->AddInput(new PentaP1Input(VzStokesEnum,vzstokes));
    8882         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8883         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8878        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8879        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8880        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
     8881        this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
     8882        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8883        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    88848884
    88858885        /*Free ressources:*/
    … …  
    89428942
    89438943        /*Add vx and vy as inputs to the tria element: */
    8944         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    8945         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    8946         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    8947         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     8944        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     8945        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     8946        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     8947        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    89488948
    89498949        /*Free ressources:*/
    … …  
    89998999                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    90009000                if (vzstokes_input){
    9001                         if (vzstokes_input->ObjectEnum()!=PentaP1InputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
     9001                        if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
    90029002                        GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    90039003                }
    … …  
    90119011                Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
    90129012                if (vzstokes_input){
    9013                         if (vzstokes_input->ObjectEnum()!=PentaP1InputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
     9013                        if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
    90149014                        GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
    90159015                }
    … …  
    90399039        if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){
    90409040                this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    9041                 this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     9041                this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    90429042        }
    90439043        else if(approximation==PattynStokesApproximationEnum){
    9044                 this->inputs->AddInput(new PentaP1Input(VzPattynEnum,vzpattyn));
     9044                this->inputs->AddInput(new PentaInput(VzPattynEnum,vzpattyn,P1Enum));
    90459045        }
    90469046        else if(approximation==MacAyealStokesApproximationEnum){
    9047                 this->inputs->AddInput(new PentaP1Input(VzMacAyealEnum,vzmacayeal));
    9048         }
    9049         this->inputs->AddInput(new PentaP1Input(VzEnum,vz));
    9050         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
     9047                this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzmacayeal,P1Enum));
     9048        }
     9049        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
     9050        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    90519051
    90529052        /*Free ressources:*/
    … …  
    91059105
    91069106        /*Add vx and vy as inputs to the tria element: */
    9107         this->inputs->AddInput(new PentaP1Input(VxEnum,vx));
    9108         this->inputs->AddInput(new PentaP1Input(VyEnum,vy));
    9109         this->inputs->AddInput(new PentaP1Input(VzEnum,vz));
    9110         this->inputs->AddInput(new PentaP1Input(VelEnum,vel));
    9111         this->inputs->AddInput(new PentaP1Input(PressureEnum,pressure));
     9107        this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
     9108        this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
     9109        this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
     9110        this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
     9111        this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    91129112
    91139113        /*Free ressources:*/
    … …  
    93509350        for(;;){
    93519351                /*Add input to the element: */
    9352                 penta->inputs->AddInput(new PentaP1Input(SedimentHeadEnum,values));
    9353                 penta->inputs->AddInput(new PentaP1Input(SedimentHeadResidualEnum,residual));
     9352                penta->inputs->AddInput(new PentaInput(SedimentHeadEnum,values,P1Enum));
     9353                penta->inputs->AddInput(new PentaInput(SedimentHeadResidualEnum,residual,P1Enum));
    93549354
    93559355                /*Stop if we have reached the surface*/
    … …  
    94519451        if(!this->IsFloating() && floatingelement==true){
    94529452                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
    9453                 this->inputs->AddInput(new PentaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
     9453                this->inputs->AddInput(new PentaInput(BasalforcingsMeltingRateEnum,&melting[0],P1Enum));
    94549454        }
    94559455
    94569456        /*Update inputs*/
    9457         this->inputs->AddInput(new PentaP1Input(SurfaceEnum,&s[0]));
    9458         this->inputs->AddInput(new PentaP1Input(BedEnum,&b[0]));
     9457        this->inputs->AddInput(new PentaInput(SurfaceEnum,&s[0],P1Enum));
     9458        this->inputs->AddInput(new PentaInput(BedEnum,&b[0],P1Enum));
    94599459   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
    94609460
    … …  
    94629462        if(migration_style==SubelementMigrationEnum){
    94639463                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
    9464                 this->inputs->AddInput(new PentaP1Input(GLlevelsetEnum,&phi[0]));
     9464                this->inputs->AddInput(new PentaInput(GLlevelsetEnum,&phi[0],P1Enum));
    94659465                this->InputExtrude(GLlevelsetEnum,ElementEnum);
    94669466        }

  • issm/trunk-jpl/src/c/classes/Inputs/CMakeLists.txt

    r14284 r15417  
    1313# }}}
    1414# THREED_SOURCES {{{
    15 set(THREED_SOURCES $ENV{ISSM_DIR}/src/c/classes/objects/Inputs/PentaP1Input.cpp PARENT_SCOPE)
     15set(THREED_SOURCES $ENV{ISSM_DIR}/src/c/classes/objects/Inputs/PentaInput.cpp PARENT_SCOPE)
    1616# }}}

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

    r15300 r15417  
    3636                        maxvalues  =new TriaInput(enum_type,pmax,P1Enum);
    3737                        break;
    38                 case PentaP1InputEnum:
    39                         values     =new PentaP1Input(enum_type,pvalues);
    40                         savedvalues=new PentaP1Input(enum_type,pvalues);
    41                         minvalues  =new PentaP1Input(enum_type,pmin);
    42                         maxvalues  =new PentaP1Input(enum_type,pmax);
     38                case PentaInputEnum:
     39                        values     =new PentaInput(enum_type,pvalues,P1Enum);
     40                        savedvalues=new PentaInput(enum_type,pvalues,P1Enum);
     41                        minvalues  =new PentaInput(enum_type,pmin,P1Enum);
     42                        maxvalues  =new PentaInput(enum_type,pmax,P1Enum);
    4343                        break;
    4444                default:
    … …  
    115115
    116116/*Object functions*/
     117/*FUNCTION ControlInput::AXPY(){{{*/
     118void ControlInput::AXPY(Input* xinput,IssmDouble scalar){
     119        values->AXPY(xinput,scalar);
     120}/*}}}*/
    117121/*FUNCTION ControlInput::Constrain(){{{*/
    118122void ControlInput::Constrain(void){

  • issm/trunk-jpl/src/c/classes/Inputs/ControlInput.h

    r15130 r15417  
    7474                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    7575                void ArtificialNoise(IssmDouble min,IssmDouble max){_error_("not implemented yet");};
    76                 void AXPY(Input* xinput,IssmDouble scalar){_error_("not implemented yet");};
     76                void AXPY(Input* xinput,IssmDouble scalar);
    7777                void Constrain(void);
    7878                void Constrain(IssmDouble min,IssmDouble max);

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

    r15130 r15417  
    263263        switch(thickness_input->ObjectEnum()){
    264264
    265                 case PentaP1InputEnum:
     265                case PentaInputEnum:
    266266                        thickness_input->GetInputAverage(&thickness_value);
    267267                        this->value=this->value*thickness_value;

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

    r15409 r15417  
    1 /*!\file PentaP1Input.c
    2  * \brief: implementation of the PentaP1Input object
     1/*!\file PentaInput.c
     2 * \brief: implementation of the PentaInput object
    33 */
    44
    … …  
    1212#include "../../shared/shared.h"
    1313
    14 /*PentaP1Input constructors and destructor*/
    15 /*FUNCTION PentaP1Input::PentaP1Input(){{{*/
    16 PentaP1Input::PentaP1Input(){
    17         return;
    18 }
    19 /*}}}*/
    20 /*FUNCTION PentaP1Input::PentaP1Input(int in_enum_type,IssmDouble* values){{{*/
    21 PentaP1Input::PentaP1Input(int in_enum_type,IssmDouble* in_values)
     14/*PentaInput constructors and destructor*/
     15/*FUNCTION PentaInput::PentaInput(){{{*/
     16PentaInput::PentaInput(){
     17        values = NULL;
     18}
     19/*}}}*/
     20/*FUNCTION PentaInput::PentaInput(int in_enum_type,IssmDouble* values,int element_type){{{*/
     21PentaInput::PentaInput(int in_enum_type,IssmDouble* in_values,int element_type_in)
    2222                :PentaRef(1)
    2323{
    2424
    2525        /*Set PentaRef*/
    26         this->SetElementType(P1Enum,0);
    27         this->element_type=P1Enum;
    28 
     26        this->SetElementType(element_type_in,0);
     27        this->element_type=element_type_in;
     28
     29        /*Set Enum*/
    2930        enum_type=in_enum_type;
    30         values[0]=in_values[0];
    31         values[1]=in_values[1];
    32         values[2]=in_values[2];
    33         values[3]=in_values[3];
    34         values[4]=in_values[4];
    35         values[5]=in_values[5];
    36 }
    37 /*}}}*/
    38 /*FUNCTION PentaP1Input::~PentaP1Input(){{{*/
    39 PentaP1Input::~PentaP1Input(){
    40         return;
     31
     32        /*Set values*/
     33        this->values=xNew<IssmDouble>(this->NumberofNodes());
     34        for(int i=0;i<this->NumberofNodes();i++) values[i]=in_values[i];
     35}
     36/*}}}*/
     37/*FUNCTION PentaInput::~PentaInput(){{{*/
     38PentaInput::~PentaInput(){
     39        xDelete<IssmDouble>(this->values);
    4140}
    4241/*}}}*/
    4342
    4443/*Object virtual functions definitions:*/
    45 /*FUNCTION PentaP1Input::Echo {{{*/
    46 void PentaP1Input::Echo(void){
     44/*FUNCTION PentaInput::Echo {{{*/
     45void PentaInput::Echo(void){
    4746        this->DeepEcho();
    4847}
    4948/*}}}*/
    50 /*FUNCTION PentaP1Input::DeepEcho{{{*/
    51 void PentaP1Input::DeepEcho(void){
    52 
    53         _printf_("PentaP1Input:\n");
     49/*FUNCTION PentaInput::DeepEcho{{{*/
     50void PentaInput::DeepEcho(void){
     51
     52        _printf_("PentaInput:\n");
    5453        _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
    55         _printf_("   values: [" << this->values[0] << " " << this->values[1] << " " << this->values[2] << " " << this->values[3] << " " << this->values[4] << " " << this->values[5] << "]\n");
    56 }
    57 /*}}}*/
    58 /*FUNCTION PentaP1Input::Id{{{*/
    59 int    PentaP1Input::Id(void){ return -1; }
    60 /*}}}*/
    61 /*FUNCTION PentaP1Input::ObjectEnum{{{*/
    62 int PentaP1Input::ObjectEnum(void){
    63 
    64         return PentaP1InputEnum;
    65 
    66 }
    67 /*}}}*/
    68 
    69 /*PentaP1Input management*/
    70 /*FUNCTION PentaP1Input::copy{{{*/
    71 Object* PentaP1Input::copy() {
    72 
    73         return new PentaP1Input(this->enum_type,this->values);
    74 
    75 }
    76 /*}}}*/
    77 /*FUNCTION PentaP1Input::InstanceEnum{{{*/
    78 int PentaP1Input::InstanceEnum(void){
     54        _printf_("   values: [");
     55        for(int i=0;i<this->NumberofNodes();i++) _printf_(" "<<this->values[i]);
     56        _printf_("]\n");
     57}
     58/*}}}*/
     59/*FUNCTION PentaInput::Id{{{*/
     60int    PentaInput::Id(void){ return -1; }
     61/*}}}*/
     62/*FUNCTION PentaInput::ObjectEnum{{{*/
     63int PentaInput::ObjectEnum(void){
     64
     65        return PentaInputEnum;
     66
     67}
     68/*}}}*/
     69/*FUNCTION PentaInput::copy{{{*/
     70Object* PentaInput::copy() {
     71
     72        return new PentaInput(this->enum_type,this->values,this->element_type);
     73
     74}
     75/*}}}*/
     76
     77/*PentaInput management*/
     78/*FUNCTION PentaInput::InstanceEnum{{{*/
     79int PentaInput::InstanceEnum(void){
    7980
    8081        return this->enum_type;
    … …  
    8283}
    8384/*}}}*/
    84 /*FUNCTION PentaP1Input::SpawnTriaInput{{{*/
    85 Input* PentaP1Input::SpawnTriaInput(int* indices){
     85/*FUNCTION PentaInput::SpawnTriaInput{{{*/
     86Input* PentaInput::SpawnTriaInput(int* indices){
    8687
    8788        /*output*/
    8889        TriaInput* outinput=NULL;
    89         IssmDouble newvalues[3];
     90        IssmDouble newvalues[3]; //Assume P1 interpolation only for now
    9091
    9192        /*Loop over the new indices*/
    … …  
    107108}
    108109/*}}}*/
    109 /*FUNCTION PentaP1Input::SpawnResult{{{*/
    110 ElementResult* PentaP1Input::SpawnResult(int step, IssmDouble time){
     110/*FUNCTION PentaInput::SpawnResult{{{*/
     111ElementResult* PentaInput::SpawnResult(int step, IssmDouble time){
    111112
    112113        return new PentaP1ElementResult(this->enum_type,this->values,step,time);
    … …  
    116117
    117118/*Object functions*/
    118 /*FUNCTION PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
    119 void PentaP1Input::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){
     119/*FUNCTION PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
     120void PentaInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){
    120121
    121122        /*Call PentaRef function*/
    … …  
    124125}
    125126/*}}}*/
    126 /*FUNCTION PentaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
    127 void PentaP1Input::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){
     127/*FUNCTION PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){{{*/
     128void PentaInput::GetInputDerivativeValue(IssmDouble* p, IssmDouble* xyz_list, GaussPenta* gauss){
    128129
    129130        /*Call PentaRef function*/
    … …  
    131132}
    132133/*}}}*/
    133 /*FUNCTION PentaP1Input::GetVxStrainRate3d{{{*/
    134 void PentaP1Input::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
     134/*FUNCTION PentaInput::GetVxStrainRate3d{{{*/
     135void PentaInput::GetVxStrainRate3d(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
    135136        int i,j;
    136137
    … …  
    140141        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    141142        IssmDouble velocity[numnodes][DOFVELOCITY];
     143
     144        _assert_(this->NumberofNodes()==6); //Check Tria too
     145
     146        /*Get B matrix: */
     147        GetBStokes(&B[0][0], xyz_list, gauss);
     148
     149        /*Create a reduced matrix of B to get rid of pressure */
     150        for (i=0;i<6;i++){
     151                for (j=0;j<3;j++){
     152                        B_reduced[i][j]=B[i][j];
     153                }
     154                for (j=4;j<7;j++){
     155                        B_reduced[i][j-1]=B[i][j];
     156                }
     157                for (j=8;j<11;j++){
     158                        B_reduced[i][j-2]=B[i][j];
     159                }
     160                for (j=12;j<15;j++){
     161                        B_reduced[i][j-3]=B[i][j];
     162                }
     163                for (j=16;j<19;j++){
     164                        B_reduced[i][j-4]=B[i][j];
     165                }
     166                for (j=20;j<23;j++){
     167                        B_reduced[i][j-5]=B[i][j];
     168                }
     169        }
     170
     171        /*Here, we are computing the strain rate of (vx,0,0)*/
     172        for(i=0;i<numnodes;i++){
     173                velocity[i][0]=this->values[i];
     174                velocity[i][1]=0.0;
     175                velocity[i][2]=0.0;
     176        }
     177        /*Multiply B by velocity, to get strain rate: */
     178        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0);
     179
     180}
     181/*}}}*/
     182/*FUNCTION PentaInput::GetVyStrainRate3d{{{*/
     183void PentaInput::GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
     184        int i,j;
     185
     186        const int numnodes=6;
     187        const int DOFVELOCITY=3;
     188        IssmDouble B[8][27];
     189        IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
     190        IssmDouble velocity[numnodes][DOFVELOCITY];
     191
     192        _assert_(this->NumberofNodes()==6); //Check Tria too
    142193
    143194        /*Get B matrix: */
    … …  
    165216        }
    166217
    167         /*Here, we are computing the strain rate of (vx,0,0)*/
     218        /*Here, we are computing the strain rate of (0,vy,0)*/
    168219        for(i=0;i<numnodes;i++){
    169                 velocity[i][0]=this->values[i];
    170                 velocity[i][1]=0.0;
     220                velocity[i][0]=0.0;
     221                velocity[i][1]=this->values[i];
    171222                velocity[i][2]=0.0;
    172223        }
    173224        /*Multiply B by velocity, to get strain rate: */
    174         MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvx,0);
    175 
    176 }
    177 /*}}}*/
    178 /*FUNCTION PentaP1Input::GetVyStrainRate3d{{{*/
    179 void PentaP1Input::GetVyStrainRate3d(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
     225        MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0);
     226
     227}
     228/*}}}*/
     229/*FUNCTION PentaInput::GetVzStrainRate3d{{{*/
     230void PentaInput::GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){
    180231        int i,j;
    181232
    … …  
    188239        /*Get B matrix: */
    189240        GetBStokes(&B[0][0], xyz_list, gauss);
    190         /*Create a reduced matrix of B to get rid of pressure */
    191         for (i=0;i<6;i++){
    192                 for (j=0;j<3;j++){
    193                         B_reduced[i][j]=B[i][j];
    194                 }
    195                 for (j=4;j<7;j++){
    196                         B_reduced[i][j-1]=B[i][j];
    197                 }
    198                 for (j=8;j<11;j++){
    199                         B_reduced[i][j-2]=B[i][j];
    200                 }
    201                 for (j=12;j<15;j++){
    202                         B_reduced[i][j-3]=B[i][j];
    203                 }
    204                 for (j=16;j<19;j++){
    205                         B_reduced[i][j-4]=B[i][j];
    206                 }
    207                 for (j=20;j<23;j++){
    208                         B_reduced[i][j-5]=B[i][j];
    209                 }
    210         }
    211 
    212         /*Here, we are computing the strain rate of (0,vy,0)*/
    213         for(i=0;i<numnodes;i++){
    214                 velocity[i][0]=0.0;
    215                 velocity[i][1]=this->values[i];
    216                 velocity[i][2]=0.0;
    217         }
    218         /*Multiply B by velocity, to get strain rate: */
    219         MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numnodes,0,&velocity[0][0],DOFVELOCITY*numnodes,1,0,epsilonvy,0);
    220 
    221 }
    222 /*}}}*/
    223 /*FUNCTION PentaP1Input::GetVzStrainRate3d{{{*/
    224 void PentaP1Input::GetVzStrainRate3d(IssmDouble* epsilonvz,IssmDouble* xyz_list, GaussPenta* gauss){
    225         int i,j;
    226 
    227         const int numnodes=6;
    228         const int DOFVELOCITY=3;
    229         IssmDouble B[8][27];
    230         IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
    231         IssmDouble velocity[numnodes][DOFVELOCITY];
    232 
    233         /*Get B matrix: */
    234         GetBStokes(&B[0][0], xyz_list, gauss);
     241
     242        _assert_(this->NumberofNodes()==6); //Check Tria too
     243
    235244        /*Create a reduced matrix of B to get rid of pressure */
    236245        for (i=0;i<6;i++){
    … …  
    267276}
    268277/*}}}*/
    269 /*FUNCTION PentaP1Input::GetVxStrainRate3dPattyn{{{*/
    270 void PentaP1Input::GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
     278/*FUNCTION PentaInput::GetVxStrainRate3dPattyn{{{*/
     279void PentaInput::GetVxStrainRate3dPattyn(IssmDouble* epsilonvx,IssmDouble* xyz_list, GaussPenta* gauss){
    271280
    272281        int i;
    … …  
    278287        GetBPattyn(&B[0][0], xyz_list, gauss);
    279288
     289        _assert_(this->NumberofNodes()==6); //Check Tria too
     290
     291
    280292        /*Here, we are computing the strain rate of (vx,0)*/
    281293        for(i=0;i<numnodes;i++){
    … …  
    291303}
    292304/*}}}*/
    293 /*FUNCTION PentaP1Input::GetVyStrainRate3dPattyn{{{*/
    294 void PentaP1Input::GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
     305/*FUNCTION PentaInput::GetVyStrainRate3dPattyn{{{*/
     306void PentaInput::GetVyStrainRate3dPattyn(IssmDouble* epsilonvy,IssmDouble* xyz_list, GaussPenta* gauss){
    295307
    296308        int i;
    … …  
    302314        GetBPattyn(&B[0][0], xyz_list, gauss);
    303315
     316        _assert_(this->NumberofNodes()==6); //Check Tria too
     317
     318
    304319        /*Here, we are computing the strain rate of (0,vy)*/
    305320        for(i=0;i<numnodes;i++){
    … …  
    315330}
    316331/*}}}*/
    317 /*FUNCTION PentaP1Input::ChangeEnum{{{*/
    318 void PentaP1Input::ChangeEnum(int newenumtype){
     332/*FUNCTION PentaInput::ChangeEnum{{{*/
     333void PentaInput::ChangeEnum(int newenumtype){
    319334        this->enum_type=newenumtype;
    320335}
    321336/*}}}*/
    322 /*FUNCTION PentaP1Input::GetInputAverage{{{*/
    323 void PentaP1Input::GetInputAverage(IssmDouble* pvalue){
    324         *pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
     337/*FUNCTION PentaInput::GetInputAverage{{{*/
     338void PentaInput::GetInputAverage(IssmDouble* pvalue){
     339
     340        int        numnodes  = this->NumberofNodes();
     341        IssmDouble numnodesd = reCast<int,IssmDouble>(numnodes);
     342        IssmDouble value     = 0.;
     343
     344        for(int i=0;i<numnodes;i++) value+=values[i];
     345        value = value/numnodesd;
     346
     347        *pvalue=value;
    325348}
    326349/*}}}*/
    327350
    328351/*Intermediary*/
    329 /*FUNCTION PentaP1Input::SquareMin{{{*/
    330 void PentaP1Input::SquareMin(IssmDouble* psquaremin,Parameters* parameters){
    331 
    332         int i;
    333         const int numnodes=6;
    334         IssmDouble valuescopy[numnodes];
     352/*FUNCTION PentaInput::SquareMin{{{*/
     353void PentaInput::SquareMin(IssmDouble* psquaremin,Parameters* parameters){
     354
     355        int        numnodes=this->NumberofNodes();
    335356        IssmDouble squaremin;
    336357
    337         /*First,  copy values, to process units if requested: */
    338         for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];
    339 
    340358        /*Now, figure out minimum of valuescopy: */
    341         squaremin=pow(valuescopy[0],2);
    342         for(i=1;i<numnodes;i++){
    343                 if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
     359        squaremin=pow(this->values[0],2);
     360        for(int i=1;i<numnodes;i++){
     361                if(pow(this->values[i],2)<squaremin)squaremin=pow(this->values[i],2);
    344362        }
    345363        /*Assign output pointers:*/
    … …  
    347365}
    348366/*}}}*/
    349 /*FUNCTION PentaP1Input::ConstrainMin{{{*/
    350 void PentaP1Input::ConstrainMin(IssmDouble minimum){
    351 
    352         int i;
    353         const int numnodes=6;
    354 
    355         for(i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
    356 }
    357 /*}}}*/
    358 /*FUNCTION PentaP1Input::InfinityNorm{{{*/
    359 IssmDouble PentaP1Input::InfinityNorm(void){
     367/*FUNCTION PentaInput::ConstrainMin{{{*/
     368void PentaInput::ConstrainMin(IssmDouble minimum){
     369
     370        int numnodes = this->NumberofNodes();
     371        for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
     372}
     373/*}}}*/
     374/*FUNCTION PentaInput::InfinityNorm{{{*/
     375IssmDouble PentaInput::InfinityNorm(void){
    360376
    361377        /*Output*/
    362         const int numnodes=6;
    363         IssmDouble norm=0;
     378        IssmDouble norm=0.;
     379        int numnodes=this->NumberofNodes();
    364380
    365381        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
    … …  
    367383}
    368384/*}}}*/
    369 /*FUNCTION PentaP1Input::Max{{{*/
    370 IssmDouble PentaP1Input::Max(void){
    371 
    372         const int numnodes=6;
    373         IssmDouble    max=values[0];
     385/*FUNCTION PentaInput::Max{{{*/
     386IssmDouble PentaInput::Max(void){
     387
     388        int  numnodes=this->NumberofNodes();
     389        IssmDouble max=values[0];
    374390
    375391        for(int i=1;i<numnodes;i++){
    … …  
    379395}
    380396/*}}}*/
    381 /*FUNCTION PentaP1Input::MaxAbs{{{*/
    382 IssmDouble PentaP1Input::MaxAbs(void){
    383 
    384         const int numnodes=6;
    385         IssmDouble    max=fabs(values[0]);
     397/*FUNCTION PentaInput::MaxAbs{{{*/
     398IssmDouble PentaInput::MaxAbs(void){
     399
     400        int  numnodes=this->NumberofNodes();
     401        IssmDouble max=fabs(values[0]);
    386402
    387403        for(int i=1;i<numnodes;i++){
    … …  
    391407}
    392408/*}}}*/
    393 /*FUNCTION PentaP1Input::Min{{{*/
    394 IssmDouble PentaP1Input::Min(void){
    395 
    396         const int numnodes=6;
    397         IssmDouble    min=values[0];
     409/*FUNCTION PentaInput::Min{{{*/
     410IssmDouble PentaInput::Min(void){
     411
     412        const int  numnodes=this->NumberofNodes();
     413        IssmDouble min=values[0];
    398414
    399415        for(int i=1;i<numnodes;i++){
    … …  
    403419}
    404420/*}}}*/
    405 /*FUNCTION PentaP1Input::MinAbs{{{*/
    406 IssmDouble PentaP1Input::MinAbs(void){
    407 
    408         const int numnodes=6;
    409         IssmDouble    min=fabs(values[0]);
     421/*FUNCTION PentaInput::MinAbs{{{*/
     422IssmDouble PentaInput::MinAbs(void){
     423
     424        const int  numnodes=this->NumberofNodes();
     425        IssmDouble min=fabs(values[0]);
    410426
    411427        for(int i=1;i<numnodes;i++){
    … …  
    415431}
    416432/*}}}*/
    417 /*FUNCTION PentaP1Input::Scale{{{*/
    418 void PentaP1Input::Scale(IssmDouble scale_factor){
     433/*FUNCTION PentaInput::Scale{{{*/
     434void PentaInput::Scale(IssmDouble scale_factor){
     435
     436        const int numnodes=this->NumberofNodes();
     437        for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
     438}
     439/*}}}*/
     440/*FUNCTION PentaInput::AXPY{{{*/
     441void PentaInput::AXPY(Input* xinput,IssmDouble scalar){
     442
     443        const int numnodes=this->NumberofNodes();
     444        PentaInput* xpentainput=NULL;
     445
     446        /*If xinput is a ControlInput, take its values directly*/
     447        if(xinput->ObjectEnum()==ControlInputEnum){
     448                xinput=((ControlInput*)xinput)->values;
     449        }
     450
     451        /*xinput is of the same type, so cast it: */
     452        if(xinput->ObjectEnum()!=PentaInputEnum)
     453          _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     454        xpentainput=(PentaInput*)xinput;
     455        if(xpentainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xpentainput->element_type));
     456
     457        /*Carry out the AXPY operation depending on type:*/
     458        for(int i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*xpentainput->values[i];
     459
     460}
     461/*}}}*/
     462/*FUNCTION PentaInput::Constrain{{{*/
     463void PentaInput::Constrain(IssmDouble cm_min, IssmDouble cm_max){
    419464
    420465        int i;
    421         const int numnodes=6;
    422 
    423         for(i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
    424 }
    425 /*}}}*/
    426 /*FUNCTION PentaP1Input::AXPY{{{*/
    427 void PentaP1Input::AXPY(Input* xinput,IssmDouble scalar){
    428 
    429         int i;
    430         const int numnodes=6;
    431 
    432         /*xinput is of the same type, so cast it: */
    433 
    434         /*Carry out the AXPY operation depending on type:*/
    435         switch(xinput->ObjectEnum()){
    436 
    437                 case PentaP1InputEnum:{
    438                         PentaP1Input* cast_input=(PentaP1Input*)xinput;
    439                         for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
    440                         return;
    441                 case ControlInputEnum:{
    442                         ControlInput* cont_input=(ControlInput*)xinput;
    443                         if(cont_input->values->ObjectEnum()!=PentaP1InputEnum) _error_("not supported yet");
    444                         PentaP1Input* cast_input=(PentaP1Input*)cont_input->values;
    445                         for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*(cast_input->values[i]);}
    446                         return;
    447                 default:
    448                         _error_("not implemented yet");
    449         }
    450 
    451 }
    452 /*}}}*/
    453 /*FUNCTION PentaP1Input::Constrain{{{*/
    454 void PentaP1Input::Constrain(IssmDouble cm_min, IssmDouble cm_max){
    455 
    456         int i;
    457         const int numnodes=6;
     466        const int numnodes=this->NumberofNodes();
    458467
    459468        if(!xIsNan<IssmDouble>(cm_min)) for(i=0;i<numnodes;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
    … …  
    462471}
    463472/*}}}*/
    464 /*FUNCTION PentaP1Input::Extrude{{{*/
    465 void PentaP1Input::Extrude(void){
    466 
    467         int i;
    468 
    469         /*First 3 values copied on 3 last values*/
    470         for(i=0;i<3;i++) this->values[3+i]=this->values[i];
    471 }
    472 /*}}}*/
    473 /*FUNCTION PentaP1Input::VerticallyIntegrate{{{*/
    474 void PentaP1Input::VerticallyIntegrate(Input* thickness_input){
     473/*FUNCTION PentaInput::Extrude{{{*/
     474void PentaInput::Extrude(void){
     475
     476        switch(this->element_type){
     477                case P1Enum:
     478                        for(int i=0;i<3;i++) this->values[3+i]=this->values[i];
     479                        break;
     480                default:
     481                        _error_("not supported yet for type "<<EnumToStringx(this->element_type));
     482        }
     483}
     484/*}}}*/
     485/*FUNCTION PentaInput::VerticallyIntegrate{{{*/
     486void PentaInput::VerticallyIntegrate(Input* thickness_input){
    475487
    476488        IssmDouble thickness;
    … …  
    480492
    481493        /*vertically integrate depending on type:*/
    482         switch(thickness_input->ObjectEnum()){
    483 
    484                 case PentaP1InputEnum:{
     494        switch(this->element_type){
     495                case P1Enum:{
    485496                        GaussPenta *gauss=new GaussPenta();
    486497                        for(int iv=0;iv<3;iv++){
    … …  
    492503                        delete gauss;
    493504                        return; }
    494 
    495505                default:
    496                         _error_("not implemented yet");
    497         }
    498 }
    499 /*}}}*/
    500 /*FUNCTION PentaP1Input::PointwiseDivide{{{*/
    501 Input* PentaP1Input::PointwiseDivide(Input* inputB){
     506                        _error_("not supported yet for type "<<EnumToStringx(this->element_type));
     507        }
     508}
     509/*}}}*/
     510/*FUNCTION PentaInput::PointwiseDivide{{{*/
     511Input* PentaInput::PointwiseDivide(Input* inputB){
    502512
    503513        /*Ouput*/
    504         PentaP1Input* outinput=NULL;
     514        PentaInput* outinput=NULL;
    505515
    506516        /*Intermediaries*/
    507         PentaP1Input *xinputB  = NULL;
    508         const int     numnodes = 6;
    509         IssmDouble    AdotBvalues[numnodes];
     517        PentaInput *xinputB  = NULL;
     518        const int   numnodes = this->NumberofNodes();
    510519
    511520        /*Check that inputB is of the same type*/
    512         if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    513         xinputB=(PentaP1Input*)inputB;
     521        if(inputB->ObjectEnum()!=PentaInputEnum)     _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     522        xinputB=(PentaInput*)inputB;
     523        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
     524
     525        /*Allocate intermediary*/
     526        IssmDouble* AdotBvalues=xNew<IssmDouble>(numnodes);
    514527
    515528        /*Create point wise sum*/
    … …  
    520533
    521534        /*Create new Penta vertex input (copy of current input)*/
    522         outinput=new PentaP1Input(this->enum_type,&AdotBvalues[0]);
     535        outinput=new PentaInput(this->enum_type,AdotBvalues,this->element_type);
    523536
    524537        /*Return output pointer*/
     538        xDelete<IssmDouble>(AdotBvalues);
    525539        return outinput;
    526540
    527541}
    528542/*}}}*/
    529 /*FUNCTION PentaP1Input::PointwiseMin{{{*/
    530 Input* PentaP1Input::PointwiseMin(Input* inputB){
     543/*FUNCTION PentaInput::PointwiseMin{{{*/
     544Input* PentaInput::PointwiseMin(Input* inputB){
    531545
    532546        /*Ouput*/
    533         PentaP1Input* outinput=NULL;
     547        PentaInput* outinput=NULL;
    534548
    535549        /*Intermediaries*/
    536         int               i;
    537         PentaP1Input *xinputB     = NULL;
    538         const int         numnodes    = 6;
    539         IssmDouble            minvalues[numnodes];
     550        int         i;
     551        PentaInput  *xinputB   = NULL;
     552        const int   numnodes  = this->NumberofNodes();
     553        IssmDouble *minvalues = xNew<IssmDouble>(numnodes);
    540554
    541555        /*Check that inputB is of the same type*/
    542         if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    543         xinputB=(PentaP1Input*)inputB;
     556        if(inputB->ObjectEnum()!=PentaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     557        xinputB=(PentaInput*)inputB;
     558        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
    544559
    545560        /*Create point wise min*/
    … …  
    550565
    551566        /*Create new Penta vertex input (copy of current input)*/
    552         outinput=new PentaP1Input(this->enum_type,&minvalues[0]);
     567        outinput=new PentaInput(this->enum_type,&minvalues[0],this->element_type);
    553568
    554569        /*Return output pointer*/
     570        xDelete<IssmDouble>(minvalues);
    555571        return outinput;
    556 
    557 }
    558 /*}}}*/
    559 /*FUNCTION PentaP1Input::PointwiseMax{{{*/
    560 Input* PentaP1Input::PointwiseMax(Input* inputB){
     572}
     573/*}}}*/
     574/*FUNCTION PentaInput::PointwiseMax{{{*/
     575Input* PentaInput::PointwiseMax(Input* inputB){
    561576
    562577        /*Ouput*/
    563         PentaP1Input* outinput=NULL;
     578        PentaInput* outinput=NULL;
    564579
    565580        /*Intermediaries*/
    566         int               i;
    567         PentaP1Input *xinputB     = NULL;
    568         const int         numnodes    = 6;
    569         IssmDouble            maxvalues[numnodes];
     581        int         i;
     582        PentaInput  *xinputB   = NULL;
     583        const int   numnodes  = this->NumberofNodes();
     584        IssmDouble *maxvalues = xNew<IssmDouble>(numnodes);
    570585
    571586        /*Check that inputB is of the same type*/
    572         if (inputB->ObjectEnum()!=PentaP1InputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    573         xinputB=(PentaP1Input*)inputB;
     587        if(inputB->ObjectEnum()!=PentaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     588        xinputB=(PentaInput*)inputB;
     589        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
    574590
    575591        /*Create point wise max*/
    … …  
    580596
    581597        /*Create new Penta vertex input (copy of current input)*/
    582         outinput=new PentaP1Input(this->enum_type,&maxvalues[0]);
     598        outinput=new PentaInput(this->enum_type,&maxvalues[0],this->element_type);
    583599
    584600        /*Return output pointer*/
     601        xDelete<IssmDouble>(maxvalues);
    585602        return outinput;
    586 
    587 }
    588 /*}}}*/
    589 /*FUNCTION PentaP1Input::GetVectorFromInputs{{{*/
    590 void PentaP1Input::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist){
    591 
    592         const int numvertices=6;
    593         vector->SetValues(numvertices,doflist,this->values,INS_VAL);
     603}
     604/*}}}*/
     605/*FUNCTION PentaInput::GetVectorFromInputs{{{*/
     606void PentaInput::GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist){
     607
     608        const int numnodes=this->NumberofNodes();
     609        vector->SetValues(numnodes,doflist,this->values,INS_VAL);
    594610
    595611} /*}}}*/
    596 /*FUNCTION PentaP1Input::Configure{{{*/
    597 void PentaP1Input::Configure(Parameters* parameters){
     612/*FUNCTION PentaInput::Configure{{{*/
     613void PentaInput::Configure(Parameters* parameters){
    598614        /*do nothing: */
    599615}

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

    r15409 r15417  
    1 /*! \file PentaP1Input.h
    2  *  \brief: header file for PentaP1Input object
     1/*! \file PentaInput.h
     2 *  \brief: header file for PentaInput object
    33 */
    44
    5 #ifndef _PENTAP1INPUT_H_
    6 #define _PENTAP1INPUT_H_
     5#ifndef _PENTAINPUT_H_
     6#define _PENTAINPUT_H_
    77
    88/*Headers:*/
    … …  
    1414/*}}}*/
    1515
    16 class PentaP1Input: public Input, public PentaRef{
     16class PentaInput: public Input, public PentaRef{
    1717
    1818        public:
    19                 /*just hold 6 values for 6 vertices: */
    2019                int        enum_type;
    21                 IssmDouble values[6];
     20                IssmDouble* values;
    2221
    23                 /*PentaP1Input constructors, destructors: {{{*/
    24                 PentaP1Input();
    25                 PentaP1Input(int enum_type,IssmDouble* values);
    26                 ~PentaP1Input();
    27                 /*}}}*/
    28                 /*Object virtual functions definitions:{{{ */
     22                /*PentaInput constructors, destructors*/
     23                PentaInput();
     24                PentaInput(int enum_type,IssmDouble* values,int element_type_in);
     25                ~PentaInput();
     26
     27                /*Object virtual functions definitions */
    2928                void    Echo();
    3029                void    DeepEcho();
    … …  
    3231                int     ObjectEnum();
    3332                Object *copy();
    34                 /*}}}*/
    35                 /*PentaP1Input management: {{{*/
     33
     34                /*PentaInput management*/
    3635                int   InstanceEnum();
    3736                Input* SpawnTriaInput(int* indices);
    … …  
    4241                void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};
    4342                void Configure(Parameters* parameters);
    44                 /*}}}*/
    45                 /*numerics: {{{*/
     43                /*numerics*/
    4644                void GetInputValue(bool* pvalue){_error_("not implemented yet");};
    4745                void GetInputValue(int* pvalue){_error_("not implemented yet");};
    … …  
    8280                void VerticallyIntegrate(Input* thickness_input);
    8381                void GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist);
    84                 /*}}}*/
    8582
    8683};
    87 #endif  /* _PENTAP1INPUT_H */
     84#endif  /* _PENTAINPUT_H */

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

    r15298 r15417  
    3838TriaInput::~TriaInput(){
    3939        xDelete<IssmDouble>(this->values);
    40         return;
    4140}
    4241/*}}}*/
    … …  
    252251void TriaInput::ConstrainMin(IssmDouble minimum){
    253252
    254         const int numnodes = this->NumberofNodes();
     253        int numnodes = this->NumberofNodes();
    255254        for(int i=0;i<numnodes;i++) if (values[i]<minimum) values[i]=minimum;
    256255}
    … …  
    260259
    261260        /*Output*/
    262         IssmDouble norm=0;
    263         const int numnodes=this->NumberofNodes();
     261        IssmDouble norm=0.;
     262        int numnodes=this->NumberofNodes();
    264263
    265264        for(int i=0;i<numnodes;i++) if(fabs(values[i])>norm) norm=fabs(values[i]);
    … …  
    270269IssmDouble TriaInput::Max(void){
    271270
    272         const int  numnodes=this->NumberofNodes();
     271        int  numnodes=this->NumberofNodes();
    273272        IssmDouble max=values[0];
    274273
    … …  
    282281IssmDouble TriaInput::MaxAbs(void){
    283282
    284         const int  numnodes=this->NumberofNodes();
     283        int  numnodes=this->NumberofNodes();
    285284        IssmDouble max=fabs(values[0]);
    286285
    … …  
    319318
    320319        const int numnodes=this->NumberofNodes();
    321 
    322320        for(int i=0;i<numnodes;i++)values[i]=values[i]*scale_factor;
    323321}
    … …  
    342340void TriaInput::AXPY(Input* xinput,IssmDouble scalar){
    343341
    344         int i;
    345342        const int numnodes=this->NumberofNodes();
    346         TriaInput*  xtriavertexinput=NULL;
     343        TriaInput*  xtriainput=NULL;
    347344
    348345        /*xinput is of the same type, so cast it: */
    349346        if(xinput->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    350         xtriavertexinput=(TriaInput*)xinput;
    351         if(xtriavertexinput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
     347        xtriainput=(TriaInput*)xinput;
     348        if(xtriainput->element_type!=this->element_type) _error_("Operation not permitted because xinput is of type " << EnumToStringx(xinput->ObjectEnum()));
    352349
    353350        /*Carry out the AXPY operation depending on type:*/
    354         for(i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*xtriavertexinput->values[i];
     351        for(int i=0;i<numnodes;i++)this->values[i]=this->values[i]+scalar*xtriainput->values[i];
    355352
    356353}
    … …  
    387384
    388385        /*Check that inputB is of the same type*/
    389         if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     386        if(inputB->ObjectEnum()!=TriaInputEnum)       _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    390387        xinputB=(TriaInput*)inputB;
    391         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     388        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
    392389
    393390        /*Create point wise min*/
    … …  
    421418        if(inputB->ObjectEnum()!=TriaInputEnum) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
    422419        xinputB=(TriaInput*)inputB;
    423         if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(inputB->ObjectEnum()));
     420        if(xinputB->element_type!=this->element_type) _error_("Operation not permitted because inputB is of type " << EnumToStringx(xinputB->element_type));
    424421
    425422        /*Create point wise max*/

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r15300 r15417  
    8080                void VerticallyIntegrate(Input* thickness_input){_error_("not supported yet");};
    8181                void GetVectorFromInputs(Vector<IssmDouble>* vector,int* doflist);
    82                 /*}}}*/
    8382
    8483};

  • issm/trunk-jpl/src/c/classes/Materials/Matdamageice.cpp

    r15300 r15417  
    817817                if (iomodel->Data(MaterialsRheologyBEnum)) {
    818818                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    819                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs));
     819                        this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum));
    820820                }
    821821
    … …  
    823823                if (iomodel->Data(MaterialsRheologyNEnum)) {
    824824                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    825                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
     825                        this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum));
    826826                }
    827827
    … …  
    829829                if (iomodel->Data(MaterialsRheologyZEnum)) {
    830830                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    831                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
     831                        this->inputs->AddInput(new PentaInput(MaterialsRheologyZEnum,nodeinputs,P1Enum));
    832832                }
    833833
    … …  
    843843                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    844844                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    845                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     845                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    846846                                                }
    847847                                                break;
    … …  
    852852                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    853853                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    854                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     854                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    855855                                                }
    856856                                                break;

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

    r15384 r15417  
    1414#include "../Inputs/Inputs.h"
    1515#include "../Inputs/TriaInput.h"
    16 #include "../Inputs/PentaP1Input.h"
     16#include "../Inputs/PentaInput.h"
    1717#include "../Inputs/ControlInput.h"
    1818#include "../Elements/Element.h"
    … …  
    571571                                        IssmDouble valuesp[6];
    572572                                        for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector
    573                                         this->inputs->AddInput(new PentaP1Input(name,valuesp));
     573                                        this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum));
    574574                                        return;
    575575                                }
    … …  
    637637                                        IssmDouble valuesp[6];
    638638                                        for (int i=0;i<6;i++) valuesp[i]=vector[((Penta*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector
    639                                         this->inputs->AddInput(new PentaP1Input(name,valuesp));
     639                                        this->inputs->AddInput(new PentaInput(name,valuesp,P1Enum));
    640640                                        return;
    641641                                }
    … …  
    751751                if (iomodel->Data(MaterialsRheologyBEnum)) {
    752752                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    753                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs));
     753                        this->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,nodeinputs,P1Enum));
    754754                }
    755755
    … …  
    757757                if (iomodel->Data(MaterialsRheologyNEnum)) {
    758758                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    759                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
     759                        this->inputs->AddInput(new PentaInput(MaterialsRheologyNEnum,nodeinputs,P1Enum));
    760760                }
    761761
    … …  
    771771                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    772772                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
    773                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     773                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    774774                                                }
    775775                                                break;

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

    r15372 r15417  
    5656#include "./Inputs/DoubleInput.h"
    5757#include "./Inputs/IntInput.h"
    58 #include "./Inputs/PentaP1Input.h"
     58#include "./Inputs/PentaInput.h"
    5959#include "./Inputs/TriaInput.h"
    6060#include "./Inputs/ControlInput.h"

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

    r15401 r15417  
    370370        PenpairEnum,
    371371        PentaEnum,
    372         PentaP1InputEnum,
     372        PentaInputEnum,
    373373        ProfilerEnum,
    374374        MatrixParamEnum,

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

    r15401 r15417  
    370370                case PenpairEnum : return "Penpair";
    371371                case PentaEnum : return "Penta";
    372                 case PentaP1InputEnum : return "PentaP1Input";
     372                case PentaInputEnum : return "PentaInput";
    373373                case ProfilerEnum : return "Profiler";
    374374                case MatrixParamEnum : return "MatrixParam";

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

    r15401 r15417  
    376376              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
    377377              else if (strcmp(name,"Penta")==0) return PentaEnum;
    378               else if (strcmp(name,"PentaP1Input")==0) return PentaP1InputEnum;
     378              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    379379              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    380380              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
Note: See TracChangeset for help on using the changeset viewer.