Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26291) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26292) @@ -21,4 +21,5 @@ //#define LATERALFRICTION 1 //#define DISCSLOPE 1 //testing for SSA +#define NOSPCSHEARVEL 1 //MLHO /*Model processing*/ @@ -194,9 +195,15 @@ } } - else{//MLHO itapopo FIXME base and shear velocities are constrained + else{//MLHO + #ifdef NOSPCSHEARVEL + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); + #else + /*Default: apply spcs to shear vx and shear vy*/ IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1); IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3); + #endif } } @@ -781,5 +788,6 @@ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.); - if(isMLHO){//itapopo FIXME shear and base velocities should be initialized correctly + if(isMLHO){ + /*itapopo FIXME applying the same initialization for shear vx and shear vy for now*/ iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.); @@ -2907,6 +2915,8 @@ Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy + Input* vxshear_input = element->GetInput(VxShearEnum); _assert_(vxshear_input); + Input* vyshear_input = element->GetInput(VyShearEnum); _assert_(vyshear_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); @@ -2923,6 +2933,6 @@ thickness_input->GetInputValue(&thickness, gauss); n_input->GetInputValue(&n,gauss); - //FIXME testing with L1L2-type viscosity - element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); + IssmDouble dim=2;//itapopo + element->material->ViscosityMLHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input); //viscosity=10e13;//itapopo Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26291) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26292) @@ -4212,4 +4212,39 @@ }/*}}}*/ +void Element::StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){/*{{{*/ + /*Compute the 2d Blatter/MLHO Strain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + /*Intermediaries*/ + IssmDouble dvx[2]; // [dvx/dx, dvx/dy] + IssmDouble dvy[2]; // [dvy/dx, dvy/dy] + IssmDouble vxshear,vyshear,thickness; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vxshear_input || !vyshear_input || !thickness_input){ + _error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vxshear: " << vxshear_input << ", vyshear: " << vyshear_input << ", thickness: " << thickness_input << "\n"); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss); + vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss); + thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0); + vxshear_input->GetInputValue(&vxshear,gauss); + vyshear_input->GetInputValue(&vyshear,gauss); + + epsilon[0] = dvx[0]; + epsilon[1] = dvy[1]; + epsilon[2] = 0.5*(dvx[1] + dvy[0]); + epsilon[3] = 0.5*vxshear/thickness; + epsilon[4] = 0.5*vyshear/thickness; + +}/*}}}*/ void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26291) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26292) @@ -181,4 +181,5 @@ void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void StrainRateMLHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input); void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input); Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 26291) +++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 26292) @@ -355,4 +355,14 @@ default: _error_("Approximation "<StrainRateMLHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vxshear_input,vyshear_input,thickness_input); + eps_eff = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[2]*epsilon[2] + epsilon[3]*epsilon[3] + epsilon[4]*epsilon[4] + epsilon[0]*epsilon[1]); + + /*Get viscosity*/ + this->GetViscosityBar(&viscosity,eps_eff,gauss); + + /*Assign output pointer*/ + *pviscosity=viscosity; +}/*}}}*/ void Matice::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/ /*Compute the L1L2 viscosity Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 26291) +++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h (revision 26292) @@ -78,4 +78,5 @@ void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); + void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input); void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf); void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h (revision 26291) +++ /issm/trunk-jpl/src/c/classes/Materials/Matlitho.h (revision 26292) @@ -72,4 +72,5 @@ void ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");}; void ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");}; + void ViscosityMLHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vxshear_input,Input* vyshear_input,Input* thickness_input){_error_("not supported");}; void ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");}; void ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};