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

[walter at geodynamics.org]

Author: walter
Date: 2008-12-13 03:57:58 -0800 (Sat, 13 Dec 2008)
New Revision: 13680

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
   long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
Log:
 r2435 at dante:  boo | 2008-12-13 03:55:54 -0800
 Make MohrCoulomb use HydrostaticTerm



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

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-12-13 11:57:56 UTC (rev 13679)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-12-13 11:57:58 UTC (rev 13680)
@@ -116,7 +116,8 @@
 		double                                             cohesionAfterSoftening,
 		double                                             frictionCoefficient,
 		double                                             frictionCoefficientAfterSoftening,
-		double                                             minimumYieldStress)
+		double                                             minimumYieldStress,
+                HydrostaticTerm*                                   hydrostaticTerm)
 {
 	self->materialPointsSwarm     = materialPointsSwarm;
 	self->pressureField           = pressureField;
@@ -134,6 +135,7 @@
 	self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
 
 	self->minimumYieldStress = minimumYieldStress;
+        self->hydrostaticTerm=hydrostaticTerm;
 }
 
 void* _MohrCoulomb_DefaultNew( Name name ) {
@@ -189,7 +191,9 @@
 			Stg_ComponentFactory_GetDouble( cf, self->name, "cohesionAfterSoftening", 0.0 ),
 			Stg_ComponentFactory_GetDouble( cf, self->name, "frictionCoefficient", 0.0 ),
 			Stg_ComponentFactory_GetDouble( cf, self->name, "frictionCoefficientAfterSoftening", 0.0 ),
-			Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ) );
+			Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ),
+                        Stg_ComponentFactory_ConstructByKey( cf, self->name, 
+                                                             "HydrostaticTerm", HydrostaticTerm, False, data ) );
 }
 
 void _MohrCoulomb_Build( void* rheology, void* data ) {
@@ -274,52 +278,6 @@
 	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,
                                       MaterialPointsSwarm* materialPointsSwarm,
@@ -343,10 +301,13 @@
   
   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);
+  if(self->hydrostaticTerm)
+    {
+      cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
+                                             lElement_I );
+      FeMesh_CoordLocalToGlobal(pressureField->feMesh, cell_I, xi, coord);
+      pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,coord);
+    }
   
   /* Calculate frictional strength.  We modify the friction and
      cohesion because we have grouped terms from the normal

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h	2008-12-13 11:57:56 UTC (rev 13679)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h	2008-12-13 11:57:58 UTC (rev 13680)
@@ -74,6 +74,7 @@
 		double                                              currentPressure;                       \
 		double                                              strainRateSecondInvariant;                   \
 		FeVariable*                                         strainRateField;                      \
+                HydrostaticTerm*                                    hydrostaticTerm;
 
 	
 	struct MohrCoulomb { __MohrCoulomb };