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

walter at geodynamics.org walter at geodynamics.org
Thu Feb 26 12:44:51 PST 2009


Author: walter
Date: 2009-02-26 12:44:51 -0800 (Thu, 26 Feb 2009)
New Revision: 14160

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:
 r2508 at dante:  boo | 2009-02-26 12:43:36 -0800
 Implement transforming boundary layers



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

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2009-02-26 17:48:54 UTC (rev 14159)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2009-02-26 20:44:51 UTC (rev 14160)
@@ -108,16 +108,21 @@
 }
 
 void _MohrCoulomb_Init(
-		MohrCoulomb*                        self,
-		FeVariable*                                        pressureField,
-		FeVariable*                                        velocityGradientsField,
-		MaterialPointsSwarm*                               materialPointsSwarm,
-		double                                             cohesion,
-		double                                             cohesionAfterSoftening,
-		double                                             frictionCoefficient,
-		double                                             frictionCoefficientAfterSoftening,
-		double                                             minimumYieldStress,
-                HydrostaticTerm*                                   hydrostaticTerm)
+                       MohrCoulomb* self,
+                       FeVariable* pressureField,
+                       FeVariable* velocityGradientsField,
+                       MaterialPointsSwarm* materialPointsSwarm,
+                       double cohesion,
+                       double cohesionAfterSoftening,
+                       double frictionCoefficient,
+                       double frictionCoefficientAfterSoftening,
+                       double boundaryCohesion,
+                       double boundaryCohesionAfterSoftening,
+                       double boundaryFrictionCoefficient,
+                       double boundaryFrictionCoefficientAfterSoftening,
+                       Stg_Shape *boundaryShape,
+                       double minimumYieldStress,
+                       HydrostaticTerm* hydrostaticTerm)
 {
 	self->materialPointsSwarm     = materialPointsSwarm;
 	self->pressureField           = pressureField;
@@ -126,13 +131,16 @@
 	self->cohesion = cohesion;
 	self->frictionCoefficient = frictionCoefficient;
 	
-	/* Strain softening of Cohesion - (linear weakening is assumed) */
-	/* needs a softening factor between +0 and 1 and a reference strain > 0 */
+	/* Strain softening of Cohesion and friction - (linear
+           weakening is assumed) needs a softening factor between +0
+           and 1 and a reference strain > 0 */
 	self->cohesionAfterSoftening = cohesionAfterSoftening;
-	
-	/* Strain softening of Friction - (linear weakening is assumed) */
-	/* needs a softening factor between +0 and 1 and a reference strain > 0 */
 	self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
+	self->boundaryCohesion = boundaryCohesion;
+	self->boundaryFrictionCoefficient = boundaryFrictionCoefficient;
+	self->boundaryCohesionAfterSoftening = boundaryCohesionAfterSoftening;
+	self->boundaryFrictionCoefficientAfterSoftening = boundaryFrictionCoefficientAfterSoftening;
+        self->boundaryShape=boundaryShape;
 
 	self->minimumYieldStress = minimumYieldStress;
         self->hydrostaticTerm=hydrostaticTerm;
@@ -191,6 +199,11 @@
 			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, "boundaryCohesion", 0.0 ),
+			Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryCohesionAfterSoftening", 0.0 ),
+			Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryFrictionCoefficient", 0.0 ),
+			Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryFrictionCoefficientAfterSoftening", 0.0 ),
+                        Stg_ComponentFactory_ConstructByKey(  cf,  self->name,  "boundaryShape", Stg_Shape,  False, data  ),
 			Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ),
                         Stg_ComponentFactory_ConstructByKey( cf, self->name, 
                                                              "HydrostaticTerm", HydrostaticTerm, False, data ) );
@@ -294,18 +307,17 @@
   
   FeVariable*                       pressureField  = self->pressureField;  
   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);
   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);
     }
   
@@ -313,11 +325,24 @@
      cohesion because we have grouped terms from the normal
      stresses and moved it to the yield indicator. */
   effectiveFrictionCoefficient =
-    _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint );
-  effectiveCohesion = _MohrCoulomb_EffectiveCohesion( self, materialPoint );
+    _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint, coord );
+  effectiveCohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint,coord);
   
-  theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
-  factor=sin(2*theta);
+  if(self->boundaryShape
+     && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
+    {
+      factor=1;
+/*       printf("strength %g %g %g %g %g %g %g\n",coord[0],coord[1], */
+/*              frictionalStrength, */
+/*              factor,effectiveFrictionCoefficient,pressure, */
+/*              effectiveCohesion); */
+    }
+  else
+    {
+      double theta;
+      theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
+      factor=sin(2*theta);
+    }
   
   if(dim==2)
     {
@@ -346,6 +371,7 @@
   /* Make sure frictionalStrength is above the minimum */
   if ( frictionalStrength < minimumYieldStress*factor) 
     frictionalStrength = minimumYieldStress*factor;
+
   return frictionalStrength;
 }
 
@@ -370,26 +396,48 @@
   ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
 }
 
-double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) {
+double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint,
+                                       double *coord ) {
 	MohrCoulomb*    self                          = (MohrCoulomb*) rheology;
 	double                         effectiveCohesion;
 	
         double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
 	
-        effectiveCohesion =  self->cohesion * (1.0 - strainWeakeningRatio) + 
-          self->cohesionAfterSoftening * strainWeakeningRatio;
+        if(self->boundaryShape
+           && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
+          {
+            effectiveCohesion=self->boundaryCohesion*(1.0-strainWeakeningRatio)+
+              self->boundaryCohesionAfterSoftening*strainWeakeningRatio;
+          }
+        else
+          {
+            effectiveCohesion = self->cohesion * (1.0 - strainWeakeningRatio) + 
+              self->cohesionAfterSoftening * strainWeakeningRatio;
+          }
 
 	return effectiveCohesion;
 }
 
-double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) {
+double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, double *coord ) {
 	MohrCoulomb*    self                          = (MohrCoulomb*) rheology;
-	static double                  effectiveFrictionCoefficient  = 0.0;
+	double                  effectiveFrictionCoefficient  = 0.0;
 
         double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
 	
-        effectiveFrictionCoefficient = self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
-          self->frictionCoefficientAfterSoftening * strainWeakeningRatio;	
+        if(self->boundaryShape
+           && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
+          {
+            effectiveFrictionCoefficient =
+              self->boundaryFrictionCoefficient * (1.0 - strainWeakeningRatio) +
+              self->boundaryFrictionCoefficientAfterSoftening
+              *strainWeakeningRatio;	
+          }
+        else
+          {
+            effectiveFrictionCoefficient =
+              self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
+              self->frictionCoefficientAfterSoftening * strainWeakeningRatio;	
+          }
 
 	return effectiveFrictionCoefficient;
 }

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h	2009-02-26 17:48:54 UTC (rev 14159)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h	2009-02-26 20:44:51 UTC (rev 14160)
@@ -66,6 +66,11 @@
 		double                                              cohesionAfterSoftening;                \
 		double                                              frictionCoefficient;                   \
 		double                                              frictionCoefficientAfterSoftening;     \
+		double                                              boundaryCohesion;                              \
+		double                                              boundaryCohesionAfterSoftening;                \
+		double                                              boundaryFrictionCoefficient;                   \
+		double                                              boundaryFrictionCoefficientAfterSoftening;     \
+                Stg_Shape*                                          boundaryShape;                          \
 		double                                              minimumYieldStress;                    \
 		/* Stored values that are calculated once for each particle */ \
 		Eigenvector                                         currentEigenvectorList[3];             \
@@ -140,8 +145,8 @@
 		double                                             yieldCriterion,
 		double                                             yieldIndicator );
 	
-	double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) ;
-	double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) ;
+	double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint, double *coord ) ;
+	double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, double *coord ) ;
 	double _MohrCoulomb_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) ;
 
 	void _MohrCoulomb_StoreCurrentParameters( 



More information about the CIG-COMMITS mailing list