Changeset 11462


Ignore:
Timestamp:
02/16/12 17:48:24 (13 years ago)
Author:
cborstad
Message:

patched in changes from previous branch

Location:
issm/branches/trunk-jpl-damage/src
Files:
3 added
18 edited

Legend:

Unmodified
Added
Removed

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

    r11347 r11462  
    100100        MaterialsRheologyLawEnum,
    101101        MaterialsRheologyNEnum,
     102        MaterialsRheologyZEnum,
     103        MaterialsRheologyZbarEnum,
    102104        MaterialsRhoIceEnum,
    103105        MaterialsRhoWaterEnum,

  • issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumToModelField.cpp

    r11199 r11462  
    1717                case MaterialsRheologyBEnum : return "rheology_B";
    1818                case MaterialsRheologyBbarEnum : return "rheology_B";
     19                case MaterialsRheologyZEnum : return "rheology_Z";
     20                case MaterialsRheologyZbarEnum : return "rheology_Z";
    1921                case BalancethicknessThickeningRateEnum: return "dhdt";
    2022                case VxEnum : return "vx";

  • issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r11347 r11462  
    104104                case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
    105105                case MaterialsRheologyNEnum : return "MaterialsRheologyN";
     106                case MaterialsRheologyZEnum : return "MaterialsRheologyZ";
     107                case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar";
    106108                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    107109                case MaterialsRhoWaterEnum : return "MaterialsRhoWater";

  • issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

    r10970 r11462  
    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 %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
    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/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

    r11001 r11462  
    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_THREED_
    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
    7273        /*Add new constrant material property tgo materials, at the end: */

  • issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

    r11000 r11462  
    5959        iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
    6060        iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
     61        iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);
    6162        iomodel->FetchDataToInput(elements,VxEnum);
    6263        iomodel->FetchDataToInput(elements,VyEnum);

  • issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r11407 r11462  
    105105              else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
    106106              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
     107              else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum;
     108              else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum;
    107109              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    108110              else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
    … …  
    135137              else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
    136138              else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
    137               else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
    138               else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
     142              if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
     143              else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
     144              else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
    143145              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
    144146              else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum;
    … …  
    258260              else if (strcmp(name,"Node")==0) return NodeEnum;
    259261              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    260               else if (strcmp(name,"Param")==0) return ParamEnum;
    261               else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"Pengrid")==0) return PengridEnum;
     265              if (strcmp(name,"Param")==0) return ParamEnum;
     266              else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
     267              else if (strcmp(name,"Pengrid")==0) return PengridEnum;
    266268              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
    267269              else if (strcmp(name,"Penta")==0) return PentaEnum;
    … …  
    381383              else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum;
    382384              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
    383               else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
    384               else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     388              if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
     389              else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum;
     390              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    389391              else if (strcmp(name,"J")==0) return JEnum;
    390392              else if (strcmp(name,"Patch")==0) return PatchEnum;

  • issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp

    r11332 r11462  
    14731473
    14741474        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1475         if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
     1475        if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type);
    14761476        else input=this->inputs->GetInput(enum_type);
    14771477        //if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type));
    … …  
    15791579                                        break;
    15801580                                case MaterialsRheologyBbarEnum:
     1581                                case MaterialsRheologyZbarEnum:
    15811582                                        /*Matice will take care of it*/ break;
    15821583                                default:
    … …  
    17671768
    17681769                        /*update input*/
    1769                         if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
     1770                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    17701771                                matice->inputs->AddInput(new TriaP1Input(name,values));
    17711772                        }
    … …  
    27682769                        *presponse=this->matice->GetBbar();
    27692770                        break;
     2771                case MaterialsRheologyZbarEnum:
     2772                        *presponse=this->matice->GetZbar();
     2773                        break;
    27702774                case VelEnum:
    27712775
    … …  
    33463350        for(int i=0;i<num_controls;i++){
    33473351
    3348                 if(control_type[i]==MaterialsRheologyBbarEnum){
     3352                if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
    33493353                        input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
    33503354                }
    … …  
    33733377        Input* input=NULL;
    33743378
    3375         if(enum_type==MaterialsRheologyBbarEnum){
     3379        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    33763380                input=(Input*)matice->inputs->GetInput(enum_type);
    33773381        }
    … …  
    33913395        Input* input=NULL;
    33923396
    3393         if(enum_type==MaterialsRheologyBbarEnum){
     3397        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    33943398                input=(Input*)matice->inputs->GetInput(enum_type);
    33953399        }
    … …  
    34103414        Input* input=NULL;
    34113415
    3412         if(enum_type==MaterialsRheologyBbarEnum){
     3416        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34133417                input=(Input*)matice->inputs->GetInput(enum_type);
    34143418        }
    … …  
    34313435
    34323436        /*If on water, grad = 0: */
    3433         if(IsOnWater()) return;
     3437        if(IsOnWater()) return; 
    34343438
    34353439        /*First deal with ∂/∂alpha(KU-F)*/
    … …  
    34403444                case MaterialsRheologyBbarEnum:
    34413445                        GradjBMacAyeal(gradient,control_index);
     3446                        break;
     3447                case MaterialsRheologyZbarEnum:
     3448                        GradjZMacAyeal(gradient,control_index);
    34423449                        break;
    34433450                case BalancethicknessThickeningRateEnum:
    … …  
    35253532}
    35263533/*}}}*/
     3534/*FUNCTION Tria::GradjZGradient{{{1*/
     3535void  Tria::GradjZGradient(Vec gradient,int weight_index,int control_index){
     3536
     3537        int        i,ig;
     3538        int        doflist1[NUMVERTICES];
     3539        double     Jdet,weight;
     3540        double     xyz_list[NUMVERTICES][3];
     3541        double     dbasis[NDOF2][NUMVERTICES];
     3542        double     dk[NDOF2];
     3543        double     grade_g[NUMVERTICES]={0.0};
     3544        GaussTria  *gauss=NULL;
     3545
     3546        /*Retrieve all inputs we will be needing: */
     3547        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3548        GradientIndexing(&doflist1[0],control_index);
     3549        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3550        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
     3551
     3552        /* Start  looping on the number of gaussian points: */
     3553        gauss=new GaussTria(2);
     3554        for (ig=gauss->begin();ig<gauss->end();ig++){
     3555
     3556                gauss->GaussPoint(ig);
     3557
     3558                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3559                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     3560                weights_input->GetInputValue(&weight,gauss,weight_index);
     3561
     3562                /*Build alpha_complement_list: */
     3563                rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     3564
     3565                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     3566                for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
     3567        }
     3568        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     3569
     3570        /*Clean up and return*/
     3571        delete gauss;
     3572}
     3573/*}}}*/
    35273574/*FUNCTION Tria::GradjBMacAyeal{{{1*/
    35283575void  Tria::GradjBMacAyeal(Vec gradient,int control_index){
    … …  
    35663613                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    35673614                matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
     3615
     3616                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3617                GetNodalFunctions(basis,gauss);
     3618
     3619                /*standard gradient dJ/dki*/
     3620                for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
     3621                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     3622                                        )*Jdet*gauss->weight*basis[i];
     3623        }
     3624
     3625        VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES);
     3626
     3627        /*clean-up*/
     3628        delete gauss;
     3629}
     3630/*}}}*/
     3631/*FUNCTION Tria::GradjZMacAyeal{{{1*/
     3632void  Tria::GradjZMacAyeal(Vec gradient,int control_index){
     3633
     3634        /*Intermediaries*/
     3635        int        i,ig;
     3636        int        doflist[NUMVERTICES];
     3637        double     vx,vy,lambda,mu,thickness,Jdet;
     3638        double     viscosity_complement;
     3639        double     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2];
     3640        double     xyz_list[NUMVERTICES][3];
     3641        double     basis[3],epsilon[3];
     3642        double     grad[NUMVERTICES]={0.0};
     3643        GaussTria *gauss = NULL;
     3644
     3645        /* Get node coordinates and dof list: */
     3646        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3647        GradientIndexing(&doflist[0],control_index);
     3648
     3649        /*Retrieve all inputs*/
     3650        Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
     3651        Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
     3652        Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
     3653        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
     3654        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
     3655        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3656
     3657        /* Start  looping on the number of gaussian points: */
     3658        gauss=new GaussTria(4);
     3659        for (ig=gauss->begin();ig<gauss->end();ig++){
     3660
     3661                gauss->GaussPoint(ig);
     3662
     3663                thickness_input->GetInputValue(&thickness,gauss);
     3664                rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
     3665                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     3666                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     3667                adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
     3668                adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
     3669
     3670                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3671                matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
    35683672
    35693673                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

  • issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.h

    r11332 r11462  
    144144                void   Gradj(Vec gradient,int control_type,int control_index);
    145145                void   GradjBGradient(Vec gradient,int weight_index,int control_index);
     146                void   GradjZGradient(Vec gradient,int weight_index,int control_index);
    146147                void   GradjBMacAyeal(Vec gradient,int control_index);
     148                void   GradjZMacAyeal(Vec gradient,int control_index);
    147149                void   GradjDragMacAyeal(Vec gradient,int control_index);
    148150                void   GradjDragStokes(Vec gradient,int control_index);

  • issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp

    r11332 r11462  
    228228}
    229229/*}}}*/
     230/*FUNCTION Matice::GetZ {{{1*/
     231double Matice::GetZ(){
     232
     233        /*Output*/
     234        double Z;
     235
     236        inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
     237        return Z;
     238}
     239/*}}}*/
     240/*FUNCTION Matice::GetZbar {{{1*/
     241double Matice::GetZbar(){
     242
     243        /*Output*/
     244        double Zbar;
     245
     246        inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
     247        return Zbar;
     248}
     249/*}}}*/
    230250/*FUNCTION Matice::GetVectorFromInputs{{{1*/
    231251void  Matice::GetVectorFromInputs(Vec vector,int input_enum){
    … …  
    282302        /*Intermediary: */
    283303        double A,e;
    284         double B,n;
     304        double Btmp,B,n,Z;
    285305
    286306        /*Get B and n*/
    287         B=GetBbar();
     307        Btmp=GetBbar();
     308        Z=GetZbar();
    288309        n=GetN();
     310        B=Z*Btmp;
    289311
    290312        if (n==1){
    … …  
    318340        if(viscosity<=0) _error_("Negative viscosity");
    319341        _assert_(B>0);
     342        _assert_(Z>0);
    320343        _assert_(n>0);
    321344
    … …  
    348371        /*Intermediaries: */
    349372        double A,e;
    350         double B,n;
     373        double Btmp,B,n,Z;
    351374
    352375        /*Get B and n*/
    353         B=GetB();
     376        Btmp=GetB();
     377        Z=GetZ();
    354378        n=GetN();
     379        B=Z*Btmp;
    355380
    356381        if (n==1){
    … …  
    389414        if(viscosity3d<=0) _error_("Negative viscosity");
    390415        _assert_(B>0);
     416        _assert_(Z>0);
    391417        _assert_(n>0);
    392418
    … …  
    418444        /*Intermediaries: */
    419445        double A,e;
    420         double B,n;
     446        double Btmp,B,n,Z;
    421447        double eps0;
    422448
    423449        /*Get B and n*/
    424450        eps0=pow((double)10,(double)-27);
    425         B=GetB();
     451        Btmp=GetB();
     452        Z=GetZ();
    426453        n=GetN();
     454        B=Z*Btmp;
    427455       
    428456        if (n==1){
    … …  
    490518
    491519        /*Get B and n*/
    492         B=GetBbar();
     520        B=GetBbar(); /* why is this needed in this function? */
    493521        n=GetN();
    494522
    … …  
    508536               
    509537                        viscosity_complement=1/(2*pow(A,e));
     538                }
     539        }
     540        else{
     541                viscosity_complement=4.5*pow((double)10,(double)17);
     542        }
     543
     544        /*Checks in debugging mode*/
     545        _assert_(B>0);
     546        _assert_(n>0);
     547        _assert_(viscosity_complement>0);
     548               
     549        /*Return: */
     550        *pviscosity_complement=viscosity_complement;
     551}
     552/*}}}*/
     553/*FUNCTION Matice::GetViscosityZComplement {{{1*/
     554void  Matice::GetViscosityZComplement(double* pviscosity_complement, double* epsilon){
     555        /*Return viscosity derivative for control method d(mu)/dZ:
     556         *
     557         *                                                                                             B
     558         * dviscosity= -------------------------------------------------------------------
     559         *                                2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
     560         *
     561         * If epsilon is NULL, it means this is the first time Gradjb is being run, and we
     562         * return mu20, initial viscosity.
     563         */
     564       
     565        /*output: */
     566        double viscosity_complement;
     567
     568        /*input strain rate: */
     569        double exx,eyy,exy;
     570
     571        /*Intermediary value A and exponent e: */
     572        double A,e;
     573        double B,n;
     574
     575        /*Get B and n*/
     576        B=GetBbar();
     577        n=GetN();
     578
     579        if(epsilon){
     580                exx=*(epsilon+0);
     581                eyy=*(epsilon+1);
     582                exy=*(epsilon+2);
     583
     584                /*Build viscosity: mu2=B/(2*A^e) */
     585                A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
     586                if(A==0){
     587                        /*Maximum viscosity_complement for 0 shear areas: */
     588                        viscosity_complement=2.25*pow((double)10,(double)17);
     589                }
     590                else{
     591                        e=(n-1)/(2*n);
     592               
     593                        viscosity_complement=B/(2*pow(A,e));
    510594                }
    511595        }
    … …  
    756840                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyNEnum,nodeinputs));
    757841                }
     842               
     843                /*Get Z*/
     844                if (iomodel->Data(MaterialsRheologyZEnum)) {
     845                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     846                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
     847                }
    758848
    759849                /*Control Inputs*/
    … …  
    771861                                                }
    772862                                                break;
     863                                        case MaterialsRheologyZbarEnum:
     864                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     865                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     866                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     867                                                        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];
     868                                                        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];
     869                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     870                                                }
     871                                                break;
    773872                                }
    774873                        }
    … …  
    797896                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    798897                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
     898                }
     899               
     900                /*Get Z*/
     901                if (iomodel->Data(MaterialsRheologyZEnum)) {
     902                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     903                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
    799904                }
    800905
    … …  
    813918                                                }
    814919                                                break;
     920                                        case MaterialsRheologyZbarEnum:
     921                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     922                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     923                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     924                                                        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];
     925                                                        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];
     926                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     927                                                }
     928                                                break;
    815929                                }
    816930                        }
    … …  
    831945                                name==MaterialsRheologyBEnum ||
    832946                                name==MaterialsRheologyBbarEnum ||
    833                                 name==MaterialsRheologyNEnum
     947                                name==MaterialsRheologyNEnum ||
     948                                name==MaterialsRheologyZEnum ||
     949                                name==MaterialsRheologyZbarEnum
    834950                ){
    835951                return true;

  • issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.h

    r11332 r11462  
    6868                void   GetViscosity3dStokes(double* pviscosity3d, double* epsilon);
    6969                void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
     70                void   GetViscosityZComplement(double* pviscosity_complement, double* pepsilon);
    7071                void   GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
    7172                void   GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
    … …  
    7374                double GetBbar();
    7475                double GetN();
     76                double GetZ();
     77                double GetZbar();
    7578                bool   IsInput(int name);
    7679                /*}}}*/

  • issm/branches/trunk-jpl-damage/src/m/classes/inversion.m

    r11276 r11462  
    8787
    8888                        checkfield(md,'inversion.iscontrol','values',[0 1]);
    89                         checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'});
     89                        checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy'});
    9090                        checkfield(md,'inversion.nsteps','numel',1,'>=',1);
    9191                        checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
    … …  
    109109                        fielddisplay(obj,'control_parameters','parameter where inverse control is carried out; ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
    110110                        fielddisplay(obj,'nsteps','number of optimization searches');
    111                         fielddisplay(obj,'cost_functions','indicate the type of response for each optimization steps');
     111                        fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step');
    112112                        fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
    113113                        fielddisplay(obj,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied');

  • issm/branches/trunk-jpl-damage/src/m/classes/materials.m

    r10969 r11462  
    1818                rheology_B   = NaN;
    1919                rheology_n   = NaN;
     20                rheology_Z   = NaN;
    2021                rheology_law = '';
    2122        end
    … …  
    7879                        checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
    7980                        checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
     81                        checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]);
    8082                        checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'});
    8183                end % }}}
    … …  
    9597                        fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]');
    9698                        fielddisplay(obj,'rheology_n','Glen''s flow law exponent');
     99                        fielddisplay(obj,'rheology_Z','rheology multiplier');
    97100                        fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius''');
    98101                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/branches/trunk-jpl-damage/src/m/classes/model/model.m

    r11139 r11462  
    192192                         if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
    193193                         if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
     194                         if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end
    194195                         if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
    195196                         if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
     197                         if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end
    196198                         if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end
    197199                         if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end

  • issm/branches/trunk-jpl-damage/src/m/enum/EnumToModelField.m

    r9681 r11462  
    1515                case MaterialsRheologyBEnum(), string='rheology_B'; return
    1616                case MaterialsRheologyBbarEnum(), string='rheology_B'; return
     17                case MaterialsRheologyZEnum(), string='rheology_Z'; return
     18                case MaterialsRheologyZbarEnum(), string='rheology_Z'; return
    1719                case BalancethicknessThickeningRateEnum: string='dhdt'; return
    1820                case VxEnum(), string='vx'; return

  • issm/branches/trunk-jpl-damage/src/m/model/collapse.m

    r11408 r11462  
    7575md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
    7676md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
     77md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z);
    7778
    7879%special for thermal modeling:

  • issm/branches/trunk-jpl-damage/src/m/model/extrude.m

    r11315 r11462  
    201201md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
    202202md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
     203md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node');
    203204
    204205%parameters

  • issm/branches/trunk-jpl-damage/src/m/model/mechanicalproperties.m

    r9734 r11462  
    11function md=mechanicalproperties(md,vx,vy)
    2 %MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity
     2%MECHANICALPROPERTIES - compute stress and strain rate for a given velocity
    33%
    44%   this routine computes the components of the stress tensor
    … …  
    4646clear vxlist vylist
    4747
    48 %compute viscosity
     48%compute viscosity (it's not clear if this section is even used, since the variables are
     49%cleared before anything is saved to the model...
    4950nu=zeros(numberofelements,1);
     51Z_bar=md.materials.rheology_Z(index)*summation/3;
    5052B_bar=md.materials.rheology_B(index)*summation/3;
    5153power=(md.materials.rheology_n-1)./(2*md.materials.rheology_n);
    … …  
    5860location=find(second_inv==0 & power==0);
    5961nu(location)=B_bar(location);
    60 clear B_bar location second_inv power
     62clear B_bar Z_bar location second_inv power
    6163
    6264%compute stress
Note: See TracChangeset for help on using the changeset viewer.