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

[walter at geodynamics.org]

Author: walter
Date: 2008-05-22 11:16:14 -0700 (Thu, 22 May 2008)
New Revision: 11996

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
Log:
 r2166 at earth:  boo | 2008-05-22 11:09:57 -0700
 Fix a bug in the factor applied to the yield criterion. Also add commented out support for averaging viscosities between nonlinear iterations.



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

Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-05-22 18:16:10 UTC (rev 11995)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c	2008-05-22 18:16:14 UTC (rev 11996)
@@ -303,7 +303,6 @@
 
 	FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
 
-		
 	/* 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. */
@@ -312,7 +311,7 @@
 
         theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
 
-        factor=1/(sin(2*theta) - effectiveFrictionCoefficient*cos(2*theta));
+        factor=sin(2*theta);
         
         frictionalStrength = effectiveFrictionCoefficient * pressure*factor + effectiveCohesion*factor;
 
@@ -325,8 +324,8 @@
           minimumYieldStress=effectiveCohesion;
 	
 	/* Make sure frictionalStrength is above the minimum */
-	if ( frictionalStrength < minimumYieldStress) 
-		frictionalStrength = minimumYieldStress;
+	if ( frictionalStrength < minimumYieldStress*factor) 
+		frictionalStrength = minimumYieldStress*factor;
 	return frictionalStrength;
 }
 
@@ -336,6 +335,10 @@
    (and thus the viscosity), while the yieldIndicator depends linearly
    on the deviatoric stresses. */
 
+typedef struct  {
+  double effVisc;
+} StoreVisc_ParticleExt;
+
 void _MohrCoulomb_HasYielded( 
 		void*                                              rheology,
 		ConstitutiveMatrix*                                constitutiveMatrix,
@@ -349,6 +352,32 @@
 
   double beta=1.0-yieldCriterion/(yieldIndicator);
 
+/*     { */
+/*       StoreVisc_ParticleExt*            particleExt; */
+      
+/*       /\* Get Parameters From Material Extension *\/ */
+/*       particleExt          = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, 2 ); */
+      
+/* /\*       if((particleExt->effVisc>viscosity*(1-beta) && particleExt->effVisc/viscosity*(1-beta)<1.0e3) *\/ */
+/* /\*          || (particleExt->effVisc<viscosity*(1-beta) && particleExt->effVisc/viscosity*(1-beta)>1.0e-3)) *\/ */
+/* /\*       if(particleExt->effVisc!=viscosity) *\/ */
+/* /\*         { *\/ */
+/* /\*           double new_viscosity; *\/ */
+/* /\*           new_viscosity=sqrt(particleExt->effVisc*viscosity*(1-beta)); *\/ */
+/* /\* /\\*           new_viscosity=(particleExt->effVisc+(viscosity*(1-beta)))/2; *\\/ *\/ */
+/* /\* /\\*           new_viscosity=pow(pow(particleExt->effVisc,1.0)*viscosity*(1-beta),1/2.0); *\\/ *\/ */
+/* /\*           beta=1.0-new_viscosity/viscosity; *\/ */
+/* /\*         } *\/ */
+
+/* /\*       if(lElement_I==340) *\/ */
+/* /\*         printf("has yielded %d %g %g %g\n",lElement_I,particleExt->effVisc,viscosity,viscosity*(1-beta)); *\/ */
+
+
+/*       if(lElement_I>680 && lElement_I<710) */
+/*         printf("has yielded %d %g %g %g %g\n",lElement_I,yieldCriterion, */
+/*                yieldIndicator,particleExt->effVisc,viscosity*(1-beta)); */
+/*     } */
+
   ConstitutiveMatrix_IsotropicCorrection( constitutiveMatrix, -viscosity * beta );
 }