Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.cpp (revision 12897) @@ -1396,5 +1396,5 @@ /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); + if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type); else input=this->inputs->GetInput(enum_type); //if (!input) _error2_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs"); @@ -1502,4 +1502,5 @@ break; case MaterialsRheologyBbarEnum: + case MaterialsRheologyZbarEnum: /*Matice will take care of it*/ break; default: @@ -1690,5 +1691,5 @@ /*update input*/ - if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){ + if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ matice->inputs->AddInput(new TriaP1Input(name,values)); } @@ -1841,4 +1842,5 @@ name==FrictionCoefficientEnum || name==MaterialsRheologyBbarEnum || + name==MaterialsRheologyZbarEnum || name==GradientEnum || name==OldGradientEnum || @@ -2797,4 +2799,7 @@ *presponse=this->matice->GetBbar(); break; + case MaterialsRheologyZbarEnum: + *presponse=this->matice->GetZbar(); + break; case VelEnum: @@ -3407,5 +3412,5 @@ for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); } @@ -3434,5 +3439,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3452,5 +3457,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3471,5 +3476,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3501,4 +3506,7 @@ case MaterialsRheologyBbarEnum: GradjBMacAyeal(gradient,control_index); + break; + case MaterialsRheologyZbarEnum: + GradjZMacAyeal(gradient,control_index); break; case BalancethicknessThickeningRateEnum: @@ -3586,4 +3594,44 @@ } /*}}}*/ +/*FUNCTION Tria::GradjZGradient{{{*/ +void Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){ + + int i,ig; + int doflist1[NUMVERTICES]; + IssmDouble Jdet,weight; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble dbasis[NDOF2][NUMVERTICES]; + IssmDouble dk[NDOF2]; + IssmDouble grade_g[NUMVERTICES]={0.0}; + GaussTria *gauss=NULL; + + /*Retrieve all inputs we will be needing: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GradientIndexing(&doflist1[0],control_index); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + weights_input->GetInputValue(&weight,gauss,weight_index); + + /*Build alpha_complement_list: */ + rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]); + } + gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ /*FUNCTION Tria::GradjBMacAyeal{{{*/ void Tria::GradjBMacAyeal(Vector* gradient,int control_index){ @@ -3627,4 +3675,61 @@ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + + /*standard gradient dJ/dki*/ + for (i=0;iweight*basis[i]; + } + + gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjZMacAyeal{{{*/ +void Tria::GradjZMacAyeal(Vector* gradient,int control_index){ + + /*Intermediaries*/ + int i,ig; + int doflist[NUMVERTICES]; + IssmDouble vx,vy,lambda,mu,thickness,Jdet; + IssmDouble viscosity_complement; + IssmDouble dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2]; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble basis[3],epsilon[3]; + IssmDouble grad[NUMVERTICES]={0.0}; + GaussTria *gauss = NULL; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GradientIndexing(&doflist[0],control_index); + + /*Retrieve all inputs*/ + Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input); + Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetInputValue(&thickness,gauss); + rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss); + vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); + adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); + + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]); GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Elements/Tria.h (revision 12897) @@ -143,5 +143,7 @@ void Gradj(Vector* gradient,int control_type,int control_index); void GradjBGradient(Vector* gradient,int weight_index,int control_index); + void GradjZGradient(Vector* gradient,int weight_index,int control_index); void GradjBMacAyeal(Vector* gradient,int control_index); + void GradjZMacAyeal(Vector* gradient,int control_index); void GradjDragMacAyeal(Vector* gradient,int control_index); void GradjDragStokes(Vector* gradient,int control_index); Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.cpp (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.cpp (revision 12897) @@ -158,4 +158,24 @@ } /*}}}*/ +/*FUNCTION Matice::GetZ {{{*/ +IssmDouble Matice::GetZ(){ + + /*Output*/ + IssmDouble Z; + + inputs->GetInputAverage(&Z,MaterialsRheologyZEnum); + return Z; +} +/*}}}*/ +/*FUNCTION Matice::GetZbar {{{*/ +IssmDouble Matice::GetZbar(){ + + /*Output*/ + IssmDouble Zbar; + + inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum); + return Zbar; +} +/*}}}*/ /*FUNCTION Matice::GetVectorFromInputs{{{*/ void Matice::GetVectorFromInputs(Vector* vector,int input_enum){ @@ -193,5 +213,5 @@ void Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){ /*From a string tensor and a material object, return viscosity, using Glen's flow law. - B + Z * B viscosity= ------------------------------------------------------------------- 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -212,9 +232,11 @@ /*Intermediary: */ IssmDouble A,e; - IssmDouble B,n; + IssmDouble Btmp,B,n,Z; /*Get B and n*/ - B=GetBbar(); + Btmp=GetBbar(); + Z=GetZbar(); n=GetN(); + B=Z*Btmp; if (n==1){ @@ -438,4 +460,60 @@ viscosity_complement=1/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17); + } + + /*Checks in debugging mode*/ + _assert_(B>0); + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosityZComplement {{{*/ +void Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){ + /*Return viscosity derivative for control method d(mu)/dZ: + * + * B + * dviscosity= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * If epsilon is NULL, it means this is the first time Gradjb is being run, and we + * return mu20, initial viscosity. + */ + + /*output: */ + IssmDouble viscosity_complement; + + /*input strain rate: */ + IssmDouble exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + IssmDouble A,e; + IssmDouble B,n; + + /*Get B and n*/ + B=GetBbar(); + n=GetN(); + + if(epsilon){ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu2=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity_complement for 0 shear areas: */ + viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=B/(2*pow(A,e)); } } @@ -687,4 +765,10 @@ } + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[reCast(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + } + /*Control Inputs*/ #ifdef _HAVE_CONTROL_ @@ -701,4 +785,14 @@ } break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + } } @@ -727,4 +821,10 @@ for(i=0;iData(MaterialsRheologyNEnum)[index]; this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs)); + } + + /*Get Z*/ + if (iomodel->Data(MaterialsRheologyZEnum)) { + for(i=0;iData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs)); } @@ -743,4 +843,13 @@ } break; + case MaterialsRheologyZbarEnum: + if (iomodel->Data(MaterialsRheologyZEnum)){ + _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); + for(j=0;jData(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)]; + for(j=0;jData(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + for(j=0;jData(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i]; + this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; } } @@ -761,5 +870,7 @@ name==MaterialsRheologyBEnum || name==MaterialsRheologyBbarEnum || - name==MaterialsRheologyNEnum + name==MaterialsRheologyNEnum || + name==MaterialsRheologyZEnum || + name==MaterialsRheologyZbarEnum ){ return true; Index: /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.h (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/c/classes/objects/Materials/Matice.h (revision 12897) @@ -63,4 +63,5 @@ void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); @@ -68,4 +69,6 @@ IssmDouble GetBbar(); IssmDouble GetN(); + IssmDouble GetZ(); + IssmDouble GetZbar(); bool IsInput(int name); /*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/m/classes/materials.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/materials.m (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/m/classes/materials.m (revision 12897) @@ -77,4 +77,5 @@ md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); end % }}} Index: /issm/branches/trunk-jpl-damage/src/m/model/EnumToModelField.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/EnumToModelField.m (revision 12896) +++ /issm/branches/trunk-jpl-damage/src/m/model/EnumToModelField.m (revision 12897) @@ -13,4 +13,6 @@ case MaterialsRheologyBEnum(), string='rheology_B'; return case MaterialsRheologyBbarEnum(), string='rheology_B'; return + case MaterialsRheologyZEnum(), string='rheology_Z'; return + case MaterialsRheologyZbarEnum(), string='rheology_Z'; return case BalancethicknessThickeningRateEnum: string='dhdt'; return case VxEnum(), string='vx'; return