Changeset 5682

Timestamp:

09/07/10 08:04:36 (15 years ago)

Author:

seroussi

Message:

fixed friction stokes

Location:

issm/trunk/src/c

Files:

• Container/Inputs.cpp (modified) (1 diff)
• Container/Inputs.h (modified) (2 diffs)
• objects/Elements/Penta.cpp (modified) (1 diff)
• objects/Loads/Friction.cpp (modified) (1 diff)
• objects/Loads/Friction.h (modified) (1 diff)

issm/trunk/src/c/Container/Inputs.cpp

@@ -99 +99 @@
 }
 /*}}}*/
+/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
+void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){
+
+        vector<Object*>::iterator object;
+        Input* input=NULL;
+        bool   found=false;
+
+        /*Go through inputs and check whether any input with the same name is already in: */
+        for ( object=objects.begin() ; object < objects.end(); object++ ){
+
+                input=(Input*)(*object); 
+                if (input->EnumType()==enum_type){
+                        found=true;
+                        break;
+                }
+        }
+
+        if (!found){
+                /*we could not find an input with the correct enum type. No defaults values were provided, 
+                 * error out: */
+                ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
+        }
+
+        /*Ok, we have an input if we made it here, request the input to return the values: */
+        input->GetParameterValue(pvalue,gauss);
+
+}
+/*}}}*/
 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){

issm/trunk/src/c/Container/Inputs.h

@@ -17 +17 @@
 class Node;
 class GaussTria;
+class GaussPenta;
 
 class Inputs: public DataSet{
@@ -49 +50 @@
                 void GetParameterValue(double* pvalue,double* gauss,int enum_type);
                 void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
+                void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);
                 void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type);

issm/trunk/src/c/objects/Elements/Penta.cpp

@@ -2879 +2879 @@
 
                         /*Calculate DL on gauss point */
-                        friction->GetAlpha2(&alpha2_gauss, gauss_tria,VxEnum,VyEnum,VzEnum);
+                        friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 
                         DLStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;

issm/trunk/src/c/objects/Loads/Friction.cpp

@@ -177 +177 @@
 }
 /*}}}*/
+/*FUNCTION Friction::GetAlpha2{{{1*/
+void Friction::GetAlpha2(double* 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: */
+        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")==0){
+                inputs->GetParameterValue(&vx, gauss,vxenum);
+                inputs->GetParameterValue(&vy, gauss,vyenum);
+                vmag=sqrt(pow(vx,2)+pow(vy,2));
+        }
+        else if (strcmp(element_type,"3d")==0){
+                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));
+        }
+        else ISSMERROR("element_type %s not supported yet",element_type);
+
+        alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}
+/*}}}*/
 /*FUNCTION Friction::GetAlphaComplement {{{1*/
 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){

issm/trunk/src/c/objects/Loads/Friction.h

@@ -29 +29 @@
                 void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
                 void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
+                void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
                 void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
                 void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);