Index: ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26462) +++ ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26463) @@ -408,12 +408,9 @@ Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input); Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input); Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); /*Build friction element, needed later: */ - Friction* friction=new Friction(element,dim); + Friction* friction=new Friction(element,3); /******** MELTING RATES ************************************//*{{{*/ element->NormalBase(&normal_base[0],xyz_list_base); @@ -469,7 +466,7 @@ /*basal friction*/ friction->GetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); vz_input->GetInputValue(&vz,gauss); + friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss); basalfriction=alpha2*(vx*vx + vy*vy + vz*vz); geothermalflux_input->GetInputValue(&geothermalflux,gauss); /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/ @@ -946,9 +943,6 @@ element->GetInputValue(&converged,ConvergedEnum); if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum; // use enthalpy from last iteration else enthalpy_enum=EnthalpyEnum; // use enthalpy from last time step - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); Input* enthalpy_input = element->GetInput(enthalpy_enum); _assert_(enthalpy_input); Input* pressure_input = element->GetInput(PressureEnum); _assert_(pressure_input); Input* watercolumn_input = element->GetInput(WatercolumnEnum); _assert_(watercolumn_input); @@ -987,9 +981,7 @@ // although one node might have become temperate. So keep heat flux switched on. geothermalflux_input->GetInputValue(&geothermalflux,gauss); friction->GetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); + friction->GetBasalSlidingSpeeds(&vx, &vy, &vz, gauss); basalfriction=alpha2*(vx*vx+vy*vy+vz*vz); heatflux=(basalfriction+geothermalflux)/(rho_ice); scalar=gauss->weight*Jdet*heatflux; Index: ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26462) +++ ../trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26463) @@ -359,8 +359,6 @@ IssmDouble g = element->FindParam(ConstantsGEnum); IssmDouble e_v = element->FindParam(HydrologyEnglacialVoidRatioEnum); Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input); - Input* vx_input = element->GetInput(VxBaseEnum);_assert_(vx_input); - Input* vy_input = element->GetInput(VyBaseEnum);_assert_(vy_input); Input* h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input); Input* H_input = element->GetInput(ThicknessEnum); _assert_(H_input); Input* b_input = element->GetInput(BedEnum); _assert_(b_input); @@ -382,8 +380,6 @@ element->NodalFunctions(basis,gauss); /*Get input values at gauss points*/ - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); h_input->GetInputValue(&h,gauss); G_input->GetInputValue(&G,gauss); B_input->GetInputValue(&B,gauss); @@ -394,6 +390,7 @@ H_input->GetInputValue(&H,gauss); /*Get basal velocity*/ + friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss); ub = sqrt(vx*vx + vy*vy); /*Compute cavity opening w*/ Index: ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26462) +++ ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26463) @@ -118,6 +118,10 @@ iomodel->FetchDataToInput(inputs,elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum); iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum); + if(iomodel->domaintype==Domain2DhorizontalEnum){ + iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum); + } iomodel->FindConstant(&frictionlaw,"md.friction.law"); /*Friction law variables*/ @@ -253,8 +257,6 @@ Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); Input* englacial_input = element->GetInput(HydrologyEnglacialInputEnum); _assert_(englacial_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); Input* headold_input = element->GetInput(HydrologyHeadOldEnum); _assert_(headold_input); @@ -286,8 +288,6 @@ englacial_input->GetInputValue(&ieb,gauss); lr_input->GetInputValue(&lr,gauss); br_input->GetInputValue(&br,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); headold_input->GetInputValue(&head_old,gauss); /*Get ice A parameter*/ @@ -303,8 +303,7 @@ /*Compute frictional heat flux*/ friction->GetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + friction->GetBasalSlidingSpeeds(&vx, &vy, gauss); frictionheat=alpha2*(vx*vx+vy*vy); /*Get water and ice pressures*/ @@ -493,8 +492,6 @@ Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); Input* B_input = element->GetInput(MaterialsRheologyBEnum); _assert_(B_input); Input* n_input = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); Input* lr_input = element->GetInput(HydrologyBumpSpacingEnum); _assert_(lr_input); Input* br_input = element->GetInput(HydrologyBumpHeightEnum); _assert_(br_input); @@ -522,8 +519,6 @@ head_input->GetInputDerivativeValue(&dh[0],xyz_list,gauss); lr_input->GetInputValue(&lr,gauss); br_input->GetInputValue(&br,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); /*Get ice A parameter*/ B_input->GetInputValue(&B,gauss); @@ -538,8 +533,7 @@ /*Compute frictional heat flux*/ friction->GetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); + friction->GetBasalSlidingSpeeds(&vx, &vy, gauss); frictionheat=alpha2*(vx*vx+vy*vy); /*Get water and ice pressures*/ Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26462) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26463) @@ -3488,7 +3488,7 @@ /*build friction object, used later on: */ /*dim=4 is special for HO, which is actually 2.5D*/ - Friction* friction=new Friction(element,dim==3?4:1); + Friction* friction=new Friction(element,dim==3?2.5:1); /*Recover portion of element that is grounded*/ if(!(friction_style==SubelementFriction2Enum)) phi=element->GetGroundedPortion(xyz_list_base); Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26462) +++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26463) @@ -17,15 +17,75 @@ /*Constructors/destructors*/ Friction::Friction(){/*{{{*/ this->element=NULL; - this->dim=0; this->law=0; + this->apply_dim = 1; + this->domaintype=-1; + this->vx_input=NULL; + this->vy_input=NULL; + this->vz_input=NULL; } /*}}}*/ -Friction::Friction(Element* element_in,int dim_in){/*{{{*/ +Friction::Friction(Element* element_in){/*{{{*/ + /* Determine the dimension according to the domain type automatically. + * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/ + this->element=element_in; + /* Load necessary parameters */ + element_in->FindParam(&this->law,FrictionLawEnum); + element_in->FindParam(&this->domaintype,DomainTypeEnum); + + /* Load VxBase and VyBase for this special case */ + switch(this->domaintype){ + case Domain2DhorizontalEnum: + this->apply_dim = 2; + this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input); + this->vz_input = NULL; + break; + case Domain2DverticalEnum: + this->apply_dim = 2; + this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); + this->vz_input = NULL; + break; + case Domain3DEnum: + this->apply_dim = 3; + this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); + this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input); + break; + default: _error_("mesh "<element=element_in; - this->dim=dim_in; + this->apply_dim = dim; + + /* Load necessary parameters */ element_in->FindParam(&this->law,FrictionLawEnum); + element_in->FindParam(&this->domaintype,DomainTypeEnum); + + /* Load VxBase and VyBase for this special case */ + switch(this->domaintype){ + case Domain2DhorizontalEnum: + this->vx_input = element_in->GetInput(VxBaseEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyBaseEnum); _assert_(this->vy_input); + this->vz_input = NULL; + break; + case Domain2DverticalEnum: + this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); + this->vz_input = NULL; + break; + case Domain3DEnum: + this->vx_input = element_in->GetInput(VxEnum); _assert_(this->vx_input); + this->vy_input = element_in->GetInput(VyEnum); _assert_(this->vy_input); + this->vz_input = element_in->GetInput(VzEnum); _assert_(this->vz_input); + break; + default: _error_("mesh "<dim<< "\n"); + _printf_(" Domain type: " << this->domaintype<< "\n"); } /*}}}*/ void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ @@ -854,33 +914,45 @@ /*diverse*/ IssmDouble vx,vy,vz,vmag; + + this->vx_input->GetInputValue(&vx, gauss); + this->vy_input->GetInputValue(&vy, gauss); + + if ((this->vz_input == NULL) || (this->apply_dim<3)) vz = 0.0; + else this->vz_input->GetInputValue(&vz, gauss); - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxBaseEnum); - element->GetInputValue(&vy,gauss,VyBaseEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - case 4: - /* This is for HO 3D case, since it requires the horizontal velocities */ - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - default: - _error_("not supported"); - } + if (this->apply_dim<2) vy = 0.0; + vmag = sqrt(vx*vx+vy*vy+vz*vz); return vmag; +}/*}}}*/ +void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/ + this->vx_input->GetInputValue(pvx, gauss); + /*Checks*/ + _assert_(!xIsNan(*pvx)); + _assert_(!xIsInf(*pvx)); }/*}}}*/ +void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/ + + this->vx_input->GetInputValue(pvx, gauss); + this->vy_input->GetInputValue(pvy, gauss); + /*Checks*/ + _assert_(!xIsNan(*pvx)); + _assert_(!xIsInf(*pvx)); + _assert_(!xIsNan(*pvy)); + _assert_(!xIsInf(*pvy)); +}/*}}}*/ +void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/ + + this->vx_input->GetInputValue(pvx, gauss); + this->vy_input->GetInputValue(pvy, gauss); + this->vz_input->GetInputValue(pvz, gauss); + /*Checks*/ + _assert_(!xIsNan(*pvx)); + _assert_(!xIsInf(*pvx)); + _assert_(!xIsNan(*pvy)); + _assert_(!xIsInf(*pvy)); + _assert_(!xIsNan(*pvz)); + _assert_(!xIsInf(*pvz)); +}/*}}}*/ Index: ../trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 26462) +++ ../trunk-jpl/src/c/classes/Loads/Friction.h (revision 26463) @@ -14,12 +14,17 @@ public: Element* element; - int dim; int law; + int domaintype; + double apply_dim; + Input* vx_input; + Input* vy_input; + Input* vz_input; /*methods: */ Friction(); - Friction(Element* element_in,int dim_in); + Friction(Element* element_in); + Friction(Element* element_in, double dim); ~Friction(); void Echo(void); @@ -46,6 +51,9 @@ IssmDouble EffectivePressure(Gauss* gauss); IssmDouble VelMag(Gauss* gauss); + void GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss); + void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss); + void GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss); }; #endif /* _FRICTION_H_ */