Changeset 12897


Ignore:
Timestamp:
08/03/12 16:22:18 (13 years ago)
Author:
cborstad
Message:

BUG: patched some changes missed during merge, NR on branch now fixed

Location:
issm/branches/trunk-jpl-damage/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r12870 r12897  
    13961396
    13971397        /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    1398         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);
    13991399        else input=this->inputs->GetInput(enum_type);
    14001400        //if (!input) _error2_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs");
    … …  
    15021502                                        break;
    15031503                                case MaterialsRheologyBbarEnum:
     1504                                case MaterialsRheologyZbarEnum:
    15041505                                        /*Matice will take care of it*/ break;
    15051506                                default:
    … …  
    16901691
    16911692                        /*update input*/
    1692                         if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
     1693                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    16931694                                matice->inputs->AddInput(new TriaP1Input(name,values));
    16941695                        }
    … …  
    18411842                                name==FrictionCoefficientEnum ||
    18421843                                name==MaterialsRheologyBbarEnum ||
     1844                                name==MaterialsRheologyZbarEnum ||
    18431845                                name==GradientEnum ||
    18441846                                name==OldGradientEnum ||
    … …  
    27972799                        *presponse=this->matice->GetBbar();
    27982800                        break;
     2801                case MaterialsRheologyZbarEnum:
     2802                        *presponse=this->matice->GetZbar();
     2803                        break;
    27992804                case VelEnum:
    28002805
    … …  
    34073412        for(int i=0;i<num_controls;i++){
    34083413
    3409                 if(control_type[i]==MaterialsRheologyBbarEnum){
     3414                if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
    34103415                        input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
    34113416                }
    … …  
    34343439        Input* input=NULL;
    34353440
    3436         if(enum_type==MaterialsRheologyBbarEnum){
     3441        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34373442                input=(Input*)matice->inputs->GetInput(enum_type);
    34383443        }
    … …  
    34523457        Input* input=NULL;
    34533458
    3454         if(enum_type==MaterialsRheologyBbarEnum){
     3459        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34553460                input=(Input*)matice->inputs->GetInput(enum_type);
    34563461        }
    … …  
    34713476        Input* input=NULL;
    34723477
    3473         if(enum_type==MaterialsRheologyBbarEnum){
     3478        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
    34743479                input=(Input*)matice->inputs->GetInput(enum_type);
    34753480        }
    … …  
    35013506                case MaterialsRheologyBbarEnum:
    35023507                        GradjBMacAyeal(gradient,control_index);
     3508                        break;
     3509                case MaterialsRheologyZbarEnum:
     3510                        GradjZMacAyeal(gradient,control_index);
    35033511                        break;
    35043512                case BalancethicknessThickeningRateEnum:
    … …  
    35863594}
    35873595/*}}}*/
     3596/*FUNCTION Tria::GradjZGradient{{{*/
     3597void  Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){
     3598
     3599        int        i,ig;
     3600        int        doflist1[NUMVERTICES];
     3601        IssmDouble     Jdet,weight;
     3602        IssmDouble     xyz_list[NUMVERTICES][3];
     3603        IssmDouble     dbasis[NDOF2][NUMVERTICES];
     3604        IssmDouble     dk[NDOF2];
     3605        IssmDouble     grade_g[NUMVERTICES]={0.0};
     3606        GaussTria  *gauss=NULL;
     3607
     3608        /*Retrieve all inputs we will be needing: */
     3609        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3610        GradientIndexing(&doflist1[0],control_index);
     3611        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3612        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
     3613
     3614        /* Start  looping on the number of gaussian points: */
     3615        gauss=new GaussTria(2);
     3616        for (ig=gauss->begin();ig<gauss->end();ig++){
     3617
     3618                gauss->GaussPoint(ig);
     3619
     3620                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3621                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
     3622                weights_input->GetInputValue(&weight,gauss,weight_index);
     3623
     3624                /*Build alpha_complement_list: */
     3625                rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     3626
     3627                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     3628                for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
     3629        }
     3630        gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
     3631
     3632        /*Clean up and return*/
     3633        delete gauss;
     3634}
     3635/*}}}*/
    35883636/*FUNCTION Tria::GradjBMacAyeal{{{*/
    35893637void  Tria::GradjBMacAyeal(Vector* gradient,int control_index){
    … …  
    36273675                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    36283676                matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
     3677
     3678                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3679                GetNodalFunctions(basis,gauss);
     3680
     3681                /*standard gradient dJ/dki*/
     3682                for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
     3683                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
     3684                                        )*Jdet*gauss->weight*basis[i];
     3685        }
     3686
     3687        gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
     3688
     3689        /*clean-up*/
     3690        delete gauss;
     3691}
     3692/*}}}*/
     3693/*FUNCTION Tria::GradjZMacAyeal{{{*/
     3694void  Tria::GradjZMacAyeal(Vector* gradient,int control_index){
     3695
     3696        /*Intermediaries*/
     3697        int        i,ig;
     3698        int        doflist[NUMVERTICES];
     3699        IssmDouble     vx,vy,lambda,mu,thickness,Jdet;
     3700        IssmDouble     viscosity_complement;
     3701        IssmDouble     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2];
     3702        IssmDouble     xyz_list[NUMVERTICES][3];
     3703        IssmDouble     basis[3],epsilon[3];
     3704        IssmDouble     grad[NUMVERTICES]={0.0};
     3705        GaussTria *gauss = NULL;
     3706
     3707        /* Get node coordinates and dof list: */
     3708        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3709        GradientIndexing(&doflist[0],control_index);
     3710
     3711        /*Retrieve all inputs*/
     3712        Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
     3713        Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
     3714        Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
     3715        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
     3716        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
     3717        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
     3718
     3719        /* Start  looping on the number of gaussian points: */
     3720        gauss=new GaussTria(4);
     3721        for (ig=gauss->begin();ig<gauss->end();ig++){
     3722
     3723                gauss->GaussPoint(ig);
     3724
     3725                thickness_input->GetInputValue(&thickness,gauss);
     3726                rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
     3727                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     3728                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     3729                adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
     3730                adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
     3731
     3732                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3733                matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
    36293734
    36303735                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

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

    r12832 r12897  
    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/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.cpp

    r12832 r12897  
    158158}
    159159/*}}}*/
     160/*FUNCTION Matice::GetZ {{{*/
     161IssmDouble Matice::GetZ(){
     162
     163        /*Output*/
     164        IssmDouble Z;
     165
     166        inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
     167        return Z;
     168}
     169/*}}}*/
     170/*FUNCTION Matice::GetZbar {{{*/
     171IssmDouble Matice::GetZbar(){
     172
     173        /*Output*/
     174        IssmDouble Zbar;
     175
     176        inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
     177        return Zbar;
     178}
     179/*}}}*/
    160180/*FUNCTION Matice::GetVectorFromInputs{{{*/
    161181void  Matice::GetVectorFromInputs(Vector* vector,int input_enum){
    … …  
    193213void  Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
    194214        /*From a string tensor and a material object, return viscosity, using Glen's flow law.
    195                                                                                                     B
     215                                                                                                   Z * B
    196216          viscosity= -------------------------------------------------------------------
    197217                                                  2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
    … …  
    212232        /*Intermediary: */
    213233        IssmDouble A,e;
    214         IssmDouble B,n;
     234        IssmDouble Btmp,B,n,Z;
    215235
    216236        /*Get B and n*/
    217         B=GetBbar();
     237        Btmp=GetBbar();
     238        Z=GetZbar();
    218239        n=GetN();
     240        B=Z*Btmp;
    219241
    220242        if (n==1){
    … …  
    438460               
    439461                        viscosity_complement=1/(2*pow(A,e));
     462                }
     463        }
     464        else{
     465                viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
     466        }
     467
     468        /*Checks in debugging mode*/
     469        _assert_(B>0);
     470        _assert_(n>0);
     471        _assert_(viscosity_complement>0);
     472               
     473        /*Return: */
     474        *pviscosity_complement=viscosity_complement;
     475}
     476/*}}}*/
     477/*FUNCTION Matice::GetViscosityZComplement {{{*/
     478void  Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
     479        /*Return viscosity derivative for control method d(mu)/dZ:
     480         *
     481         *                                                                                             B
     482         * dviscosity= -------------------------------------------------------------------
     483         *                                2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
     484         *
     485         * If epsilon is NULL, it means this is the first time Gradjb is being run, and we
     486         * return mu20, initial viscosity.
     487         */
     488       
     489        /*output: */
     490        IssmDouble viscosity_complement;
     491
     492        /*input strain rate: */
     493        IssmDouble exx,eyy,exy;
     494
     495        /*Intermediary value A and exponent e: */
     496        IssmDouble A,e;
     497        IssmDouble B,n;
     498
     499        /*Get B and n*/
     500        B=GetBbar();
     501        n=GetN();
     502
     503        if(epsilon){
     504                exx=*(epsilon+0);
     505                eyy=*(epsilon+1);
     506                exy=*(epsilon+2);
     507
     508                /*Build viscosity: mu2=B/(2*A^e) */
     509                A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
     510                if(A==0){
     511                        /*Maximum viscosity_complement for 0 shear areas: */
     512                        viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
     513                }
     514                else{
     515                        e=(n-1)/(2*n);
     516               
     517                        viscosity_complement=B/(2*pow(A,e));
    440518                }
    441519        }
    … …  
    687765                }
    688766
     767                /*Get Z*/
     768                if (iomodel->Data(MaterialsRheologyZEnum)) {
     769                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     770                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
     771                }
     772
    689773                /*Control Inputs*/
    690774                #ifdef _HAVE_CONTROL_
    … …  
    701785                                                }
    702786                                                break;
     787                                        case MaterialsRheologyZbarEnum:
     788                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     789                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     790                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     791                                                        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];
     792                                                        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];
     793                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     794                                                }
     795                                                break;
     796
    703797                                }
    704798                        }
    … …  
    727821                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
    728822                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
     823                }
     824
     825                /*Get Z*/
     826                if (iomodel->Data(MaterialsRheologyZEnum)) {
     827                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
     828                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
    729829                }
    730830
    … …  
    743843                                                }
    744844                                                break;
     845                                        case MaterialsRheologyZbarEnum:
     846                                                if (iomodel->Data(MaterialsRheologyZEnum)){
     847                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum));
     848                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
     849                                                        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];
     850                                                        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];
     851                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     852                                                }
     853                                                break;
    745854                                }
    746855                        }
    … …  
    761870                                name==MaterialsRheologyBEnum ||
    762871                                name==MaterialsRheologyBbarEnum ||
    763                                 name==MaterialsRheologyNEnum
     872                                name==MaterialsRheologyNEnum ||
     873                                name==MaterialsRheologyZEnum ||
     874                                name==MaterialsRheologyZbarEnum
    764875                ){
    765876                return true;

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

    r12472 r12897  
    6363                void   GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
    6464                void   GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
     65                void   GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
    6566                void   GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6667                void   GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    … …  
    6869                IssmDouble GetBbar();
    6970                IssmDouble GetN();
     71                IssmDouble GetZ();
     72                IssmDouble GetZbar();
    7073                bool   IsInput(int name);
    7174                /*}}}*/

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

    r12878 r12897  
    7777                        md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
    7878                        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]);
    7980                        md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'});
    8081                end % }}}

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

    r11684 r12897  
    1313                case MaterialsRheologyBEnum(), string='rheology_B'; return
    1414                case MaterialsRheologyBbarEnum(), string='rheology_B'; return
     15                case MaterialsRheologyZEnum(), string='rheology_Z'; return
     16                case MaterialsRheologyZbarEnum(), string='rheology_Z'; return
    1517                case BalancethicknessThickeningRateEnum: string='dhdt'; return
    1618                case VxEnum(), string='vx'; return
Note: See TracChangeset for help on using the changeset viewer.