Index: /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11372) +++ /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11373) @@ -218,4 +218,24 @@ } /*}}}*/ +/*FUNCTION Matice::GetZ {{{1*/ +double Matice::GetZ(){ + + /*Output*/ + double Z; + + inputs->GetInputAverage(&Z,MaterialsRheologyZEnum); + return Z; +} +/*}}}*/ +/*FUNCTION Matice::GetZbar {{{1*/ +double Matice::GetZbar(){ + + /*Output*/ + double Zbar; + + inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum); + return Zbar; +} +/*}}}*/ /*FUNCTION Matice::GetN {{{1*/ double Matice::GetN(){ @@ -263,10 +283,10 @@ void Matice::GetViscosity2d(double* pviscosity, double* 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] - where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity - vector, and n the flow law exponent. + where viscosity is the viscosity, Z the damage effect variable (1-D), + B the flow law parameter, (u,v) the velocity vector, and n the flow law exponent. If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we @@ -282,9 +302,11 @@ /*Intermediary: */ double A,e; - double B,n; - - /*Get B and n*/ - B=GetBbar(); + double Btmp,B,n,Z; + + /*Get B, Z and n*/ + Btmp=GetBbar(); + Z=GetZbar(); n=GetN(); + B=Z*Btmp; if (n==1){ @@ -318,4 +340,5 @@ if(viscosity<=0) _error_("Negative viscosity"); _assert_(B>0); + _assert_(Z>0); _assert_(n>0); @@ -329,5 +352,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * - * B + * Z * B * viscosity3d= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -348,9 +371,11 @@ /*Intermediaries: */ double A,e; - double B,n; + double Btmp,B,n,Z; /*Get B and n*/ - B=GetB(); + Btmp=GetB(); + Z=GetZ(); n=GetN(); + B=Z*Btmp; if (n==1){ @@ -389,4 +414,5 @@ if(viscosity3d<=0) _error_("Negative viscosity"); _assert_(B>0); + _assert_(Z>0); _assert_(n>0); @@ -399,5 +425,5 @@ /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * - * B + * Z * B * viscosity3d= ------------------------------------------------------------------- * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] @@ -418,11 +444,13 @@ /*Intermediaries: */ double A,e; - double B,n; + double Btmp,B,Z,n; double eps0; /*Get B and n*/ eps0=pow((double)10,(double)-27); - B=GetB(); - n=GetN(); + Btmp=GetB(); + Z=GetZ(); + n=GetB(); + B=Z*Btmp; if (n==1){ @@ -461,4 +489,5 @@ if(viscosity3d<=0) _error_("Negative viscosity"); _assert_(B>0); + _assert_(Z>0); _assert_(n>0); @@ -490,5 +519,5 @@ /*Get B and n*/ - B=GetBbar(); + B=GetBbar(); /* why is this needed in this function? */ n=GetN(); @@ -508,4 +537,60 @@ viscosity_complement=1/(2*pow(A,e)); + } + } + else{ + viscosity_complement=4.5*pow((double)10,(double)17); + } + + /*Checks in debugging mode*/ + _assert_(B>0); /* shouldn't be necessary as B is never used! */ + _assert_(n>0); + _assert_(viscosity_complement>0); + + /*Return: */ + *pviscosity_complement=viscosity_complement; +} +/*}}}*/ +/*FUNCTION Matice::GetViscosityZComplement {{{1*/ +void Matice::GetViscosityZComplement(double* pviscosity_complement, double* 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: */ + double viscosity_complement; + + /*input strain rate: */ + double exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + double A,e; + double 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((double)10,(double)17); + } + else{ + e=(n-1)/(2*n); + + viscosity_complement=B/(2*pow(A,e)); } } @@ -685,4 +770,10 @@ this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,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 TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + } /*Get n*/ @@ -706,4 +797,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(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; } } @@ -726,4 +826,10 @@ for(i=0;iData(MaterialsRheologyBEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)]; this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,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)); } @@ -748,4 +854,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,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1)); + } + break; + } } @@ -766,4 +882,6 @@ name==MaterialsRheologyBEnum || name==MaterialsRheologyBbarEnum || + name==MaterialsRheologyZEnum || + name==MaterialsRheologyZbarEnum || name==MaterialsRheologyNEnum ){ Index: /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11372) +++ /issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11373) @@ -68,6 +68,9 @@ void GetViscosity3dStokes(double* pviscosity3d, double* epsilon); void GetViscosityComplement(double* pviscosity_complement, double* pepsilon); + void GetViscosityZComplement(double* pviscosity_complement, double* pepsilon); double GetB(); double GetBbar(); + double GetZ(); + double GetZbar(); double GetN(); bool IsInput(int name);