Index: /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 11462) @@ -100,4 +100,6 @@ MaterialsRheologyLawEnum, MaterialsRheologyNEnum, + MaterialsRheologyZEnum, + MaterialsRheologyZbarEnum, MaterialsRhoIceEnum, MaterialsRhoWaterEnum, Index: /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumToModelField.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumToModelField.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumToModelField.cpp (revision 11462) @@ -17,4 +17,6 @@ case MaterialsRheologyBEnum : return "rheology_B"; case MaterialsRheologyBbarEnum : return "rheology_B"; + case MaterialsRheologyZEnum : return "rheology_Z"; + case MaterialsRheologyZbarEnum : return "rheology_Z"; case BalancethicknessThickeningRateEnum: return "dhdt"; case VxEnum : return "vx"; Index: /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 11462) @@ -104,4 +104,6 @@ case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw"; case MaterialsRheologyNEnum : return "MaterialsRheologyN"; + case MaterialsRheologyZEnum : return "MaterialsRheologyZ"; + case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar"; case MaterialsRhoIceEnum : return "MaterialsRhoIce"; case MaterialsRhoWaterEnum : return "MaterialsRhoWater"; Index: /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 11462) @@ -48,4 +48,5 @@ case FrictionCoefficientEnum: iomodel->FetchData(1,FrictionCoefficientEnum); break; case MaterialsRheologyBbarEnum: iomodel->FetchData(1,MaterialsRheologyBEnum); break; + case MaterialsRheologyZbarEnum: iomodel->FetchData(1,MaterialsRheologyZEnum); break; default: _error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i])); } @@ -66,4 +67,4 @@ /*Free data: */ - iomodel->DeleteData(1+4+5,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum); + iomodel->DeleteData(1+4+6,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,MaterialsRheologyZEnum); } Index: /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 11462) @@ -45,5 +45,5 @@ /*Fetch data needed: */ - iomodel->FetchData(4,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum); + iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum); #ifdef _HAVE_THREED_ if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum); @@ -67,6 +67,7 @@ /*Free data: */ - iomodel->DeleteData(9,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, - MaterialsRheologyBEnum,MaterialsRheologyNEnum,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum); + iomodel->DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum, + MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum, + InversionMaxParametersEnum); /*Add new constrant material property tgo materials, at the end: */ Index: /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 11462) @@ -59,4 +59,5 @@ iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum); iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum); + iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum); iomodel->FetchDataToInput(elements,VxEnum); iomodel->FetchDataToInput(elements,VyEnum); Index: /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 11462) @@ -105,4 +105,6 @@ else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum; else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum; + else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum; + else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum; else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum; else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum; @@ -135,10 +137,10 @@ else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum; else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum; - else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum; - else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum; + if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum; + else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum; + else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum; else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum; else if (strcmp(name,"QmuMassFluxSegments")==0) return QmuMassFluxSegmentsEnum; @@ -258,10 +260,10 @@ else if (strcmp(name,"Node")==0) return NodeEnum; else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; - else if (strcmp(name,"Param")==0) return ParamEnum; - else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"Pengrid")==0) return PengridEnum; + if (strcmp(name,"Param")==0) return ParamEnum; + else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum; + else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; else if (strcmp(name,"Penta")==0) return PentaEnum; @@ -381,10 +383,10 @@ else if (strcmp(name,"DoubleElementResult")==0) return DoubleElementResultEnum; else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; - else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; - else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; + if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; + else if (strcmp(name,"DoubleVecExternalResult")==0) return DoubleVecExternalResultEnum; + else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum; else if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"Patch")==0) return PatchEnum; Index: /issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp (revision 11462) @@ -1473,5 +1473,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) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type)); @@ -1579,4 +1579,5 @@ break; case MaterialsRheologyBbarEnum: + case MaterialsRheologyZbarEnum: /*Matice will take care of it*/ break; default: @@ -1767,5 +1768,5 @@ /*update input*/ - if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){ + if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){ matice->inputs->AddInput(new TriaP1Input(name,values)); } @@ -2768,4 +2769,7 @@ *presponse=this->matice->GetBbar(); break; + case MaterialsRheologyZbarEnum: + *presponse=this->matice->GetZbar(); + break; case VelEnum: @@ -3346,5 +3350,5 @@ for(int i=0;iinputs->GetInput(control_type[i]); _assert_(input); } @@ -3373,5 +3377,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3391,5 +3395,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3410,5 +3414,5 @@ Input* input=NULL; - if(enum_type==MaterialsRheologyBbarEnum){ + if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){ input=(Input*)matice->inputs->GetInput(enum_type); } @@ -3431,5 +3435,5 @@ /*If on water, grad = 0: */ - if(IsOnWater()) return; + if(IsOnWater()) return; /*First deal with ∂/∂alpha(KU-F)*/ @@ -3440,4 +3444,7 @@ case MaterialsRheologyBbarEnum: GradjBMacAyeal(gradient,control_index); + break; + case MaterialsRheologyZbarEnum: + GradjZMacAyeal(gradient,control_index); break; case BalancethicknessThickeningRateEnum: @@ -3525,4 +3532,44 @@ } /*}}}*/ +/*FUNCTION Tria::GradjZGradient{{{1*/ +void Tria::GradjZGradient(Vec gradient,int weight_index,int control_index){ + + int i,ig; + int doflist1[NUMVERTICES]; + double Jdet,weight; + double xyz_list[NUMVERTICES][3]; + double dbasis[NDOF2][NUMVERTICES]; + double dk[NDOF2]; + double 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]); + } + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ /*FUNCTION Tria::GradjBMacAyeal{{{1*/ void Tria::GradjBMacAyeal(Vec gradient,int control_index){ @@ -3566,4 +3613,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]; + } + + VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES); + + /*clean-up*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjZMacAyeal{{{1*/ +void Tria::GradjZMacAyeal(Vec gradient,int control_index){ + + /*Intermediaries*/ + int i,ig; + int doflist[NUMVERTICES]; + double vx,vy,lambda,mu,thickness,Jdet; + double viscosity_complement; + double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2]; + double xyz_list[NUMVERTICES][3]; + double basis[3],epsilon[3]; + double 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/objects/Elements/Tria.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.h (revision 11462) @@ -144,5 +144,7 @@ void Gradj(Vec gradient,int control_type,int control_index); void GradjBGradient(Vec gradient,int weight_index,int control_index); + void GradjZGradient(Vec gradient,int weight_index,int control_index); void GradjBMacAyeal(Vec gradient,int control_index); + void GradjZMacAyeal(Vec gradient,int control_index); void GradjDragMacAyeal(Vec gradient,int control_index); void GradjDragStokes(Vec gradient,int control_index); Index: /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp (revision 11462) @@ -228,4 +228,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::GetVectorFromInputs{{{1*/ void Matice::GetVectorFromInputs(Vec vector,int input_enum){ @@ -282,9 +302,11 @@ /*Intermediary: */ double A,e; - double B,n; + double Btmp,B,n,Z; /*Get B and n*/ - B=GetBbar(); + 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); @@ -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); @@ -418,11 +444,13 @@ /*Intermediaries: */ double A,e; - double B,n; + double Btmp,B,n,Z; double eps0; /*Get B and n*/ eps0=pow((double)10,(double)-27); - B=GetB(); + Btmp=GetB(); + Z=GetZ(); n=GetN(); + B=Z*Btmp; if (n==1){ @@ -490,5 +518,5 @@ /*Get B and n*/ - B=GetBbar(); + B=GetBbar(); /* why is this needed in this function? */ n=GetN(); @@ -508,4 +536,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); + _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)); } } @@ -756,4 +840,10 @@ this->inputs->AddInput(new TriaP1Input(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 TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs)); + } /*Control Inputs*/ @@ -771,4 +861,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; } } @@ -797,4 +896,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)); } @@ -813,4 +918,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; } } @@ -831,5 +945,7 @@ name==MaterialsRheologyBEnum || name==MaterialsRheologyBbarEnum || - name==MaterialsRheologyNEnum + name==MaterialsRheologyNEnum || + name==MaterialsRheologyZEnum || + name==MaterialsRheologyZbarEnum ){ return true; Index: /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/c/objects/Materials/Matice.h (revision 11462) @@ -68,4 +68,5 @@ void GetViscosity3dStokes(double* pviscosity3d, double* epsilon); void GetViscosityComplement(double* pviscosity_complement, double* pepsilon); + void GetViscosityZComplement(double* pviscosity_complement, double* pepsilon); void GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon); void GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon); @@ -73,4 +74,6 @@ double GetBbar(); double GetN(); + double GetZ(); + double GetZbar(); bool IsInput(int name); /*}}}*/ Index: /issm/branches/trunk-jpl-damage/src/m/classes/inversion.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/inversion.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/classes/inversion.m (revision 11462) @@ -87,5 +87,5 @@ checkfield(md,'inversion.iscontrol','values',[0 1]); - checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'}); + checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy'}); checkfield(md,'inversion.nsteps','numel',1,'>=',1); checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0); @@ -109,5 +109,5 @@ fielddisplay(obj,'control_parameters','parameter where inverse control is carried out; ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'); fielddisplay(obj,'nsteps','number of optimization searches'); - fielddisplay(obj,'cost_functions','indicate the type of response for each optimization steps'); + fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step'); fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'); fielddisplay(obj,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied'); Index: /issm/branches/trunk-jpl-damage/src/m/classes/materials.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/materials.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/classes/materials.m (revision 11462) @@ -18,4 +18,5 @@ rheology_B = NaN; rheology_n = NaN; + rheology_Z = NaN; rheology_law = ''; end @@ -78,4 +79,5 @@ checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]); checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); + checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]); checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'}); end % }}} @@ -95,4 +97,5 @@ fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]'); fielddisplay(obj,'rheology_n','Glen''s flow law exponent'); + fielddisplay(obj,'rheology_Z','rheology multiplier'); fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius'''); end % }}} @@ -110,4 +113,5 @@ WriteData(fid,'object',obj,'fieldname','rheology_B','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','rheology_n','format','DoubleMat','mattype',2); + WriteData(fid,'object',obj,'fieldname','rheology_Z','format','DoubleMat','mattype',1); WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum,'format','Integer'); end % }}} Index: /issm/branches/trunk-jpl-damage/src/m/classes/model/model.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/model/model.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/classes/model/model.m (revision 11462) @@ -192,6 +192,8 @@ if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end + if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end + if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end Index: /issm/branches/trunk-jpl-damage/src/m/enum/EnumToModelField.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/EnumToModelField.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/enum/EnumToModelField.m (revision 11462) @@ -15,4 +15,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 Index: /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZEnum.m (revision 11462) +++ /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZEnum.m (revision 11462) @@ -0,0 +1,11 @@ +function macro=MaterialsRheologyZEnum() +%MATERIALSRHEOLOGYZENUM - Enum of MaterialsRheologyZ +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=MaterialsRheologyZEnum() + +macro=StringToEnum('MaterialsRheologyZ'); Index: /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZbarEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZbarEnum.m (revision 11462) +++ /issm/branches/trunk-jpl-damage/src/m/enum/MaterialsRheologyZbarEnum.m (revision 11462) @@ -0,0 +1,11 @@ +function macro=MaterialsRheologyZbarEnum() +%MATERIALSRHEOLOGYZBARENUM - Enum of MaterialsRheologyZbar +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=MaterialsRheologyZbarEnum() + +macro=StringToEnum('MaterialsRheologyZbar'); Index: /issm/branches/trunk-jpl-damage/src/m/model/collapse.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/collapse.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/model/collapse.m (revision 11462) @@ -75,4 +75,5 @@ md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B); md.materials.rheology_n=project2d(md,md.materials.rheology_n,1); +md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z); %special for thermal modeling: Index: /issm/branches/trunk-jpl-damage/src/m/model/extrude.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/extrude.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/model/extrude.m (revision 11462) @@ -201,4 +201,5 @@ md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node'); md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element'); +md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node'); %parameters Index: /issm/branches/trunk-jpl-damage/src/m/model/mechanicalproperties.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/mechanicalproperties.m (revision 11461) +++ /issm/branches/trunk-jpl-damage/src/m/model/mechanicalproperties.m (revision 11462) @@ -1,4 +1,4 @@ function md=mechanicalproperties(md,vx,vy) -%MECHANICALPROPERTIES - compute stress and strain rate for a goven velocity +%MECHANICALPROPERTIES - compute stress and strain rate for a given velocity % % this routine computes the components of the stress tensor @@ -46,6 +46,8 @@ clear vxlist vylist -%compute viscosity +%compute viscosity (it's not clear if this section is even used, since the variables are +%cleared before anything is saved to the model... nu=zeros(numberofelements,1); +Z_bar=md.materials.rheology_Z(index)*summation/3; B_bar=md.materials.rheology_B(index)*summation/3; power=(md.materials.rheology_n-1)./(2*md.materials.rheology_n); @@ -58,5 +60,5 @@ location=find(second_inv==0 & power==0); nu(location)=B_bar(location); -clear B_bar location second_inv power +clear B_bar Z_bar location second_inv power %compute stress Index: /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parametercontrolZ.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parametercontrolZ.m (revision 11462) +++ /issm/branches/trunk-jpl-damage/src/m/model/parameterization/parametercontrolZ.m (revision 11462) @@ -0,0 +1,121 @@ +function md=parametercontrolZ(md,varargin), +%PARAMETERCONTROLZ - parameterization for control method on Z +% +% It is possible to specify the number of steps, values for the +% minimum and maximum values of Z, the +% kind of cm_responses to use or the the optscal. +% +% Usage: +% md=parametercontrolZ(md,varargin) +% +% Example: +% md=parametercontrolZ(md) +% md=parametercontrolZ(md,'nsteps',20,'cm_responses',0) +% md=parametercontrolZ(md,'cm_min',0,'cm_max',1,'cm_jump',0.99,'maxiter',20) +% +% See also PARAMETERCONTROLB + +%process options +options=pairoptions(varargin{:}); + +%control type +md.inversion.control_parameters={'MaterialsRheologyZbar'}; + +%weights +weights=getfieldvalue(options,'weights',ones(md.mesh.numberofvertices,1)); +if (length(weights)~=md.mesh.numberofvertices) + md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,1); +else + md.inversion.cost_functions_coefficients=weights; +end + +%nsteps +nsteps=getfieldvalue(options,'nsteps',100); +if (length(nsteps)~=1 | nsteps<=0 | floor(nsteps)~=nsteps) + md.inversion.nsteps=100; +else + md.inversion.nsteps=nsteps; +end + + +%cm_min +cm_min=getfieldvalue(options,'cm_min',10^-13*ones(md.mesh.numberofvertices,1)); +if (length(cm_min)==1) + md.inversion.min_parameters=cm_min*ones(md.mesh.numberofvertices,1); +elseif (length(cm_min)==md.mesh.numberofvertices) + md.inversion.min_parameters=cm_min; +else + md.inversion.min_parameters=cm_min; +end + +%cm_max +cm_max=getfieldvalue(options,'cm_max',ones(md.mesh.numberofvertices,1)); +if (length(cm_max)==1) + md.inversion.max_parameters=cm_max*ones(md.mesh.numberofvertices,1); +elseif (length(cm_max)==md.mesh.numberofvertices) + md.inversion.max_parameters=cm_max; +else + md.inversion.max_parameters=cm_max; +end + +%eps_cm +eps_cm=getfieldvalue(options,'eps_cm',NaN); +if (length(eps_cm)~=1 | eps_cm<0 ) + md.inversion.cost_function_threshold=NaN; +else + md.inversion.cost_function_threshold=eps_cm; +end + +%maxiter +maxiter=getfieldvalue(options,'maxiter',10*ones(md.inversion.nsteps,1)); +if (any(maxiter<0) | any(floor(maxiter)~=maxiter)) + md.inversion.maxiter_per_step=10*ones(md.inversion.nsteps,1); +else + md.inversion.maxiter_per_step=repmat(maxiter(:),md.inversion.nsteps,1); + md.inversion.maxiter_per_step(md.inversion.nsteps+1:end)=[]; +end + +%cm_jump +cm_jump=getfieldvalue(options,'cm_jump',0.9*ones(md.inversion.nsteps,1)); +if ~isreal(cm_jump) + md.inversion.step_threshold=0.9*ones(md.inversion.nsteps,1); +else + md.inversion.step_threshold=repmat(cm_jump(:),md.inversion.nsteps,1); + md.inversion.step_threshold(md.inversion.nsteps+1:end)=[]; +end + +%cm_responses +found=0; +if exist(options,'cm_responses'), + cm_responses=getfieldvalue(options,'cm_responses'); + if ~any(~ismember(cm_responses,[ 101:105])), + md.inversion.cost_functions=repmat(cm_responses(:),md.inversion.nsteps,1); + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.cost_functions=[... + 103*ones(third,1);... + 101*ones(third,1);... + repmat([101;101;103;101],third,1)... + ]; + md.inversion.cost_functions(md.inversion.nsteps+1:end)=[]; +end + +%optscal +found=0; +if exist(options,'optscal'), + optscal=getfieldvalue(options,'optscal'); + if ~any(optscal<0), + md.inversion.gradient_scaling=repmat(optscal(:),md.inversion.nsteps,1); + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; + found=1; + end +end +if ~found + third=ceil(md.inversion.nsteps/3); + md.inversion.gradient_scaling=[ones(3,1);0.8*ones(third-3,1);0.5*ones(2*third,1);]; + md.inversion.gradient_scaling(md.inversion.nsteps+1:end)=[]; +end