[cig-commits] r11823 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Thu Apr 17 10:17:38 PDT 2008
Author: walter
Date: 2008-04-17 10:17:38 -0700 (Thu, 17 Apr 2008)
New Revision: 11823
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
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta
Log:
r2096 at earth: boo | 2008-04-17 10:17:36 -0700
Rough, ugly fix for MohrCoulomb yield indicator
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2095
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2096
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-04-17 17:17:31 UTC (rev 11822)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2008-04-17 17:17:38 UTC (rev 11823)
@@ -111,6 +111,7 @@
MohrCoulomb* self,
FeVariable* pressureField,
FeVariable* velocityGradientsField,
+ FeVariable* strainRateField,
MaterialPointsSwarm* materialPointsSwarm,
double cohesion,
double cohesionAfterSoftening,
@@ -124,6 +125,7 @@
self->materialPointsSwarm = materialPointsSwarm;
self->pressureField = pressureField;
self->velocityGradientsField = velocityGradientsField;
+ self->strainRateField = strainRateField;
self->cohesion = cohesion;
self->frictionCoefficient = frictionCoefficient;
@@ -165,6 +167,7 @@
FeVariable* pressureField;
MaterialPointsSwarm* materialPointsSwarm;
FeVariable* velocityGradientsField;
+ FeVariable* strainRateField;
/* Construct Parent */
_YieldRheology_Construct( self, cf, data );
@@ -182,11 +185,14 @@
"PressureField", FeVariable, True, data );
velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
"VelocityGradientsField", FeVariable, True, data );
+ strainRateField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "StrainRateField", FeVariable, True, data );
_MohrCoulomb_Init(
self,
pressureField,
velocityGradientsField,
+ strainRateField,
materialPointsSwarm,
Stg_ComponentFactory_GetDouble( cf, self->name, "cohesion", 0.0 ),
Stg_ComponentFactory_GetDouble( cf, self->name, "cohesionAfterSoftening", 0.0 ),
@@ -266,24 +272,40 @@
Coord xi )
{
MohrCoulomb* self = (MohrCoulomb*) rheology;
- double* stress = self->currentStress;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double sigma_ns;
- Dimension_Index dim = constitutiveMatrix->dim;
- double stressMin;
- double stressMax;
- double theta;
+/* double* stress = self->currentStress; */
+/* Eigenvector* eigenvectorList = self->currentEigenvectorList; */
+/* double sigma_ns; */
+/* Dimension_Index dim = constitutiveMatrix->dim; */
+/* double stressMin; */
+/* double stressMax; */
+/* double theta; */
- /* Failure criterion from stress field - we can calculate
- * whether the material has failed at the current stress from
- * the principle stresses */
+/* /\* Failure criterion from stress field - we can calculate */
+/* * whether the material has failed at the current stress from */
+/* * the principle stresses *\/ */
- stressMin = eigenvectorList[0].eigenvalue;
- stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
+/* stressMin = eigenvectorList[0].eigenvalue; */
+/* stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue); */
- sigma_ns = 0.5 * ( stressMax - stressMin );
+/* sigma_ns = 0.5 * ( stressMax - stressMin ); */
- return fabs(sigma_ns);
+/* return fabs(sigma_ns); */
+
+ SymmetricTensor stress;
+ SymmetricTensor strainRate;
+ double stressInv;
+
+ /* Get Strain Rate */
+ FeVariable_InterpolateWithinElement( self->strainRateField, lElement_I, xi, strainRate );
+
+ /* Get Stress */
+ ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate, stress );
+ stressInv = SymmetricTensor_2ndInvariant( stress, constitutiveMatrix->dim );
+
+/* printf("stress %lf %lf %lf\n",stressInv,xi[0],xi[1]); */
+
+ return stressInv;
+
}
double _MohrCoulomb_GetYieldCriterion(
@@ -333,8 +355,16 @@
MohrCoulomb* self = (MohrCoulomb*)rheology;
Eigenvector* eigenvectorList = self->currentEigenvectorList;
Dimension_Index dim = constitutiveMatrix->dim;
- double viscosity;
+ double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
+/* printf("yielded %lf %lf %lf %lf %lf %lf\n", */
+/* viscosity, yieldCriterion, */
+/* yieldIndicator, self->strainRateSecondInvariant, */
+/* viscosity * yieldCriterion/yieldIndicator, */
+/* yieldCriterion/(2*self->strainRateSecondInvariant)); */
+
+/* viscosity = viscosity * yieldCriterion/yieldIndicator; */
+
viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2008-04-17 17:17:31 UTC (rev 11822)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2008-04-17 17:17:38 UTC (rev 11823)
@@ -73,6 +73,7 @@
SymmetricTensor currentStress; \
double currentPressure; \
double strainRateSecondInvariant; \
+ FeVariable* strainRateField; \
struct MohrCoulomb { __MohrCoulomb };
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2008-04-17 17:17:31 UTC (rev 11822)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2008-04-17 17:17:38 UTC (rev 11823)
@@ -75,6 +75,12 @@
<param name="Description">Define the weakening behaviour of material parameters. This component dependency is first defined in YieldRheology class, but as non essential. Since we want it for MC criterion, we make it essential here.</param>
</struct>
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Name">StrainRateField</param>
+ <param name="Type">FeVariable</param>
+ <param name="Description">The field that provides the $\dot\epsilon$ in the equation above.</param>
+ </struct>
</list>
<!-- Add an example XML if possible -->
More information about the cig-commits
mailing list