[cig-commits] r14450 - in long/3D/Gale/trunk: . input/benchmarks input/cookbook input/examples src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Wed Mar 25 12:32:45 PDT 2009
Author: walter
Date: 2009-03-25 12:32:43 -0700 (Wed, 25 Mar 2009)
New Revision: 14450
Added:
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta
Removed:
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta
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
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/benchmarks/geomod_2004_extension.xml
long/3D/Gale/trunk/input/benchmarks/geomod_2004_shortening.xml
long/3D/Gale/trunk/input/benchmarks/geomod_2008_brittle.xml
long/3D/Gale/trunk/input/benchmarks/geomod_2008_stable.xml
long/3D/Gale/trunk/input/benchmarks/geomod_2008_unstable.xml
long/3D/Gale/trunk/input/benchmarks/isostasy.xml
long/3D/Gale/trunk/input/benchmarks/mohr_coulomb.xml
long/3D/Gale/trunk/input/benchmarks/shear_kinetic_friction.xml
long/3D/Gale/trunk/input/benchmarks/shear_kinetic_transform.xml
long/3D/Gale/trunk/input/benchmarks/shear_static_friction.xml
long/3D/Gale/trunk/input/benchmarks/shear_static_transform.xml
long/3D/Gale/trunk/input/cookbook/yielding.xml
long/3D/Gale/trunk/input/examples/rifting.xml
long/3D/Gale/trunk/input/examples/tibet.xml
long/3D/Gale/trunk/input/examples/tibet3D.xml
long/3D/Gale/trunk/input/examples/weak_zone.xml
long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/SConscript
long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h
Log:
r2597 at dante: boo | 2009-03-25 12:32:22 -0700
Rename MohrCoulomb to DruckerPrager
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2588
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2597
Modified: long/3D/Gale/trunk/input/benchmarks/geomod_2004_extension.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/geomod_2004_extension.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/geomod_2004_extension.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -351,7 +351,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/geomod_2004_shortening.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/geomod_2004_shortening.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/geomod_2004_shortening.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -355,7 +355,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="sandYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
@@ -390,7 +390,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="beadsYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/geomod_2008_brittle.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/geomod_2008_brittle.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/geomod_2008_brittle.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -355,7 +355,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="quartzYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
@@ -390,7 +390,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="corundumYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/geomod_2008_stable.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/geomod_2008_stable.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/geomod_2008_stable.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -378,7 +378,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="quartzYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/geomod_2008_unstable.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/geomod_2008_unstable.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/geomod_2008_unstable.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -355,7 +355,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="quartzYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
@@ -390,7 +390,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="corundumYielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/isostasy.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/isostasy.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/isostasy.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -350,7 +350,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/mohr_coulomb.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/mohr_coulomb.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/mohr_coulomb.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -404,7 +404,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/shear_kinetic_friction.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/shear_kinetic_friction.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/shear_kinetic_friction.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -354,7 +354,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">NodalPressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
@@ -427,11 +427,11 @@
<param name="VelocityField">VelocityField</param>
<param name="friction_type">double</param>
<!-- Use this value for nx=64, ny=16 -->
-<!-- <param name="friction_value">0.531914893617</param> -->
+ <param name="friction_value">0.531914893617</param>
<!-- Use this value for nx=128, ny=32 -->
<!-- <param name="friction_value">0.515463917526</param> -->
<!-- Use this value for nx=ny=infinity -->
- <param name="friction_value">0.50</param>
+<!-- <param name="friction_value">0.50</param> -->
<param name="vx_type">double</param>
<param name="vx_value">0.0</param>
</struct>
@@ -524,12 +524,12 @@
<param name="type">WallVC</param>
<param name="wall">left</param>
<list name="variables">
+<!-- <struct> -->
+<!-- <param name="name">vx</param> -->
+<!-- <param name="type">func</param> -->
+<!-- <param name="value">Velocity_SimpleShear</param> -->
+<!-- </struct> -->
<struct>
- <param name="name">vx</param>
- <param name="type">func</param>
- <param name="value">Velocity_SimpleShear</param>
- </struct>
- <struct>
<param name="name">vy</param>
<param name="type">double</param>
<param name="value">0</param>
Modified: long/3D/Gale/trunk/input/benchmarks/shear_kinetic_transform.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/shear_kinetic_transform.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/shear_kinetic_transform.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -357,7 +357,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/benchmarks/shear_static_friction.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/shear_static_friction.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/shear_static_friction.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,15 +1,6 @@
<?xml version="1.0"?>
<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
<struct name="components">
- <struct name="mesh-constant">
- <param name="Type">FeMesh</param>
- <param name="elementType">constant</param>
- </struct>
- <struct name="constantMesh-generator">
- <param name="Type">C0Generator</param>
- <param name="mesh">mesh-constant</param>
- <param name="elementMesh">mesh-linear</param>
- </struct>
<struct name="mesh-linear">
<param name="Type">FeMesh</param>
<param name="elementType">linear</param>
@@ -105,44 +96,24 @@
</struct>
<struct name="pressure">
<param name="Type">MeshVariable</param>
- <param name="mesh">mesh-constant</param>
+ <param name="mesh">mesh-linear</param>
<param name="Rank">Scalar</param>
<param name="DataType">Double</param>
</struct>
<struct name="pressureDofLayout">
<param name="Type">DofLayout</param>
- <param name="mesh">mesh-constant</param>
+ <param name="mesh">mesh-linear</param>
<list name="BaseVariables">
<param>pressure</param>
</list>
</struct>
<struct name="PressureField">
<param name="Type">FeVariable</param>
- <param name="FEMesh">mesh-constant</param>
+ <param name="FEMesh">mesh-linear</param>
<param name="DofLayout">pressureDofLayout</param>
<param name="LinkedDofInfo">pressureLinkedDofs</param>
</struct>
- <struct name="nodalPressure">
- <param name="Type">MeshVariable</param>
- <param name="mesh">mesh-linear</param>
- <param name="Rank">Scalar</param>
- <param name="DataType">Double</param>
- </struct>
- <struct name="nodalPressureDofLayout">
- <param name="Type">DofLayout</param>
- <param name="mesh">mesh-linear</param>
- <list name="BaseVariables">
- <param>nodalPressure</param>
- </list>
- </struct>
- <struct name="NodalPressureField">
- <param name="Type">FeVariable</param>
- <param name="FEMesh">mesh-linear</param>
- <param name="DofLayout">nodalPressureDofLayout</param>
- <param name="LinkedDofInfo">nodalPressureLinkedDofs</param>
- </struct>
-
<struct name="StressField">
<param name="Type">StressField</param>
<param name="StrainRateField">StrainRateField</param>
@@ -329,6 +300,20 @@
<param name="nonLinearTolerance">nonLinearTolerance</param>
<param name="makeConvergenceFile">false</param>
</struct>
+ <struct name="c_matrix">
+ <param name="Type">StiffnessMatrix</param>
+ <param name="RowVariable">PressureField</param>
+ <param name="ColumnVariable">PressureField</param>
+ <param name="RHS">cont_force</param>
+ <param name="allowZeroElementContributions">True</param>
+ </struct>
+ <struct name="mixedStabiliser">
+ <param name="Type">MixedStabiliserTerm</param>
+ <param name="Swarm">gaussSwarm</param>
+ <param name="picSwarm">picIntegrationPoints</param>
+ <param name="storeVisc">storeViscosity</param>
+ <param name="StiffnessMatrix">c_matrix</param>
+ </struct>
<struct name="buoyancyForceTerm">
<param name="Type">BuoyancyForceTerm</param>
<param name="ForceVector">mom_force</param>
@@ -338,6 +323,24 @@
<struct name="backgroundShape">
<param name="Type">Everywhere</param>
</struct>
+
+ <struct name="boundaryShape">
+ <param name="Type">Box</param>
+ <param name="startX">minX</param>
+ <param name="endX">maxX</param>
+ <param name="startY">minY</param>
+ <param name="endY">boundaryUpper</param>
+ <param name="startZ">minZ</param>
+ <param name="endZ">maxZ</param>
+ </struct>
+ <struct name="pdmsShape">
+ <param name="Type">Intersection</param>
+ <list name="shapes">
+ <param>backgroundShape</param>
+ <param>!boundaryShape</param>
+ </list>
+ </struct>
+
<struct name="pdmsViscosity">
<param name="Type">MaterialViscosity</param>
<param name="eta0">1.0</param>
@@ -354,8 +357,8 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
- <param name="PressureField">NodalPressureField</param>
+ <param name="Type">DruckerPrager</param>
+ <param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
<param name="Context">context</param>
@@ -375,7 +378,7 @@
</struct>
<struct name="pdms">
<param name="Type">RheologyMaterial</param>
- <param name="Shape">backgroundShape</param>
+ <param name="Shape">pdmsShape</param>
<param name="density">1.0</param>
<list name="Rheology">
<param>pdmsViscosity</param>
@@ -384,6 +387,36 @@
<param>storeStress</param>
</list>
</struct>
+
+ <struct name="boundaryStrainWeakening">
+ <param name="Type">StrainWeakening</param>
+ <param name="TimeIntegrator">timeIntegrator</param>
+ <param name="MaterialPointsSwarm">materialSwarm</param>
+ <param name="initialSofteningStrain">0.5</param>
+ <param name="finalSofteningStrain">1.0</param>
+ <param name="initialDamageFraction">0.0</param>
+ <param name="initialDamageWavenumber">0.5</param>
+ <param name="initialDamageFactor">0.5</param>
+ <param name="healingRate">0.0</param>
+ </struct>
+ <struct name="boundaryYielding">
+ <param name="Type">DruckerPrager</param>
+ <param name="PressureField">PressureField</param>
+ <param name="VelocityGradientsField">VelocityGradientsField</param>
+ <param name="MaterialPointsSwarm">materialSwarm</param>
+ <param name="Context">context</param>
+ <param name="StrainWeakening">boundaryStrainWeakening</param>
+ <param name="cohesion">1.0</param>
+ <param name="cohesionAfterSoftening">1.0</param>
+ <param name="frictionCoefficient">0.286745385759</param>
+ <param name="frictionCoefficientAfterSoftening">0.249328002843</param>
+ </struct>
+ <struct name="boundaryViscosity">
+ <param name="Type">MaterialViscosity</param>
+ <param name="eta0">1.0e4</param>
+ </struct>
+
+
<struct name="escapedRoutine">
<param name="Type">EscapedRoutine</param>
<param name="idealParticleCount">0</param>
@@ -398,13 +431,6 @@
<param>true</param>
</list>
</struct>
- <struct name="pressureRemesher">
- <param name="Type">CellRemesher</param>
- <param name="mesh">mesh-constant</param>
- <param name="meshType">regular</param>
- <param name="dim">3</param>
- <param name="cellMesh">mesh-linear</param>
- </struct>
<struct name="matrixSolver">
<param name="Type"> PETScMGSolver </param>
<param name="levels"> mgLevels </param>
@@ -448,17 +474,6 @@
<struct name="EulerDeform">
<list name="systems">
<struct>
- <param name="mesh">mesh-constant</param>
- <param name="remesher">pressureRemesher</param>
- <param name="velocityField">VelocityField</param>
- <list name="fields">
- <struct>
- <param name="field">PressureField</param>
- <param name="variable">pressure</param>
- </struct>
- </list>
- </struct>
- <struct>
<param name="mesh">mesh-linear</param>
<param name="remesher">velocityRemesher</param>
<param name="velocityField">VelocityField</param>
@@ -472,6 +487,10 @@
<param name="field">VelocityField</param>
<param name="variable">velocity</param>
</struct>
+ <struct>
+ <param name="field">PressureField</param>
+ <param name="variable">pressure</param>
+ </struct>
</list>
</struct>
</list>
@@ -564,41 +583,9 @@
</list>
</struct>
- <struct>
- <param name="type">StaticFrictionVC</param>
- <param name="wall">bottom</param>
- <param name="StaticFriction">2.0643</param>
- <param name="StickOnFirstStep">True</param>
- <param name="includeUpperX">False</param>
- <param name="StressField">StressField</param>
- <param name="PressureField">NodalPressureField</param>
- <list name="variables">
- <struct>
- <param name="name">vx</param>
- <param name="type">double</param>
- <param name="value">0.0</param>
- </struct>
- </list>
- </struct>
</list>
</struct>
- <struct name="stressICs">
- <param name="type">CompositeVC</param>
- <list name="vcList">
- <struct>
- <param name="type">AllNodesVC</param>
- <list name="variables">
- <struct>
- <param name="name">stress</param>
- <param name="type">double</param>
- <param name="value">0.0</param>
- </struct>
- </list>
- </struct>
- </list>
- </struct>
-
<param name="SimpleShearFactor">1.0</param>
<param name="SimpleShearCentreY">0</param>
<!-- <param name="checkpointEvery">1</param> -->
Modified: long/3D/Gale/trunk/input/benchmarks/shear_static_transform.xml
===================================================================
--- long/3D/Gale/trunk/input/benchmarks/shear_static_transform.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/benchmarks/shear_static_transform.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -357,7 +357,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/cookbook/yielding.xml
===================================================================
--- long/3D/Gale/trunk/input/cookbook/yielding.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/cookbook/yielding.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -403,7 +403,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/examples/rifting.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/rifting.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/examples/rifting.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -400,7 +400,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/examples/tibet.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/tibet.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/examples/tibet.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -453,7 +453,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/examples/tibet3D.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/tibet3D.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/examples/tibet3D.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -438,7 +438,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Modified: long/3D/Gale/trunk/input/examples/weak_zone.xml
===================================================================
--- long/3D/Gale/trunk/input/examples/weak_zone.xml 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/input/examples/weak_zone.xml 2009-03-25 19:32:43 UTC (rev 14450)
@@ -451,7 +451,7 @@
<param name="healingRate">0.0</param>
</struct>
<struct name="yielding">
- <param name="Type">MohrCoulomb</param>
+ <param name="Type">DruckerPrager</param>
<param name="PressureField">PressureField</param>
<param name="VelocityGradientsField">VelocityGradientsField</param>
<param name="MaterialPointsSwarm">materialSwarm</param>
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,444 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-** Copyright (c) 2005, Monash Cluster Computing
-** All rights reserved.
-** Redistribution and use in source and binary forms, with or without modification,
-** are permitted provided that the following conditions are met:
-**
-** * Redistributions of source code must retain the above copyright notice,
-** this list of conditions and the following disclaimer.
-** * Redistributions in binary form must reproduce the above copyright
-** notice, this list of conditions and the following disclaimer in the
-** documentation and/or other materials provided with the distribution.
-** * Neither the name of the Monash University nor the names of its contributors
-** may be used to endorse or promote products derived from this software
-** without specific prior written permission.
-**
-** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
-** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
-** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-**
-**
-** Contact:
-*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
-*%
-** Contributors:
-*+ Robert Turnbull
-*+ Vincent Lemiale
-*+ Louis Moresi
-*+ David May
-*+ David Stegman
-*+ Mirko Velic
-*+ Patrick Sunter
-*+ Julian Giordani
-*+
-** $Id$
-**
-**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
-#include <mpi.h>
-#include <StGermain/StGermain.h>
-#include <StgFEM/StgFEM.h>
-#include <PICellerator/PICellerator.h>
-
-#include "types.h"
-#include "RheologyClass.h"
-#include "StrainWeakening.h"
-#include "YieldRheology.h"
-#include "VonMises.h"
-#include "DruckerPrager.h"
-#include "ConstitutiveMatrix.h"
-
-#include <assert.h>
-
-/* Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
-const Type DruckerPrager_Type = "DruckerPrager";
-
-/* Private Constructor: This will accept all the virtual functions for this class as arguments. */
-DruckerPrager* _DruckerPrager_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
- YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
- YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
- YieldRheology_HasYieldedFunction* _hasYielded,
- Name name )
-{
- DruckerPrager* self;
-
- /* Call private constructor of parent - this will set virtual functions of parent and continue up the hierarchy tree. At the beginning of the tree it will allocate memory of the size of object and initialise all the memory to zero. */
- assert( sizeOfSelf >= sizeof(DruckerPrager) );
- self = (DruckerPrager*) _VonMises_New(
- sizeOfSelf,
- type,
- _delete,
- _print,
- _copy,
- _defaultConstructor,
- _construct,
- _build,
- _initialise,
- _execute,
- _destroy,
- _modifyConstitutiveMatrix,
- _getYieldCriterion,
- _getYieldIndicator,
- _hasYielded,
- name );
-
- /* Function pointers for this class that are not on the parent class should be set here */
-
- return self;
-}
-
-void _DruckerPrager_Init(
- DruckerPrager* self,
- FeVariable* pressureField,
- MaterialPointsSwarm* materialPointsSwarm,
- FiniteElementContext* context,
- double frictionCoefficient,
- double frictionCoefficientAfterSoftening )
-{
- DruckerPrager_Particle* particleExt;
- StandardParticle materialPoint;
-
- self->particleExtHandle = ExtensionManager_Add( materialPointsSwarm->particleExtensionMgr,
- DruckerPrager_Type, sizeof(DruckerPrager_Particle) );
-
- /* Assign Pointers */
- self->pressureField = pressureField;
-
- self->frictionCoefficient = frictionCoefficient;
-
- /* Strain softening of Friction - (linear weakening is assumed) */
- /* needs a softening factor between +0 and 1 and a reference strain > 0 */
- self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
-
- /* Update Drawing Parameters */
- EP_PrependClassHook( Context_GetEntryPoint( context, AbstractContext_EP_DumpClass ),
- _DruckerPrager_UpdateDrawParameters, self );
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, &materialPoint, self->particleExtHandle );
-
- /* Setup Variables for Visualisation */
- self->brightness = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "DruckerPragerBrightness",
- (ArithPointer) &particleExt->brightness - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- self->opacity = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "DruckerPragerOpacity",
- (ArithPointer) &particleExt->opacity - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- self->diameter = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "DruckerPragerDiameter",
- (ArithPointer) &particleExt->diameter - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- /* The tensileFailure variable allows to check whether a materialPoint has failed in tensile mode or not */
- self->tensileFailure = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "DruckerPragerTensileFailure",
- (ArithPointer) &particleExt->tensileFailure - (ArithPointer) &materialPoint,
- Variable_DataType_Char );
-
-}
-
-void* _DruckerPrager_DefaultNew( Name name ) {
- return (void*) _DruckerPrager_New(
- sizeof(DruckerPrager),
- DruckerPrager_Type,
- _YieldRheology_Delete,
- _YieldRheology_Print,
- _YieldRheology_Copy,
- _DruckerPrager_DefaultNew,
- _DruckerPrager_Construct,
- _DruckerPrager_Build,
- _DruckerPrager_Initialise,
- _YieldRheology_Execute,
- _YieldRheology_Destroy,
- _YieldRheology_ModifyConstitutiveMatrix,
- _DruckerPrager_GetYieldCriterion,
- _VonMises_GetYieldIndicator,
- _VonMises_HasYielded,
- name );
-}
-
-void _DruckerPrager_Construct( void* druckerPrager, Stg_ComponentFactory* cf, void* data ){
- DruckerPrager* self = (DruckerPrager*)druckerPrager;
- FeVariable* pressureField;
- MaterialPointsSwarm* materialPointsSwarm;
- FiniteElementContext* context;
-
- /* Construct Parent */
- _VonMises_Construct( self, cf, data );
-
- /* TODO: KeyFallback soon to be deprecated/updated */
- pressureField = Stg_ComponentFactory_ConstructByNameWithKeyFallback(
- cf, self->name, "PressureField", "PressureField", FeVariable, True, data );
- /*pressureField = Stg_ComponentFactory_ConstructByKey(
- cf, self->name, "PressureField", FeVariable, True );*/
- materialPointsSwarm = (MaterialPointsSwarm*)
- Stg_ComponentFactory_ConstructByKey( cf, self->name, "MaterialPointsSwarm", MaterialPointsSwarm, True, data );
-
- context = (FiniteElementContext*)
- Stg_ComponentFactory_ConstructByName( cf, "context", FiniteElementContext, True, data );
-
- _DruckerPrager_Init( self,
- pressureField,
- materialPointsSwarm,
- context,
- Stg_ComponentFactory_GetDouble( cf, self->name, "frictionCoefficient", 0.0 ),
- Stg_ComponentFactory_GetDouble( cf, self->name, "frictionCoefficientAfterSoftening", 0.0 ) );
-}
-
-void _DruckerPrager_Build( void* rheology, void* data ) {
- DruckerPrager* self = (DruckerPrager*) rheology;
-
- /* Build parent */
- _YieldRheology_Build( self, data );
-
- Stg_Component_Build( self->brightness, data, False );
- Stg_Component_Build( self->opacity, data, False );
- Stg_Component_Build( self->diameter, data, False );
- Stg_Component_Build( self->tensileFailure, data, False );
-
-}
-
-
-void _DruckerPrager_Initialise( void* rheology, void* data ) {
- DruckerPrager* self = (DruckerPrager*) rheology;
- Particle_Index lParticle_I;
- Particle_Index particleLocalCount;
- AbstractContext* context = (AbstractContext*)data;
-
- _YieldRheology_Initialise( self, data );
-
- /* Initialise variables that I've created - (mainly just SwarmVariables)
- * This will run a Variable_Update for us */
- Stg_Component_Initialise( self->brightness, data, False );
- Stg_Component_Initialise( self->opacity, data, False );
- Stg_Component_Initialise( self->diameter, data, False );
- Stg_Component_Initialise( self->tensileFailure, data, False );
-
- /* We should only set initial conditions if in regular non-restart mode. If in restart mode, then
- the particle-based variables will be set correcty when we re-load the Swarm. */
- if ( !(context && (True == context->loadFromCheckPoint)) ) {
- /* We don't need to Initialise hasYieldedVariable because it's a parent variable and _YieldRheology_Initialise
- * has already been called */
- particleLocalCount = self->hasYieldedVariable->variable->arraySize;
-
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
-
- Variable_SetValueFloat( self->brightness->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->diameter->variable, lParticle_I, 0.0 );
-
- Variable_SetValueChar( self->hasYieldedVariable->variable, lParticle_I, False );
- Variable_SetValueChar( self->tensileFailure->variable,lParticle_I, False );
- }
- }
-}
-
-double _DruckerPrager_GetYieldCriterion(
- void* druckerPrager,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi )
-{
- DruckerPrager* self = (DruckerPrager*) druckerPrager;
- double cohesion;
- double cohesionAfterSoftening;
- double frictionCoefficient;
- double frictionCoefficientAfterSoftening;
- double minimumYieldStress;
- double strainWeakeningRatio;
- double effectiveCohesion;
- double effectiveFrictionCoefficient;
- double phi;
- double sinphi;
- double oneOverDenominator;
- double dpFrictionCoefficient;
- double dpCohesion;
- double frictionalStrength;
- FeVariable* pressureField = self->pressureField;
- double pressure;
- DruckerPrager_Particle* particleExt;
-
- /* Get Parameters From Material Extension */
- cohesion = self->cohesion;
- cohesionAfterSoftening = self->cohesionAfterSoftening;
- frictionCoefficient = self->frictionCoefficient;
- frictionCoefficientAfterSoftening = self->frictionCoefficientAfterSoftening;
- minimumYieldStress = self->minimumYieldStress;
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
-
- FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
-
- /* Strain softening of yield stress - if the strain weakening is not defined then this function returns
- * a 0 value */
- strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- effectiveCohesion = cohesion * (1.0 - strainWeakeningRatio) +
- cohesionAfterSoftening * strainWeakeningRatio;
-
- effectiveFrictionCoefficient = frictionCoefficient * (1.0 - strainWeakeningRatio) +
- frictionCoefficientAfterSoftening * strainWeakeningRatio;
-
- phi = atan(effectiveFrictionCoefficient);
- sinphi = sin(phi);
- oneOverDenominator = 1.0 / (sqrt(3.0) * (3.0 - sinphi));
-
- dpFrictionCoefficient = 6.0 * sinphi * oneOverDenominator;
- dpCohesion = 6.0 * effectiveCohesion * cos(phi) * oneOverDenominator;
-
- frictionalStrength = dpFrictionCoefficient * pressure + dpCohesion ;
-
- particleExt->tensileFailure = (frictionalStrength <= 0.0);
-
- if(minimumYieldStress==0)
- minimumYieldStress=dpCohesion;
-
- if ( frictionalStrength < minimumYieldStress)
- frictionalStrength = minimumYieldStress;
-
- return frictionalStrength;
-}
-
-void _DruckerPrager_UpdateDrawParameters( void* rheology ) {
- DruckerPrager* self = (DruckerPrager*) rheology;
- Particle_Index lParticle_I;
- Particle_Index particleLocalCount;
- StrainWeakening* strainWeakening = self->strainWeakening;
- MaterialPoint* materialPoint;
-
- double length;
- double brightness;
- double opacity;
- double strainWeakeningRatio;
- double localMaxStrainIncrement;
- double localMeanStrainIncrement;
- Particle_Index localFailed;
-
- double globalMaxStrainIncrement;
- double globalMeanStrainIncrement;
- Particle_Index globalFailed;
-
- double averagedGlobalMaxStrainIncrement = 0.0;
-
- double oneOverGlobalMaxStrainIncrement;
- double postFailureWeakeningIncrement;
-
- /* Note : this function defines some drawing parameters (brightness, opacity, diameter) as
- * functions of the strain weakening - this needs to be improved since most of the parameters
- * that define this dependency are hard coded here. We need to have a more flexible way
- * to construct the viz parameters as functions of material parameters */
-
- /* We should only update the drawing parameters if the strain weakening is defined */
- if (strainWeakening==NULL)
- return;
-
- localMaxStrainIncrement = 0.0;
- localMeanStrainIncrement = 0.0;
- localFailed = 0;
-
- /* Update all variables */
- Variable_Update( self->hasYieldedVariable->variable );
- Variable_Update( self->brightness->variable );
- Variable_Update( self->opacity->variable );
- Variable_Update( self->diameter->variable );
- Variable_Update( strainWeakening->postFailureWeakeningIncrement->variable );
-
- particleLocalCount = self->hasYieldedVariable->variable->arraySize;
-
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
- if ( Variable_GetValueChar( self->hasYieldedVariable->variable, lParticle_I )) {
- localFailed++;
-
- postFailureWeakeningIncrement =
- Variable_GetValueDouble( strainWeakening->postFailureWeakeningIncrement->variable, lParticle_I );
-
- localMeanStrainIncrement += postFailureWeakeningIncrement;
-
- if(localMaxStrainIncrement < postFailureWeakeningIncrement)
- localMaxStrainIncrement = postFailureWeakeningIncrement;
- }
- }
-
- MPI_Allreduce( &localMaxStrainIncrement, &globalMaxStrainIncrement, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
- MPI_Allreduce( &localMeanStrainIncrement, &globalMeanStrainIncrement, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
- MPI_Allreduce( &localFailed, &globalFailed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
-
- if(globalFailed == 0)
- return;
-
- globalMeanStrainIncrement /= (double) globalFailed;
-
- averagedGlobalMaxStrainIncrement =
- 0.5 * averagedGlobalMaxStrainIncrement +
- 0.25 * globalMeanStrainIncrement +
- 0.25 * globalMaxStrainIncrement;
-
- /* Let's simply assume that twice the mean is a good place to truncate these values */
- oneOverGlobalMaxStrainIncrement = 1.0 / averagedGlobalMaxStrainIncrement;
-
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
- materialPoint = (MaterialPoint*)Swarm_ParticleAt( strainWeakening->swarm, lParticle_I );
-
- if ( Variable_GetValueChar( self->hasYieldedVariable->variable, lParticle_I ) == False ||
- StrainWeakening_GetPostFailureWeakening( strainWeakening, materialPoint ) < 0.0 )
- {
- Variable_SetValueFloat( self->brightness->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->diameter->variable, lParticle_I, 0.0 );
- continue;
- }
-
- postFailureWeakeningIncrement =
- Variable_GetValueDouble( strainWeakening->postFailureWeakeningIncrement->variable, lParticle_I );
-
- strainWeakeningRatio = StrainWeakening_CalcRatio( strainWeakening, materialPoint );
-
- length = 0.001 + 0.003 * strainWeakeningRatio;
- brightness = strainWeakeningRatio * postFailureWeakeningIncrement * oneOverGlobalMaxStrainIncrement;
-
- opacity = 0.5 * brightness;
-
- if( brightness > 1.0 )
- brightness = 1.0;
-
- if( opacity > 0.5 )
- opacity = 0.5;
-
- if( opacity < 0.1 )
- opacity = 0.0;
-
- Variable_SetValueFloat( self->brightness->variable, lParticle_I, brightness );
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, opacity );
- Variable_SetValueFloat( self->diameter->variable, lParticle_I, (float) length );
- }
-
-}
Copied: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c (from rev 14446, long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c)
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c (rev 0)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.c 2009-03-25 19:32:43 UTC (rev 14450)
@@ -0,0 +1,509 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** Copyright (c) 2005, Monash Cluster Computing
+** All rights reserved.
+** Redistribution and use in source and binary forms, with or without modification,
+** are permitted provided that the following conditions are met:
+**
+** * Redistributions of source code must retain the above copyright notice,
+** this list of conditions and the following disclaimer.
+** * Redistributions in binary form must reproduce the above copyright
+** notice, this list of conditions and the following disclaimer in the
+** documentation and/or other materials provided with the distribution.
+** * Neither the name of the Monash University nor the names of its contributors
+** may be used to endorse or promote products derived from this software
+** without specific prior written permission.
+**
+** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+**
+**
+** Contact:
+*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
+*%
+** Contributors:
+*+ Robert Turnbull
+*+ Vincent Lemiale
+*+ Louis Moresi
+*+ David May
+*+ David Stegman
+*+ Mirko Velic
+*+ Patrick Sunter
+*+ Julian Giordani
+*+
+** $Id$
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
+
+#include "types.h"
+#include "RheologyClass.h"
+#include "StrainWeakening.h"
+#include "YieldRheology.h"
+#include "DruckerPrager.h"
+#include "ConstitutiveMatrix.h"
+
+#include <assert.h>
+#include <string.h>
+
+/* Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
+const Type DruckerPrager_Type = "DruckerPrager";
+
+/* Private Constructor: This will accept all the virtual functions for this class as arguments. */
+DruckerPrager* _DruckerPrager_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
+ YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
+ YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
+ YieldRheology_HasYieldedFunction* _hasYielded,
+ Name name )
+{
+ DruckerPrager* self;
+
+ /* Call private constructor of parent - this will set virtual functions of parent and continue up the hierarchy tree. At the beginning of the tree it will allocate memory of the size of object and initialise all the memory to zero. */
+ assert( sizeOfSelf >= sizeof(DruckerPrager) );
+ self = (DruckerPrager*) _YieldRheology_New(
+ sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _modifyConstitutiveMatrix,
+ _getYieldCriterion,
+ _getYieldIndicator,
+ _hasYielded,
+ name );
+
+ /* Function pointers for this class that are not on the parent class should be set here */
+
+ return self;
+}
+
+void _DruckerPrager_Init(
+ DruckerPrager* self,
+ FeVariable* pressureField,
+ FeVariable* velocityGradientsField,
+ MaterialPointsSwarm* materialPointsSwarm,
+ double cohesion,
+ double cohesionAfterSoftening,
+ double frictionCoefficient,
+ double frictionCoefficientAfterSoftening,
+ double boundaryCohesion,
+ double boundaryCohesionAfterSoftening,
+ double boundaryFrictionCoefficient,
+ double boundaryFrictionCoefficientAfterSoftening,
+ Bool boundaryBottom,
+ Bool boundaryTop,
+ Bool boundaryLeft,
+ Bool boundaryRight,
+ Bool boundaryFront,
+ Bool boundaryBack,
+ double minimumYieldStress,
+ HydrostaticTerm* hydrostaticTerm)
+{
+ self->materialPointsSwarm = materialPointsSwarm;
+ self->pressureField = pressureField;
+ self->velocityGradientsField = velocityGradientsField;
+
+ self->cohesion = cohesion;
+ self->frictionCoefficient = frictionCoefficient;
+
+ /* 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;
+ self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
+ self->boundaryCohesion = boundaryCohesion;
+ self->boundaryFrictionCoefficient = boundaryFrictionCoefficient;
+ self->boundaryCohesionAfterSoftening = boundaryCohesionAfterSoftening;
+ self->boundaryFrictionCoefficientAfterSoftening = boundaryFrictionCoefficientAfterSoftening;
+ self->boundaryBottom=boundaryBottom;
+ self->boundaryTop=boundaryTop;
+ self->boundaryLeft=boundaryLeft;
+ self->boundaryRight=boundaryRight;
+ self->boundaryFront=boundaryFront;
+ self->boundaryBack=boundaryBack;
+
+ self->minimumYieldStress = minimumYieldStress;
+ self->hydrostaticTerm=hydrostaticTerm;
+}
+
+void* _DruckerPrager_DefaultNew( Name name ) {
+ return (void*) _DruckerPrager_New(
+ sizeof(DruckerPrager),
+ DruckerPrager_Type,
+ _YieldRheology_Delete,
+ _YieldRheology_Print,
+ _YieldRheology_Copy,
+ _DruckerPrager_DefaultNew,
+ _DruckerPrager_Construct,
+ _DruckerPrager_Build,
+ _DruckerPrager_Initialise,
+ _YieldRheology_Execute,
+ _YieldRheology_Destroy,
+ _DruckerPrager_ModifyConstitutiveMatrix,
+ _DruckerPrager_GetYieldCriterion,
+ _DruckerPrager_GetYieldIndicator,
+ _DruckerPrager_HasYielded,
+ name );
+}
+
+void _DruckerPrager_Construct( void* rheology, Stg_ComponentFactory* cf,
+ void *data ){
+ DruckerPrager* self = (DruckerPrager*)rheology;
+ FeVariable* pressureField;
+ MaterialPointsSwarm* materialPointsSwarm;
+ FeVariable* velocityGradientsField;
+
+ /* Construct Parent */
+ _YieldRheology_Construct( self, cf, data );
+
+ /* Make sure that there is strain weakening */
+ Journal_Firewall(
+ self->strainWeakening != NULL,
+ Journal_Register( Error_Type, self->type ),
+ "Error in func '%s' for %s '%s': DruckerPrager rheology needs strain weakening.\n",
+ __func__, self->type, self->name );
+
+ materialPointsSwarm = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "MaterialPointsSwarm", MaterialPointsSwarm, True, data );
+ pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "PressureField", FeVariable, True, data );
+ velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "VelocityGradientsField", FeVariable, True, data );
+
+ _DruckerPrager_Init(
+ self,
+ pressureField,
+ velocityGradientsField,
+ materialPointsSwarm,
+ Stg_ComponentFactory_GetDouble( cf, self->name, "cohesion", 0.0 ),
+ 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_GetBool( cf, self->name, "boundaryBottom", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryTop", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryLeft", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryRight", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryFront", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryBack", False ),
+ Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ),
+ Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "HydrostaticTerm", HydrostaticTerm, False, data ) );
+}
+
+void _DruckerPrager_Build( void* rheology, void* data ) {
+ DruckerPrager* self = (DruckerPrager*) rheology;
+
+ /* Build parent */
+ _YieldRheology_Build( self, data );
+}
+
+void _DruckerPrager_Initialise( void* rheology, void* data ) {
+ DruckerPrager* self = (DruckerPrager*) rheology;
+ Particle_Index lParticle_I;
+ Particle_Index particleLocalCount;
+ AbstractContext* context = (AbstractContext*)data;
+
+ _YieldRheology_Initialise( self, data );
+
+ /* Initialise these components now just in case this function is called before their own _Initialise function */
+ Stg_Component_Initialise( self->materialPointsSwarm, data, False );
+ Stg_Component_Initialise( self->strainWeakening, data, False );
+
+ /* We don't need to Initialise hasYieldedVariable because it's a parent variable and _YieldRheology_Initialise
+ * has already been called */
+ particleLocalCount = self->hasYieldedVariable->variable->arraySize;
+
+ /* If restarting from checkpoint, don't change the parameters on the particles */
+ if ( !(context && (True == context->loadFromCheckPoint) ) ) {
+ for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
+ Variable_SetValueChar( self->hasYieldedVariable->variable, lParticle_I, False );
+ }
+ }
+}
+
+void _DruckerPrager_ModifyConstitutiveMatrix(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ DruckerPrager* self = (DruckerPrager*) rheology;
+
+ /* Don't want to yield on the first ever solve */
+ if ( constitutiveMatrix->previousSolutionExists == False ) {
+ return;
+ }
+
+ /* Calculate and store values of stress, pressure, velocity
+ gradients, inavriants so they are only calculated once */
+ _DruckerPrager_StoreCurrentParameters( self, constitutiveMatrix,
+ materialPointsSwarm, lElement_I,
+ materialPoint, xi );
+
+ _YieldRheology_ModifyConstitutiveMatrix(self, constitutiveMatrix,
+ materialPointsSwarm,lElement_I,
+ materialPoint, xi );
+}
+
+double _DruckerPrager_GetYieldIndicator(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ DruckerPrager* self = (DruckerPrager*) rheology;
+
+ return SymmetricTensor_2ndInvariant( self->stress,
+ constitutiveMatrix->dim );
+}
+
+double _DruckerPrager_GetYieldCriterion(void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi)
+{
+ DruckerPrager* self = (DruckerPrager*) rheology;
+ double minimumYieldStress;
+ double effectiveCohesion;
+ double effectiveFrictionCoefficient;
+ double frictionalStrength;
+ Dimension_Index dim = constitutiveMatrix->dim;
+
+ FeVariable* pressureField = self->pressureField;
+ double pressure;
+ double factor;
+ Cell_Index cell_I;
+ Coord coord;
+ Element_GlobalIndex element_gI = 0;
+ unsigned inds[3];
+ Grid* elGrid;
+ Bool inside_boundary;
+ 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)
+ {
+ pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,coord);
+ }
+
+ /* Normally we add the average of the trace of the stress. With
+ compressible material, you have to do it. But with stabilized
+ linear pressure elements, the non-zero trace is a numerical
+ artifact. So we do not add it. */
+
+/* pressure+=self->trace/dim; */
+
+ /* 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. */
+
+
+ /* Big song and dance to see if we are at a boundary that we care about */
+ elGrid = *(Grid**)ExtensionManager_Get(pressureField->feMesh->info,
+ pressureField->feMesh,
+ ExtensionManager_GetHandle(pressureField->feMesh->info, "elementGrid" ) );
+
+ element_gI = FeMesh_ElementDomainToGlobal( pressureField->feMesh, lElement_I );
+ RegularMeshUtils_Element_1DTo3D( pressureField->feMesh, element_gI, inds );
+
+ inside_boundary=(self->boundaryBottom && inds[1]==0)
+ || (self->boundaryTop && inds[1]==elGrid->sizes[1]-1)
+ || (self->boundaryLeft && inds[0]==0)
+ || (self->boundaryRight && inds[0]==elGrid->sizes[0]-1)
+ || (dim==3 && ((self->boundaryBack && inds[2]==0)
+ || (self->boundaryFront && inds[2]==elGrid->sizes[2]-1)));
+
+ effectiveFrictionCoefficient =
+ _DruckerPrager_EffectiveFrictionCoefficient( self, materialPoint,
+ inside_boundary );
+ effectiveCohesion = _DruckerPrager_EffectiveCohesion(self,materialPoint,
+ inside_boundary);
+
+
+ if(dim==2)
+ {
+ /* effectiveFrictionCoefficient=tan(phi). If
+ factor=sin(atan(1/tan(phi))) =>
+ factor=cos(phi)=1/sqrt(1+tan(phi)**2) */
+ factor=1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
+ frictionalStrength = effectiveFrictionCoefficient*pressure*factor
+ + effectiveCohesion*factor;
+ }
+ else
+ {
+ double cos_phi, sin_phi;
+ /* cos(phi)=1/sqrt(1+tan(phi)**2) */
+ cos_phi=
+ 1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
+ sin_phi=effectiveFrictionCoefficient*cos_phi;
+ factor=2*cos_phi/(sqrt(3.0)*(3-sin_phi));
+
+ /* The full expression is
+
+ sqrt(J2)=p*2*sin(phi)/(sqrt(3)*(3-sin(phi)))
+ + C*6*cos(phi)/(sqrt(3)*(3-sin(phi)))
+
+ Note the extra factor of 3 for cohesion */
+
+ frictionalStrength = effectiveFrictionCoefficient*factor*pressure
+ + effectiveCohesion*3*factor;
+ }
+
+ /* If the minimumYieldStress is not set, then use the
+ effective cohesion. Maybe it should be the modified
+ effective cohesion, though that probably should not matter
+ much. */
+ minimumYieldStress = self->minimumYieldStress;
+ if(minimumYieldStress==0.0)
+ minimumYieldStress=effectiveCohesion;
+
+ /* Make sure frictionalStrength is above the minimum */
+ if ( frictionalStrength < minimumYieldStress*factor)
+ frictionalStrength = minimumYieldStress*factor;
+
+ return frictionalStrength;
+}
+
+/* This scales the viscosity so that the stress is the yield stress.
+ We can simply scale the viscosity by yieldCriterion/yieldIndicator
+ because yieldCriterion does not depend on the deviatoric stress
+ (and thus the viscosity), while the yieldIndicator depends linearly
+ on the deviatoric stresses. */
+
+void _DruckerPrager_HasYielded(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ double yieldCriterion,
+ double yieldIndicator )
+{
+ DruckerPrager* self = (DruckerPrager*)rheology;
+ double viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
+
+ ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
+}
+
+double _DruckerPrager_EffectiveCohesion( void* rheology, void* materialPoint,
+ Bool inside_boundary ) {
+ DruckerPrager* self = (DruckerPrager*) rheology;
+ double effectiveCohesion;
+
+ double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
+
+ if(inside_boundary)
+ {
+ effectiveCohesion=self->boundaryCohesion*(1.0-strainWeakeningRatio)+
+ self->boundaryCohesionAfterSoftening*strainWeakeningRatio;
+ }
+ else
+ {
+ effectiveCohesion = self->cohesion * (1.0 - strainWeakeningRatio) +
+ self->cohesionAfterSoftening * strainWeakeningRatio;
+ }
+
+ return effectiveCohesion;
+}
+
+double _DruckerPrager_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, Bool inside_boundary ) {
+ DruckerPrager* self = (DruckerPrager*) rheology;
+ double effectiveFrictionCoefficient = 0.0;
+
+ double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
+
+ if(inside_boundary)
+ {
+ effectiveFrictionCoefficient =
+ self->boundaryFrictionCoefficient * (1.0 - strainWeakeningRatio) +
+ self->boundaryFrictionCoefficientAfterSoftening
+ *strainWeakeningRatio;
+ }
+ else
+ {
+ effectiveFrictionCoefficient =
+ self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
+ self->frictionCoefficientAfterSoftening * strainWeakeningRatio;
+ }
+
+ return effectiveFrictionCoefficient;
+}
+
+void _DruckerPrager_StoreCurrentParameters(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ DruckerPrager* self = (DruckerPrager*) rheology;
+ SymmetricTensor strainRate;
+ Dimension_Index dim = constitutiveMatrix->dim;
+ double strainRateTrace;
+ int i;
+
+ FeVariable_InterpolateWithinElement( self->pressureField, lElement_I, xi,
+ &self->currentPressure );
+ FeVariable_InterpolateWithinElement(self->velocityGradientsField,lElement_I,
+ xi, self->currentVelocityGradient );
+ TensorArray_GetSymmetricPart(self->currentVelocityGradient,dim,strainRate);
+
+ ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate,
+ self->stress );
+ SymmetricTensor_GetTrace(self->stress, dim, &(self->trace));
+ SymmetricTensor_GetTrace(strainRate, dim, &strainRateTrace);
+
+ /* Subtract the trace so that the 2nd invariant is not polluted by
+ the trace */
+
+ for(i=0;i<dim;++i)
+ {
+ strainRate[TensorMapST3D[i][i]]-=strainRateTrace/dim;
+ self->stress[TensorMapST3D[i][i]]-=self->trace/dim;
+ }
+
+ self->strainRateSecondInvariant=SymmetricTensor_2ndInvariant(strainRate,dim);
+}
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,115 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-** Copyright (c) 2005, Monash Cluster Computing
-** All rights reserved.
-** Redistribution and use in source and binary forms, with or without modification,
-** are permitted provided that the following conditions are met:
-**
-** * Redistributions of source code must retain the above copyright notice,
-** this list of conditions and the following disclaimer.
-** * Redistributions in binary form must reproduce the above copyright
-** notice, this list of conditions and the following disclaimer in the
-** documentation and/or other materials provided with the distribution.
-** * Neither the name of the Monash University nor the names of its contributors
-** may be used to endorse or promote products derived from this software
-** without specific prior written permission.
-**
-** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
-** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
-** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-**
-**
-** Contact:
-*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
-*%
-** Contributors:
-*+ Robert Turnbull
-*+ Vincent Lemiale
-*+ Louis Moresi
-*+ David May
-*+ David Stegman
-*+ Mirko Velic
-*+ Patrick Sunter
-*+ Julian Giordani
-*+
-** $Id$
-**
-**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
-#ifndef __Underworld_Rheology_DruckerPrager_h__
-#define __Underworld_Rheology_DruckerPrager_h__
-
- /** Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
- extern const Type DruckerPrager_Type;
-
- typedef struct {
- float brightness;
- float opacity;
- float diameter;
- Particle_Bool tensileFailure;
- } DruckerPrager_Particle;
-
- /** Rheology class contents - this is defined as a macro so that sub-classes of this class can use this macro at the start of the definition of their struct */
- #define __DruckerPrager \
- /* Parent info */ \
- __VonMises \
- /* Virtual functions go here */ \
- /* General Info */\
- SwarmVariable* brightness; \
- SwarmVariable* opacity; \
- SwarmVariable* diameter; \
- SwarmVariable* tensileFailure; \
- ExtensionInfo_Index particleExtHandle; \
- /* Param passed in */\
- FeVariable* pressureField; \
- double frictionCoefficient; \
- double frictionCoefficientAfterSoftening;
-
-
- struct DruckerPrager { __DruckerPrager };
-
- /** Private Constructor: This will accept all the virtual functions for this class as arguments. */
- DruckerPrager* _DruckerPrager_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
- YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
- YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
- YieldRheology_HasYieldedFunction* _hasYielded,
- Name name ) ;
-
- /* 'Stg_Component' implementations */
- void* _DruckerPrager_DefaultNew( Name name ) ;
- void _DruckerPrager_Construct( void* rheology, Stg_ComponentFactory* cf, void* data );
-
- void _DruckerPrager_Build( void* rheology, void* data );
- void _DruckerPrager_Initialise( void* rheology, void* data ) ;
-
- /* 'YieldRheology' implementations */
- double _DruckerPrager_GetYieldCriterion(
- void* druckerPrager,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi );
-
- void _DruckerPrager_UpdateDrawParameters( void* rheology ) ;
-
-#endif
Copied: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h (from rev 14446, long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h)
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h (rev 0)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.h 2009-03-25 19:32:43 UTC (rev 14450)
@@ -0,0 +1,165 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** Copyright (c) 2005, Monash Cluster Computing
+** All rights reserved.
+** Redistribution and use in source and binary forms, with or without modification,
+** are permitted provided that the following conditions are met:
+**
+** * Redistributions of source code must retain the above copyright notice,
+** this list of conditions and the following disclaimer.
+** * Redistributions in binary form must reproduce the above copyright
+** notice, this list of conditions and the following disclaimer in the
+** documentation and/or other materials provided with the distribution.
+** * Neither the name of the Monash University nor the names of its contributors
+** may be used to endorse or promote products derived from this software
+** without specific prior written permission.
+**
+** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+**
+**
+** Contact:
+*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
+*%
+** Contributors:
+*+ Robert Turnbull
+*+ Vincent Lemiale
+*+ Louis Moresi
+*+ David May
+*+ David Stegman
+*+ Mirko Velic
+*+ Patrick Sunter
+*+ Julian Giordani
+*+
+** $Id$
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#ifndef __Underworld_Rheology_DruckerPrager_h__
+#define __Underworld_Rheology_DruckerPrager_h__
+
+ /** Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
+ extern const Type DruckerPrager_Type;
+
+ /** Rheology class contents - this is defined as a macro so that sub-classes of this class can use this macro at the start of the definition of their struct */
+ #define __DruckerPrager \
+ /* Parent info */ \
+ __YieldRheology \
+ /* Virtual functions go here */ \
+ /* Material Parameters */\
+ ExtensionInfo_Index particleExtHandle; \
+ /* Param passed in */ \
+ MaterialPointsSwarm* materialPointsSwarm; \
+ FeVariable* pressureField; \
+ FeVariable* velocityGradientsField; \
+ /* Director component is used to update the normal */\
+ Director* director; \
+ double cohesion; \
+ double cohesionAfterSoftening; \
+ double frictionCoefficient; \
+ double frictionCoefficientAfterSoftening; \
+ double boundaryCohesion; \
+ double boundaryCohesionAfterSoftening; \
+ double boundaryFrictionCoefficient; \
+ double boundaryFrictionCoefficientAfterSoftening; \
+ Bool boundaryTop; \
+ Bool boundaryBottom; \
+ Bool boundaryLeft; \
+ Bool boundaryRight; \
+ Bool boundaryFront; \
+ Bool boundaryBack; \
+ double minimumYieldStress; \
+ /* Stored values that are calculated once for each particle */ \
+ double trace; \
+ TensorArray currentVelocityGradient; \
+ SymmetricTensor stress; \
+ double currentPressure; \
+ double strainRateSecondInvariant; \
+ FeVariable* strainRateField; \
+ HydrostaticTerm* hydrostaticTerm;
+
+
+ struct DruckerPrager { __DruckerPrager };
+
+ /** Private Constructor: This will accept all the virtual functions for this class as arguments. */
+ DruckerPrager* _DruckerPrager_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
+ YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
+ YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
+ YieldRheology_HasYieldedFunction* _hasYielded,
+ Name name ) ;
+
+ /* 'Stg_Component' implementations */
+ void* _DruckerPrager_DefaultNew( Name name ) ;
+ void _DruckerPrager_Construct( void* rheology, Stg_ComponentFactory* cf,
+ void *data );
+
+ void _DruckerPrager_Build( void* rheology, void* data );
+ void _DruckerPrager_Initialise( void* rheology, void* data ) ;
+
+ /* 'YieldRheology' implementations */
+ void _DruckerPrager_ModifyConstitutiveMatrix(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ double _DruckerPrager_GetYieldCriterion(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ double _DruckerPrager_GetYieldIndicator(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ void _DruckerPrager_HasYielded(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ double yieldCriterion,
+ double yieldIndicator );
+
+ double _DruckerPrager_EffectiveCohesion( void* rheology, void* materialPoint, Bool inside_boundary ) ;
+ double _DruckerPrager_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, Bool inside_boundary ) ;
+ double _DruckerPrager_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) ;
+
+ void _DruckerPrager_StoreCurrentParameters(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi ) ;
+
+#endif
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,65 +0,0 @@
-<?xml version="1.0"?>
-<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
-<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
-
-<param name="Name">DruckerPrager</param>
-<param name="Author">...</param>
-<param name="Organisation">MCC</param>
-<param name="Project">Underworld</param>
-<param name="Location">./Underworld/Rheology/src/</param>
-<param name="Project Web">http://mcc.monash.edu.au/Underworld</param>
-<param name="Copyright">Copyright (c) 2005, Monash Cluster Computing</param>
-<param name="License">http://www.opensource.org/licenses/bsd-license.php</param>
-<param name="Parent">VonMises</param>
-<param name="Reference">...</param>
-<param name="Summary">...</param>
-<param name="Description">Implements a Drucker-Prager yield Rheology. The parameters of this criterion are calibrated so that the yield surface coincides with the Mohr Coulomb one for the condition of a compression test. </param>
-<param name="Equation">Yield Surface: $f = \sqrt{3} \bar{\tau} - \alpha p - k \le 0, \alpha = \frac{6 \sin\varphi}{3 - \sin\varphi}, k = \frac{6 c \cos\varphi}{3 - \sin\varphi}$</param>
-
-<!--Now the interesting stuff-->
-
-
-<list name="Params">
- <struct>
- <param name="Name">frictionCoefficient</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">Defines the dependance of the criterion to the pressure.</param>
- </struct>
- <struct>
- <param name="Name">frictionCoefficientAfterSoftening</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">This is the value of the frictionCoefficient after the material has weakened</param>
- </struct>
-
-</list>
-
-<list name="Dependencies">
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">FeVariable</param>
- <param name="Fallback Key">PressureField</param>
- <param name="Description">Essential component since the DruckerPrager criterion is pressure dependant</param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">IntegrationPointsSwarm</param>
- <param name="Fallback Key">MaterialPoints</param>
- <param name="Description">Defines the set of material points. It is needed because infos stored by each material point are added for viz purpose. </param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">FiniteElementContext</param>
- <param name="Fallback Key">Context</param>
- <param name="Description">The context is essential because the DruckerPrager_UpdateDrawParameters function is prepended to an existing Entry Point. </param>
- </struct>
-
-</list>
-<!-- Add an example XML if possible -->
-<param name="Example">Underworld/InputFiles/ExtensionDP.xml</param>
-
-</StGermainData>
Copied: long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta (from rev 11823, long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta)
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta (rev 0)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DruckerPrager.meta 2009-03-25 19:32:43 UTC (rev 14450)
@@ -0,0 +1,90 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">DruckerPrager</param>
+<param name="Organisation">MCC</param>
+<param name="Project">Underworld</param>
+<param name="Location">./Underworld/Rheology/src/</param>
+<param name="Project Web">http://mcc.monash.edu.au/Underworld</param>
+<param name="Copyright">Copyright (c) 2005, Monash Cluster Computing</param>
+<param name="License">http://www.opensource.org/licenses/bsd-license.php</param>
+<param name="Parent">YieldRheology</param>
+<param name="Description">Implements faulting behaviour based on Mohr-Coulomb failure criterion. See following reference for details : Anisotropic viscous models of large-deformation Mohr–Coulomb failure, Moresi, L.; Mühlhaus, H.-B., Philosophical Magazine, Volume 86, Numbers 21-22, -22/21 July–1 August 2006, pp. 3287-3305</param>
+
+<!--Now the interesting stuff-->
+
+<list name="Params">
+ <struct>
+ <param name="Name">cohesion</param>
+ <param name="Type">Double</param>
+ <param name="Default">0.0</param>
+ <param name="Description">Cohesion for a pristine material</param>
+ </struct>
+ <struct>
+ <param name="Name">cohesionAfterSoftening</param>
+ <param name="Type">Double</param>
+ <param name="Default">0.0</param>
+ <param name="Description">Value of the cohesion after the material has weakened</param>
+ </struct>
+ <struct>
+ <param name="Name">frictionCoefficient</param>
+ <param name="Type">Double</param>
+ <param name="Default">0.0</param>
+ <param name="Description">Defines the dependance of the criterion to the normal stress</param>
+ </struct>
+ <struct>
+ <param name="Name">frictionCoefficientAfterSoftening</param>
+ <param name="Type">Double</param>
+ <param name="Default">0.0</param>
+ <param name="Description">This is the value of the frictionCoefficient after the material has weakened</param>
+ </struct>
+ <struct>
+ <param name="Name">minimumYieldStress</param>
+ <param name="Type">Double</param>
+ <param name="Default">0.0</param>
+ <param name="Description">The value of the failure criterion (defining the yield surface) can not be smaller than this parameter</param>
+ </struct>
+</list>
+
+<list name="Dependencies">
+
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Type">FeVariable</param>
+ <param name="Fallback Key">PressureField</param>
+ <param name="Description">Essential since the DruckerPrager failure criterion is pressure dependant</param>
+ </struct>
+
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Type">IntegrationPointsSwarm</param>
+ <param name="Description">Defines the set of material points</param>
+ </struct>
+
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Type">FeVariable</param>
+ <param name="Fallback Key">VelocityGradientsFieldField</param>
+ <param name="Description">The resolved component of this tensor wrt the failure plane are mainly used to determine the most favourable oriented plane for failure to occur </param>
+ </struct>
+
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Type">StrainWeakening</param>
+ <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 -->
+<param name="Example">Underworld/InputFiles/ExtensionFMM.xml</param>
+<param name="Example">Underworld/InputFiles/ExtensionFMM3D.xml</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c 2009-03-25 19:32:43 UTC (rev 14450)
@@ -73,7 +73,6 @@
Stg_ComponentRegister_Add( componentRegister, Byerlee_Type, "0", _Byerlee_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, DruckerPrager_Type, "0", _DruckerPrager_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, FaultingMoresiMuhlhaus2006_Type, "0", _FaultingMoresiMuhlhaus2006_DefaultNew );
- Stg_ComponentRegister_Add( componentRegister, MohrCoulomb_Type, "0", _MohrCoulomb_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, StrainWeakening_Type, "0", _StrainWeakening_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, BuiterStrainWeakening_Type, "0", _BuiterStrainWeakening_DefaultNew );
@@ -103,7 +102,6 @@
RegisterParent( Byerlee_Type, VonMises_Type );
RegisterParent( DruckerPrager_Type, VonMises_Type );
RegisterParent( FaultingMoresiMuhlhaus2006_Type, YieldRheology_Type );
- RegisterParent( MohrCoulomb_Type, YieldRheology_Type );
RegisterParent( StrainWeakening_Type, TimeIntegratee_Type );
RegisterParent( BuiterStrainWeakening_Type, StrainWeakening_Type );
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,509 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-** Copyright (c) 2005, Monash Cluster Computing
-** All rights reserved.
-** Redistribution and use in source and binary forms, with or without modification,
-** are permitted provided that the following conditions are met:
-**
-** * Redistributions of source code must retain the above copyright notice,
-** this list of conditions and the following disclaimer.
-** * Redistributions in binary form must reproduce the above copyright
-** notice, this list of conditions and the following disclaimer in the
-** documentation and/or other materials provided with the distribution.
-** * Neither the name of the Monash University nor the names of its contributors
-** may be used to endorse or promote products derived from this software
-** without specific prior written permission.
-**
-** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
-** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
-** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-**
-**
-** Contact:
-*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
-*%
-** Contributors:
-*+ Robert Turnbull
-*+ Vincent Lemiale
-*+ Louis Moresi
-*+ David May
-*+ David Stegman
-*+ Mirko Velic
-*+ Patrick Sunter
-*+ Julian Giordani
-*+
-** $Id$
-**
-**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
-#include <mpi.h>
-#include <StGermain/StGermain.h>
-#include <StgFEM/StgFEM.h>
-#include <PICellerator/PICellerator.h>
-
-#include "types.h"
-#include "RheologyClass.h"
-#include "StrainWeakening.h"
-#include "YieldRheology.h"
-#include "MohrCoulomb.h"
-#include "ConstitutiveMatrix.h"
-
-#include <assert.h>
-#include <string.h>
-
-/* Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
-const Type MohrCoulomb_Type = "MohrCoulomb";
-
-/* Private Constructor: This will accept all the virtual functions for this class as arguments. */
-MohrCoulomb* _MohrCoulomb_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
- YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
- YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
- YieldRheology_HasYieldedFunction* _hasYielded,
- Name name )
-{
- MohrCoulomb* self;
-
- /* Call private constructor of parent - this will set virtual functions of parent and continue up the hierarchy tree. At the beginning of the tree it will allocate memory of the size of object and initialise all the memory to zero. */
- assert( sizeOfSelf >= sizeof(MohrCoulomb) );
- self = (MohrCoulomb*) _YieldRheology_New(
- sizeOfSelf,
- type,
- _delete,
- _print,
- _copy,
- _defaultConstructor,
- _construct,
- _build,
- _initialise,
- _execute,
- _destroy,
- _modifyConstitutiveMatrix,
- _getYieldCriterion,
- _getYieldIndicator,
- _hasYielded,
- name );
-
- /* Function pointers for this class that are not on the parent class should be set here */
-
- return self;
-}
-
-void _MohrCoulomb_Init(
- MohrCoulomb* self,
- FeVariable* pressureField,
- FeVariable* velocityGradientsField,
- MaterialPointsSwarm* materialPointsSwarm,
- double cohesion,
- double cohesionAfterSoftening,
- double frictionCoefficient,
- double frictionCoefficientAfterSoftening,
- double boundaryCohesion,
- double boundaryCohesionAfterSoftening,
- double boundaryFrictionCoefficient,
- double boundaryFrictionCoefficientAfterSoftening,
- Bool boundaryBottom,
- Bool boundaryTop,
- Bool boundaryLeft,
- Bool boundaryRight,
- Bool boundaryFront,
- Bool boundaryBack,
- double minimumYieldStress,
- HydrostaticTerm* hydrostaticTerm)
-{
- self->materialPointsSwarm = materialPointsSwarm;
- self->pressureField = pressureField;
- self->velocityGradientsField = velocityGradientsField;
-
- self->cohesion = cohesion;
- self->frictionCoefficient = frictionCoefficient;
-
- /* 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;
- self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
- self->boundaryCohesion = boundaryCohesion;
- self->boundaryFrictionCoefficient = boundaryFrictionCoefficient;
- self->boundaryCohesionAfterSoftening = boundaryCohesionAfterSoftening;
- self->boundaryFrictionCoefficientAfterSoftening = boundaryFrictionCoefficientAfterSoftening;
- self->boundaryBottom=boundaryBottom;
- self->boundaryTop=boundaryTop;
- self->boundaryLeft=boundaryLeft;
- self->boundaryRight=boundaryRight;
- self->boundaryFront=boundaryFront;
- self->boundaryBack=boundaryBack;
-
- self->minimumYieldStress = minimumYieldStress;
- self->hydrostaticTerm=hydrostaticTerm;
-}
-
-void* _MohrCoulomb_DefaultNew( Name name ) {
- return (void*) _MohrCoulomb_New(
- sizeof(MohrCoulomb),
- MohrCoulomb_Type,
- _YieldRheology_Delete,
- _YieldRheology_Print,
- _YieldRheology_Copy,
- _MohrCoulomb_DefaultNew,
- _MohrCoulomb_Construct,
- _MohrCoulomb_Build,
- _MohrCoulomb_Initialise,
- _YieldRheology_Execute,
- _YieldRheology_Destroy,
- _MohrCoulomb_ModifyConstitutiveMatrix,
- _MohrCoulomb_GetYieldCriterion,
- _MohrCoulomb_GetYieldIndicator,
- _MohrCoulomb_HasYielded,
- name );
-}
-
-void _MohrCoulomb_Construct( void* rheology, Stg_ComponentFactory* cf,
- void *data ){
- MohrCoulomb* self = (MohrCoulomb*)rheology;
- FeVariable* pressureField;
- MaterialPointsSwarm* materialPointsSwarm;
- FeVariable* velocityGradientsField;
-
- /* Construct Parent */
- _YieldRheology_Construct( self, cf, data );
-
- /* Make sure that there is strain weakening */
- Journal_Firewall(
- self->strainWeakening != NULL,
- Journal_Register( Error_Type, self->type ),
- "Error in func '%s' for %s '%s': MohrCoulomb rheology needs strain weakening.\n",
- __func__, self->type, self->name );
-
- materialPointsSwarm = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "MaterialPointsSwarm", MaterialPointsSwarm, True, data );
- pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "PressureField", FeVariable, True, data );
- velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "VelocityGradientsField", FeVariable, True, data );
-
- _MohrCoulomb_Init(
- self,
- pressureField,
- velocityGradientsField,
- materialPointsSwarm,
- Stg_ComponentFactory_GetDouble( cf, self->name, "cohesion", 0.0 ),
- 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_GetBool( cf, self->name, "boundaryBottom", False ),
- Stg_ComponentFactory_GetBool( cf, self->name, "boundaryTop", False ),
- Stg_ComponentFactory_GetBool( cf, self->name, "boundaryLeft", False ),
- Stg_ComponentFactory_GetBool( cf, self->name, "boundaryRight", False ),
- Stg_ComponentFactory_GetBool( cf, self->name, "boundaryFront", False ),
- Stg_ComponentFactory_GetBool( cf, self->name, "boundaryBack", False ),
- Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ),
- Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "HydrostaticTerm", HydrostaticTerm, False, data ) );
-}
-
-void _MohrCoulomb_Build( void* rheology, void* data ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
-
- /* Build parent */
- _YieldRheology_Build( self, data );
-}
-
-void _MohrCoulomb_Initialise( void* rheology, void* data ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- Particle_Index lParticle_I;
- Particle_Index particleLocalCount;
- AbstractContext* context = (AbstractContext*)data;
-
- _YieldRheology_Initialise( self, data );
-
- /* Initialise these components now just in case this function is called before their own _Initialise function */
- Stg_Component_Initialise( self->materialPointsSwarm, data, False );
- Stg_Component_Initialise( self->strainWeakening, data, False );
-
- /* We don't need to Initialise hasYieldedVariable because it's a parent variable and _YieldRheology_Initialise
- * has already been called */
- particleLocalCount = self->hasYieldedVariable->variable->arraySize;
-
- /* If restarting from checkpoint, don't change the parameters on the particles */
- if ( !(context && (True == context->loadFromCheckPoint) ) ) {
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
- Variable_SetValueChar( self->hasYieldedVariable->variable, lParticle_I, False );
- }
- }
-}
-
-void _MohrCoulomb_ModifyConstitutiveMatrix(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi )
-{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
-
- /* Don't want to yield on the first ever solve */
- if ( constitutiveMatrix->previousSolutionExists == False ) {
- return;
- }
-
- /* Calculate and store values of stress, pressure, velocity
- gradients, inavriants so they are only calculated once */
- _MohrCoulomb_StoreCurrentParameters( self, constitutiveMatrix,
- materialPointsSwarm, lElement_I,
- materialPoint, xi );
-
- _YieldRheology_ModifyConstitutiveMatrix(self, constitutiveMatrix,
- materialPointsSwarm,lElement_I,
- materialPoint, xi );
-}
-
-double _MohrCoulomb_GetYieldIndicator(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi )
-{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
-
- return SymmetricTensor_2ndInvariant( self->stress,
- constitutiveMatrix->dim );
-}
-
-double _MohrCoulomb_GetYieldCriterion(void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi)
-{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double minimumYieldStress;
- double effectiveCohesion;
- double effectiveFrictionCoefficient;
- double frictionalStrength;
- Dimension_Index dim = constitutiveMatrix->dim;
-
- FeVariable* pressureField = self->pressureField;
- double pressure;
- double factor;
- Cell_Index cell_I;
- Coord coord;
- Element_GlobalIndex element_gI = 0;
- unsigned inds[3];
- Grid* elGrid;
- Bool inside_boundary;
- 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)
- {
- pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,coord);
- }
-
- /* Normally we add the average of the trace of the stress. With
- compressible material, you have to do it. But with stabilized
- linear pressure elements, the non-zero trace is a numerical
- artifact. So we do not add it. */
-
-/* pressure+=self->trace/dim; */
-
- /* 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. */
-
-
- /* Big song and dance to see if we are at a boundary that we care about */
- elGrid = *(Grid**)ExtensionManager_Get(pressureField->feMesh->info,
- pressureField->feMesh,
- ExtensionManager_GetHandle(pressureField->feMesh->info, "elementGrid" ) );
-
- element_gI = FeMesh_ElementDomainToGlobal( pressureField->feMesh, lElement_I );
- RegularMeshUtils_Element_1DTo3D( pressureField->feMesh, element_gI, inds );
-
- inside_boundary=(self->boundaryBottom && inds[1]==0)
- || (self->boundaryTop && inds[1]==elGrid->sizes[1]-1)
- || (self->boundaryLeft && inds[0]==0)
- || (self->boundaryRight && inds[0]==elGrid->sizes[0]-1)
- || (dim==3 && ((self->boundaryBack && inds[2]==0)
- || (self->boundaryFront && inds[2]==elGrid->sizes[2]-1)));
-
- effectiveFrictionCoefficient =
- _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint,
- inside_boundary );
- effectiveCohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint,
- inside_boundary);
-
-
- if(dim==2)
- {
- /* effectiveFrictionCoefficient=tan(phi). If
- factor=sin(atan(1/tan(phi))) =>
- factor=cos(phi)=1/sqrt(1+tan(phi)**2) */
- factor=1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
- frictionalStrength = effectiveFrictionCoefficient*pressure*factor
- + effectiveCohesion*factor;
- }
- else
- {
- double cos_phi, sin_phi;
- /* cos(phi)=1/sqrt(1+tan(phi)**2) */
- cos_phi=
- 1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
- sin_phi=effectiveFrictionCoefficient*cos_phi;
- factor=2*cos_phi/(sqrt(3.0)*(3-sin_phi));
-
- /* The full expression is
-
- sqrt(J2)=p*2*sin(phi)/(sqrt(3)*(3-sin(phi)))
- + C*6*cos(phi)/(sqrt(3)*(3-sin(phi)))
-
- Note the extra factor of 3 for cohesion */
-
- frictionalStrength = effectiveFrictionCoefficient*factor*pressure
- + effectiveCohesion*3*factor;
- }
-
- /* If the minimumYieldStress is not set, then use the
- effective cohesion. Maybe it should be the modified
- effective cohesion, though that probably should not matter
- much. */
- minimumYieldStress = self->minimumYieldStress;
- if(minimumYieldStress==0.0)
- minimumYieldStress=effectiveCohesion;
-
- /* Make sure frictionalStrength is above the minimum */
- if ( frictionalStrength < minimumYieldStress*factor)
- frictionalStrength = minimumYieldStress*factor;
-
- return frictionalStrength;
-}
-
-/* This scales the viscosity so that the stress is the yield stress.
- We can simply scale the viscosity by yieldCriterion/yieldIndicator
- because yieldCriterion does not depend on the deviatoric stress
- (and thus the viscosity), while the yieldIndicator depends linearly
- on the deviatoric stresses. */
-
-void _MohrCoulomb_HasYielded(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- double yieldCriterion,
- double yieldIndicator )
-{
- MohrCoulomb* self = (MohrCoulomb*)rheology;
- double viscosity = yieldCriterion/(2*self->strainRateSecondInvariant);
-
- ConstitutiveMatrix_SetIsotropicViscosity( constitutiveMatrix, viscosity );
-}
-
-double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint,
- Bool inside_boundary ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double effectiveCohesion;
-
- double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- if(inside_boundary)
- {
- 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, Bool inside_boundary ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double effectiveFrictionCoefficient = 0.0;
-
- double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- if(inside_boundary)
- {
- effectiveFrictionCoefficient =
- self->boundaryFrictionCoefficient * (1.0 - strainWeakeningRatio) +
- self->boundaryFrictionCoefficientAfterSoftening
- *strainWeakeningRatio;
- }
- else
- {
- effectiveFrictionCoefficient =
- self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
- self->frictionCoefficientAfterSoftening * strainWeakeningRatio;
- }
-
- return effectiveFrictionCoefficient;
-}
-
-void _MohrCoulomb_StoreCurrentParameters(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi )
-{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- SymmetricTensor strainRate;
- Dimension_Index dim = constitutiveMatrix->dim;
- double strainRateTrace;
- int i;
-
- FeVariable_InterpolateWithinElement( self->pressureField, lElement_I, xi,
- &self->currentPressure );
- FeVariable_InterpolateWithinElement(self->velocityGradientsField,lElement_I,
- xi, self->currentVelocityGradient );
- TensorArray_GetSymmetricPart(self->currentVelocityGradient,dim,strainRate);
-
- ConstitutiveMatrix_CalculateStress( constitutiveMatrix, strainRate,
- self->stress );
- SymmetricTensor_GetTrace(self->stress, dim, &(self->trace));
- SymmetricTensor_GetTrace(strainRate, dim, &strainRateTrace);
-
- /* Subtract the trace so that the 2nd invariant is not polluted by
- the trace */
-
- for(i=0;i<dim;++i)
- {
- strainRate[TensorMapST3D[i][i]]-=strainRateTrace/dim;
- self->stress[TensorMapST3D[i][i]]-=self->trace/dim;
- }
-
- self->strainRateSecondInvariant=SymmetricTensor_2ndInvariant(strainRate,dim);
-}
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,165 +0,0 @@
-/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-** Copyright (c) 2005, Monash Cluster Computing
-** All rights reserved.
-** Redistribution and use in source and binary forms, with or without modification,
-** are permitted provided that the following conditions are met:
-**
-** * Redistributions of source code must retain the above copyright notice,
-** this list of conditions and the following disclaimer.
-** * Redistributions in binary form must reproduce the above copyright
-** notice, this list of conditions and the following disclaimer in the
-** documentation and/or other materials provided with the distribution.
-** * Neither the name of the Monash University nor the names of its contributors
-** may be used to endorse or promote products derived from this software
-** without specific prior written permission.
-**
-** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
-** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
-** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
-** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
-** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
-** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-**
-**
-** Contact:
-*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
-*%
-** Contributors:
-*+ Robert Turnbull
-*+ Vincent Lemiale
-*+ Louis Moresi
-*+ David May
-*+ David Stegman
-*+ Mirko Velic
-*+ Patrick Sunter
-*+ Julian Giordani
-*+
-** $Id$
-**
-**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
-#ifndef __Underworld_Rheology_MohrCoulomb_h__
-#define __Underworld_Rheology_MohrCoulomb_h__
-
- /** Textual name of this class - This is a global pointer which is used for times when you need to refer to class and not a particular instance of a class */
- extern const Type MohrCoulomb_Type;
-
- /** Rheology class contents - this is defined as a macro so that sub-classes of this class can use this macro at the start of the definition of their struct */
- #define __MohrCoulomb \
- /* Parent info */ \
- __YieldRheology \
- /* Virtual functions go here */ \
- /* Material Parameters */\
- ExtensionInfo_Index particleExtHandle; \
- /* Param passed in */ \
- MaterialPointsSwarm* materialPointsSwarm; \
- FeVariable* pressureField; \
- FeVariable* velocityGradientsField; \
- /* Director component is used to update the normal */\
- Director* director; \
- double cohesion; \
- double cohesionAfterSoftening; \
- double frictionCoefficient; \
- double frictionCoefficientAfterSoftening; \
- double boundaryCohesion; \
- double boundaryCohesionAfterSoftening; \
- double boundaryFrictionCoefficient; \
- double boundaryFrictionCoefficientAfterSoftening; \
- Bool boundaryTop; \
- Bool boundaryBottom; \
- Bool boundaryLeft; \
- Bool boundaryRight; \
- Bool boundaryFront; \
- Bool boundaryBack; \
- double minimumYieldStress; \
- /* Stored values that are calculated once for each particle */ \
- double trace; \
- TensorArray currentVelocityGradient; \
- SymmetricTensor stress; \
- double currentPressure; \
- double strainRateSecondInvariant; \
- FeVariable* strainRateField; \
- HydrostaticTerm* hydrostaticTerm;
-
-
- struct MohrCoulomb { __MohrCoulomb };
-
- /** Private Constructor: This will accept all the virtual functions for this class as arguments. */
- MohrCoulomb* _MohrCoulomb_New(
- SizeT sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- Rheology_ModifyConstitutiveMatrixFunction* _modifyConstitutiveMatrix,
- YieldRheology_GetYieldCriterionFunction* _getYieldCriterion,
- YieldRheology_GetYieldIndicatorFunction* _getYieldIndicator,
- YieldRheology_HasYieldedFunction* _hasYielded,
- Name name ) ;
-
- /* 'Stg_Component' implementations */
- void* _MohrCoulomb_DefaultNew( Name name ) ;
- void _MohrCoulomb_Construct( void* rheology, Stg_ComponentFactory* cf,
- void *data );
-
- void _MohrCoulomb_Build( void* rheology, void* data );
- void _MohrCoulomb_Initialise( void* rheology, void* data ) ;
-
- /* 'YieldRheology' implementations */
- void _MohrCoulomb_ModifyConstitutiveMatrix(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi );
-
- double _MohrCoulomb_GetYieldCriterion(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi );
-
- double _MohrCoulomb_GetYieldIndicator(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi );
-
- void _MohrCoulomb_HasYielded(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- double yieldCriterion,
- double yieldIndicator );
-
- double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint, Bool inside_boundary ) ;
- double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, Bool inside_boundary ) ;
- double _MohrCoulomb_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) ;
-
- void _MohrCoulomb_StoreCurrentParameters(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- Coord xi ) ;
-
-#endif
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2009-03-25 19:32:43 UTC (rev 14450)
@@ -1,90 +0,0 @@
-<?xml version="1.0"?>
-<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
-<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
-
-<param name="Name">MohrCoulomb</param>
-<param name="Organisation">MCC</param>
-<param name="Project">Underworld</param>
-<param name="Location">./Underworld/Rheology/src/</param>
-<param name="Project Web">http://mcc.monash.edu.au/Underworld</param>
-<param name="Copyright">Copyright (c) 2005, Monash Cluster Computing</param>
-<param name="License">http://www.opensource.org/licenses/bsd-license.php</param>
-<param name="Parent">YieldRheology</param>
-<param name="Description">Implements faulting behaviour based on Mohr-Coulomb failure criterion. See following reference for details : Anisotropic viscous models of large-deformation Mohr–Coulomb failure, Moresi, L.; Mühlhaus, H.-B., Philosophical Magazine, Volume 86, Numbers 21-22, -22/21 July–1 August 2006, pp. 3287-3305</param>
-
-<!--Now the interesting stuff-->
-
-<list name="Params">
- <struct>
- <param name="Name">cohesion</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">Cohesion for a pristine material</param>
- </struct>
- <struct>
- <param name="Name">cohesionAfterSoftening</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">Value of the cohesion after the material has weakened</param>
- </struct>
- <struct>
- <param name="Name">frictionCoefficient</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">Defines the dependance of the criterion to the normal stress</param>
- </struct>
- <struct>
- <param name="Name">frictionCoefficientAfterSoftening</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">This is the value of the frictionCoefficient after the material has weakened</param>
- </struct>
- <struct>
- <param name="Name">minimumYieldStress</param>
- <param name="Type">Double</param>
- <param name="Default">0.0</param>
- <param name="Description">The value of the failure criterion (defining the yield surface) can not be smaller than this parameter</param>
- </struct>
-</list>
-
-<list name="Dependencies">
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">FeVariable</param>
- <param name="Fallback Key">PressureField</param>
- <param name="Description">Essential since the MohrCoulomb failure criterion is pressure dependant</param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">IntegrationPointsSwarm</param>
- <param name="Description">Defines the set of material points</param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">FeVariable</param>
- <param name="Fallback Key">VelocityGradientsFieldField</param>
- <param name="Description">The resolved component of this tensor wrt the failure plane are mainly used to determine the most favourable oriented plane for failure to occur </param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">StrainWeakening</param>
- <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 -->
-<param name="Example">Underworld/InputFiles/ExtensionFMM.xml</param>
-<param name="Example">Underworld/InputFiles/ExtensionFMM3D.xml</param>
-
-</StGermainData>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h 2009-03-25 19:32:43 UTC (rev 14450)
@@ -67,7 +67,6 @@
#include "Byerlee.h"
#include "DruckerPrager.h"
#include "FaultingMoresiMuhlhaus2006.h"
- #include "MohrCoulomb.h"
#include "StrainWeakening.h"
#include "BuiterStrainWeakening.h"
#include "Director.h"
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/SConscript
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/SConscript 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/SConscript 2009-03-25 19:32:43 UTC (rev 14450)
@@ -35,7 +35,6 @@
FrankKamenetskii.h
Init.h
MaterialViscosity.h
-MohrCoulomb.h
MultiRheologyMaterial.h
NonNewtonian.h
OrthotropicAligned.h
@@ -68,7 +67,6 @@
FrankKamenetskii.c
Init.c
MaterialViscosity.c
-MohrCoulomb.c
MultiRheologyMaterial.c
NonNewtonian.c
OrthotropicAligned.c
@@ -97,7 +95,6 @@
FaultingMoresiMuhlhaus2006.meta
FrankKamenetskii.meta
MaterialViscosity.meta
-MohrCoulomb.meta
MultiRheologyMaterial.meta
NonNewtonian.meta
OrthotropicAligned.meta
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h 2009-03-25 19:19:07 UTC (rev 14449)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h 2009-03-25 19:32:43 UTC (rev 14450)
@@ -69,7 +69,6 @@
typedef struct DruckerPrager DruckerPrager;
typedef struct FaultingMoresiMuhlhaus2006 FaultingMoresiMuhlhaus2006;
- typedef struct MohrCoulomb MohrCoulomb;
typedef struct StrainWeakening StrainWeakening;
typedef struct BuiterStrainWeakening BuiterStrainWeakening;
More information about the CIG-COMMITS
mailing list