Changeset 13119


Ignore:
Timestamp:
08/22/12 07:36:45 (13 years ago)
Author:
cborstad
Message:

NEW: damage inversion capability merged from branches/trunk-jpl-damage into trunk-jpl. This affects EVERYONE. You will need to initialize a new rheology parameter md.materials.rheology_Z in all your model runs. Check the updated .par files in test/Par for syntax.

Location:
issm/trunk-jpl
Files:
40 edited
9 copied

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src

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

    r13083 r13119  
    100100        MaterialsRheologyLawEnum,
    101101        MaterialsRheologyNEnum,
     102        MaterialsRheologyZEnum,
     103        MaterialsRheologyZbarEnum,
    102104        MaterialsRhoIceEnum,
    103105        MaterialsRhoWaterEnum,
     
    361363        TemperatureOldEnum,
    362364        TemperaturePicardEnum,
     365        TemperatureSurfaceEnum,
     366        TemperatureBasalEnum,
    363367        ThicknessAbsMisfitEnum,
    364368        TypeEnum,
  • issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

    r13092 r13119  
    15551555        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    15561556        if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyBEnum);
     1557        else if (enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyZEnum);
    15571558        else input=this->inputs->GetInput(enum_type);
    15581559        //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still
     
    16671668                                        break;
    16681669                                case MaterialsRheologyBbarEnum:
     1670                                case MaterialsRheologyZbarEnum:
    16691671                                        /*Matice will take care of it*/ break;
    16701672                                default:
     
    19621964
    19631965                                /*update input*/
    1964                                 this->inputs->AddInput(new PentaP1Input(name,values));
     1966                                if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
     1967                                        matice->inputs->AddInput(new PentaP1Input(name,values));
     1968                                }
     1969                                else{
     1970                                        this->inputs->AddInput(new PentaP1Input(name,values));
     1971                                }
    19651972                                return;
    19661973                                break;
     
    32303237                        *presponse=this->matice->GetBbar();
    32313238                        break;
     3239                case MaterialsRheologyZbarEnum:
     3240                        *presponse=this->matice->GetZbar();
     3241                        break;
    32323242                case VelEnum:
    32333243                        {
     
    44604470                input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
    44614471        }
     4472        else if(enum_type==MaterialsRheologyZbarEnum){
     4473                if(!IsOnBed()) return;
     4474                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
     4475        }
     4476               
    44624477        else{
    44634478                input=inputs->GetInput(enum_type);
     
    44774492        if(enum_type==MaterialsRheologyBbarEnum){
    44784493                input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
     4494        }
     4495        else if(enum_type==MaterialsRheologyZbarEnum){
     4496                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
    44794497        }
    44804498        else{
     
    44964514        if(enum_type==MaterialsRheologyBbarEnum){
    44974515                input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
     4516        }
     4517        else if(enum_type==MaterialsRheologyZbarEnum){
     4518                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
    44984519        }
    44994520        else{
     
    45404561        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    45414562
     4563        /*Depth Averaging Z*/
     4564        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     4565
    45424566        /*Call Tria function*/
    45434567        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     
    45474571        /*Delete B averaged*/
    45484572        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     4573
     4574        /*Delete Z averaged*/
     4575        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
    45494576
    45504577        /*clean up and return*/
     
    51115138                        input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
    51125139                }
     5140                else if(control_type[i]==MaterialsRheologyZbarEnum){
     5141                        if (!IsOnBed()) goto cleanup_and_return;
     5142                        input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input);
     5143                }
    51135144                else{
    51145145                        input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
     
    51235154                if(control_type[i]==MaterialsRheologyBbarEnum){
    51245155                        this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum);
     5156                }
     5157                else if(control_type[i]==MaterialsRheologyZbarEnum){
     5158                        this->InputExtrude(MaterialsRheologyZEnum,MaterialsEnum);
    51255159                }
    51265160        }
     
    63046338        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    63056339
     6340        /*Depth Averaging Z*/
     6341        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     6342
    63066343        /*Call Tria function*/
    63076344        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     
    63116348        /*Delete B averaged*/
    63126349        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     6350
     6351        /*Delete Z averaged*/
     6352        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
    63136353
    63146354        /*clean up and return*/
     
    78607900        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    78617901
     7902        /*Depth Averaging Z*/
     7903        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     7904
    78627905        /*Call Tria function*/
    78637906        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     
    78677910        /*Delete B averaged*/
    78687911        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     7912
     7913        /*Delete Z averaged*/
     7914        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
    78697915
    78707916        /*clean up and return*/
  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp

    r13092 r13119  
    491491        /*Initialize Element matrix*/
    492492        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     493
    493494        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
    494495
     
    13951396
    13961397        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1397         if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
     1398        if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type);
    13981399        else input=this->inputs->GetInput(enum_type);
    13991400        //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs");
     
    15011502                                        break;
    15021503                                case MaterialsRheologyBbarEnum:
     1504                                case MaterialsRheologyZbarEnum:
    15031505                                        /*Matice will take care of it*/ break;
    15041506                                default:
     
    16891691
    16901692                        /*update input*/
    1691                         if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
     1693                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    16921694                                matice->inputs->AddInput(new TriaP1Input(name,values));
    16931695                        }
     
    18401842                                name==FrictionCoefficientEnum ||
    18411843                                name==MaterialsRheologyBbarEnum ||
     1844                                name==MaterialsRheologyZbarEnum ||
    18421845                                name==GradientEnum ||
    18431846                                name==OldGradientEnum ||
     
    27912794                        *presponse=this->matice->GetBbar();
    27922795                        break;
     2796                case MaterialsRheologyZbarEnum:
     2797                        *presponse=this->matice->GetZbar();
     2798                        break;
    27932799                case VelEnum:{
    27942800
     
    34023408        for(int i=0;i<num_controls;i++){
    34033409
    3404                 if(control_type[i]==MaterialsRheologyBbarEnum){
     3410                if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
    34053411                        input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
    34063412                }
     
    34293435        Input* input=NULL;
    34303436
    3431         if(enum_type==MaterialsRheologyBbarEnum){
     3437        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34323438                input=(Input*)matice->inputs->GetInput(enum_type);
    34333439        }
     
    34473453        Input* input=NULL;
    34483454
    3449         if(enum_type==MaterialsRheologyBbarEnum){
     3455        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34503456                input=(Input*)matice->inputs->GetInput(enum_type);
    34513457        }
     
    34663472        Input* input=NULL;
    34673473
    3468         if(enum_type==MaterialsRheologyBbarEnum){
     3474        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34693475                input=(Input*)matice->inputs->GetInput(enum_type);
    34703476        }
     
    34963502                case MaterialsRheologyBbarEnum:
    34973503                        GradjBMacAyeal(gradient,control_index);
     3504                        break;
     3505                case MaterialsRheologyZbarEnum:
     3506                        GradjZMacAyeal(gradient,control_index);
    34983507                        break;
    34993508                case BalancethicknessThickeningRateEnum:
     
    35833592}
    35843593/*}}}*/
     3594/*FUNCTION Tria::GradjZGradient{{{*/
     3595void  Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){
     3596
     3597        int        i,ig;
     3598        int        doflist1[NUMVERTICES];
     3599        IssmDouble     Jdet,weight;
     3600        IssmDouble     xyz_list[NUMVERTICES][3];
     3601        IssmDouble     dbasis[NDOF2][NUMVERTICES];
     3602        IssmDouble     dk[NDOF2];
     3603        IssmDouble     grade_g[NUMVERTICES]={0.0};
     3604        GaussTria  *gauss=NULL;
     3605
     3606        /*Retrieve all inputs we will be needing: */
     3607        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3608        GradientIndexing(&doflist1[0],control_index);
     3609        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3610        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
     3611
     3612        /* Start  looping on the number of gaussian points: */
     3613        gauss=new GaussTria(2);
     3614        for (ig=gauss->begin();ig<gauss->end();ig++){
     3615
     3616                gauss->GaussPoint(ig);
     3617
     3618                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3619                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     3620                weights_input->GetInputValue(&weight,gauss,weight_index);
     3621
     3622                /*Build alpha_complement_list: */
     3623                rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     3624
     3625                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     3626                for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
     3627        }
     3628        gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
     3629
     3630        /*Clean up and return*/
     3631        delete gauss;
     3632}
     3633/*}}}*/
    35853634/*FUNCTION Tria::GradjBMacAyeal{{{*/
    35863635void  Tria::GradjBMacAyeal(Vector* gradient,int control_index){
     
    36243673                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    36253674                matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
     3675
     3676                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3677                GetNodalFunctions(basis,gauss);
     3678
     3679                /*standard gradient dJ/dki*/
     3680                for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
     3681                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     3682                                        )*Jdet*gauss->weight*basis[i];
     3683        }
     3684
     3685        gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
     3686
     3687        /*clean-up*/
     3688        delete gauss;
     3689}
     3690/*}}}*/
     3691/*FUNCTION Tria::GradjZMacAyeal{{{*/
     3692void  Tria::GradjZMacAyeal(Vector* gradient,int control_index){
     3693
     3694        /*Intermediaries*/
     3695        int        i,ig;
     3696        int        doflist[NUMVERTICES];
     3697        IssmDouble     vx,vy,lambda,mu,thickness,Jdet;
     3698        IssmDouble     viscosity_complement;
     3699        IssmDouble     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2];
     3700        IssmDouble     xyz_list[NUMVERTICES][3];
     3701        IssmDouble     basis[3],epsilon[3];
     3702        IssmDouble     grad[NUMVERTICES]={0.0};
     3703        GaussTria *gauss = NULL;
     3704
     3705        /* Get node coordinates and dof list: */
     3706        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3707        GradientIndexing(&doflist[0],control_index);
     3708
     3709        /*Retrieve all inputs*/
     3710        Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
     3711        Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
     3712        Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
     3713        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
     3714        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
     3715        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3716
     3717        /* Start  looping on the number of gaussian points: */
     3718        gauss=new GaussTria(4);
     3719        for (ig=gauss->begin();ig<gauss->end();ig++){
     3720
     3721                gauss->GaussPoint(ig);
     3722
     3723                thickness_input->GetInputValue(&thickness,gauss);
     3724                rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
     3725                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     3726                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     3727                adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
     3728                adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
     3729
     3730                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3731                matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
    36263732
    36273733                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     
    51095215
    51105216        /*Get input (either in element or material)*/
    5111         if(control_enum==MaterialsRheologyBbarEnum){
     5217        if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){
    51125218                input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
    51135219        }
     
    51445250        new_input = new TriaP1Input(control_enum,values);
    51455251
    5146         if(control_enum==MaterialsRheologyBbarEnum){
     5252        if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){
    51475253                input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
    51485254        }
  • issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h

    r13092 r13119  
    143143                void   Gradj(Vector* gradient,int control_type,int control_index);
    144144                void   GradjBGradient(Vector* gradient,int weight_index,int control_index);
     145                void   GradjZGradient(Vector* gradient,int weight_index,int control_index);
    145146                void   GradjBMacAyeal(Vector* gradient,int control_index);
     147                void   GradjZMacAyeal(Vector* gradient,int control_index);
    146148                void   GradjDragMacAyeal(Vector* gradient,int control_index);
    147149                void   GradjDragStokes(Vector* gradient,int control_index);
  • issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp

    r13073 r13119  
    172172}
    173173/*}}}*/
     174/*FUNCTION Matice::GetZ {{{*/
     175IssmDouble Matice::GetZ(){
     176
     177        /*Output*/
     178        IssmDouble Z;
     179
     180        inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
     181        return Z;
     182}
     183/*}}}*/
     184/*FUNCTION Matice::GetZbar {{{*/
     185IssmDouble Matice::GetZbar(){
     186
     187        /*Output*/
     188        IssmDouble Zbar;
     189
     190        inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
     191        return Zbar;
     192}
     193/*}}}*/
    174194/*FUNCTION Matice::GetVectorFromInputs{{{*/
    175195void  Matice::GetVectorFromInputs(Vector* vector,int input_enum){
     
    207227void  Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
    208228        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    209                                                                                                     B
     229                                                                                                   Z * B
    210230          viscosity= -------------------------------------------------------------------
    211231                                                  2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
     
    226246        /*Intermediary: */
    227247        IssmDouble A,e;
    228         IssmDouble B,n;
     248        IssmDouble Btmp,B,n,Z;
    229249
    230250        /*Get B and n*/
    231         B=GetBbar();
     251        Btmp=GetBbar();
     252        Z=GetZbar();
    232253        n=GetN();
     254        B=Z*Btmp;
    233255
    234256        if (n==1){
     
    292314        /*Intermediaries: */
    293315        IssmDouble A,e;
    294         IssmDouble B,n;
    295 
    296         /*Get B and n*/
    297         B=GetB();
     316        IssmDouble B,n,Z;
     317
     318        /*Get B, Z and n*/
    298319        n=GetN();
     320        Z=GetZ();
     321        B=Z*GetB();
    299322
    300323        if (n==1){
     
    362385        /*Intermediaries: */
    363386        IssmDouble A,e;
    364         IssmDouble B,n;
     387        IssmDouble B,n,Z;
    365388        IssmDouble eps0;
    366389
    367390        /*Get B and n*/
    368391        eps0=pow((IssmDouble)10,(IssmDouble)-27);
    369         B=GetB();
    370392        n=GetN();
     393        Z=GetZ();
     394        B=Z*GetB();
    371395       
    372396        if (n==1){
     
    452476               
    453477                        viscosity_complement=1/(2*pow(A,e));
     478                }
     479        }
     480        else{
     481                viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
     482        }
     483
     484        /*Checks in debugging mode*/
     485        _assert_(B>0);
     486        _assert_(n>0);
     487        _assert_(viscosity_complement>0);
     488               
     489        /*Return: */
     490        *pviscosity_complement=viscosity_complement;
     491}
     492/*}}}*/
     493/*FUNCTION Matice::GetViscosityZComplement {{{*/
     494void  Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
     495        /*Return viscosity derivative for control method d(mu)/dZ:
     496         *
     497         *                                                                                             B
     498         * dviscosity= -------------------------------------------------------------------
     499         *                                2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
     500         *
     501         * If epsilon is NULL, it means this is the first time Gradjb is being run, and we
     502         * return mu20, initial viscosity.
     503         */
     504       
     505        /*output: */
     506        IssmDouble viscosity_complement;
     507
     508        /*input strain rate: */
     509        IssmDouble exx,eyy,exy;
     510
     511        /*Intermediary value A and exponent e: */
     512        IssmDouble A,e;
     513        IssmDouble B,n;
     514
     515        /*Get B and n*/
     516        B=GetBbar();
     517        n=GetN();
     518
     519        if(epsilon){
     520                exx=*(epsilon+0);
     521                eyy=*(epsilon+1);
     522                exy=*(epsilon+2);
     523
     524                /*Build viscosity: mu2=B/(2*A^e) */
     525                A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
     526                if(A==0){
     527                        /*Maximum viscosity_complement for 0 shear areas: */
     528                        viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
     529                }
     530                else{
     531                        e=(n-1)/(2*n);
     532               
     533                        viscosity_complement=B/(2*pow(A,e));
    454534                }
    455535        }
     
    701781                }
    702782
     783                /*Get Z*/
     784                if (iomodel->Data(MaterialsRheologyZEnum)) {
     785                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     786                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
     787                }
     788
    703789                /*Control Inputs*/
    704790                #ifdef _HAVE_CONTROL_
     
    715801                                                }
    716802                                                break;
     803                                        case MaterialsRheologyZbarEnum:
     804                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     805                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     806                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     807                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
     808                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
     809                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     810                                                }
     811                                                break;
     812
    717813                                }
    718814                        }
     
    741837                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    742838                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
     839                }
     840
     841                /*Get Z*/
     842                if (iomodel->Data(MaterialsRheologyZEnum)) {
     843                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     844                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
    743845                }
    744846
     
    757859                                                }
    758860                                                break;
     861                                        case MaterialsRheologyZbarEnum:
     862                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     863                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     864                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     865                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
     866                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
     867                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     868                                                }
     869                                                break;
    759870                                }
    760871                        }
     
    775886                                name==MaterialsRheologyBEnum ||
    776887                                name==MaterialsRheologyBbarEnum ||
    777                                 name==MaterialsRheologyNEnum
     888                                name==MaterialsRheologyNEnum ||
     889                                name==MaterialsRheologyZEnum ||
     890                                name==MaterialsRheologyZbarEnum
    778891                ){
    779892                return true;
  • issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h

    r13051 r13119  
    5858                /*}}}*/
    5959                /*Matice Numerics: {{{*/
    60                 void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    61                 void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
    62                 void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
    63                 void       GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
    64                 void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
    65                 void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    66                 void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     60                void   SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin);
     61                void   GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
     62                void   GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
     63                void   GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
     64                void   GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
     65                void   GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
     66                void   GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
     67                void   GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6768                IssmDouble GetA();
    6869                IssmDouble GetB();
    6970                IssmDouble GetBbar();
    7071                IssmDouble GetN();
    71                 bool       IsInput(int name);
     72                IssmDouble GetZ();
     73                IssmDouble GetZbar();
     74                bool   IsInput(int name);
    7275                /*}}}*/
    7376};
  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r13083 r13119  
    105105                case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
    106106                case MaterialsRheologyNEnum : return "MaterialsRheologyN";
     107                case MaterialsRheologyZEnum : return "MaterialsRheologyZ";
     108                case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar";
    107109                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    108110                case MaterialsRhoWaterEnum : return "MaterialsRhoWater";
     
    352354                case TemperatureOldEnum : return "TemperatureOld";
    353355                case TemperaturePicardEnum : return "TemperaturePicard";
     356                case TemperatureSurfaceEnum : return "TemperatureSurface";
     357                case TemperatureBasalEnum : return "TemperatureBasal";
    354358                case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit";
    355359                case TypeEnum : return "Type";
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r13073 r13119  
    4848                        case FrictionCoefficientEnum: iomodel->FetchData(1,FrictionCoefficientEnum); break;
    4949                        case MaterialsRheologyBbarEnum:    iomodel->FetchData(1,MaterialsRheologyBEnum); break;
     50                        case MaterialsRheologyZbarEnum:    iomodel->FetchData(1,MaterialsRheologyZEnum); break;
    5051                        default: _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
    5152                }
     
    6667       
    6768        /*Free data: */
    68         iomodel->DeleteData(1+4+5,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum);
     69        iomodel->DeleteData(1+4+6,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,MaterialsRheologyZEnum);
    6970}
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r12832 r13119  
    4545       
    4646        /*Fetch data needed: */
    47         iomodel->FetchData(4,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum);
     47        iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum);
    4848        #ifdef _HAVE_3D_
    4949        if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum);
     
    6767       
    6868        /*Free data: */
    69         iomodel->DeleteData(9,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum,
    70                                 MaterialsRheologyBEnum,MaterialsRheologyNEnum,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum);
     69        iomodel->DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum,
     70                                MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum,
     71                                InversionMaxParametersEnum);
    7172
    72         /*Add new constrant material property tgo materials, at the end: */
     73        /*Add new constant material property to materials, at the end: */
    7374        materials->AddObject(new Matpar(numberofelements+1,iomodel));//put it at the end of the materials
    7475       
     
    8182        for (i=0;i<numberofvertices;i++){
    8283
    83                 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
     84                /*vertices and nodes (same number, as we are running continuous galerkin formulation): */
    8485                if(iomodel->my_vertices[i]){
    8586                       
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r13051 r13119  
    6161        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
    6262        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     63        iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);
    6364        iomodel->FetchDataToInput(elements,VxEnum);
    6465        iomodel->FetchDataToInput(elements,VyEnum);
  • issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp

    r13056 r13119  
    4747        }
    4848
    49         /*Assign output pointers: */
     49        /*Assign output pointers:*/
    5050        *puf=uf;
    5151}
  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r13083 r13119  
    106106              else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
    107107              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
     108              else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum;
     109              else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum;
    108110              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    109111              else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
     
    136138              else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum;
    137139              else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
    138               else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
    139               else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
    140140         else stage=2;
    141141   }
    142142   if(stage==2){
    143               if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
     143              if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
     144              else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
     145              else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
    144146              else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
    145147              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
     
    259261              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    260262              else if (strcmp(name,"Contour")==0) return ContourEnum;
    261               else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    262               else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    263263         else stage=3;
    264264   }
    265265   if(stage==3){
    266               if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
     266              if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
     267              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
     268              else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
    267269              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    268270              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
     
    359361              else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
    360362              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
     363              else if (strcmp(name,"TemperatureSurface")==0) return TemperatureSurfaceEnum;
     364              else if (strcmp(name,"TemperatureBasal")==0) return TemperatureBasalEnum;
    361365              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    362366              else if (strcmp(name,"Type")==0) return TypeEnum;
     
    380384              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    381385              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    382               else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
     386         else stage=4;
     387   }
     388   if(stage==4){
     389              if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    383390              else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    384391              else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
    385392              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    386          else stage=4;
    387    }
    388    if(stage==4){
    389               if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
     393              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    390394              else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
    391395              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
  • issm/trunk-jpl/src/m/classes/initialization.m

    r13020 r13119  
    1212                pressure      = NaN;
    1313                temperature   = NaN;
     14                surfacetemp   = NaN;
     15                basaltemp     = NaN;
    1416                watercolumn   = NaN;
    1517                waterfraction = NaN;
     
    6870                        fielddisplay(obj,'pressure','pressure field');
    6971                        fielddisplay(obj,'temperature','temperature in Kelvins');
     72                        fielddisplay(obj,'surfacetemp','surface temperature in Kelvins');
     73                        fielddisplay(obj,'basaltemp','basal temperature in Kelvins');
    7074                        fielddisplay(obj,'watercolumn','thickness of subglacial water');
    7175                        fielddisplay(obj,'waterfraction','fraction of water in the ice');
     
    7377                end % }}}
    7478                function marshall(obj,fid) % {{{
    75                         WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum());
    76                         WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum());
    77                         WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum());
    78                         WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum());
    79                         WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum());
    80                         WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum());
    81                         WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum());
     79                        WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum);
     80                        WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum);
     81                        WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum);
     82                        WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum);
     83                        WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum);
     84                        WriteData(fid,'data',obj.surfacetemp,'format','DoubleMat','mattype',1,'enum',TemperatureSurfaceEnum);
     85                        WriteData(fid,'data',obj.basaltemp,'format','DoubleMat','mattype',1,'enum',TemperatureBasalEnum);
     86                        WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);
     87                        WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);
    8288                end % }}}
    8389        end
  • issm/trunk-jpl/src/m/classes/inversion.m

    r13059 r13119  
    8383                        num_costfunc=size(md.inversion.cost_functions,2);
    8484
    85                         md = checkfield(md,'inversion.iscontrol','values',[0 1]);
    86                         md = checkfield(md,'inversion.tao','values',[0 1]);
    87                         md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
    88                         md = checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'});
    89                         md = checkfield(md,'inversion.nsteps','numel',[1],'>=',1);
    90                         md = checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
    91                         md = checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]);
    92                         md = checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:505]);
    93                         md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
    94                         md = checkfield(md,'inversion.gradient_only','values',[0 1]);
    95                         md = checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]);
    96                         md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
    97                         md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
     85                        checkfield(md,'inversion.iscontrol','values',[0 1]);
     86                        checkfield(md,'inversion.tao','values',[0 1]);
     87                        checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
     88                        checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy'});
     89                        checkfield(md,'inversion.nsteps','numel',1,'>=',1);
     90                        checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
     91                        checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]);
     92                        checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:503]);
     93                        checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
     94                        checkfield(md,'inversion.gradient_only','values',[0 1]);
     95                        checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]);
     96                        checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
     97                        checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
    9898
    9999                        if solution==BalancethicknessSolutionEnum()
  • issm/trunk-jpl/src/m/classes/materials.m

    r13094 r13119  
    1919                rheology_B   = NaN;
    2020                rheology_n   = NaN;
     21                rheology_Z   = NaN;
    2122                rheology_law = '';
    2223        end
     
    7677                        md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
    7778                        md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
     79                        md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]);
    7880                        md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'});
    7981                end % }}}
     
    9496                        fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]');
    9597                        fielddisplay(obj,'rheology_n','Glen''s flow law exponent');
     98                        fielddisplay(obj,'rheology_Z','rheology multiplier');
    9699                        fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius''');
    97100                end % }}}
     
    110113                        WriteData(fid,'object',obj,'fieldname','rheology_B','format','DoubleMat','mattype',1);
    111114                        WriteData(fid,'object',obj,'fieldname','rheology_n','format','DoubleMat','mattype',2);
     115                        WriteData(fid,'object',obj,'fieldname','rheology_Z','format','DoubleMat','mattype',1);
    112116                        WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer');
    113117                end % }}}
  • issm/trunk-jpl/src/m/classes/model/model.m

    r13043 r13119  
    168168                         md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
    169169                         md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
     170                         md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z);
    170171
    171172                         %special for thermal modeling:
     
    716717                         md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
    717718                         md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
     719                         md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node');
    718720
    719721                         %parameters
     
    829831                         if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
    830832                         if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
     833                         if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end
    831834                         if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
    832835                         if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
     836                         if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end
    833837                         if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end
    834838                         if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end
  • issm/trunk-jpl/src/m/contrib/hack/tres.m

    r13021 r13119  
    105105                case MaterialsRheologyBEnum(), string='rheology_B'; return
    106106                case MaterialsRheologyBbarEnum(), string='rheology_B'; return
     107                case MaterialsRheologyZEnum(), string='rheology_Z'; return
     108                case MaterialsRheologyZbarEnum(), string='rheology_Z'; return
    107109                case BalancethicknessThickeningRateEnum(), string='dhdt'; return
    108110                case VxEnum(), string='vx'; return
  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r13083 r13119  
    889889        return StringToEnum('MaterialsRheologyN')[0]
    890890
     891def MaterialsRheologyZEnum():
     892        """
     893        MATERIALSRHEOLOGYZENUM - Enum of MaterialsRheologyZ
     894
     895           Usage:
     896              macro=MaterialsRheologyZEnum()
     897        """
     898
     899        return StringToEnum('MaterialsRheologyZ')[0]
     900
     901def MaterialsRheologyZbarEnum():
     902        """
     903        MATERIALSRHEOLOGYZBARENUM - Enum of MaterialsRheologyZbar
     904
     905           Usage:
     906              macro=MaterialsRheologyZbarEnum()
     907        """
     908
     909        return StringToEnum('MaterialsRheologyZbar')[0]
     910
    891911def MaterialsRhoIceEnum():
    892912        """
     
    33593379        return StringToEnum('TemperaturePicard')[0]
    33603380
     3381def TemperatureSurfaceEnum():
     3382        """
     3383        TEMPERATURESURFACEENUM - Enum of TemperatureSurface
     3384
     3385           Usage:
     3386              macro=TemperatureSurfaceEnum()
     3387        """
     3388
     3389        return StringToEnum('TemperatureSurface')[0]
     3390
     3391def TemperatureBasalEnum():
     3392        """
     3393        TEMPERATUREBASALENUM - Enum of TemperatureBasal
     3394
     3395           Usage:
     3396              macro=TemperatureBasalEnum()
     3397        """
     3398
     3399        return StringToEnum('TemperatureBasal')[0]
     3400
    33613401def ThicknessAbsMisfitEnum():
    33623402        """
     
    46074647        """
    46084648
    4609         return 459
    4610 
     4649        return 463
     4650
  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r13083 r13119  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=459;
     11macro=463;
  • issm/trunk-jpl/test

  • issm/trunk-jpl/test/NightlyRun/IdToName.m

    r12883 r13119  
    6565        case 236, name='SquareShelfTranIspddIsdeltaM2d';
    6666        case 237, name='SquareShelfTranIspddIsdeltaM3d';
     67        case 270, name='SquareShelfDiagM2dDamage';
     68        case 272, name='SquareShelfCMZM2dDamage';
     69        case 274, name='SquareShelfDiagM2dDamageRift';
    6770        case 301, name='SquareSheetConstrainedDiagM2d';
    6871        case 302, name='SquareSheetConstrainedDiagH2d';
  • issm/trunk-jpl/test/NightlyRun/README

    r5100 r13119  
    31316: 79North
    3232Add the id and testname in IdToName.m (incresing order)
    33 We try not to create to many .par and .exp files, so try to use the existing ones as much as possible.
     33We try not to create too many .par and .exp files, so try to use the existing ones as much as possible.
    3434To modify some characteristics, do it in the testxxx.m file.
    3535The testxxx_nightly.m is used to define the parameters you want to check in the nightlyruns.
  • issm/trunk-jpl/test/Par/79North.par

    r9734 r13119  
    1313md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
    1414md.materials.rheology_B=paterson(md.initialization.temperature);
     15md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1516md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    1617
  • issm/trunk-jpl/test/Par/ISMIPA.par

    r9734 r13119  
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    1818md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     19md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1920
    2021disp('      boundary conditions for diagnostic model');
  • issm/trunk-jpl/test/Par/ISMIPB.par

    r9734 r13119  
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    1818md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     19md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1920
    2021disp('      boundary conditions for diagnostic model');
  • issm/trunk-jpl/test/Par/ISMIPC.par

    r9734 r13119  
    1818md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    1919md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     20md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2021
    2122disp('      boundary conditions for diagnostic model: ');
  • issm/trunk-jpl/test/Par/ISMIPD.par

    r9734 r13119  
    1717md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    1818md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     19md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1920
    2021disp('      boundary conditions for diagnostic model: ');
  • issm/trunk-jpl/test/Par/ISMIPE.par

    r9734 r13119  
    2626md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
    2727md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     28md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2829
    2930disp('      boundary conditions for diagnostic model: ');
  • issm/trunk-jpl/test/Par/ISMIPF.par

    r9764 r13119  
    1616md.materials.rheology_B=1.4734*10^14*ones(md.mesh.numberofvertices,1);
    1717md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
     18md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1819md.materials.rheology_law='None';
    1920
  • issm/trunk-jpl/test/Par/Pig.par

    r9734 r13119  
    1717md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
    1818md.materials.rheology_B=paterson(md.initialization.temperature);
     19md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    1920md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    2021md.initialization.temperature=md.initialization.temperature;
  • issm/trunk-jpl/test/Par/RoundSheetEISMINT.par

    r9734 r13119  
    2020md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution
    2121md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     22md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2223
    2324disp('      creating surface mass balance');
  • issm/trunk-jpl/test/Par/RoundSheetShelf.par

    r10304 r13119  
    6060md.materials.rheology_B=paterson(md.initialization.temperature);
    6161md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     62md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    6263
    6364%Surface mass balance and basal melting
  • issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par

    r9734 r13119  
    2525md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution
    2626md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     27md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2728
    2829disp('      creating surface mass balance');
  • issm/trunk-jpl/test/Par/SquareEISMINT.par

    r9864 r13119  
    2727md.materials.rheology_B=1.7687*10^8*ones(md.mesh.numberofvertices,1);
    2828md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     29md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2930
    3031disp('      creating surface mass balance');
  • issm/trunk-jpl/test/Par/SquareSheetConstrained.par

    r9734 r13119  
    2222md.materials.rheology_B=paterson(md.initialization.temperature);
    2323md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     24md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2425
    2526%Friction
  • issm/trunk-jpl/test/Par/SquareSheetShelf.par

    r9734 r13119  
    2525md.materials.rheology_B=paterson(md.initialization.temperature);
    2626md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     27md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2728
    2829%Accumulation and melting
  • issm/trunk-jpl/test/Par/SquareShelf.par

    r9734 r13119  
    2222md.materials.rheology_B=paterson(md.initialization.temperature);
    2323md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     24md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2425
    2526%Friction
  • issm/trunk-jpl/test/Par/SquareShelfConstrained.par

    r12853 r13119  
    2222md.materials.rheology_B=paterson(md.initialization.temperature);
    2323md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     24md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    2425
    2526%Surface mass balance and basal melting
  • issm/trunk-jpl/test/Par/SquareThermal.par

    r9733 r13119  
    2828md.materials.rheology_B=paterson(md.initialization.temperature);
    2929md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     30md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
    3031
    3132disp('      creating surface mass balance');
Note: See TracChangeset for help on using the changeset viewer.