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

[walter at geodynamics.org]

Author: walter
Date: 2009-03-07 19:22:52 -0800 (Sat, 07 Mar 2009)
New Revision: 14251

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
 r2538 at dante:  boo | 2009-03-07 19:22:28 -0800
 Add the trace of the stress to the pressure in Mohr Coulomb.  This is needed since we are using linear elements



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

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2009-03-08 01:27:34 UTC (rev 14250)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2009-03-08 03:22:52 UTC (rev 14251)
@@ -329,7 +329,7 @@
   unsigned		inds[3];
   Grid*			elGrid;
   Bool inside_boundary;
-  
+  Eigenvector* eigenvectorList  = self->currentEigenvectorList;
   FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
   
   cell_I=CellLayout_MapElementIdToCellId(materialPointsSwarm->cellLayout,
@@ -340,6 +340,19 @@
       pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,coord);
     }
   
+  /* Add the average of the trace, since with linear pressure elements
+     we are not guaranteed to have a stress field with zero stress */
+
+  if(dim==2)
+    {
+      pressure=(eigenvectorList[0].eigenvalue+eigenvectorList[1].eigenvalue)/2;
+    }
+  else
+    {
+      pressure=(eigenvectorList[0].eigenvalue+eigenvectorList[1].eigenvalue
+                +eigenvectorList[2].eigenvalue)/2;
+    }
+
   /* 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. */