Changeset 15041

Timestamp:

05/17/13 11:30:17 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: removed automated penalty when surface slopes > 0.06% (MOUNTAINEXPONENT)

Location:

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

Files:

• Penta.cpp (modified) (5 diffs)
• Tria.cpp (modified) (2 diffs)

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

@@ -6067 +6067 @@
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
-        IssmDouble slope[3]={0.0,0.0,0.0};
-        IssmDouble MAXSLOPE=.06; // 6 %
-        IssmDouble MOUNTAINKEXPONENT=10;
         IssmDouble L[2][numdof];
         IssmDouble DL[2][2]                  ={{ 0,0 },{0,0}}; //for basal drag
@@ -6099 +6096 @@
         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
         Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
@@ -6115 +6111 @@
                 /*Friction: */
                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
-
-                // 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: */
-                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
-                if (slope_magnitude>MAXSLOPE){
-                        alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
-                }
 
                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
@@ -6905 +6892 @@
         IssmDouble slope_magnitude,alpha2,Jdet;
         IssmDouble phi=1.0;
-        IssmDouble slope[3]={0.0,0.0,0.0};
-        IssmDouble MAXSLOPE=.06; // 6 %
-        IssmDouble MOUNTAINKEXPONENT=10;
         IssmDouble L[2][numdof];
         IssmDouble DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
@@ -6944 +6928 @@
                 GetL(&L[0][0], gauss,NDOF2);
 
-                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
                 friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
-                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
-
-                // 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: */
-                if (slope_magnitude>MAXSLOPE){
-                        alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
-                }
 
                 DL_scalar=alpha2*gauss->weight*Jdet;

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

@@ -3319 +3319 @@
         int        analysis_type,migration_style;
         int        gausspoints=2;
-        IssmDouble MAXSLOPE  = .06; // 6 %
-        IssmDouble MOUNTAINKEXPONENT = 10;
-        IssmDouble slope_magnitude,alpha2;
+        IssmDouble alpha2;
         IssmDouble Jdet;
         IssmDouble phi=1.0;
@@ -3363 +3361 @@
                 gauss->GaussPoint(ig);
 
-                // 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: */
-                surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-                slope_magnitude=sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
-                if(slope_magnitude>MAXSLOPE) alpha2=pow((IssmDouble)10,MOUNTAINKEXPONENT);
-                else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
+                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
                 if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
                 if(migration_style==SubelementMigration2Enum){