[cig-commits] r13650 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src

[walter at geodynamics.org]

Author: walter
Date: 2008-12-11 03:21:21 -0800 (Thu, 11 Dec 2008)
New Revision: 13650

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
 r2414 at dante:  boo | 2008-12-11 03:15:53 -0800
 Add in a hard coded hack for hydrostatic pressure in MohrCoulomb for tibet



Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
   - 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2413
   + 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2414

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-12-11 11:21:16 UTC (rev 13649)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-12-11 11:21:21 UTC (rev 13650)
@@ -282,6 +282,52 @@
 	return fabs(sigma_ns);
 }
 
+double Pressure_Analytic(double density, double gravity, double h,
+                         double alpha, double T_0, double A,
+                         double B, double C)
+{
+  return density*gravity*h*(1-alpha*(T_0 + A*h/2 + B*(1+exp(-C*h)/(C*h))));
+}
+
+double PressureFunction(Coord coord)
+{
+  double T_0=273;
+  double A=0.0096667;
+  double B=93.333;
+  double C=1.0e-4;
+  double y_max=100000;
+  double material_boundary=68000;
+  double gravity=9.81;
+
+  double h=y_max-coord[1];
+  double upper_alpha=3.0e-5;
+  double lower_alpha=3.0e-5;
+  double upper_density=2800;
+  double lower_density=3300;
+  double T, p;
+
+  if(coord[1]>y_max)
+    {
+      p=0;
+    }
+  else if(coord[1]>material_boundary)
+    {
+      p=Pressure_Analytic(upper_density, gravity, h,
+                          upper_alpha, T_0, A, B, C);
+    }
+  else
+    {
+      p=Pressure_Analytic(upper_density, gravity, y_max-material_boundary,
+                          upper_alpha, T_0, A, B, C)
+        - Pressure_Analytic(lower_density, gravity,y_max-material_boundary,
+                            lower_alpha, T_0, A, B, C)
+        + Pressure_Analytic(lower_density, gravity,h,
+                            lower_alpha, T_0, A, B, C);
+    }
+  return p;
+}
+
+
 double _MohrCoulomb_GetYieldCriterion( 
 		void*                                              rheology,
 		ConstitutiveMatrix*                                constitutiveMatrix,
@@ -300,9 +346,16 @@
 	double                            pressure;
         double theta;
         double factor;
+        Cell_Index                       cell_I;
+        Coord coord;
 
 	FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
 
+        cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
+                                               lElement_I );
+        FeMesh_CoordLocalToGlobal(pressureField->feMesh, cell_I, xi, coord);
+        pressure+=PressureFunction(coord);
+
 	/* Calculate frictional strength.  We modify the friction and
            cohesion because we have grouped terms from the normal
            stresses and moved it to the yield indicator. */