Changeset 16823

Timestamp:

11/18/13 16:03:02 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: adding other sliding laws in SIA

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (3 diffs)
• Tria.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -8447 +8447 @@
         IssmDouble  z_list[NUMVERTICES];
         IssmDouble  slope[2];
-        IssmDouble  slope2,constant_part;
+        IssmDouble  slope2,constant_part,drag;
         IssmDouble  rho_ice,gravity,n,B;
         IssmDouble  ub,vb,z_g,surface,thickness;
@@ -8468 +8468 @@
         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
+        Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
         for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
 
@@ -8512 +8513 @@
                 /*Deal with lower surface*/
                 if(IsOnBed()){
-                        constant_part=-1.58*1.e-10*rho_ice*gravity*thickness;
-                        ub=constant_part*slope[0];
-                        vb=constant_part*slope[1];
+                        drag_input->GetInputValue(&drag,gauss);
+
+                        /*No ref...*/
+                        ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+                        vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+                        ///*Ritz et al. 1996*/
+                        //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+                        //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+                        ///*Rutt et al. 2009*/
+                        //ub=-drag*rho_ice*gravity*thickness*slope[0];
+                        //vb=-drag*rho_ice*gravity*thickness*slope[1];
 
                         pe->values[2*node0+0]+=ub/connectivity[0];

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -4220 +4220 @@
 
         /*Intermediaries */
-        IssmDouble constant_part,ub,vb;
-        IssmDouble rho_ice,gravity,n,B;
+        IssmDouble ub,vb;
+        IssmDouble rho_ice,gravity,n,B,drag;
         IssmDouble slope2,thickness,connectivity;
         IssmDouble slope[2];
@@ -4239 +4239 @@
         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
         Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
+        Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
 
         /*Spawn 3 sing elements: */
@@ -4249 +4250 @@
 
                 thickness_input->GetInputValue(&thickness,gauss);
+                drag_input->GetInputValue(&drag,gauss);
                 slopex_input->GetInputValue(&slope[0],gauss);
                 slopey_input->GetInputValue(&slope[1],gauss);
-                slope2=pow(slope[0],2)+pow(slope[1],2);
-
-                constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
-
+                slope2=slope[0]*slope[0]+slope[1]*slope[1];
+
+                /*No ref...*/
                 ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
                 vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+                ///*Ritz et al. 1996*/
+                //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+                //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+                ///*Rutt et al. 2009*/
+                //ub=-drag*rho_ice*gravity*thickness*slope[0];
+                //vb=-drag*rho_ice*gravity*thickness*slope[1];
 
                 pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;