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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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                        }
Note: See TracChangeset for help on using the changeset viewer.