[cig-commits] commit: Subtract out the trace when computing the trace of the stress and strain rate, because Q1-Q1 elements add some in. Also, use a more robust relation to set the viscosity when yielding.

[Mercurial]

changeset:   817:0e3fa2f27435
tag:         tip
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Feb 10 07:05:49 2010 -0800
files:       Rheology/src/DruckerPrager.c Rheology/src/VonMises.c Rheology/src/VonMises.h
description:
Subtract out the trace when computing the trace of the stress and strain rate, because Q1-Q1 elements add some in.  Also, use a more robust relation to set the viscosity when yielding.


diff -r af3a9bc2f150 -r 0e3fa2f27435 Rheology/src/DruckerPrager.c
--- a/Rheology/src/DruckerPrager.c	Mon Feb 08 12:28:54 2010 -0800
+++ b/Rheology/src/DruckerPrager.c	Wed Feb 10 07:05:49 2010 -0800
@@ -345,14 +345,8 @@ double _DruckerPrager_GetYieldCriterion(
 	double                            frictionCoefficient;
 	double                            frictionCoefficientAfterSoftening;
 	double                            minimumYieldStress;
-	double                            strainWeakeningRatio;
 	double                            effectiveCohesion;
 	double                            effectiveFrictionCoefficient;
-	double                            phi;
-	double                            sinphi;
-	double                            oneOverDenominator;
-	double                            dpFrictionCoefficient;
-	double                            dpCohesion;
 	double                            frictionalStrength;
 	double                            pressure;
 	DruckerPrager_Particle*           particleExt;
@@ -479,15 +473,14 @@ void _DruckerPrager_HasYielded(
 		double                           yieldIndicator )
 {
    DruckerPrager* self = (DruckerPrager*)rheology;
-   double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
+   double old_viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
+   double viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
 
-   _VonMises_HasYielded(
-      self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint,
-      yieldCriterion, yieldIndicator );
+   ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
 
    if( constitutiveMatrix->sle && constitutiveMatrix->sle->nlFormJacobian ) {
 
-      constitutiveMatrix->derivs[8] += viscosity * self->curFrictionCoef / yieldIndicator;
+      constitutiveMatrix->derivs[8] += old_viscosity * self->curFrictionCoef / yieldIndicator;
 
    }
 
diff -r af3a9bc2f150 -r 0e3fa2f27435 Rheology/src/VonMises.c
--- a/Rheology/src/VonMises.c	Mon Feb 08 12:28:54 2010 -0800
+++ b/Rheology/src/VonMises.c	Wed Feb 10 07:05:49 2010 -0800
@@ -230,6 +230,8 @@ double _VonMises_GetYieldIndicator(
 {
 	VonMises*                         self             = (VonMises*) rheology;
 	SymmetricTensor                   strainRate;
+        int i;
+        double stressTrace, strainRateTrace;
 	
 	/* Get Strain Rate */
 	if( self->strainRateField ) {
@@ -243,8 +245,27 @@ double _VonMises_GetYieldIndicator(
 	}
 
 	/* Get Stress */
-	ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, self->stress );
-	self->stressInv = SymmetricTensor_2ndInvariant( self->stress, constitutiveMatrix->dim );
+	ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate,
+                                            self->stress );
+        SymmetricTensor_GetTrace(self->stress, constitutiveMatrix->dim,
+                                 &stressTrace);
+        SymmetricTensor_GetTrace(strainRate, constitutiveMatrix->dim,
+                                 &strainRateTrace);
+
+        for(i=0;i<constitutiveMatrix->dim;++i)
+          {
+            strainRate[TensorMapST3D[i][i]]-=
+              strainRateTrace/constitutiveMatrix->dim;
+            self->stress[TensorMapST3D[i][i]]-=
+              stressTrace/constitutiveMatrix->dim;
+          }
+
+        /* Save the strain rate invariant */
+        self->strainRateSecondInvariant=
+          SymmetricTensor_2ndInvariant( strainRate, constitutiveMatrix->dim );
+
+	self->stressInv = 
+          SymmetricTensor_2ndInvariant( self->stress, constitutiveMatrix->dim );
 
         self->yieldIndicator = self->stressInv;
 
diff -r af3a9bc2f150 -r 0e3fa2f27435 Rheology/src/VonMises.h
--- a/Rheology/src/VonMises.h	Mon Feb 08 12:28:54 2010 -0800
+++ b/Rheology/src/VonMises.h	Wed Feb 10 07:05:49 2010 -0800
@@ -61,7 +61,8 @@
 		double                                              cohesionAfterSoftening;               \
 		Bool                                                strainRateSoftening;                  \
                 SymmetricTensor                                     stress;\
-                double                                              stressInv;
+                double                                              stressInv; \
+                double                                              strainRateSecondInvariant;
 
 	struct VonMises { __VonMises };