Changeset 15298


Ignore:
Timestamp:
06/21/13 08:41:36 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: preparing quadratic elements for MacAyeal (Still some work to do...). Removed TriaP1Input and renamed TriaInput, which can have multiple interpolations

Location:
issm/trunk-jpl/src
Files:
3 added
2 deleted
20 edited

Legend:

Unmodified
Added
Removed

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

    r15174 r15298  
    8787                                        ./classes/Inputs/Input.h\
    8888                                        ./classes/Inputs/InputLocal.h\
    89                                         ./classes/Inputs/TriaP1Input.h\
    90                                         ./classes/Inputs/TriaP1Input.cpp\
     89                                        ./classes/Inputs/TriaInput.h\
     90                                        ./classes/Inputs/TriaInput.cpp\
    9191                                        ./classes/Inputs/BoolInput.h\
    9292                                        ./classes/Inputs/BoolInput.cpp\

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

    r15012 r15298  
    12571257}
    12581258/*}}}*/
     1259/*FUNCTION PentaRef::NumberofNodes{{{*/
     1260int PentaRef::NumberofNodes(void){
     1261
     1262        switch(this->element_type){
     1263                case P1Enum:   return NUMNODESP1;
     1264                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1265        }
     1266
     1267        return -1;
     1268}
     1269/*}}}*/

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

    r14951 r15298  
    5959                void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){_error_("only PentaGauss are supported");};
    6060
     61                int  NumberofNodes(void);
    6162};
    6263#endif

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

    r15286 r15298  
    887887
    888888        /*Add Stress tensor components into inputs*/
    889         this->inputs->AddInput(new TriaP1Input(StressTensorxxEnum,&sigma_xx[0]));
    890         this->inputs->AddInput(new TriaP1Input(StressTensorxyEnum,&sigma_xy[0]));
    891         this->inputs->AddInput(new TriaP1Input(StressTensorxzEnum,&sigma_xz[0]));
    892         this->inputs->AddInput(new TriaP1Input(StressTensoryyEnum,&sigma_yy[0]));
    893         this->inputs->AddInput(new TriaP1Input(StressTensoryzEnum,&sigma_yz[0]));
    894         this->inputs->AddInput(new TriaP1Input(StressTensorzzEnum,&sigma_zz[0]));
     889        this->inputs->AddInput(new TriaInput(StressTensorxxEnum,&sigma_xx[0]));
     890        this->inputs->AddInput(new TriaInput(StressTensorxyEnum,&sigma_xy[0]));
     891        this->inputs->AddInput(new TriaInput(StressTensorxzEnum,&sigma_xz[0]));
     892        this->inputs->AddInput(new TriaInput(StressTensoryyEnum,&sigma_yy[0]));
     893        this->inputs->AddInput(new TriaInput(StressTensoryzEnum,&sigma_yz[0]));
     894        this->inputs->AddInput(new TriaInput(StressTensorzzEnum,&sigma_zz[0]));
    895895
    896896        /*Clean up and return*/
    … …  
    10281028        for (int imonth=0;imonth<12;imonth++) {
    10291029                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    1030                 TriaP1Input* newmonthinput1 = new TriaP1Input(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0]);
     1030                TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0]);
    10311031                NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    10321032
    10331033                for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    1034                 TriaP1Input* newmonthinput2 = new TriaP1Input(SurfaceforcingsPrecipitationEnum,&tmp[0]);
     1034                TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0]);
    10351035                NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    10361036        }
    … …  
    13941394                for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
    13951395        }
     1396}
     1397/*}}}*/
     1398/*FUNCTION Tria::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/
     1399void Tria::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){
     1400
     1401        _assert_(pvalue);
     1402
     1403        Input *input    = inputs->GetInput(enumtype);
     1404        int    numnodes = this->NumberofNodes();
     1405
     1406        /* Start looping on the number of vertices: */
     1407        if (input){
     1408                GaussTria* gauss=new GaussTria();
     1409                for (int iv=0;iv<this->NumberofNodes();iv++){
     1410                        gauss->GaussNode(numnodes,iv);
     1411                        input->GetInputValue(&pvalue[iv],gauss);
     1412                }
     1413                delete gauss;
     1414        }
     1415        else{
     1416                for (int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
     1417        }
     1418}
     1419/*}}}*/
     1420/*FUNCTION Tria::GetInputListOnNodes(IssmDouble* pvalue,int enumtype) {{{*/
     1421void Tria::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){
     1422
     1423        _assert_(pvalue);
     1424
     1425        /*Recover input*/
     1426        Input* input=inputs->GetInput(enumtype);
     1427        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
     1428        int    numnodes = this->NumberofNodes();
     1429
     1430        /* Start looping on the number of vertices: */
     1431        GaussTria* gauss=new GaussTria();
     1432        for (int iv=0;iv<this->NumberofNodes();iv++){
     1433                gauss->GaussNode(numnodes,iv);
     1434                input->GetInputValue(&pvalue[iv],gauss);
     1435        }
     1436        delete gauss;
    13961437}
    13971438/*}}}*/
    … …  
    16871728                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    16881729                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1689                                                 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1730                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    16901731                                        }
    16911732                                        break;
    … …  
    16951736                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    16961737                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1697                                                 this->inputs->AddInput(new ControlInput(VxEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1738                                                this->inputs->AddInput(new ControlInput(VxEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    16981739                                        }
    16991740                                        break;
    … …  
    17031744                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    17041745                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
    1705                                                 this->inputs->AddInput(new ControlInput(VyEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1746                                                this->inputs->AddInput(new ControlInput(VyEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17061747                                        }
    17071748                                        break;
    … …  
    17111752                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    17121753                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    1713                                                 this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1754                                                this->inputs->AddInput(new ControlInput(ThicknessEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17141755                                        }
    17151756                                        break;
    … …  
    17191760                                                for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    17201761                                                for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
    1721                                                 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1762                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    17221763                                        }
    17231764                                        break;
    … …  
    17391780                for(i=0;i<num_cm_responses;i++){
    17401781                        for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
    1741                         datasetinput->inputs->AddObject(new TriaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs));
     1782                        datasetinput->inputs->AddObject(new TriaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs));
    17421783                }
    17431784
    … …  
    18251866
    18261867        /*Add input to the element: */
    1827         this->inputs->AddInput(new TriaP1Input(enum_type,values));
     1868        this->inputs->AddInput(new TriaInput(enum_type,values));
    18281869
    18291870        /*Free ressources:*/
    … …  
    18901931
    18911932        /*Add input to the element: */
    1892         this->inputs->AddInput(new TriaP1Input(ThicknessEnum,newthickness));
    1893         this->inputs->AddInput(new TriaP1Input(SurfaceEnum,newsurface));
    1894         this->inputs->AddInput(new TriaP1Input(BedEnum,newbed));
     1933        this->inputs->AddInput(new TriaInput(ThicknessEnum,newthickness));
     1934        this->inputs->AddInput(new TriaInput(SurfaceEnum,newsurface));
     1935        this->inputs->AddInput(new TriaInput(BedEnum,newbed));
    18951936
    18961937        /*Free ressources:*/
    … …  
    19151956                /*update input*/
    19161957                if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    1917                         material->inputs->AddInput(new TriaP1Input(name,values));
     1958                        material->inputs->AddInput(new TriaInput(name,values));
    19181959                }
    19191960                else{
    1920                         this->inputs->AddInput(new TriaP1Input(name,values));
     1961                        this->inputs->AddInput(new TriaInput(name,values));
    19211962                }
    19221963                return;
    … …  
    19281969                /*update input*/
    19291970                if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    1930                         material->inputs->AddInput(new TriaP1Input(name,values));
     1971                        material->inputs->AddInput(new TriaInput(name,values));
    19311972                }
    19321973                else{
    1933                         this->inputs->AddInput(new TriaP1Input(name,values));
     1974                        this->inputs->AddInput(new TriaInput(name,values));
    19341975                }
    19351976                return;
    … …  
    19451986                }
    19461987                /*Add input to the element: */
    1947                 this->inputs->AddInput(new TriaP1Input(name,values));
     1988                this->inputs->AddInput(new TriaInput(name,values));
    19481989
    19491990                /*Free ressources:*/
    … …  
    19571998                }
    19581999                /*Add input to the element: */
    1959                 this->inputs->AddInput(new TriaP1Input(name,values));
     2000                this->inputs->AddInput(new TriaInput(name,values));
    19602001
    19612002                /*Free ressources:*/
    … …  
    20332074
    20342075                        /*create static input: */
    2035                         this->inputs->AddInput(new TriaP1Input(vector_enum,nodeinputs));
     2076                        this->inputs->AddInput(new TriaInput(vector_enum,nodeinputs));
    20362077                }
    20372078                else if(M==numberofvertices+1){
    … …  
    20492090
    20502091                                if(t==0) transientinput=new TransientInput(vector_enum);
    2051                                 transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time);
     2092                                transientinput->AddTimeInput(new TriaInput(vector_enum,nodeinputs),time);
    20522093                        }
    20532094                        this->inputs->AddInput(transientinput);
    … …  
    23162357        if(!this->IsFloating() && floatingelement==true){
    23172358                for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
    2318                 this->inputs->AddInput(new TriaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
     2359                this->inputs->AddInput(new TriaInput(BasalforcingsMeltingRateEnum,&melting[0]));
    23192360        }
    23202361
    23212362        /*Update inputs*/
    23222363   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,floatingelement));
    2323         this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0]));
    2324         this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0]));
     2364        this->inputs->AddInput(new TriaInput(SurfaceEnum,&s[0]));
     2365        this->inputs->AddInput(new TriaInput(BedEnum,&b[0]));
    23252366
    23262367        /*Recalculate phi*/
    23272368        if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){
    23282369                for(i=0;i<NUMVERTICES;i++) phi[i]=h[i]+ba[i]/density;
    2329                 this->inputs->AddInput(new TriaP1Input(GLlevelsetEnum,&phi[0]));
     2370                this->inputs->AddInput(new TriaInput(GLlevelsetEnum,&phi[0]));
    23302371        }
    23312372}
    … …  
    24862527
    24872528   /*Update inputs*/   
    2488    this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
     2529   this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0]));
    24892530   // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    24902531
    … …  
    25762617                }  //end of the loop over the vertices
    25772618          /*Update inputs*/
    2578           this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0]));
     2619          this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&smb[0]));
    25792620}
    25802621/*}}}*/
    … …  
    27602801                        if(!iomodel->Data(VxEnum)){
    27612802                                for(i=0;i<3;i++)nodeinputs[i]=0;
    2762                                 this->inputs->AddInput(new TriaP1Input(VxEnum,nodeinputs));
    2763                                 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVxEnum,nodeinputs));
     2803                                this->inputs->AddInput(new TriaInput(VxEnum,nodeinputs));
     2804                                if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVxEnum,nodeinputs));
    27642805                        }
    27652806                        if(!iomodel->Data(VyEnum)){
    27662807                                for(i=0;i<3;i++)nodeinputs[i]=0;
    2767                                 this->inputs->AddInput(new TriaP1Input(VyEnum,nodeinputs));
    2768                                 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVyEnum,nodeinputs));
     2808                                this->inputs->AddInput(new TriaInput(VyEnum,nodeinputs));
     2809                                if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVyEnum,nodeinputs));
    27692810                        }
    27702811                        if(!iomodel->Data(VzEnum)){
    27712812                                for(i=0;i<3;i++)nodeinputs[i]=0;
    2772                                 this->inputs->AddInput(new TriaP1Input(VzEnum,nodeinputs));
    2773                                 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVzEnum,nodeinputs));
     2813                                this->inputs->AddInput(new TriaInput(VzEnum,nodeinputs));
     2814                                if(dakota_analysis) this->inputs->AddInput(new TriaInput(QmuVzEnum,nodeinputs));
    27742815                        }
    27752816                        if(!iomodel->Data(PressureEnum)){
    27762817                                for(i=0;i<3;i++)nodeinputs[i]=0;
    27772818                                if(dakota_analysis){
    2778                                         this->inputs->AddInput(new TriaP1Input(PressureEnum,nodeinputs));
    2779                                         this->inputs->AddInput(new TriaP1Input(QmuPressureEnum,nodeinputs));
     2819                                        this->inputs->AddInput(new TriaInput(PressureEnum,nodeinputs));
     2820                                        this->inputs->AddInput(new TriaInput(QmuPressureEnum,nodeinputs));
    27802821                                }
    27812822                        }
    … …  
    32363277ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
    32373278
    3238         /*Constants*/
    3239         const int  numdof=NDOF2*NUMVERTICES;
    3240 
    32413279        /*Intermediaries*/
    32423280        int        i,j;
    … …  
    32453283        IssmDouble viscosity_overshoot,thickness,Jdet;
    32463284        IssmDouble epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
    3247         IssmDouble B[3][numdof];
    3248         IssmDouble Bprime[3][numdof];
    3249         IssmDouble D[3][3]   = {0.0};
    32503285        IssmDouble D_scalar;
    3251         GaussTria *gauss = NULL;
    3252 
    3253         /*Initialize Element matrix*/
    3254         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     3286
     3287        /*Fetch number of nodes and dof for this finite element*/
     3288        int numnodes = this->NumberofNodes();
     3289        int numdof   = numnodes*NDOF2;
     3290
     3291        /*Initialize Element matrix, vectors and Gaussian points*/
     3292        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum);
     3293        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
     3294        IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
     3295        IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
    32553296
    32563297        /*Retrieve all inputs and parameters*/
    … …  
    32643305
    32653306        /* Start  looping on the number of gaussian points: */
    3266         gauss=new GaussTria(2);
     3307        GaussTria*     gauss  = new GaussTria(2);
    32673308        for(int ig=gauss->begin();ig<gauss->end();ig++){
    32683309
    … …  
    32703311
    32713312                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3272                 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
    3273                 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
     3313                GetBMacAyeal(&B[0], &xyz_list[0][0], gauss);
     3314                GetBprimeMacAyeal(&Bprime[0], &xyz_list[0][0], gauss);
    32743315
    32753316                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    … …  
    32803321
    32813322                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    3282                 D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
    3283                 for (i=0;i<3;i++) D[i][i]=D_scalar;
    3284 
    3285                 TripleMultiply(&B[0][0],3,numdof,1,
    3286                                         &D[0][0],3,3,0,
    3287                                         &Bprime[0][0],3,numdof,0,
     3323                D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
     3324                for(i=0;i<3;i++) D[i*3+i]=D_scalar;
     3325
     3326                TripleMultiply(B,3,numdof,1,
     3327                                        D,3,3,0,
     3328                                        Bprime,3,numdof,0,
    32883329                                        &Ke->values[0],1);
    32893330        }
    … …  
    32943335        /*Clean up and return*/
    32953336        delete gauss;
     3337        xDelete<IssmDouble>(B);
     3338        xDelete<IssmDouble>(D);
     3339        xDelete<IssmDouble>(Bprime);
    32963340        return Ke;
    32973341}
    … …  
    32993343/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{*/
    33003344ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
    3301 
    3302         /*Constants*/
    3303         const int  numdof=NDOF2*NUMVERTICES;
    33043345
    33053346        /*Intermediaries*/
    … …  
    33123353        IssmDouble fraction1,fraction2;
    33133354        IssmDouble phi=1.0;
    3314         IssmDouble L[2][numdof];
    3315         IssmDouble DL[2][2]  = {{ 0,0 },{0,0}};
    3316         IssmDouble DL_scalar;
    3317         IssmDouble slope[2]  = {0.0,0.0};
     3355        IssmDouble D_scalar;
    33183356        IssmDouble gllevelset;
    33193357        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    33213359        GaussTria *gauss    = NULL;
    33223360
    3323         /*Initialize Element matrix and return if necessary*/
     3361        /*Return if element is inactive*/
    33243362        if(IsFloating()) return NULL;
    3325         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     3363
     3364        /*Fetch number of nodes and dof for this finite element*/
     3365        int numnodes = this->NumberofNodes();
     3366        int numdof   = numnodes*NDOF2;
     3367
     3368        /*Initialize Element matrix and vectors*/
     3369        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum);
     3370        IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
     3371        IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
    33263372
    33273373        /*Retrieve all inputs and parameters*/
    … …  
    33613407                }
    33623408
    3363                 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
     3409                GetL(&B[0], &xyz_list[0][0], gauss,NDOF2);
    33643410                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3365                 DL_scalar=alpha2*gauss->weight*Jdet;
    3366                 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    3367 
    3368                 TripleMultiply( &L[0][0],2,numdof,1,
    3369                                         &DL[0][0],2,2,0,
    3370                                         &L[0][0],2,numdof,0,
     3411                D_scalar=alpha2*gauss->weight*Jdet;
     3412                for(i=0;i<2;i++) D[i*2+i]=D_scalar;
     3413
     3414                TripleMultiply(B,2,numdof,1,
     3415                                        D,2,2,0,
     3416                                        B,2,numdof,0,
    33713417                                        &Ke->values[0],1);
    33723418        }
    … …  
    33783424        delete gauss;
    33793425        delete friction;
     3426        xDelete<IssmDouble>(B);
     3427        xDelete<IssmDouble>(D);
    33803428        return Ke;
    33813429}
    … …  
    34043452/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{*/
    34053453ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
    3406 
    3407         /*Constants*/
    3408         const int    numdof=NDOF2*NUMVERTICES;
    34093454
    34103455        /*Intermediaries */
    … …  
    34143459        IssmDouble     xyz_list[NUMVERTICES][3];
    34153460        IssmDouble     slope[2];
    3416         IssmDouble     basis[3];
    3417         IssmDouble     pe_g_gaussian[numdof];
    3418         GaussTria*     gauss=NULL;
    3419 
    3420         /*Initialize Element vector*/
     3461
     3462        /*Fetch number of nodes and dof for this finite element*/
     3463        int numnodes = this->NumberofNodes();
     3464        int numdof   = numnodes*NDOF2;
     3465
     3466        /*Initialize Element vector and vectors*/
    34213467        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
     3468        GaussTria*     gauss  = new GaussTria(2);
     3469        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    34223470
    34233471        /*Retrieve all inputs and parameters*/
    … …  
    34283476
    34293477        /* Start  looping on the number of gaussian points: */
    3430         gauss=new GaussTria(2);
    34313478        for(int ig=gauss->begin();ig<gauss->end();ig++){
    34323479
    … …  
    34403487                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
    34413488
    3442                 /*Build pe_g_gaussian vector: */
    3443                 for (i=0;i<NUMVERTICES;i++){
     3489                /*Build load vector: */
     3490                for (i=0;i<numnodes;i++){
    34443491                        for (j=0;j<NDOF2;j++){
    34453492                                pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    … …  
    34533500        /*Clean up and return*/
    34543501        delete gauss;
     3502        xDelete<IssmDouble>(basis);
    34553503        return pe;
    34563504}
    … …  
    35083556/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{*/
    35093557ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
    3510 
    3511         /*Constants*/
    3512         const int    numdof=NDOF2*NUMVERTICES;
    35133558
    35143559        /*Intermediaries */
    … …  
    35213566        IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];*/
    35223567        IssmDouble eps1[2],eps2[2];
    3523         IssmDouble phi[NUMVERTICES];
    3524         IssmDouble dphi[2][NUMVERTICES];
    3525         GaussTria *gauss=NULL;
    3526 
    3527         /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
    3528         ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
     3568        GaussTria* gauss = NULL;
     3569
     3570        /*Fetch number of nodes and dof for this finite element*/
     3571        int numnodes = this->NumberofNodes();
     3572        int numdof   = numnodes*NDOF2;
     3573
     3574        /*Initialize Element matrix, vectors and Gaussian points*/
     3575        ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); //Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)
     3576        IssmDouble*    dphi   = xNew<IssmDouble>(2*numnodes);
    35293577
    35303578        /*Retrieve all inputs and parameters*/
    … …  
    35413589
    35423590                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3543                 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     3591                GetNodalFunctionsDerivatives(dphi,&xyz_list[0][0],gauss);
    35443592
    35453593                thickness_input->GetInputValue(&thickness, gauss);
    … …  
    35493597                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    35503598
    3551                 for(i=0;i<3;i++){
    3552                         for(j=0;j<3;j++){
    3553                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
    3554                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
    3555                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
    3556                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
    3557 
    3558                                 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
    3559                                 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
    3560                                 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
    3561                                 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
     3599                for(i=0;i<numnodes;i++){
     3600                        for(j=0;j<numnodes;j++){
     3601                                eps1dotdphii=eps1[0]*dphi[2*i+0]+eps1[1]*dphi[2*i+1];
     3602                                eps1dotdphij=eps1[0]*dphi[2*j+0]+eps1[1]*dphi[2*j+1];
     3603                                eps2dotdphii=eps2[0]*dphi[2*i+0]+eps2[1]*dphi[2*i+1];
     3604                                eps2dotdphij=eps2[0]*dphi[2*j+0]+eps2[1]*dphi[2*j+1];
     3605
     3606                                Ke->values[2*numnodes*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
     3607                                Ke->values[2*numnodes*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
     3608                                Ke->values[2*numnodes*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
     3609                                Ke->values[2*numnodes*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
    35623610                        }
    35633611                }
    … …  
    35683616
    35693617        /*Clean up and return*/
     3618        xDelete<IssmDouble>(dphi);
    35703619        delete gauss;
    35713620        return Ke;
    … …  
    35753624void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solution){
    35763625
    3577         const int    numdof=NDOF2*NUMVERTICES;
    3578 
    3579         int          i;
    3580         int*         doflist=NULL;
    3581         IssmDouble       vx,vy;
    3582         IssmDouble       values[numdof];
    3583         GaussTria*   gauss=NULL;
    3584 
    3585         /*Get dof list: */
     3626        IssmDouble   vx,vy;
     3627        int*         doflist = NULL;
     3628        GaussTria*   gauss   = NULL;
     3629
     3630        /*Fetch number of nodes and dof for this finite element*/
     3631        int numnodes = this->NumberofNodes();
     3632        int numdof   = numnodes*NDOF2;
     3633
     3634        /*Fetch dof list and allocate solution vector*/
    35863635        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3636        IssmDouble* values = xNew<IssmDouble>(numdof);
    35873637
    35883638        /*Get inputs*/
    … …  
    35933643        /*P1 element only for now*/
    35943644        gauss=new GaussTria();
    3595         for(i=0;i<NUMVERTICES;i++){
    3596 
    3597                 gauss->GaussVertex(i);
     3645        for(int i=0;i<numnodes;i++){
     3646
     3647                gauss->GaussNode(numnodes,i);
    35983648
    35993649                /*Recover vx and vy*/
    … …  
    36093659        delete gauss;
    36103660        xDelete<int>(doflist);
     3661        xDelete<IssmDouble>(values);
    36113662}
    36123663/*}}}*/
    … …  
    36533704void  Tria::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){
    36543705
    3655         const int numdof=NDOF2*NUMVERTICES;
    3656 
    3657         int       i;
    3658         int*      doflist=NULL;
    3659         IssmDouble    rho_ice,g;
    3660         IssmDouble    values[numdof];
    3661         IssmDouble    vx[NUMVERTICES];
    3662         IssmDouble    vy[NUMVERTICES];
    3663         IssmDouble    vz[NUMVERTICES];
    3664         IssmDouble    vel[NUMVERTICES];
    3665         IssmDouble    pressure[NUMVERTICES];
    3666         IssmDouble    thickness[NUMVERTICES];
    3667 
    3668         /*Get dof list: */
     3706        int        i;
     3707        IssmDouble rho_ice,g;
     3708        int*       doflist=NULL;
     3709
     3710        /*Fetch number of nodes and dof for this finite element*/
     3711        int numnodes = this->NumberofNodes();
     3712        int numdof   = numnodes*NDOF2;
     3713
     3714        /*Fetch dof list and allocate solution vectors*/
    36693715        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3716        IssmDouble* values    = xNew<IssmDouble>(numdof);
     3717        IssmDouble* vx        = xNew<IssmDouble>(numdof);
     3718        IssmDouble* vy        = xNew<IssmDouble>(numdof);
     3719        IssmDouble* vz        = xNew<IssmDouble>(numdof);
     3720        IssmDouble* vel       = xNew<IssmDouble>(numdof);
     3721        IssmDouble* pressure  = xNew<IssmDouble>(numdof);
     3722        IssmDouble* thickness = xNew<IssmDouble>(numdof);
    36703723
    36713724        /*Use the dof list to index into the solution vector: */
    … …  
    36733726
    36743727        /*Transform solution in Cartesian Space*/
    3675         TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
     3728        TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    36763729
    36773730        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    3678         for(i=0;i<NUMVERTICES;i++){
     3731        for(i=0;i<numnodes;i++){
    36793732                vx[i]=values[i*NDOF2+0];
    36803733                vy[i]=values[i*NDOF2+1];
    … …  
    36863739
    36873740        /*Get Vz and compute vel*/
    3688         GetInputListOnVertices(&vz[0],VzEnum,0);
    3689         for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     3741        GetInputListOnNodes(&vz[0],VzEnum,0.);
     3742        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    36903743
    36913744        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    … …  
    36933746        rho_ice=matpar->GetRhoIce();
    36943747        g=matpar->GetG();
    3695         GetInputListOnVertices(&thickness[0],ThicknessEnum);
    3696         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
     3748        GetInputListOnNodes(&thickness[0],ThicknessEnum);
     3749        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
    36973750
    36983751        /*Now, we have to move the previous Vx and Vy inputs  to old
    … …  
    37033756
    37043757        /*Add vx and vy as inputs to the tria element: */
    3705         this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
    3706         this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
    3707         this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
    3708         this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
     3758        this->inputs->AddInput(new TriaInput(VxEnum,vx));
     3759        this->inputs->AddInput(new TriaInput(VyEnum,vy));
     3760        this->inputs->AddInput(new TriaInput(VelEnum,vel));
     3761        this->inputs->AddInput(new TriaInput(PressureEnum,pressure));
    37093762
    37103763        /*Free ressources:*/
    37113764        xDelete<int>(doflist);
    3712 
     3765        xDelete<IssmDouble>(values);
     3766        xDelete<IssmDouble>(vx);
     3767        xDelete<IssmDouble>(vy);
     3768        xDelete<IssmDouble>(vz);
     3769        xDelete<IssmDouble>(vel);
     3770        xDelete<IssmDouble>(pressure);
     3771        xDelete<IssmDouble>(thickness);
    37133772}
    37143773/*}}}*/
    … …  
    37163775void  Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){
    37173776
    3718         const int numdof=NDOF2*NUMVERTICES;
    3719 
    3720         int       i;
    3721         int*      doflist=NULL;
    3722         IssmDouble    rho_ice,g;
    3723         IssmDouble    values[numdof];
    3724         IssmDouble    vx[NUMVERTICES];
    3725         IssmDouble    vy[NUMVERTICES];
    3726         IssmDouble    vz[NUMVERTICES];
    3727         IssmDouble    vel[NUMVERTICES];
    3728         IssmDouble    pressure[NUMVERTICES];
    3729         IssmDouble    thickness[NUMVERTICES];
    3730 
    3731         /*Get dof list: */
     3777        int        i;
     3778        IssmDouble rho_ice,g;
     3779        int*       doflist=NULL;
     3780
     3781        /*Fetch number of nodes and dof for this finite element*/
     3782        int numnodes = this->NumberofNodes();
     3783        int numdof   = numnodes*NDOF2;
     3784
     3785        /*Fetch dof list and allocate solution vectors*/
    37323786        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3787        IssmDouble* values    = xNew<IssmDouble>(numdof);
     3788        IssmDouble* vx        = xNew<IssmDouble>(numdof);
     3789        IssmDouble* vy        = xNew<IssmDouble>(numdof);
     3790        IssmDouble* vz        = xNew<IssmDouble>(numdof);
     3791        IssmDouble* vel       = xNew<IssmDouble>(numdof);
     3792        IssmDouble* pressure  = xNew<IssmDouble>(numdof);
     3793        IssmDouble* thickness = xNew<IssmDouble>(numdof);
    37333794
    37343795        /*Use the dof list to index into the solution vector: */
    37353796        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    37363797
     3798        /*Transform solution in Cartesian Space*/
     3799        TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
     3800
    37373801        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    3738         for(i=0;i<NUMVERTICES;i++){
     3802        for(i=0;i<numnodes;i++){
    37393803                vx[i]=values[i*NDOF2+0];
    37403804                vy[i]=values[i*NDOF2+1];
    … …  
    37453809        }
    37463810
    3747         /*Now Compute vel*/
    3748         GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    3749         for(i=0;i<NUMVERTICES;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     3811        /*Get Vz and compute vel*/
     3812        GetInputListOnNodes(&vz[0],VzEnum,0.);
     3813        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    37503814
    37513815        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    … …  
    37533817        rho_ice=matpar->GetRhoIce();
    37543818        g=matpar->GetG();
    3755         GetInputListOnVertices(&thickness[0],ThicknessEnum);
    3756         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
     3819        GetInputListOnNodes(&thickness[0],ThicknessEnum);
     3820        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
    37573821
    37583822        /*Now, we have to move the previous Vx and Vy inputs  to old
    … …  
    37633827
    37643828        /*Add vx and vy as inputs to the tria element: */
    3765         this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
    3766         this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
    3767         this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
    3768         this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
     3829        this->inputs->AddInput(new TriaInput(VxEnum,vx));
     3830        this->inputs->AddInput(new TriaInput(VyEnum,vy));
     3831        this->inputs->AddInput(new TriaInput(VelEnum,vel));
     3832        this->inputs->AddInput(new TriaInput(PressureEnum,pressure));
    37693833
    37703834        /*Free ressources:*/
    37713835        xDelete<int>(doflist);
     3836        xDelete<IssmDouble>(values);
     3837        xDelete<IssmDouble>(vx);
     3838        xDelete<IssmDouble>(vy);
     3839        xDelete<IssmDouble>(vz);
     3840        xDelete<IssmDouble>(vel);
     3841        xDelete<IssmDouble>(pressure);
     3842        xDelete<IssmDouble>(thickness);
    37723843}
    37733844/*}}}*/
    … …  
    39223993        GradientIndexing(&vertexpidlist[0],control_index);
    39233994        for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
    3924         grad_input=new TriaP1Input(GradientEnum,grad_list);
     3995        grad_input=new TriaInput(GradientEnum,grad_list);
    39253996
    39263997        ((ControlInput*)input)->SetGradient(grad_input);
    … …  
    56895760
    56905761        /*Add vx and vy as inputs to the tria element: */
    5691         this->inputs->AddInput(new TriaP1Input(AdjointxEnum,lambdax));
    5692         this->inputs->AddInput(new TriaP1Input(AdjointyEnum,lambday));
     5762        this->inputs->AddInput(new TriaInput(AdjointxEnum,lambdax));
     5763        this->inputs->AddInput(new TriaInput(AdjointyEnum,lambday));
    56935764
    56945765        /*Free ressources:*/
    … …  
    57195790
    57205791        /*Add vx and vy as inputs to the tria element: */
    5721         this->inputs->AddInput(new TriaP1Input(AdjointEnum,lambda));
     5792        this->inputs->AddInput(new TriaInput(AdjointEnum,lambda));
    57225793
    57235794        /*Free ressources:*/
    … …  
    57715842                values[i]=vector[vertexpidlist[i]];
    57725843        }
    5773         new_input = new TriaP1Input(control_enum,values);
     5844        new_input = new TriaInput(control_enum,values);
    57745845
    57755846        if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){
    … …  
    58655936
    58665937        /*Add to inputs*/
    5867         this->inputs->AddInput(new TriaP1Input(HydrologyWaterVxEnum,vx));
    5868         this->inputs->AddInput(new TriaP1Input(HydrologyWaterVyEnum,vy));
     5938        this->inputs->AddInput(new TriaInput(HydrologyWaterVxEnum,vx));
     5939        this->inputs->AddInput(new TriaInput(HydrologyWaterVyEnum,vy));
    58695940}
    58705941/*}}}*/
    … …  
    63996470
    64006471        /*Add input to the element: */
    6401         this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));
     6472        this->inputs->AddInput(new TriaInput(WatercolumnEnum,values));
    64026473
    64036474        /*Free ressources:*/
    … …  
    64536524
    64546525        /*Add input to the element: */
    6455         this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
    6456         this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
     6526        this->inputs->AddInput(new TriaInput(SedimentHeadEnum,values));
     6527        this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual));
    64576528
    64586529        /*Free ressources:*/
    … …  
    66576728                case VertexEnum:
    66586729
    6659                         /*New TriaP1Input*/
     6730                        /*New TriaInput*/
    66606731                        IssmDouble values[3];
    66616732
    … …  
    67186789
    67196790                                        /*Add new inputs: */
    6720                                         this->inputs->AddInput(new TriaP1Input(ThicknessEnum,thickness));
    6721                                         this->inputs->AddInput(new TriaP1Input(BedEnum,bed));
    6722                                         this->inputs->AddInput(new TriaP1Input(SurfaceEnum,surface));
     6791                                        this->inputs->AddInput(new TriaInput(ThicknessEnum,thickness));
     6792                                        this->inputs->AddInput(new TriaInput(BedEnum,bed));
     6793                                        this->inputs->AddInput(new TriaInput(SurfaceEnum,surface));
    67236794
    67246795                                        break;
    67256796                                default:
    6726                                         this->inputs->AddInput(new TriaP1Input(name,values));
     6797                                        this->inputs->AddInput(new TriaInput(name,values));
    67276798                        }
    67286799                        break;
    … …  
    67766847
    67776848                                if(t==0) transientinput=new TransientInput(name);
    6778                                 transientinput->AddTimeInput(new TriaP1Input(name,values),time);
     6849                                transientinput->AddTimeInput(new TriaInput(name,values),time);
    67796850                                transientinput->Configure(parameters);
    67806851                        }

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

    r15249 r15298  
    206206                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    207207                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue,int index); //TO BE REMOVED
     208                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
     209                void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
    208210                void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    209211                void           GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);

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

    r15282 r15298  
    2626/*}}}*/
    2727/*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{*/
    28 
    2928TriaRef::TriaRef(const int nummodels){
    3029
    … …  
    222221        /*Build B': */
    223222        for (int i=0;i<NUMNODESP1;i++){
    224                 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)=2*dbasis[0][i];
    225                 *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
    226                 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
    227                 *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2*dbasis[1][i];
    228                 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i];
    229                 *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i];
     223                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)   = 2.*dbasis[0][i];
     224                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1) = dbasis[1][i];
     225                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)   = dbasis[0][i];
     226                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1) = 2.*dbasis[1][i];
     227                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)   = dbasis[1][i];
     228                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1) = dbasis[0][i];
    230229        }
    231230}
    … …  
    395394        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    396395
    397         basis[0]=gauss->coord1;
    398         basis[1]=gauss->coord2;
    399         basis[2]=gauss->coord3;
    400 
    401 }
    402 /*}}}*/
    403 /*FUNCTION TriaRef::GetNodalFunctionsP2{{{*/
    404 void TriaRef::GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss){
    405         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    406 
    407         /*Allocate*/
    408         IssmDouble* basis =xNew<IssmDouble>(NUMNODESP2);
    409 
    410         basis[0]=gauss->coord1*(2.*gauss->coord1-1.);
    411         basis[1]=gauss->coord2*(2.*gauss->coord2-1.);
    412         basis[2]=gauss->coord3*(2.*gauss->coord3-1.);
    413         basis[3]=4.*gauss->coord3*gauss->coord2;
    414         basis[4]=4.*gauss->coord3*gauss->coord1;
    415         basis[5]=4.*gauss->coord1*gauss->coord2;
    416 
    417         /*Assign output pointer*/
    418         *pbasis = basis;
    419 
     396        switch(this->element_type){
     397                case P1Enum:
     398                case P1DGEnum:
     399                        basis[0]=gauss->coord1;
     400                        basis[1]=gauss->coord2;
     401                        basis[2]=gauss->coord3;
     402                        return;
     403                case P2Enum:
     404                        basis[0]=gauss->coord1*(2.*gauss->coord1-1.);
     405                        basis[1]=gauss->coord2*(2.*gauss->coord2-1.);
     406                        basis[2]=gauss->coord3*(2.*gauss->coord3-1.);
     407                        basis[3]=4.*gauss->coord3*gauss->coord2;
     408                        basis[4]=4.*gauss->coord3*gauss->coord1;
     409                        basis[5]=4.*gauss->coord1*gauss->coord2;
     410                        return;
     411                default:
     412                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     413        }
    420414}
    421415/*}}}*/
    … …  
    424418        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    425419
    426         IssmDouble BasisFunctions[3];
    427 
    428         GetNodalFunctions(&BasisFunctions[0],gauss);
    429 
    430420        _assert_(index1>=0 && index1<3);
    431421        _assert_(index2>=0 && index2<3);
    432         basis[0]=BasisFunctions[index1];
    433         basis[1]=BasisFunctions[index2];
     422
     423        /*Fetch number of nodes for this finite element*/
     424        int numnodes = this->NumberofNodes();
     425
     426        /*Get nodal functions*/
     427        IssmDouble* BasisFunctions=xNew<IssmDouble>(numnodes);
     428        GetNodalFunctions(BasisFunctions,gauss);
     429
     430        switch(this->element_type){
     431                case P1Enum:
     432                case P1DGEnum:
     433                        basis[0]=BasisFunctions[index1];
     434                        basis[1]=BasisFunctions[index2];
     435                        return;
     436                case P2Enum:
     437                        _assert_(index2<index1);
     438                        basis[0]=BasisFunctions[index1];
     439                        basis[1]=BasisFunctions[index2];
     440                        basis[2]=BasisFunctions[3+index2-1];
     441                        return;
     442                default:
     443                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     444        }
     445
     446        xDelete<IssmDouble>(BasisFunctions);
    434447}
    435448/*}}}*/
    … …  
    439452        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    440453         * actual coordinate system): */
    441         int       i;
    442         IssmDouble    dbasis_ref[NDOF2][NUMNODESP1];
    443454        IssmDouble    Jinv[NDOF2][NDOF2];
    444455
    445         /*Get derivative values with respect to parametric coordinate system: */
    446         GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss);
     456        /*Fetch number of nodes for this finite element*/
     457        int numnodes = this->NumberofNodes();
     458
     459        /*Get nodal functions derivatives in reference triangle*/
     460        IssmDouble* dbasis_ref=xNew<IssmDouble>(2*numnodes);
     461        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss);
    447462
    448463        /*Get Jacobian invert: */
    … …  
    450465
    451466        /*Build dbasis:
    452          *
    453467         * [dhi/dx]= Jinv*[dhi/dr]
    454468         * [dhi/dy]       [dhi/ds]
    455469         */
    456         for (i=0;i<NUMNODESP1;i++){
    457                 dbasis[NUMNODESP1*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
    458                 dbasis[NUMNODESP1*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
     470        for(int i=0;i<numnodes;i++){
     471                dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
     472                dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
    459473        }
    460474
    … …  
    462476/*}}}*/
    463477/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
    464 void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss){
     478void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){
    465479        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    466480         * natural coordinate system) at the gaussian point. */
    467481
    468         /*First nodal function: */
    469         *(dl1dl3+NUMNODESP1*0+0)=-0.5;
    470         *(dl1dl3+NUMNODESP1*1+0)=-1.0/(2.0*SQRT3);
    471 
    472         /*Second nodal function: */
    473         *(dl1dl3+NUMNODESP1*0+1)=0.5;
    474         *(dl1dl3+NUMNODESP1*1+1)=-1.0/(2.0*SQRT3);
    475 
    476         /*Third nodal function: */
    477         *(dl1dl3+NUMNODESP1*0+2)=0;
    478         *(dl1dl3+NUMNODESP1*1+2)=1.0/SQRT3;
    479 
    480 }
    481 /*}}}*/
    482 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesP2Reference{{{*/
    483 void TriaRef::GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss){
    484         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    485          * natural coordinate system) at the gaussian point. */
    486 
    487         /*Allocate*/
    488         IssmDouble* dbasis =xNew<IssmDouble>(NUMNODESP2*2);
    489 
    490         /*Nodal function 1*/
    491         dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
    492         dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
    493         /*Nodal function 2*/
    494         dbasis[NUMNODESP2*0+0]=+2.*gauss->coord2 + 0.5;
    495         dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
    496         /*Nodal function 3*/
    497         dbasis[NUMNODESP2*0+0]=0.;
    498         dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
    499         /*Nodal function 4*/
    500         dbasis[NUMNODESP2*0+0]=+2.*gauss->coord3;
    501         dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
    502         /*Nodal function 5*/
    503         dbasis[NUMNODESP2*0+0]=-2.*gauss->coord3;
    504         dbasis[NUMNODESP2*0+1]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
    505         /*Nodal function 6*/
    506         dbasis[NUMNODESP2*0+0]=2.*(gauss->coord1-gauss->coord2);
    507         dbasis[NUMNODESP2*0+1]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
    508 
    509         /*Assign output pointer*/
    510         *pdbasis = dbasis;
     482        switch(this->element_type){
     483                case P1Enum: case P1DGEnum:
     484                        /*Nodal function 1*/
     485                        dbasis[NUMNODESP1*0+0]=-0.5;
     486                        dbasis[NUMNODESP1*1+0]=-1.0/(2.0*SQRT3);
     487                        /*Nodal function 2*/
     488                        dbasis[NUMNODESP1*0+1]=0.5;
     489                        dbasis[NUMNODESP1*1+1]=-1.0/(2.0*SQRT3);
     490                        /*Nodal function 3*/
     491                        dbasis[NUMNODESP1*0+2]=0;
     492                        dbasis[NUMNODESP1*1+2]=1.0/SQRT3;
     493                        return;
     494                case P2Enum:
     495                        /*Nodal function 1*/
     496                        dbasis[NUMNODESP2*0+0]=-2.*gauss->coord1 + 0.5;
     497                        dbasis[NUMNODESP2*1+0]=-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
     498                        /*Nodal function 2*/
     499                        dbasis[NUMNODESP2*0+1]=+2.*gauss->coord2 + 0.5;
     500                        dbasis[NUMNODESP2*1+1]=-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
     501                        /*Nodal function 3*/
     502                        dbasis[NUMNODESP2*0+2]=0.;
     503                        dbasis[NUMNODESP2*1+2]=+4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
     504                        /*Nodal function 4*/
     505                        dbasis[NUMNODESP2*0+3]=+2.*gauss->coord3;
     506                        dbasis[NUMNODESP2*1+3]=+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
     507                        /*Nodal function 5*/
     508                        dbasis[NUMNODESP2*0+4]=-2.*gauss->coord3;
     509                        dbasis[NUMNODESP2*1+4]=+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
     510                        /*Nodal function 6*/
     511                        dbasis[NUMNODESP2*0+5]=2.*(gauss->coord1-gauss->coord2);
     512                        dbasis[NUMNODESP2*1+5]=-2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
     513                        return;
     514                default:
     515                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     516        }
     517
    511518}
    512519/*}}}*/
    … …  
    519526         *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    520527         *
    521          * p is a vector of size 2x1 already allocated.
    522          */
    523 
    524         /*Nodal Derivatives*/
    525         IssmDouble dbasis[2][3]; //nodal derivative functions in actual coordinate system.
    526 
    527         /*Get dh1dh2dh3 in actual coordinate system: */
    528         GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss);
     528         * p is a vector already allocated.
     529         */
     530
     531        /*Output*/
     532        IssmDouble dpx=0.;
     533        IssmDouble dpy=0.;
     534
     535        /*Fetch number of nodes for this finite element*/
     536        int numnodes = this->NumberofNodes();
     537
     538        /*Get nodal functions derivatives*/
     539        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     540        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     541
     542        /*Calculate parameter for this Gauss point*/
     543        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
     544        for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
    529545
    530546        /*Assign values*/
    531         *(p+0)=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2];
    532         *(p+1)=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2];
     547        *(p+0)=dpx;
     548        *(p+1)=dpy;
    533549
    534550}
    … …  
    537553void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){
    538554
    539         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
    540          * point specifie by gauss: */
    541 
    542         /*nodal functions annd output: */
    543         IssmDouble basis[3];
     555        /*Output*/
     556        IssmDouble value =0.;
     557
     558        /*Fetch number of nodes for this finite element*/
     559        int numnodes = this->NumberofNodes();
    544560
    545561        /*Get nodal functions*/
     562        IssmDouble* basis=xNew<IssmDouble>(numnodes);
    546563        GetNodalFunctions(basis, gauss);
    547564
    548         /*Get parameter*/
    549         *p=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2];
     565        /*Calculate parameter for this Gauss point*/
     566        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
     567
     568        /*Assign output pointer*/
     569        *p = value;
    550570}
    551571/*}}}*/
    … …  
    557577                case P1DGEnum: return NUMNODESP1;
    558578                case P2Enum:   return NUMNODESP2;
    559                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     579                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    560580        }
    561581

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

    r15282 r15298  
    3636                void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
    3737                void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
    38                 void GetNodalFunctionsP2(IssmDouble** pbasis,GaussTria* gauss);
    3938                void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
    4039                void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
    … …  
    4241                void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss);
    4342                void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss);
    44                 void GetNodalFunctionsDerivativesP2Reference(IssmDouble** pdbasis,GaussTria* gauss);
    4543                void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
    4644                void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
    4745
    4846                int  NumberofNodes(void);
    49 
    5047};
    5148#endif

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

    r15104 r15298  
    3030
    3131        switch(enum_input){
    32                 case TriaP1InputEnum:
    33                         values     =new TriaP1Input(enum_type,pvalues);
    34                         savedvalues=new TriaP1Input(enum_type,pvalues);
    35                         minvalues  =new TriaP1Input(enum_type,pmin);
    36                         maxvalues  =new TriaP1Input(enum_type,pmax);
     32                case TriaInputEnum:
     33                        values     =new TriaInput(enum_type,pvalues);
     34                        savedvalues=new TriaInput(enum_type,pvalues);
     35                        minvalues  =new TriaInput(enum_type,pmin);
     36                        maxvalues  =new TriaInput(enum_type,pmax);
    3737                        break;
    3838                case PentaP1InputEnum:

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

    r15130 r15298  
    8686
    8787        /*output*/
    88         TriaP1Input* outinput=NULL;
     88        TriaInput* outinput=NULL;
    8989        IssmDouble newvalues[3];
    9090
    … …  
    100100
    101101        /*Create new Tria input*/
    102         outinput=new TriaP1Input(this->enum_type,&newvalues[0]);
     102        outinput=new TriaInput(this->enum_type,&newvalues[0]);
    103103
    104104        /*Assign output*/

  • TabularUnified issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp

    r15104 r15298  
    4141        int segment_width;
    4242        int element;
    43         int num_nodes;
     43        int numnodes;
     44        int numvertices;
    4445        int dim;
    4546        int numberofelements;
    … …  
    8889        else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");
    8990
    90         if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum)
    91          num_nodes=4;
    92         else
    93          num_nodes=2;
     91        if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
     92                numnodes=4;
     93                numvertices=4;
     94        }
     95        else{
     96                numnodes=2;
     97                numvertices=2;
     98        }
    9499
    95100        /*Fill*/
    … …  
    101106
    102107        /*Hooks: */
    103         this->hnodes=new Hook(icefront_node_ids,num_nodes);
    104         this->hvertices=new Hook(icefront_vertex_ids,num_nodes);
     108        this->hnodes=new Hook(icefront_node_ids,numnodes);
     109        this->hvertices=new Hook(icefront_vertex_ids,numvertices);
    105110        this->helement=new Hook(&icefront_eid,1);
    106111        this->hmatpar=new Hook(&icefront_mparid,1);
    … …  
    246251
    247252        /*Checks in debugging mode*/
    248         /*{{{*/
    249253        _assert_(nodes);
    250254        _assert_(element);
    251255        _assert_(matpar);
    252         /*}}}*/
    253256
    254257        /*Retrieve parameters: */
    … …  
    496499ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){
    497500
    498         /*Constants*/
    499         const int numnodes= NUMVERTICESSEG;
    500         const int numdofs = numnodes *NDOF2;
    501 
    502501        /*Intermediary*/
    503         int        ig,index1,index2,fill;
    504         IssmDouble     Jdet;
    505         IssmDouble     thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
    506         IssmDouble     water_pressure,air_pressure,surface_under_water,base_under_water;
    507         IssmDouble     xyz_list[numnodes][3];
    508         IssmDouble     normal[2];
    509         IssmDouble     L[2];
     502        int         ig,index1,index2,fill;
     503        IssmDouble  Jdet;
     504        IssmDouble  thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
     505        IssmDouble  water_pressure,air_pressure,surface_under_water,base_under_water;
     506        IssmDouble  xyz_list[NUMVERTICESSEG][3];
     507        IssmDouble  normal[2];
    510508        GaussTria *gauss;
    511509
     510        /*return of element is on water*/
    512511        Tria* tria=((Tria*)element);
    513 
    514         /*Initialize Element vector and return if necessary*/
    515512        if(tria->IsOnWater()) return NULL;
     513
     514        /*Fetch number of nodes and dof for this finite element*/
     515        //int numnodes = this->NumberofNodes();
     516        int numnodes = 2;
     517        int numdof   = numnodes*NDOF2;
     518
     519        /*Initialize Element vector and vectors*/
    516520        ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);
     521        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    517522
    518523        /*Retrieve all inputs and parameters*/
    … …  
    558563
    559564                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    560                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     565                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2);
    561566
    562567                for (int i=0;i<numnodes;i++){
    563                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i];
    564                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i];
     568                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     569                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    565570                }
    566571        }
    … …  
    570575
    571576        /*Clean up and return*/
     577        xDelete<IssmDouble>(basis);
    572578        delete gauss;
    573579        return pe;

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

    r15104 r15298  
    630630                                        IssmDouble values[3];
    631631                                        for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()];
    632                                         this->inputs->AddInput(new TriaP1Input(name,values));
     632                                        this->inputs->AddInput(new TriaInput(name,values));
    633633                                        return;
    634634                                }
    … …  
    672672                                        IssmDouble values[3];
    673673                                        for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //use sid list, to index into serial oriented vector
    674                                         this->inputs->AddInput(new TriaP1Input(name,values));
     674                                        this->inputs->AddInput(new TriaInput(name,values));
    675675                                        /*Special case for rheology B in 2D: Pourave land for this solution{{{*/
    676676                                        if(name==MaterialsRheologyBEnum){
    … …  
    685685                                                if(dim==2){
    686686                                                        /*Dupliacte rheology input: */
    687                                                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values));
     687                                                        this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values));
    688688                                                }
    689689                                        }
    … …  
    759759                if (iomodel->Data(MaterialsRheologyBEnum)) {
    760760                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    761                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs));
     761                        this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs));
    762762                }
    763763
    … …  
    765765                if (iomodel->Data(MaterialsRheologyNEnum)) {
    766766                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    767                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs));
     767                        this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs));
    768768                }
    769769
    … …  
    771771                if (iomodel->Data(MaterialsRheologyZEnum)) {
    772772                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    773                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
     773                        this->inputs->AddInput(new TriaInput(MaterialsRheologyZbarEnum,nodeinputs));
    774774                }
    775775
    … …  
    785785                                                        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];
    786786                                                        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];
    787                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     787                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    788788                                                }
    789789                                                break;
    … …  
    794794                                                        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];
    795795                                                        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];
    796                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     796                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    797797                                                }
    798798                                                break;

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

    r15293 r15298  
    1313#include "../Inputs/Input.h"
    1414#include "../Inputs/Inputs.h"
    15 #include "../Inputs/TriaP1Input.h"
     15#include "../Inputs/TriaInput.h"
    1616#include "../Inputs/PentaP1Input.h"
    1717#include "../Inputs/ControlInput.h"
    … …  
    564564                                        IssmDouble values[3];
    565565                                        for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Pid()];
    566                                         this->inputs->AddInput(new TriaP1Input(name,values));
     566                                        this->inputs->AddInput(new TriaInput(name,values));
    567567                                        return;
    568568                                }
    … …  
    606606                                        IssmDouble values[3];
    607607                                        for (int i=0;i<3;i++) values[i]=vector[((Tria*)element)->vertices[i]->Sid()]; //index into serial oriented vector
    608                                         this->inputs->AddInput(new TriaP1Input(name,values));
     608                                        this->inputs->AddInput(new TriaInput(name,values));
    609609                                        /*Special case for rheology B in 2D: Pourave land for this solution{{{*/
    610610                                        if(name==MaterialsRheologyBEnum){
    … …  
    619619                                                if(dim==2){
    620620                                                        /*Dupliacte rheology input: */
    621                                                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,values));
     621                                                        this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,values));
    622622                                                }
    623623                                        }
    … …  
    699699                if (iomodel->Data(MaterialsRheologyBEnum)) {
    700700                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
    701                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs));
     701                        this->inputs->AddInput(new TriaInput(MaterialsRheologyBbarEnum,nodeinputs));
    702702                }
    703703
    … …  
    705705                if (iomodel->Data(MaterialsRheologyNEnum)) {
    706706                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    707                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs));
     707                        this->inputs->AddInput(new TriaInput(MaterialsRheologyNEnum,nodeinputs));
    708708                }
    709709
    … …  
    719719                                                        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];
    720720                                                        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];
    721                                                         this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     721                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyBbarEnum,TriaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    722722                                                }
    723723                                                break;

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

    r15106 r15298  
    5858#include "./Inputs/IntInput.h"
    5959#include "./Inputs/PentaP1Input.h"
    60 #include "./Inputs/TriaP1Input.h"
     60#include "./Inputs/TriaInput.h"
    6161#include "./Inputs/ControlInput.h"
    6262#include "./Inputs/DatasetInput.h"

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

    r15104 r15298  
    388388        /*update static arrays*/
    389389        switch(iv){
    390                 case 0:
    391                         coord1=1; coord2=0; coord3=0;
     390                case 0: coord1=1.; coord2=0.; coord3=0.; break;
     391                case 1: coord1=0.; coord2=1.; coord3=0.; break;
     392                case 2: coord1=0.; coord2=0.; coord3=1.; break;
     393                default: _error_("vertex index should be in [0 2]");
     394        }
     395
     396}
     397/*}}}*/
     398/*FUNCTION GaussTria::GaussNode{{{*/
     399void GaussTria::GaussNode(int numnodes,int iv){
     400
     401        /*in debugging mode: check that the default constructor has been called*/
     402        _assert_(numgauss==-1);
     403
     404        /*update static arrays*/
     405        switch(numnodes){
     406                case 3: //P1 Lagrange element
     407                        switch(iv){
     408                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
     409                                case 1: coord1=0.; coord2=1.; coord3=0.; break;
     410                                case 2: coord1=0.; coord2=0.; coord3=1.; break;
     411                                default: _error_("node index should be in [0 2]");
     412                        }
    392413                        break;
    393                 case 1:
    394                         coord1=0; coord2=1; coord3=0;
     414                case 6: //P2 Lagrange element
     415                        switch(iv){
     416                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
     417                                case 1: coord1=0.; coord2=1.; coord3=0.; break;
     418                                case 2: coord1=0.; coord2=0.; coord3=1.; break;
     419                                case 3: coord1=0.; coord2=.5; coord3=.5; break;
     420                                case 4: coord1=.5; coord2=0.; coord3=.5; break;
     421                                case 5: coord1=.5; coord2=.5; coord3=0.; break;
     422                                default: _error_("node index should be in [0 5]");
     423                        }
    395424                        break;
    396                 case 2:
    397                         coord1=0; coord2=0; coord3=1;
    398                         break;
    399                 default:
    400                         _error_("vertex index should be in [0 2]");
    401 
     425                default: _error_("supported number of nodes are 3 and 6");
    402426        }
    403427

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

    r15071 r15298  
    4040                void GaussPoint(int ig);
    4141                void GaussVertex(int iv);
     42                void GaussNode(int numnodes,int iv);
    4243                void GaussCenter(void);
    4344                void GaussEdgeCenter(int index1,int index2);

  • TabularUnified issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp

    r15000 r15298  
    5151        }
    5252
     53        if(false){
     54                /*Quadratic element*/
     55                int numberofedges;
     56                int  element1,element2;
     57                bool my_edge;
     58
     59                iomodel->Constant(&numberofedges,MeshNumberofedgesEnum);
     60
     61                for (i=0;i<numberofedges;i++){
     62
     63                        /*Get left and right elements*/
     64                        element1=reCast<int>(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2]
     65                        element2=reCast<int>(iomodel->Data(MeshEdgesEnum)[4*i+3])-1; //edges are [node1 node2 elem1 elem2]
     66
     67                        /*Check whether we should include this edge (element2 is -2 for boundary edges)*/
     68                        my_edge = iomodel->my_elements[element1];
     69                        if(!my_edge && element2>=0){
     70                                my_edge = iomodel->my_elements[element2];
     71                        }
     72
     73                        /*Add node on edge*/
     74                        if(my_edge){
     75                                nodes->AddObject(new Node(iomodel->nodecounter+numberofvertices+i+1,numberofvertices+i,numberofvertices+i+1,i,iomodel,DiagnosticHorizAnalysisEnum));
     76                        }
     77                }
     78        }
     79
    5380        /*Clean fetched data: */
    5481        iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBordermacayealEnum,FlowequationBorderstokesEnum,

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

    r15282 r15298  
    384384        StringParamEnum,
    385385        TriaEnum,
    386         TriaP1InputEnum,
     386        TriaInputEnum,
    387387        VertexEnum,
    388388        VertexPIdEnum,

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

    r15282 r15298  
    384384                case StringParamEnum : return "StringParam";
    385385                case TriaEnum : return "Tria";
    386                 case TriaP1InputEnum : return "TriaP1Input";
     386                case TriaInputEnum : return "TriaInput";
    387387                case VertexEnum : return "Vertex";
    388388                case VertexPIdEnum : return "VertexPId";

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

    r15282 r15298  
    393393              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    394394              else if (strcmp(name,"Tria")==0) return TriaEnum;
    395               else if (strcmp(name,"TriaP1Input")==0) return TriaP1InputEnum;
     395              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    396396              else if (strcmp(name,"Vertex")==0) return VertexEnum;
    397397              else if (strcmp(name,"VertexPId")==0) return VertexPIdEnum;

  • TabularUnified issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r15282 r15298  
    51475147        return StringToEnum('Tria')[0]
    51485148
    5149 def TriaP1InputEnum():
    5150         """
    5151         TRIAP1INPUTENUM - Enum of TriaP1Input
    5152 
    5153         WARNING: DO NOT MODIFY THIS FILE
    5154                                 this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
    5155                                 Please read src/c/shared/Enum/README for more information
    5156 
    5157            Usage:
    5158               macro=TriaP1InputEnum()
    5159         """
    5160 
    5161         return StringToEnum('TriaP1Input')[0]
     5149def TriaInputEnum():
     5150        """
     5151        TRIAINPUTENUM - Enum of TriaInput
     5152
     5153        WARNING: DO NOT MODIFY THIS FILE
     5154                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     5155                                Please read src/c/shared/Enum/README for more information
     5156
     5157           Usage:
     5158              macro=TriaInputEnum()
     5159        """
     5160
     5161        return StringToEnum('TriaInput')[0]
    51625162
    51635163def VertexEnum():
Note: See TracChangeset for help on using the changeset viewer.