Changeset 5682
- Timestamp:
- 09/07/10 08:04:36 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Inputs.cpp
r5660 r5682 99 99 } 100 100 /*}}}*/ 101 /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/ 102 void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){ 103 104 vector<Object*>::iterator object; 105 Input* input=NULL; 106 bool found=false; 107 108 /*Go through inputs and check whether any input with the same name is already in: */ 109 for ( object=objects.begin() ; object < objects.end(); object++ ){ 110 111 input=(Input*)(*object); 112 if (input->EnumType()==enum_type){ 113 found=true; 114 break; 115 } 116 } 117 118 if (!found){ 119 /*we could not find an input with the correct enum type. No defaults values were provided, 120 * error out: */ 121 ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type)); 122 } 123 124 /*Ok, we have an input if we made it here, request the input to return the values: */ 125 input->GetParameterValue(pvalue,gauss); 126 127 } 128 /*}}}*/ 101 129 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/ 102 130 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){ -
issm/trunk/src/c/Container/Inputs.h
r5660 r5682 17 17 class Node; 18 18 class GaussTria; 19 class GaussPenta; 19 20 20 21 class Inputs: public DataSet{ … … 49 50 void GetParameterValue(double* pvalue,double* gauss,int enum_type); 50 51 void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type); 52 void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type); 51 53 void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue); 52 54 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); -
issm/trunk/src/c/objects/Elements/Penta.cpp
r5678 r5682 2879 2879 2880 2880 /*Calculate DL on gauss point */ 2881 friction->GetAlpha2(&alpha2_gauss, gauss _tria,VxEnum,VyEnum,VzEnum);2881 friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum); 2882 2882 2883 2883 DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d; -
issm/trunk/src/c/objects/Loads/Friction.cpp
r5640 r5682 177 177 } 178 178 /*}}}*/ 179 /*FUNCTION Friction::GetAlpha2{{{1*/ 180 void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){ 181 182 /*This routine calculates the basal friction coefficient 183 alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ 184 185 /*diverse: */ 186 double r,s; 187 double drag_p, drag_q; 188 double gravity,rho_ice,rho_water; 189 double Neff; 190 double thickness,bed; 191 double vx,vy,vz,vmag; 192 double drag_coefficient; 193 double alpha2; 194 195 196 /*Recover parameters: */ 197 inputs->GetParameterValue(&drag_p,DragPEnum); 198 inputs->GetParameterValue(&drag_q,DragQEnum); 199 inputs->GetParameterValue(&thickness, gauss,ThicknessEnum); 200 inputs->GetParameterValue(&bed, gauss,BedEnum); 201 inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum); 202 203 /*Get material parameters: */ 204 gravity=matpar->GetG(); 205 rho_ice=matpar->GetRhoIce(); 206 rho_water=matpar->GetRhoWater(); 207 208 //compute r and q coefficients: */ 209 r=drag_q/drag_p; 210 s=1./drag_p; 211 212 //From bed and thickness, compute effective pressure when drag is viscous: 213 Neff=gravity*(rho_ice*thickness+rho_water*bed); 214 215 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 216 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 217 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 218 replace it by Neff=0 (ie, equival it to an ice shelf)*/ 219 if (Neff<0)Neff=0; 220 221 if(strcmp(element_type,"2d")==0){ 222 inputs->GetParameterValue(&vx, gauss,vxenum); 223 inputs->GetParameterValue(&vy, gauss,vyenum); 224 vmag=sqrt(pow(vx,2)+pow(vy,2)); 225 } 226 else if (strcmp(element_type,"3d")==0){ 227 inputs->GetParameterValue(&vx, gauss,vxenum); 228 inputs->GetParameterValue(&vy, gauss,vyenum); 229 inputs->GetParameterValue(&vz, gauss,vzenum); 230 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); 231 } 232 else ISSMERROR("element_type %s not supported yet",element_type); 233 234 alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1)); 235 236 /*Assign output pointers:*/ 237 *palpha2=alpha2; 238 } 239 /*}}}*/ 179 240 /*FUNCTION Friction::GetAlphaComplement {{{1*/ 180 241 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){ -
issm/trunk/src/c/objects/Loads/Friction.h
r5640 r5682 29 29 void GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum); 30 30 void GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum); 31 void GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum); 31 32 void GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum); 32 33 void GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
Note:
See TracChangeset
for help on using the changeset viewer.