Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3831) +++ /issm/trunk/src/c/Makefile.am (revision 3832) @@ -50,4 +50,6 @@ ./objects/Loads/Friction.h\ ./objects/Loads/Friction.cpp\ + ./objects/Loads/Friction2.h\ + ./objects/Loads/Friction2.cpp\ ./objects/DakotaPlugin.h\ ./objects/DakotaPlugin.cpp\ @@ -497,4 +499,6 @@ ./objects/Loads/Friction.h\ ./objects/Loads/Friction.cpp\ + ./objects/Loads/Friction2.h\ + ./objects/Loads/Friction2.cpp\ ./objects/DakotaPlugin.h\ ./objects/DakotaPlugin.cpp\ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3831) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3832) @@ -1282,5 +1282,4 @@ double gauss_weight; double gauss_l1l2l3[3]; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* matrices: */ @@ -1299,14 +1298,6 @@ double slope_magnitude; - /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]; - double vy_list[numgrids]; - double thickness_list[numgrids]; - double bed_list[numgrids]; - double dragcoefficient_list[numgrids]; - double drag_p,drag_q; - /*friction: */ - double alpha2_list[numgrids]={0.0,0.0,0.0}; + Friction2* friction2=NULL; double alpha2; @@ -1331,38 +1322,7 @@ } + /*build friction object, used later on: */ if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - - /*Recover inputs: */ - inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); - inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); - inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); - inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); - inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); - inputs->GetParameterValue(&drag_p,DragPEnum); - inputs->GetParameterValue(&drag_q,DragQEnum); - - /*Build alpha2_list used by drag stiffness matrix*/ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ - friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); - strcpy(friction->element_type,"2d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&dragcoefficient_list[0]; - friction->bed=&bed_list[0]; - friction->thickness=&thickness_list[0]; - friction->vx=&vx_list[0]; - friction->vy=&vy_list[0]; - friction->p=drag_p; - friction->q=drag_q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); + friction2=new Friction2("2d",inputs,matpar); /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ @@ -1378,4 +1338,7 @@ + /*Friction: */ + friction2->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum); + // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ @@ -1384,7 +1347,5 @@ if (slope_magnitude>MAXSLOPE){ - alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); + alpha2=pow((double)10,MOUNTAINKEXPONENT); } @@ -1395,7 +1356,5 @@ GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - /*Now, take care of the basal friction if there is any: */ - GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); - + DL_scalar=alpha2*gauss_weight*Jdet; for (i=0;i<2;i++){ @@ -1421,4 +1380,5 @@ xfree((void**)&third_gauss_area_coord); xfree((void**)&gauss_weights); + delete friction2; } @@ -2940,16 +2900,9 @@ double pressure_list[3]; double pressure; - double alpha2_list[3]; - double basalfriction_list[3]; + int drag_type; double basalfriction; + Friction2* friction2=NULL; + double alpha2,vx,vy; double geothermalflux_value; - - double vx_list[numgrids]; - double vy_list[numgrids]; - double thickness_list[numgrids]; - double bed_list[numgrids]; - double dragcoefficient_list[numgrids]; - double drag_p,drag_q; - int drag_type; /* gaussian points: */ @@ -2961,5 +2914,4 @@ double gauss_weight; double gauss_coord[3]; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /*matrices: */ @@ -2969,7 +2921,4 @@ double scalar; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&drag_type,DragTypeEnum); /* Get node coordinates and dof list: */ @@ -2984,44 +2933,8 @@ this->parameters->FindParam(&dt,DtEnum); - - /*Recover inputs: */ - inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); - inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); - inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); - inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); - inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); - inputs->GetParameterValue(&drag_p,DragPEnum); - inputs->GetParameterValue(&drag_q,DragQEnum); - - /*Build alpha2_list used by drag stiffness matrix*/ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - - friction->element_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char)); - strcpy(friction->element_type,"3d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&dragcoefficient_list[0]; - friction->bed=&bed_list[0]; - friction->thickness=&thickness_list[0]; - friction->vx=&vx_list[0]; - friction->vy=&vy_list[0]; - friction->p=drag_p; - friction->q=drag_q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - - /* Compute basal friction */ - for(i=0;iGetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); - GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord); - + + friction2->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum); + inputs->GetParameterValue(&vx, &gauss_coord[0],VxAverageEnum); + inputs->GetParameterValue(&vy, &gauss_coord[0],VyAverageEnum); + basalfriction= alpha2*(pow(vx,(double)2.0)+pow(vy,(double)2.0)); + /*Calculate scalar parameter*/ scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); @@ -3064,4 +2981,5 @@ xfree((void**)&third_gauss_area_coord); xfree((void**)&gauss_weights); + delete friction2; } Index: /issm/trunk/src/c/objects/Loads/Friction2.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction2.cpp (revision 3832) +++ /issm/trunk/src/c/objects/Loads/Friction2.cpp (revision 3832) @@ -0,0 +1,169 @@ +/*!\file Friction2.c + * \brief: implementation of the Friction2 object + */ + +/*Headers:*/ +/*{{{1*/ +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include +#include "../objects.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +/*}}}*/ + +/*methods: */ +/*FUNCTION Friction2::Friction2() {{{1*/ +Friction2::Friction2(){ + this->element_type=NULL; + this->inputs=NULL; + this->matpar=NULL; +} +/*}}}*/ +/*FUNCTION Friction2::Friction2(char* element_type, Inputs* inputs,Matpar* matpar){{{1*/ +Friction2::Friction2(char* element_type_in,Inputs* inputs_in,Matpar* matpar_in){ + + this->inputs=inputs_in; + this->element_type=(char*)xmalloc((strlen(element_type_in)+1)*sizeof(char)); + strcpy(this->element_type,element_type); + this->matpar=matpar_in; + +} +/*}}}*/ +/*FUNCTION Friction2::~Friction2() {{{1*/ +Friction2::~Friction2(){ + xfree((void**)&element_type); +} +/*}}}*/ +/*FUNCTION Friction2::Echo {{{1*/ +void Friction2::Echo(void){ + printf("Friction2:\n"); + printf(" element_type: %s\n",this->element_type); + inputs->Echo(); + matpar->Echo(); +} +/*}}}*/ +/*FUNCTION Friction2::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){{{1*/ +void Friction2::GetAlpha2(double* palpha2, double* 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: */ + double r,s; + double drag_p, drag_q; + double gravity,rho_ice,rho_water; + double Neff; + double thickness,bed; + double vx,vy,vz,vmag; + double drag_coefficient; + double alpha2; + + + /*Recover parameters: */ + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); + inputs->GetParameterValue(&bed, gauss,BedEnum); + inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); + + /*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")){ + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)); + } + else{ + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + inputs->GetParameterValue(&vz, gauss,vzenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); + } + + alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +} +/*}}}*/ +/*FUNCTION Friction2::GetAlphaComplement {{{1*/ +void Friction2::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){ + + /* 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: */ + int i; + double Neff; + double r,s; + double vx; + double vy; + double vmag; + double drag_p,drag_q; + double drag_coefficient; + double bed,thickness; + double gravity,rho_ice,rho_water; + double alpha_complement; + + /*Recover parameters: */ + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); + inputs->GetParameterValue(&bed, gauss,BedEnum); + inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); + + /*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: + inputs->GetParameterValue(&vx, gauss,vxenum); + inputs->GetParameterValue(&vy, gauss,vyenum); + vmag=sqrt(pow(vx,2)+pow(vy,2)); + + alpha_complement=pow(Neff,r)*pow(vmag,(s-1)); + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Friction2.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Friction2.h (revision 3832) +++ /issm/trunk/src/c/objects/Loads/Friction2.h (revision 3832) @@ -0,0 +1,33 @@ +/*!\file Friction2.h + * \brief: header file for friction2 object + */ + +#ifndef _FRICTION2_H_ +#define _FRICTION2_H_ + +/*Headers:*/ +/*{{{1*/ +class Inputs; +class Matpar; +/*}}}*/ + +class Friction2{ + + public: + + char* element_type; + Inputs* inputs; + Matpar* matpar; + + /*methods: */ + Friction2(); + Friction2(char* element_type, Inputs* inputs,Matpar* matpar); + ~Friction2(); + + void Echo(void); + void GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum); + void GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum); + +}; + +#endif /* _FRICTION2_H_ */ Index: /issm/trunk/src/c/objects/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 3831) +++ /issm/trunk/src/c/objects/objects.h (revision 3832) @@ -26,4 +26,5 @@ /*Loads: */ #include "./Loads/Friction.h" +#include "./Loads/Friction2.h" #include "./Loads/Icefront.h" #include "./Loads/Numericalflux.h"