Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26462) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 26463) @@ -409,10 +409,7 @@ 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 ************************************//*{{{*/ @@ -470,5 +467,5 @@ /*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); @@ -947,7 +944,4 @@ 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); @@ -988,7 +982,5 @@ 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); Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26462) +++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (revision 26463) @@ -360,6 +360,4 @@ 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); @@ -383,6 +381,4 @@ /*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); @@ -395,4 +391,5 @@ /*Get basal velocity*/ + friction->GetBasalSlidingSpeeds(&vx, &vy ,gauss); ub = sqrt(vx*vx + vy*vy); Index: /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26462) +++ /issm/trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp (revision 26463) @@ -119,4 +119,8 @@ 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"); @@ -254,6 +258,4 @@ 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); @@ -287,6 +289,4 @@ 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); @@ -304,6 +304,5 @@ /*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); @@ -494,6 +493,4 @@ 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); @@ -523,6 +520,4 @@ lr_input->GetInputValue(&lr,gauss); br_input->GetInputValue(&br,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); /*Get ice A parameter*/ @@ -539,6 +534,5 @@ /*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); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26462) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26463) @@ -3489,5 +3489,5 @@ /*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*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26462) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 26463) @@ -18,13 +18,73 @@ 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; - this->dim=dim_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->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"); } /*}}}*/ @@ -855,32 +915,44 @@ /*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: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 26462) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 26463) @@ -15,10 +15,15 @@ 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(); @@ -47,4 +52,7 @@ 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); };