Changeset 21388

Timestamp:

11/17/16 16:13:30 (8 years ago)

Author:

hongjuy

Message:

CHG: Use the right coordinate system and bed slope for shelf dampending

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -3045 +3045 @@
         IssmDouble  rho_water     = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
         IssmDouble  gravity       = element->GetMaterialParameter(ConstantsGEnum);
-        Input*      surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
+        Input*      base_input = element->GetInput(BaseEnum); _assert_(base_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -3052 +3052 @@
                 gauss->GaussPoint(ig);
 
-                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+                base_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
                 element->NodalFunctionsVelocity(vbasis,gauss);
                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
@@ -3942 +3942 @@
 
                 for(i=0;i<vnumnodes;i++){
-                        pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
-                        pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
-                        if(dim==3){
-                                pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
-                        }
-                }
-        }
-
-        /*Transform coordinate system*/
-        element->TransformLoadVectorCoord(pe,cs_list);
+                                pe->values[i*dim+(dim-1)]+=-water_pressure*gauss->weight*Jdet*vbasis[i];
+                }
+        }
 
         /* shelf dampening*/
@@ -3964 +3957 @@
                         element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
                         element->NodalFunctionsVelocity(vbasis,gauss);
-                        element->NormalBase(&normal[0],xyz_list_base);
-                        if (dim==2) normal_b=normal[1];
-                        else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
                         mb_input->GetInputValue(&mb, gauss);
                         for(i=0;i<vnumnodes;i++){
-                                pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
+                                pe->values[i*dim+(dim-1)] += -dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i];
                         }
                 }
         }
+
+        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
 
         /*Clean up and return*/