Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16872) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16873) @@ -118,4 +118,5 @@ virtual void GetInputValue(int* pvalue,int enum_type)=0; virtual void GetInputValue(IssmDouble* pvalue,int enum_type)=0; + virtual void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type)=0; virtual void GetVerticesCoordinates(IssmDouble** xyz_list)=0; virtual void GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16872) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16873) @@ -193,5 +193,5 @@ /*Build friction element, needed later: */ - friction=new Friction("3d",inputs,matpar,StressbalanceAnalysisEnum); + friction=new Friction(this,3); /* Start looping on the number of gauss 2d (nodes on the bedrock) */ @@ -202,5 +202,5 @@ gauss->GaussPoint(ig); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); @@ -1445,4 +1445,12 @@ if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); input->GetInputValue(pvalue); + +}/*}}}*/ +/*FUNCTION Penta::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum) {{{*/ +void Penta::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){ + + Input* input=inputs->GetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue,gauss); }/*}}}*/ @@ -4580,5 +4588,5 @@ /*Build friction element, needed later: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gauss 2d (nodes on the bedrock) */ @@ -4603,5 +4611,5 @@ geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); @@ -4839,5 +4847,5 @@ /*Build frictoin element, needed later: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gauss 2d (nodes on the bedrock) */ @@ -4851,5 +4859,5 @@ geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); vx_input->GetInputValue(&vx,gauss); vy_input->GetInputValue(&vy,gauss); @@ -4985,4 +4993,7 @@ GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum); Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -4991,5 +5002,5 @@ /*Build friction element, needed later: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /******** MELTING RATES ************************************/ @@ -5032,5 +5043,5 @@ /*basal friction*/ - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0)); @@ -6028,5 +6039,5 @@ /*Build frictoin element, needed later: */ - friction=new Friction("2d",inputs,matpar,analysis_type); + friction=new Friction(this,2); /* Start looping on the number of gaussian points: */ @@ -6040,5 +6051,5 @@ /*Build alpha_complement_list: */ - friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,NULL); dragcoefficient_input->GetInputValue(&drag, gauss); @@ -6103,5 +6114,5 @@ /*Build frictoin element, needed later: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gaussian points: */ @@ -6112,5 +6123,5 @@ /*Recover alpha_complement and drag: */ - friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,vz_input); drag_input->GetInputValue(&drag,gauss); @@ -6905,5 +6916,5 @@ /*build friction object, used later on: */ - friction=new Friction("2d",inputs,matpar,analysis_type); + friction=new Friction(this,2); /* Start looping on the number of gaussian points: */ @@ -6914,5 +6925,5 @@ /*Friction: */ - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss); @@ -7139,5 +7150,5 @@ /*build friction object, used later on: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gaussian points: */ @@ -7157,5 +7168,5 @@ NormalBase(&bed_normal[0],&xyz_list_tria[0][0]); - friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d; @@ -7762,5 +7773,5 @@ /*build friction object, used later on: */ - friction=new Friction("2d",inputs,matpar,analysis_type); + friction=new Friction(this,2); /*Recover portion of element that is grounded*/ @@ -7784,5 +7795,5 @@ GetBHOFriction(&B[0],gauss); - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; if(migration_style==SubelementMigration2Enum){ @@ -8066,5 +8077,5 @@ /*build friction object, used later on: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gaussian points: */ @@ -8077,5 +8088,5 @@ GetLFS(BFriction,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input); D[0*2+0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx @@ -8324,5 +8335,5 @@ /*build friction object, used later on: */ - friction=new Friction("3d",inputs,matpar,analysis_type); + friction=new Friction(this,3); /* Start looping on the number of gauss 2d (nodes on the bedrock) */ @@ -8341,5 +8352,5 @@ this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); - friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); for(i=0;iGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); material->GetViscosity3dFS(&viscosity,&epsilon[0]); - friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input); for(i=0;iGetInput(inputenum); + if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element"); + input->GetInputValue(pvalue,gauss); + +}/*}}}*/ /*FUNCTION Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/ void Tria::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){ @@ -3805,5 +3813,5 @@ /*build friction object, used later on: */ - friction=new Friction("1d",inputs,matpar,analysis_type); + friction=new Friction(this,1); /* Start looping on the number of gauss 1d (nodes on the bedrock) */ @@ -3816,5 +3824,5 @@ GetLFS(BFriction,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,NULL); D = +alpha2*gauss->weight*Jdet; //taub_x = -alpha2 vx @@ -3954,5 +3962,5 @@ /*build friction object, used later on: */ - friction=new Friction("2d",inputs,matpar,analysis_type); + friction=new Friction(this,2); /*Recover portion of element that is grounded*/ @@ -3972,5 +3980,5 @@ gauss->GaussPoint(ig); - friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlpha2(&alpha2, gauss,vx_input,vy_input,NULL); if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2; if(migration_style==SubelementMigration2Enum){ @@ -4961,5 +4969,5 @@ /*Build frictoin element, needed later: */ - friction=new Friction("2d",inputs,matpar,analysis_type); + friction=new Friction(this,2); /*Retrieve all inputs we will be needing: */ @@ -4980,5 +4988,5 @@ /*Build alpha_complement_list: */ - friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum); + friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,NULL); dragcoefficient_input->GetInputValue(&drag, gauss); @@ -5005,5 +5013,5 @@ //for (int iv=0;ivGaussVertex(iv); - // friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum); + // friction->GetAlphaComplement(&alpha_complement,gauss,vx_input,vy_input,vz_input); // dragcoefficient_input->GetInputValue(&drag, gauss); // adjointx_input->GetInputValue(&lambda, gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16872) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16873) @@ -283,4 +283,5 @@ void GetInputValue(int* pvalue,int enum_type); void GetInputValue(IssmDouble* pvalue,int enum_type); + void GetInputValue(IssmDouble* pvalue,Gauss* gauss,int enum_type); void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype); void GetStrainRate2d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input); Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 16872) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 16873) @@ -18,23 +18,17 @@ /*FUNCTION Friction::Friction() {{{*/ Friction::Friction(){ - this->element_type=NULL; - this->inputs=NULL; - this->matpar=NULL; + this->element=NULL; + this->dim=0; } /*}}}*/ -/*FUNCTION Friction::Friction(const char* element_type, Inputs* inputs,Matpar* matpar,int analysis_type){{{*/ -Friction::Friction(const char* element_type_in,Inputs* inputs_in,Matpar* matpar_in, int in_analysis_type){ +/*FUNCTION Friction::Friction(Element* element,int dim){{{*/ +Friction::Friction(Element* element_in,int dim_in){ - this->analysis_type=in_analysis_type; - this->inputs=inputs_in; - this->element_type=xNew(strlen(element_type_in)+1); - xMemCpy(this->element_type,element_type_in,(strlen(element_type_in)+1)); - - this->matpar=matpar_in; + this->element=element_in; + this->dim=dim_in; } /*}}}*/ /*FUNCTION Friction::~Friction() {{{*/ Friction::~Friction(){ - xDelete(element_type); } /*}}}*/ @@ -44,12 +38,9 @@ void Friction::Echo(void){ _printf_("Friction:\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - _printf_(" element_type: " << this->element_type << "\n"); - inputs->Echo(); - matpar->Echo(); + _printf_(" dim: " << this->dim<< "\n"); } /*}}}*/ -/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{*/ -void Friction::GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ +/*FUNCTION Friction::GetAlpha2{{{*/ +void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ /*This routine calculates the basal friction coefficient @@ -59,5 +50,4 @@ IssmDouble r,s; IssmDouble drag_p, drag_q; - IssmDouble gravity,rho_ice,rho_water; IssmDouble Neff; IssmDouble thickness,bed; @@ -67,14 +57,12 @@ /*Recover parameters: */ - inputs->GetInputValue(&drag_p,FrictionPEnum); - inputs->GetInputValue(&drag_q,FrictionQEnum); - this->GetInputValue(&thickness, gauss,ThicknessEnum); - this->GetInputValue(&bed, gauss,BedEnum); - this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); + element->GetInputValue(&drag_p,FrictionPEnum); + element->GetInputValue(&drag_q,FrictionQEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&bed, gauss,BedEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); //compute r and q coefficients: */ @@ -84,32 +72,30 @@ //From bed and thickness, compute effective pressure when drag is viscous: Neff=gravity*(rho_ice*thickness+rho_water*bed); + if(Neff<0)Neff=0; - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; - - if(strcmp(element_type,"1d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - vmag=sqrt(vx*vx); + switch(dim){ + case 1: + vx_input->GetInputValue(&vx,gauss); + vmag=sqrt(vx*vx); + break; + case 2: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vmag=sqrt(vx*vx+vy*vy); + break; + case 3: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); + vmag=sqrt(vx*vx+vy*vy+vz*vz); + break; + default: + _error_("not supported"); } - else if(strcmp(element_type,"2d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); - } - else if (strcmp(element_type,"3d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - this->GetInputValue(&vz, gauss,vzenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - else _error_("element_type "<< element_type << " not supported yet"); /*Checks that s-1>0 if v=0*/ - if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha2 is Inf"); + if(vmag==0. && (s-1.)<0.) _error_("velocity is 0 and (s-1)=" << (s-1.) << "<0, alpha2 is Inf"); - alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); + alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.)); _assert_(!xIsNan(alpha2)); @@ -118,70 +104,6 @@ } /*}}}*/ -/*FUNCTION Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{*/ -void Friction::GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ - - /*This routine calculates the basal friction coefficient - alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ - - /*diverse: */ - IssmDouble r,s; - IssmDouble drag_p, drag_q; - IssmDouble gravity,rho_ice,rho_water; - IssmDouble Neff; - IssmDouble thickness,bed; - IssmDouble vx,vy,vz,vmag; - IssmDouble drag_coefficient; - IssmDouble alpha2; - - /*Recover parameters: */ - inputs->GetInputValue(&drag_p,FrictionPEnum); - inputs->GetInputValue(&drag_q,FrictionQEnum); - this->GetInputValue(&thickness, gauss,ThicknessEnum); - this->GetInputValue(&bed, gauss,BedEnum); - this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - //compute r and q coefficients: */ - r=drag_q/drag_p; - s=1./drag_p; - - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); - - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; - - if(strcmp(element_type,"2d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); - } - else if (strcmp(element_type,"3d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - this->GetInputValue(&vz, gauss,vzenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - else _error_("element_type "<< element_type << " not supported yet"); - - /*Checks that s-1>0 if v=0*/ - if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf"); - - alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); - _assert_(!xIsNan(alpha2)); - - /*Assign output pointers:*/ - *palpha2=alpha2; -} -/*}}}*/ -/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum) {{{*/ -void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum){ +/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss,int vxenum,int vyenum,int vzenum) {{{*/ +void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){ /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. @@ -196,18 +118,15 @@ IssmDouble drag_coefficient; IssmDouble bed,thickness; - IssmDouble gravity,rho_ice,rho_water; IssmDouble alpha_complement; /*Recover parameters: */ - inputs->GetInputValue(&drag_p,FrictionPEnum); - inputs->GetInputValue(&drag_q,FrictionQEnum); - this->GetInputValue(&thickness, gauss,ThicknessEnum); - this->GetInputValue(&bed, gauss,BedEnum); - this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); + element->GetInputValue(&drag_p,FrictionPEnum); + element->GetInputValue(&drag_q,FrictionQEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&bed, gauss,BedEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); //compute r and q coefficients: */ @@ -217,29 +136,31 @@ //From bed and thickness, compute effective pressure when drag is viscous: Neff=gravity*(rho_ice*thickness+rho_water*bed); - - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; + if(Neff<0)Neff=0; //We need the velocity magnitude to evaluate the basal stress: - if(strcmp(element_type,"2d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); + switch(dim){ + case 1: + vx_input->GetInputValue(&vx,gauss); + vmag=sqrt(vx*vx); + break; + case 2: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vmag=sqrt(vx*vx+vy*vy); + break; + case 3: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vz_input->GetInputValue(&vz,gauss); + vmag=sqrt(vx*vx+vy*vy+vz*vz); + break; + default: + _error_("not supported"); } - else if (strcmp(element_type,"3d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - this->GetInputValue(&vz, gauss,vzenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - else _error_("element_type "<< element_type << " not supported yet"); /*Checks that s-1>0 if v=0*/ - if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf"); + if(vmag==0. && (s-1.)<0.) _error_("velocity is 0 and (s-1)=" << (s-1.) << "<0, alpha2 is Inf"); - alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); _assert_(!xIsNan(alpha_complement)); + alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan(alpha_complement)); /*Assign output pointers:*/ @@ -247,85 +168,2 @@ } /*}}}*/ -/*FUNCTION Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum) {{{*/ -void Friction::GetAlphaComplement(IssmDouble* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ - - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. - * FrictionGetAlphaComplement is used in control methods on drag, and it computes: - * alpha_complement= Neff ^r * vel ^s*/ - - /*diverse: */ - IssmDouble r,s; - IssmDouble vx,vy,vz,vmag; - IssmDouble drag_p,drag_q; - IssmDouble Neff; - IssmDouble drag_coefficient; - IssmDouble bed,thickness; - IssmDouble gravity,rho_ice,rho_water; - IssmDouble alpha_complement; - - /*Recover parameters: */ - inputs->GetInputValue(&drag_p,FrictionPEnum); - inputs->GetInputValue(&drag_q,FrictionQEnum); - this->GetInputValue(&thickness, gauss,ThicknessEnum); - this->GetInputValue(&bed, gauss,BedEnum); - this->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); - - /*Get material parameters: */ - gravity=matpar->GetG(); - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - //compute r and q coefficients: */ - r=drag_q/drag_p; - s=1./drag_p; - - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); - - /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - replace it by Neff=0 (ie, equival it to an ice shelf)*/ - if (Neff<0)Neff=0; - - //We need the velocity magnitude to evaluate the basal stress: - if(strcmp(element_type,"2d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)); - } - else if (strcmp(element_type,"3d")==0){ - this->GetInputValue(&vx, gauss,vxenum); - this->GetInputValue(&vy, gauss,vyenum); - this->GetInputValue(&vz, gauss,vzenum); - vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - else _error_("element_type "<< element_type << " not supported yet"); - - /*Checks that s-1>0 if v=0*/ - if(vmag==0 && (s-1)<0) _error_("velocity is 0 and (s-1)=" << (s-1) << "<0, alpha_complement is Inf"); - - alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); _assert_(!xIsNan(alpha_complement)); - - /*Assign output pointers:*/ - *palpha_complement=alpha_complement; -} -/*}}}*/ -/*FUNCTION Friction::GetInputValue{{{*/ -void Friction::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type){ - - Input* input=inputs->GetInput(enum_type); - if(!input) _error_("input " << EnumToStringx(enum_type) << " not found"); - input->GetInputValue(pvalue,gauss); - -} -/*}}}*/ -/*FUNCTION Friction::GetInputValue{{{*/ -void Friction::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type){ - - Input* input=inputs->GetInput(enum_type); - if(!input) _error_("input " << EnumToStringx(enum_type) << " not found"); - input->GetInputValue(pvalue,gauss); - -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 16872) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 16873) @@ -19,21 +19,15 @@ int analysis_type; - char* element_type; - Inputs* inputs; - Matpar* matpar; + Element* element; + int dim; /*methods: */ Friction(); - Friction(const char* element_type, Inputs* inputs,Matpar* matpar, int analysis_type); + Friction(Element* element_in,int dim_in); ~Friction(); void Echo(void); - void GetAlpha2(IssmDouble* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); - void GetAlpha2(IssmDouble* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); - void GetAlphaComplement(IssmDouble* alpha_complement, GaussTria* gauss,int vxenum,int vyenum,int vzenum); - void GetAlphaComplement(IssmDouble* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); - void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,int enum_type); - void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int enum_type); - + void GetAlpha2(IssmDouble* palpha2, Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); + void GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input); };