[cig-commits] r4353 - in long/3D/Gale/trunk/src/Underworld: .
Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Thu Aug 17 17:18:08 PDT 2006
Author: walter
Date: 2006-08-17 17:18:06 -0700 (Thu, 17 Aug 2006)
New Revision: 4353
Added:
long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.meta
Removed:
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/src/Underworld/
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/types.h
Log:
r491 at earth: boo | 2006-08-17 17:15:27 -0700
r465 at earth (orig r296): VincentLemiale | 2006-08-03 22:55:09 -0700
*** For Mohr-Coulomb rheology users - please read this commit ***
Louis and myself decided to rename 'Mohr-Coulomb' component to 'FaultingMoresiMuhlhaus2006',
because we both agreed the name was inappropriate.
'Mohr-Coulomb' actually refer to a global failure criterion that can be represented as a yield
envelope in the principle stress space.
Although this yielding component makes use of a Mohr-Coulomb type failure criterion,
it also implies other physical assumptions. Details about the formulation
can be found in the following reference :
Anisotropic viscous models of large-deformation Mohr?\226?\128?\147Coulomb failure
Moresi, L.; M?\195?\188hlhaus, H.-B.
Philosophical Magazine, Volume 86, Numbers 21-22, -22/21 July?\226?\128?\1471 August 2006,
pp. 3287-3305
The current implementation is consistent with this paper.
Basically, a Mohr-Coulomb criterion is used on each lagrangian points (particles) to
determine if a failure plane is (re)activated at this particular point. If
yes, the preferred plane of sliding (assumed to be the one that is more
closely aligned with the local sense of shear strain rate) defines a direction of anisotropy.
The material properties can then be softened to model the localization of
deformation. Once a plane of failure is activated, it is assumed to be a
preferred oriention for subsequent sliding, but in case it is no longer
active, new orientations are tested.
As one can see, this formulation does not simply model a Mohr-Coulomb failure criterion
in the strict sense. It is actually a way of incoporating brittle (or
semi-brittle) material behaviour in a viscous code. If anyone wants to use it, it is
important to bear in mind the underlying assumptions of this model.
Thanks,
Vincent
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:490
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:295
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:491
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:296
Added: long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.c 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.c 2006-08-18 00:18:06 UTC (rev 4353)
@@ -0,0 +1,937 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** 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 <StG_FEM/StG_FEM.h>
+#include <PICellerator/PICellerator.h>
+
+#include "types.h"
+#include "RheologyClass.h"
+#include "StrainWeakening.h"
+#include "YieldRheology.h"
+#include "FaultingMoresiMuhlhaus2006.h"
+#include "ConstitutiveMatrix.h"
+#include "Director.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 FaultingMoresiMuhlhaus2006_Type = "FaultingMoresiMuhlhaus2006";
+
+/* Private Constructor: This will accept all the virtual functions for this class as arguments. */
+FaultingMoresiMuhlhaus2006* _FaultingMoresiMuhlhaus2006_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 )
+{
+ FaultingMoresiMuhlhaus2006* 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(FaultingMoresiMuhlhaus2006) );
+ self = (FaultingMoresiMuhlhaus2006*) _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 _FaultingMoresiMuhlhaus2006_Init(
+ FaultingMoresiMuhlhaus2006* self,
+ FeVariable* pressureField,
+ FeVariable* velocityGradientsField,
+ MaterialPointsSwarm* materialPointsSwarm,
+ FiniteElementContext* context,
+ Director* director,
+ double cohesion,
+ double cohesionAfterSoftening,
+ double frictionCoefficient,
+ double frictionCoefficientAfterSoftening,
+ double minimumYieldStress,
+ Bool ignoreOldOrientation )
+{
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+ StandardParticle materialPoint;
+ Dimension_Index dim = materialPointsSwarm->dim;
+
+ self->materialPointsSwarm = materialPointsSwarm;
+ self->pressureField = pressureField;
+ self->velocityGradientsField = velocityGradientsField;
+
+ self->director = director;
+
+ self->particleExtHandle = ExtensionManager_Add( materialPointsSwarm->particleExtensionMgr,
+ FaultingMoresiMuhlhaus2006_Type, sizeof(FaultingMoresiMuhlhaus2006_Particle) );
+
+ self->cohesion = cohesion;
+ self->frictionCoefficient = frictionCoefficient;
+
+ /* Strain softening of Cohesion - (linear weakening is assumed) */
+ /* needs a softening factor between +0 and 1 and a reference strain > 0 */
+ self->cohesionAfterSoftening = cohesionAfterSoftening;
+
+ /* Strain softening of Friction - (linear weakening is assumed) */
+ /* needs a softening factor between +0 and 1 and a reference strain > 0 */
+ self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
+
+ self->minimumYieldStress = minimumYieldStress;
+
+ /* Should orientation from previous timestep be tested ? */
+ self->ignoreOldOrientation = ignoreOldOrientation;
+
+ /* Update Drawing Parameters */
+ EP_PrependClassHook( Context_GetEntryPoint( context, AbstractContext_EP_DumpClass ),
+ _FaultingMoresiMuhlhaus2006_UpdateDrawParameters, self );
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, &materialPoint, self->particleExtHandle );
+
+ /* Add variables for viz purpose */
+
+ self->slipRate = Swarm_NewScalarVariable(
+ materialPointsSwarm,
+ "SlipRate",
+ (ArithPointer) &particleExt->slipRate - (ArithPointer) &materialPoint,
+ Variable_DataType_Double );
+
+ /* slip vector gives the orientation of the failure plane */
+ self->slip = Swarm_NewVectorVariable(
+ materialPointsSwarm,
+ "SlipVector",
+ (ArithPointer) &particleExt->slip - (ArithPointer) &materialPoint,
+ Variable_DataType_Double,
+ dim,
+ "SlipVectorX",
+ "SlipVectorY",
+ "SlipVectorZ" );
+
+ /* Some visualisation parameters (brightness, opacity, length, thickness) */
+ self->brightness = Swarm_NewScalarVariable(
+ materialPointsSwarm,
+ "FaultingMoresiMuhlhaus2006Brightness",
+ (ArithPointer) &particleExt->brightness - (ArithPointer) &materialPoint,
+ Variable_DataType_Float );
+
+ self->opacity = Swarm_NewScalarVariable(
+ materialPointsSwarm,
+ "FaultingMoresiMuhlhaus2006Opacity",
+ (ArithPointer) &particleExt->opacity - (ArithPointer) &materialPoint,
+ Variable_DataType_Float );
+
+ self->length = Swarm_NewScalarVariable(
+ materialPointsSwarm,
+ "FaultingMoresiMuhlhaus2006Length",
+ (ArithPointer) &particleExt->length - (ArithPointer) &materialPoint,
+ Variable_DataType_Float );
+
+ self->thickness = Swarm_NewScalarVariable(
+ materialPointsSwarm,
+ "FaultingMoresiMuhlhaus2006Thickness",
+ (ArithPointer) &particleExt->thickness - (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,
+ "FaultingMoresiMuhlhaus2006TensileFailure",
+ (ArithPointer) &particleExt->tensileFailure - (ArithPointer) &materialPoint,
+ Variable_DataType_Char );
+}
+
+void* _FaultingMoresiMuhlhaus2006_DefaultNew( Name name ) {
+ return (void*) _FaultingMoresiMuhlhaus2006_New(
+ sizeof(FaultingMoresiMuhlhaus2006),
+ FaultingMoresiMuhlhaus2006_Type,
+ _YieldRheology_Delete,
+ _YieldRheology_Print,
+ _YieldRheology_Copy,
+ _FaultingMoresiMuhlhaus2006_DefaultNew,
+ _FaultingMoresiMuhlhaus2006_Construct,
+ _FaultingMoresiMuhlhaus2006_Build,
+ _FaultingMoresiMuhlhaus2006_Initialise,
+ _YieldRheology_Execute,
+ _YieldRheology_Destroy,
+ _FaultingMoresiMuhlhaus2006_ModifyConstitutiveMatrix,
+ _FaultingMoresiMuhlhaus2006_GetYieldCriterion,
+ _FaultingMoresiMuhlhaus2006_GetYieldIndicator,
+ _FaultingMoresiMuhlhaus2006_HasYielded,
+ name );
+}
+
+void _FaultingMoresiMuhlhaus2006_Construct( void* rheology, Stg_ComponentFactory* cf ){
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*)rheology;
+ FeVariable* pressureField;
+ MaterialPointsSwarm* materialPointsSwarm;
+ FeVariable* velocityGradientsField;
+ FiniteElementContext* context;
+ Director* director;
+
+ /* Construct Parent */
+ _YieldRheology_Construct( self, cf );
+
+ /* 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': FaultingMoresiMuhlhaus2006 rheology needs strain weakening.\n",
+ __func__, self->type, self->name );
+
+ context = Stg_ComponentFactory_ConstructByName( cf,
+ "context", FiniteElementContext, True );
+ materialPointsSwarm = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "MaterialPointsSwarm", MaterialPointsSwarm, True );
+ pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "PressureField", FeVariable, True );
+ velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
+ "VelocityGradientsField", FeVariable, True );
+ director = Stg_ComponentFactory_ConstructByKey( cf, self->name, "Director", Director, True );
+
+ _FaultingMoresiMuhlhaus2006_Init(
+ self,
+ pressureField,
+ velocityGradientsField,
+ materialPointsSwarm,
+ context,
+ director,
+ 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, "minimumYieldStress", 0.0 ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "ignoreOldOrientation", False ) );
+}
+
+void _FaultingMoresiMuhlhaus2006_Build( void* rheology, void* data ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+
+ /* Build parent */
+ _YieldRheology_Build( self, data );
+
+ Stg_Component_Build( self->slipRate, data, False );
+ Stg_Component_Build( self->slip, data, False );
+ Stg_Component_Build( self->brightness, data, False );
+ Stg_Component_Build( self->opacity, data, False );
+ Stg_Component_Build( self->length, data, False );
+ Stg_Component_Build( self->thickness, data, False );
+ Stg_Component_Build( self->tensileFailure, data, False );
+
+}
+
+void _FaultingMoresiMuhlhaus2006_Initialise( void* rheology, void* data ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ XYZ normal = { 0.0,1.0,0.0 };
+ XYZ slip = { 1.0,0.0,0.0 };
+ double* ptr;
+ Particle_Index lParticle_I;
+ Particle_Index particleLocalCount;
+ double normalLength2;
+ double invNormalLength;
+ double initialDamageFraction;
+ double* normalDirector;
+ MaterialPoint* materialPoint;
+ Index dof_I;
+ Dimension_Index dim = self->materialPointsSwarm->dim;
+ 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->director, data, False );
+ Stg_Component_Initialise( self->strainWeakening, data, False );
+
+ /* Initialise variables that I've created - (mainly just SwarmVariables)
+ * This will run a Variable_Update for us */
+ Stg_Component_Initialise( self->slipRate, data, False );
+ Stg_Component_Initialise( self->slip, data, False );
+ Stg_Component_Initialise( self->brightness, data, False );
+ Stg_Component_Initialise( self->opacity, data, False );
+ Stg_Component_Initialise( self->length, data, False );
+ Stg_Component_Initialise( self->thickness, data, False );
+ Stg_Component_Initialise( self->tensileFailure, 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 );
+
+ Variable_SetValueDouble( self->slipRate->variable, lParticle_I, 0.0 );
+
+ ptr = Variable_GetPtrDouble( self->slip->variable, lParticle_I );
+
+ materialPoint = (MaterialPoint*)Swarm_ParticleAt( self->materialPointsSwarm, lParticle_I );
+
+ normalDirector = Director_GetNormalPtr( self->director, materialPoint);
+
+ initialDamageFraction = StrainWeakening_GetInitialDamageFraction( self->strainWeakening, materialPoint );
+
+ if (drand48() < initialDamageFraction) {
+ normalLength2 = 0.0;
+
+ for( dof_I=0; dof_I < dim ; dof_I++) {
+
+ normal[dof_I] = 1.0 - 2.0 * drand48();
+ normalLength2 += normal[dof_I] * normal[dof_I];
+ }
+
+ invNormalLength = 1.0/sqrt(normalLength2);
+
+ for( dof_I=0; dof_I < dim ; dof_I++){
+ normal[dof_I] *= invNormalLength;
+ }
+
+ /*TODO : improve this initialisation (is it really needed ?)*/
+ slip[0] = - normal[1];
+ slip[1] = normal[0];
+ slip[2] = 0.0;
+
+ }
+
+ memcpy( normalDirector, normal, sizeof(Coord) );
+ memcpy( ptr, slip, sizeof(Coord) );
+
+ Variable_SetValueDouble( self->slipRate->variable, lParticle_I, 0.0 );
+
+ Variable_SetValueFloat( self->opacity->variable, lParticle_I, 0.0 );
+ Variable_SetValueFloat( self->length->variable, lParticle_I, 0.0 );
+ Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
+ Variable_SetValueChar( self->tensileFailure->variable,lParticle_I, False );
+ }
+ }
+}
+
+void _FaultingMoresiMuhlhaus2006_ModifyConstitutiveMatrix(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double postFailureWeakening;
+ double yieldCriterion;
+ double yieldIndicator; /* A materialPoint will yield if yieldCriterion < yieldIndicator */
+
+ /* Don't want to yield on the first ever solve */
+ if ( constitutiveMatrix->previousSolutionExists == False ) {
+ return;
+ }
+
+ /* Calculate and store values of stress, pressure, velocity gradients, eigenvectors so they are only calculated once */
+ _FaultingMoresiMuhlhaus2006_StoreCurrentParameters( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint, xi );
+
+ postFailureWeakening = StrainWeakening_GetPostFailureWeakening( self->strainWeakening, materialPoint );
+
+ /* First part : treatment of the existing weakened directions */
+ if ( !self->ignoreOldOrientation && postFailureWeakening > 0.0 ) {
+
+ /* tryingOldOrientation is a flag used to know where we are. This is necessary because some parts are treated differently
+ * if we are dealing with old orientation or with pristine material (see for example _FaultingMoresiMuhlhaus2006_GetYieldIndicator) */
+ self->tryingOldOrientation = True;
+
+ if ( _FaultingMoresiMuhlhaus2006_OldOrientationStillSoftening( self, materialPointsSwarm, materialPoint, constitutiveMatrix->dim ) ) {
+ yieldCriterion = _FaultingMoresiMuhlhaus2006_GetYieldCriterion( self, constitutiveMatrix, materialPointsSwarm,
+ lElement_I, materialPoint, xi );
+ yieldIndicator = _FaultingMoresiMuhlhaus2006_GetYieldIndicator( self, constitutiveMatrix, materialPointsSwarm,
+ lElement_I, materialPoint, xi );
+
+ /* Set a bool to TRUE or FLAG depending on whether a materialPoint has failed or not */
+ YieldRheology_SetParticleFlag( self, materialPointsSwarm, materialPoint, (yieldCriterion < yieldIndicator) );
+
+ if ( yieldCriterion < yieldIndicator ) {
+ StrainWeakening_AssignIncrement( self->strainWeakening, constitutiveMatrix,
+ materialPointsSwarm, lElement_I, materialPoint, yieldCriterion, yieldIndicator );
+ _FaultingMoresiMuhlhaus2006_HasYielded( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint,
+ yieldCriterion, yieldIndicator );
+ return;
+ }
+ }
+ }
+
+ /* Second part : treatment of pristine material with no direction favoured */
+ self->tryingOldOrientation = False;
+ _YieldRheology_ModifyConstitutiveMatrix( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint, xi );
+}
+
+double _FaultingMoresiMuhlhaus2006_GetYieldIndicator(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double* stress = self->currentStress;
+ Eigenvector* eigenvectorList = self->currentEigenvectorList;
+ double sigma_ns;
+ Dimension_Index dim = constitutiveMatrix->dim;
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+ double stressMin;
+ double stressMax;
+ double theta;
+ XYZ normal;
+
+ Director_GetNormal( self->director, materialPoint, normal );
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+
+ if ( self->tryingOldOrientation ) {
+ sigma_ns = SymmetricTensor_MultiplyByVectors( stress, normal, particleExt->slip, dim );
+ }
+ else {
+ /* Failure criterion from stress field -
+ * we can calculate whether the material has failed at the current stress from the principle stresses
+ * (without having to compute the orientation) if it does fail,
+ * we can then calculate the orientation for the orthotropy model and history */
+ theta = 0.5 * atan( 1.0/ _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( self, materialPoint ) );
+
+ stressMin = eigenvectorList[0].eigenvalue;
+ stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
+
+ sigma_ns = 0.5 * sin( 2.0 * theta ) * ( stressMax - stressMin );
+ }
+
+ return fabs(sigma_ns);
+}
+
+double _FaultingMoresiMuhlhaus2006_GetYieldCriterion(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double minimumYieldStress;
+ double effectiveCohesion;
+ double effectiveFrictionCoefficient;
+ double frictionalStrength;
+ double sigma_nn;
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+
+ /* Calculate frictional strength */
+ effectiveFrictionCoefficient = _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( self, materialPoint );
+ effectiveCohesion = _FaultingMoresiMuhlhaus2006_EffectiveCohesion( self, materialPoint );
+ sigma_nn = _FaultingMoresiMuhlhaus2006_Sigma_nn( self, materialPoint, constitutiveMatrix->dim );
+
+ frictionalStrength = effectiveFrictionCoefficient * sigma_nn + effectiveCohesion;
+
+ /* Check if it should break in tension (strictly : frictionalStrength < tensileStrength) */
+ particleExt->tensileFailure = (frictionalStrength < 0.0);
+
+ /* Make sure frictionalStrength is above the minimum */
+ minimumYieldStress = self->minimumYieldStress;
+
+ if ( frictionalStrength < minimumYieldStress)
+ frictionalStrength = minimumYieldStress;
+
+ return frictionalStrength;
+}
+
+void _FaultingMoresiMuhlhaus2006_HasYielded(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ double yieldCriterion,
+ double yieldIndicator )
+{
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double beta;
+ double* normal;
+ double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+
+ normal = _FaultingMoresiMuhlhaus2006_UpdateNormalDirection( self, materialPointsSwarm, materialPoint, constitutiveMatrix->dim );
+
+ beta = 1.0 - yieldCriterion / yieldIndicator;
+
+ ConstitutiveMatrix_SetSecondViscosity( constitutiveMatrix, viscosity * beta, normal );
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+ particleExt->slipRate = self->storedSlipRateValue;
+
+ /* Set a flag to tell the director that particles which have failed don't need to have their normal updated */
+ Director_SetDontUpdateParticleFlag( self->director, materialPoint, True );
+}
+
+Bool _FaultingMoresiMuhlhaus2006_OldOrientationStillSoftening( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double* velocityGradient = self->currentVelocityGradient;
+ double* stress = self->currentStress;
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+ double* n;
+ double* s;
+ double tau_nn;
+ double traction[3];
+ double dVparalleldXperpendicular;
+ double dVparalleldXperpendicular1;
+ XYZ normal;
+
+ Director_GetNormal( self->director, materialPoint, normal );
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+ n = normal;
+ s = particleExt->slip;
+
+ /* If an existing weakened direction exists we test this direction for further failure
+ * The slip direction should be set to the direction of maximum shear stress and stored
+ * for plotting
+ * (it does not need to be updated as a history variable like the director ) */
+ if ( dim == 2 ) {
+ traction[0] = stress[0] * n[0] + stress[2] * n[1];
+ traction[1] = stress[2] * n[0] + stress[1] * n[1];
+
+ tau_nn = n[0] * traction[0] + n[1] * traction[1];
+
+ s[0] = traction[0] - tau_nn * n[0];
+ s[1] = traction[1] - tau_nn * n[1];
+ }
+ else {
+ traction[0] = stress[0] * n[0] + stress[3] * n[1] + stress[4] * n[2];
+ traction[1] = stress[3] * n[0] + stress[1] * n[1] + stress[5] * n[2];
+ traction[2] = stress[4] * n[0] + stress[5] * n[1] + stress[2] * n[2];
+
+ tau_nn = n[0] * traction[0] + n[1] * traction[1] + n[2] * traction[2];
+
+ s[0] = traction[0] - tau_nn * n[0];
+ s[1] = traction[1] - tau_nn * n[1];
+ s[2] = traction[2] - tau_nn * n[2];
+ }
+ StGermain_VectorNormalise( s, dim );
+
+ /* Store tau_nn so that we can use it in calculating sigma_nn in _FaultingMoresiMuhlhaus2006_Sigma_nn */
+ self->tau_nn = tau_nn;
+
+ /* Softening direction */
+ dVparalleldXperpendicular = fabs(TensorArray_MultiplyByVectors( velocityGradient, s, n, dim )) ;
+
+ /* Hardening Direction */
+ dVparalleldXperpendicular1 = fabs(TensorArray_MultiplyByVectors( velocityGradient, n, s, dim )) ;
+
+ /* Store this value on class in case we need to store it on the materialPoint if it fails */
+ self->storedSlipRateValue = dVparalleldXperpendicular;
+
+ /* We should only continue testing this plane if it is also in the
+ * softening orientation with respect to the strain-rate gradient rather than hardening */
+ return (dVparalleldXperpendicular1 < dVparalleldXperpendicular);
+}
+
+double* _FaultingMoresiMuhlhaus2006_UpdateNormalDirection( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+ double* velocityGradient = self->currentVelocityGradient;
+ Eigenvector* eigenvectorList = self->currentEigenvectorList;
+ double strainRate_ns[2];
+ XYZ normal[2];
+ XYZ slip[2];
+ double theta;
+ int favourablePlane;
+ double* normalDirector;
+ double tanPhi;
+
+ normalDirector = Director_GetNormalPtr( self->director, materialPoint);
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+
+ /* This function is specific for pristine materials --
+ * We don't need to calculate it for failure in old orientations */
+ if ( self->tryingOldOrientation )
+ return normalDirector;
+
+ if ((fabs(tanPhi = _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( self, materialPoint )))<0.000001)
+ theta = M_PI / 4.0;
+ else
+ theta = 0.5 * atan( 1.0/ tanPhi );
+
+ if (dim == 2){
+
+ /* Identify potential failure directions */
+ StGermain_RotateCoordinateAxis( slip[0], eigenvectorList[0].vector, K_AXIS, +theta);
+ StGermain_RotateCoordinateAxis( slip[1], eigenvectorList[0].vector, K_AXIS, -theta);
+ StGermain_RotateCoordinateAxis( normal[0], eigenvectorList[0].vector, K_AXIS, 0.5*M_PI + theta);
+ StGermain_RotateCoordinateAxis( normal[1], eigenvectorList[0].vector, K_AXIS, 0.5*M_PI - theta);
+ }
+ else {
+
+ /* Identify potential failure directions */
+ StGermain_RotateVector( slip[0], eigenvectorList[0].vector, eigenvectorList[1].vector, + theta );
+ StGermain_RotateVector( slip[1], eigenvectorList[0].vector, eigenvectorList[1].vector, - theta );
+ StGermain_RotateVector( normal[0], eigenvectorList[0].vector, eigenvectorList[1].vector, 0.5*M_PI + theta);
+ StGermain_RotateVector( normal[1], eigenvectorList[0].vector, eigenvectorList[1].vector, 0.5*M_PI - theta);
+ }
+
+ /* Resolve shear strain-rate for the potential failure planes */
+ strainRate_ns[0] = fabs(TensorArray_MultiplyByVectors( velocityGradient, slip[0], normal[0], dim ));
+ strainRate_ns[1] = fabs(TensorArray_MultiplyByVectors( velocityGradient, slip[1], normal[1], dim ));
+
+ /* Choose the plane which is oriented favorably for continued slip */
+ favourablePlane = strainRate_ns[0] > strainRate_ns[1] ? 0 : 1;
+
+ memcpy( normalDirector, normal[favourablePlane], dim * sizeof(double) );
+ memcpy( particleExt->slip, slip[favourablePlane], dim * sizeof(double) );
+
+ self->storedSlipRateValue = strainRate_ns[ favourablePlane ];
+
+ return normalDirector;
+}
+
+double _FaultingMoresiMuhlhaus2006_EffectiveCohesion( void* rheology, void* materialPoint ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double effectiveCohesion;
+
+ if ( self->tryingOldOrientation || self->ignoreOldOrientation ) {
+ /* If this is the old orientation or old orientations are ignored
+ * then the weakening should be considered here */
+ double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
+
+ effectiveCohesion = self->cohesion * (1.0 - strainWeakeningRatio) +
+ self->cohesionAfterSoftening * strainWeakeningRatio;
+ }
+ else {
+ /* If old orientations are given priority and this is an unbroken direction
+ * then consider material as completely pristine */
+ effectiveCohesion = self->cohesion;
+ }
+
+ return effectiveCohesion;
+}
+
+double _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ static double effectiveFrictionCoefficient = 0.0;
+ static int prevTryingOldOrientationValue = -1;
+ static void* prevParticle = NULL;
+
+ /* See if we can return cached value */
+ if ( self->tryingOldOrientation == (Bool) prevTryingOldOrientationValue && materialPoint == prevParticle )
+ return effectiveFrictionCoefficient;
+
+ /* This is a different situation to the cached case - therefore recalculate */
+ if ( self->tryingOldOrientation || self->ignoreOldOrientation ) {
+ /* If this is the old orientation or old orientations are ignored
+ * then the weakening should be considered here */
+ double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
+
+ effectiveFrictionCoefficient = self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
+ self->frictionCoefficientAfterSoftening * strainWeakeningRatio;
+ }
+ else {
+ /* If old orientations are given priority and this is an unbroken direction
+ * then consider material as completely pristine */
+ effectiveFrictionCoefficient = self->frictionCoefficient;
+ }
+
+ /* Cache values */
+ prevParticle = materialPoint;
+ prevTryingOldOrientationValue = self->tryingOldOrientation;
+
+ return effectiveFrictionCoefficient;
+}
+
+double _FaultingMoresiMuhlhaus2006_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ double pressure = self->currentPressure;
+ Eigenvector* eigenvectorList = self->currentEigenvectorList;
+ double tau_nn;
+ double sigma_nn;
+
+ if ( self->tryingOldOrientation ) {
+ /* tau_nn has already been calculated in _FaultingMoresiMuhlhaus2006_OldOrientationStillSoftening */
+ tau_nn = self->tau_nn;
+ }
+ else {
+ double stressMin = eigenvectorList[0].eigenvalue;
+ double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
+ double theta = 0.5 * atan( 1.0/ _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( self, materialPoint ) );
+
+ tau_nn = 0.5 * (( stressMax + stressMin ) + cos( 2.0 * theta ) * ( stressMax - stressMin ));
+ }
+
+ sigma_nn = pressure - tau_nn;
+
+ return sigma_nn;
+}
+
+void _FaultingMoresiMuhlhaus2006_StoreCurrentParameters(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi )
+{
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ SymmetricTensor strainRate;
+ Dimension_Index dim = constitutiveMatrix->dim;
+ FaultingMoresiMuhlhaus2006_Particle* particleExt;
+
+ particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
+ particleExt->slipRate = 0.0;
+
+ 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->currentStress );
+
+ SymmetricTensor_CalcAllEigenvectors( self->currentStress, dim, self->currentEigenvectorList );
+
+ Director_SetDontUpdateParticleFlag( self->director, materialPoint, False );
+}
+
+void _FaultingMoresiMuhlhaus2006_UpdateDrawParameters( void* rheology ) {
+ FaultingMoresiMuhlhaus2006* self = (FaultingMoresiMuhlhaus2006*) rheology;
+ Particle_Index lParticle_I;
+ Particle_Index particleLocalCount;
+ StrainWeakening* strainWeakening = self->strainWeakening;
+ StandardParticle* materialPoint;
+ double slipRate;
+
+ double length;
+ double brightness;
+ double opacity;
+ double strainWeakeningRatio;
+ double localMaxStrainIncrement;
+ double localMaxSlipRate;
+ double localMeanStrainIncrement;
+ double localMeanSlipRate;
+ Particle_Index localFailed;
+
+ double globalMaxSlipRate;
+ double globalMaxStrainIncrement;
+ double globalMeanStrainIncrement;
+ double globalMeanSlipRate;
+ Particle_Index globalFailed;
+
+ double averagedGlobalMaxSlipRate = 0.0;
+ double averagedGlobalMaxStrainIncrement = 0.0;
+
+ double oneOverGlobalMaxSlipRate;
+ double oneOverGlobalMaxStrainIncrement;
+ double postFailureWeakeningIncrement;
+
+ /* This parameter is only needed for the alternative commented out set of parameters */
+ /* double globalMaxStrainWeakeningRatio = 0.0; */
+
+ /* 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 */
+
+ localMaxStrainIncrement = 0.0;
+ localMeanStrainIncrement = 0.0;
+ localMaxSlipRate = 0.0;
+ localMeanSlipRate = 0.0;
+ localFailed = 0;
+
+ /* Update all variables */
+ Variable_Update( self->hasYieldedVariable->variable );
+ Variable_Update( self->slipRate->variable );
+ Variable_Update( self->brightness->variable );
+ Variable_Update( self->opacity->variable );
+ Variable_Update( self->length->variable );
+ Variable_Update( self->thickness->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++;
+ slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
+
+ postFailureWeakeningIncrement =
+ Variable_GetValueDouble( strainWeakening->postFailureWeakeningIncrement->variable, lParticle_I );
+
+ localMeanSlipRate += slipRate;
+ localMeanStrainIncrement += postFailureWeakeningIncrement;
+
+ if(localMaxSlipRate < slipRate)
+ localMaxSlipRate = slipRate;
+
+ if(localMaxStrainIncrement < postFailureWeakeningIncrement)
+ localMaxStrainIncrement = postFailureWeakeningIncrement;
+ }
+ }
+
+ MPI_Allreduce( &localMaxSlipRate, &globalMaxSlipRate, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
+ MPI_Allreduce( &localMaxStrainIncrement, &globalMaxStrainIncrement, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
+ MPI_Allreduce( &localMeanSlipRate, &globalMeanSlipRate, 1, MPI_DOUBLE, MPI_SUM, 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;
+ globalMeanSlipRate /= (double) globalFailed;
+
+ averagedGlobalMaxStrainIncrement =
+ 0.5 * averagedGlobalMaxStrainIncrement +
+ 0.25 * globalMeanStrainIncrement +
+ 0.25 * globalMaxStrainIncrement;
+ averagedGlobalMaxSlipRate =
+ 0.5 * averagedGlobalMaxSlipRate +
+ 0.25 * globalMaxSlipRate +
+ 0.25 * globalMeanSlipRate;
+
+ /* Let's simply assume that twice the mean is a good place to truncate these values */
+ oneOverGlobalMaxSlipRate = 1.0 / averagedGlobalMaxSlipRate;
+ oneOverGlobalMaxStrainIncrement = 1.0 / averagedGlobalMaxStrainIncrement;
+
+ for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
+ 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->length->variable, lParticle_I, 0.0 );
+ Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
+ continue;
+ }
+
+ slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
+ strainWeakeningRatio = StrainWeakening_CalcRatio( strainWeakening, materialPoint );
+
+
+ length = 0.005 + 0.005 * strainWeakeningRatio;
+ length *= 10;
+ brightness = strainWeakeningRatio * slipRate * oneOverGlobalMaxSlipRate;
+ if( brightness > 1.0 )
+ brightness = 1.0;
+ opacity = 0.5 * pow(brightness,3.0);
+
+ Variable_SetValueFloat( self->brightness->variable, lParticle_I, brightness );
+ Variable_SetValueFloat( self->opacity->variable, lParticle_I, opacity );
+ Variable_SetValueFloat( self->length->variable, lParticle_I, (float) length );
+ Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.5 + 10.0 * (float)length );
+
+ /* This parameter is only needed for the alternative commented out set of parameters */
+ /* if (strainWeakeningRatio>globalMaxStrainWeakeningRatio)
+ globalMaxStrainWeakeningRatio = strainWeakeningRatio;
+ */
+ }
+
+ /* This part is another set of parameters that has been used with shear problems.
+ * it's kept here for now, but it shows that in principle we could define as many set of parameters
+ * as we want (depending on the problem ?)*/
+ /*
+ for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
+ materialPoint = Swarm_ParticleAt( strainWeakening->materialPointsSwarm, 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->length->variable, lParticle_I, 0.0 );
+ Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
+ continue;
+ }
+
+ slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
+ strainWeakeningRatio = StrainWeakening_CalcRatio( strainWeakening, materialPoint );
+
+ length = 0.03 + 0.03 * (strainWeakeningRatio/globalMaxStrainWeakeningRatio);
+
+ brightness = 1.0-(strainWeakeningRatio/globalMaxStrainWeakeningRatio);
+
+ if (brightness > 1.0)
+ brightness = 1.0;
+
+ opacity = (slipRate/globalMaxSlipRate);
+
+ if (opacity > 0.90)
+ opacity = 1.0;//this condition is to make sure we have enough planes that will be clearly seen.
+
+ Variable_SetValueFloat( self->brightness->variable, lParticle_I, brightness );
+ Variable_SetValueFloat( self->opacity->variable, lParticle_I, opacity );
+ Variable_SetValueFloat( self->length->variable, lParticle_I, (float) length );
+ Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.5 + 10.0 * (float)length );
+ }
+ */
+}
Property changes on: long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.c
___________________________________________________________________
Name: svn:keywords
+ LastChangedDate Author Id
Name: svn:eol-style
+ native
Added: long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.h 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.h 2006-08-18 00:18:06 UTC (rev 4353)
@@ -0,0 +1,178 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** 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_FaultingMoresiMuhlhaus2006_h__
+#define __Underworld_Rheology_FaultingMoresiMuhlhaus2006_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 FaultingMoresiMuhlhaus2006_Type;
+
+ typedef struct {
+ double slipRate;
+ XYZ slip;
+ float brightness;
+ float opacity;
+ float length;
+ float thickness;
+ Particle_Bool tensileFailure;
+ } FaultingMoresiMuhlhaus2006_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 __FaultingMoresiMuhlhaus2006 \
+ /* Parent info */ \
+ __YieldRheology \
+ /* Virtual functions go here */ \
+ /* Material Parameters */\
+ SwarmVariable* slipRate; \
+ SwarmVariable* slip; \
+ SwarmVariable* brightness; \
+ SwarmVariable* opacity; \
+ SwarmVariable* length; \
+ SwarmVariable* thickness; \
+ SwarmVariable* tensileFailure; \
+ 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 minimumYieldStress; \
+ Bool ignoreOldOrientation; \
+ /* Stored values that are calculated once for each particle */ \
+ Eigenvector currentEigenvectorList[3]; \
+ TensorArray currentVelocityGradient; \
+ SymmetricTensor currentStress; \
+ double currentPressure; \
+ double tau_nn; \
+ double storedSlipRateValue; \
+ /* Flags set to tell functions what's going on */ \
+ Bool tryingOldOrientation;
+
+ struct FaultingMoresiMuhlhaus2006 { __FaultingMoresiMuhlhaus2006 };
+
+ /** Private Constructor: This will accept all the virtual functions for this class as arguments. */
+ FaultingMoresiMuhlhaus2006* _FaultingMoresiMuhlhaus2006_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* _FaultingMoresiMuhlhaus2006_DefaultNew( Name name ) ;
+ void _FaultingMoresiMuhlhaus2006_Construct( void* rheology, Stg_ComponentFactory* cf );
+
+ void _FaultingMoresiMuhlhaus2006_Build( void* rheology, void* data );
+ void _FaultingMoresiMuhlhaus2006_Initialise( void* rheology, void* data ) ;
+
+ /* 'YieldRheology' implementations */
+ void _FaultingMoresiMuhlhaus2006_ModifyConstitutiveMatrix(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ double _FaultingMoresiMuhlhaus2006_GetYieldCriterion(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ double _FaultingMoresiMuhlhaus2006_GetYieldIndicator(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi );
+
+ void _FaultingMoresiMuhlhaus2006_HasYielded(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ double yieldCriterion,
+ double yieldIndicator );
+
+ Bool _FaultingMoresiMuhlhaus2006_OldOrientationStillSoftening( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) ;
+ double* _FaultingMoresiMuhlhaus2006_UpdateNormalDirection( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) ;
+
+ double _FaultingMoresiMuhlhaus2006_EffectiveCohesion( void* rheology, void* materialPoint ) ;
+ double _FaultingMoresiMuhlhaus2006_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) ;
+ double _FaultingMoresiMuhlhaus2006_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) ;
+
+ void _FaultingMoresiMuhlhaus2006_StoreCurrentParameters(
+ void* rheology,
+ ConstitutiveMatrix* constitutiveMatrix,
+ MaterialPointsSwarm* materialPointsSwarm,
+ Element_LocalIndex lElement_I,
+ MaterialPoint* materialPoint,
+ Coord xi ) ;
+
+ void _FaultingMoresiMuhlhaus2006_UpdateDrawParameters( void* rheology ) ;
+
+#endif
Property changes on: long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.h
___________________________________________________________________
Name: svn:keywords
+ LastChangedDate Author Id
Name: svn:eol-style
+ native
Added: long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.meta 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/FaultingMoresiMuhlhaus2006.meta 2006-08-18 00:18:06 UTC (rev 4353)
@@ -0,0 +1,104 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">FaultingMoresiMuhlhaus2006</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>
+ <struct>
+ <param name="Name">ignoreOldOrientation</param>
+ <param name="Type">Bool</param>
+ <param name="Default">False</param>
+ <param name="Description">Boolean to indicate if we want to take into account the oldOrientation</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">FiniteElementContext</param>
+ <param name="Fallback Key">Context</param>
+ <param name="Description">The context is essential because the FaultingMoresiMuhlhaus2006_UpdateDrawParameters function is prepended to an existing Entry Point. </param>
+ </struct>
+
+ <struct>
+ <param name="Essential">True</param>
+ <param name="Type">Director</param>
+ <param name="Description">This component allows to update the normal associated with a failure plane wrt the deformation. At the moment, it is only possible to update all normal at the same time, no matter if they correspond to active failure planes or not</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>
+
+
+</list>
+<!-- Add an example XML if possible -->
+<param name="Example">Underworld/InputFiles/ExtensionFFM.xml</param>
+<param name="Example">Underworld/InputFiles/ExtensionFFM3D.xml</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Init.c 2006-08-18 00:18:06 UTC (rev 4353)
@@ -70,7 +70,7 @@
Stg_ComponentRegister_Add( componentRegister, VonMises_Type, "0", _VonMises_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, Byerlee_Type, "0", _Byerlee_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, DruckerPrager_Type, "0", _DruckerPrager_DefaultNew );
- Stg_ComponentRegister_Add( componentRegister, MohrCoulomb_Type, "0", _MohrCoulomb_DefaultNew );
+ Stg_ComponentRegister_Add( componentRegister, FaultingMoresiMuhlhaus2006_Type, "0", _FaultingMoresiMuhlhaus2006_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, StrainWeakening_Type, "0", _StrainWeakening_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, BuiterStrainWeakening_Type, "0", _BuiterStrainWeakening_DefaultNew );
@@ -93,11 +93,11 @@
RegisterParent( StoreStress_Type, Rheology_Type );
RegisterParent( StoreVisc_Type, Rheology_Type );
- RegisterParent( YieldRheology_Type, Rheology_Type );
- RegisterParent( VonMises_Type, YieldRheology_Type );
- RegisterParent( Byerlee_Type, VonMises_Type );
- RegisterParent( DruckerPrager_Type, VonMises_Type );
- RegisterParent( MohrCoulomb_Type, YieldRheology_Type );
+ RegisterParent( YieldRheology_Type, Rheology_Type );
+ RegisterParent( VonMises_Type, YieldRheology_Type );
+ RegisterParent( Byerlee_Type, VonMises_Type );
+ RegisterParent( DruckerPrager_Type, VonMises_Type );
+ RegisterParent( FaultingMoresiMuhlhaus2006_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 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2006-08-18 00:18:06 UTC (rev 4353)
@@ -1,937 +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 <StG_FEM/StG_FEM.h>
-#include <PICellerator/PICellerator.h>
-
-#include "types.h"
-#include "RheologyClass.h"
-#include "StrainWeakening.h"
-#include "YieldRheology.h"
-#include "MohrCoulomb.h"
-#include "ConstitutiveMatrix.h"
-#include "Director.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,
- FiniteElementContext* context,
- Director* director,
- double cohesion,
- double cohesionAfterSoftening,
- double frictionCoefficient,
- double frictionCoefficientAfterSoftening,
- double minimumYieldStress,
- Bool ignoreOldOrientation )
-{
- MohrCoulomb_Particle* particleExt;
- StandardParticle materialPoint;
- Dimension_Index dim = materialPointsSwarm->dim;
-
- self->materialPointsSwarm = materialPointsSwarm;
- self->pressureField = pressureField;
- self->velocityGradientsField = velocityGradientsField;
-
- self->director = director;
-
- self->particleExtHandle = ExtensionManager_Add( materialPointsSwarm->particleExtensionMgr,
- MohrCoulomb_Type, sizeof(MohrCoulomb_Particle) );
-
- self->cohesion = cohesion;
- self->frictionCoefficient = frictionCoefficient;
-
- /* Strain softening of Cohesion - (linear weakening is assumed) */
- /* needs a softening factor between +0 and 1 and a reference strain > 0 */
- self->cohesionAfterSoftening = cohesionAfterSoftening;
-
- /* Strain softening of Friction - (linear weakening is assumed) */
- /* needs a softening factor between +0 and 1 and a reference strain > 0 */
- self->frictionCoefficientAfterSoftening = frictionCoefficientAfterSoftening;
-
- self->minimumYieldStress = minimumYieldStress;
-
- /* Should orientation from previous timestep be tested ? */
- self->ignoreOldOrientation = ignoreOldOrientation;
-
- /* Update Drawing Parameters */
- EP_PrependClassHook( Context_GetEntryPoint( context, AbstractContext_EP_DumpClass ),
- _MohrCoulomb_UpdateDrawParameters, self );
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, &materialPoint, self->particleExtHandle );
-
- /* Add variables for viz purpose */
-
- self->slipRate = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "SlipRate",
- (ArithPointer) &particleExt->slipRate - (ArithPointer) &materialPoint,
- Variable_DataType_Double );
-
- /* slip vector gives the orientation of the failure plane */
- self->slip = Swarm_NewVectorVariable(
- materialPointsSwarm,
- "SlipVector",
- (ArithPointer) &particleExt->slip - (ArithPointer) &materialPoint,
- Variable_DataType_Double,
- dim,
- "SlipVectorX",
- "SlipVectorY",
- "SlipVectorZ" );
-
- /* Some visualisation parameters (brightness, opacity, length, thickness) */
- self->brightness = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "MohrCoulombBrightness",
- (ArithPointer) &particleExt->brightness - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- self->opacity = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "MohrCoulombOpacity",
- (ArithPointer) &particleExt->opacity - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- self->length = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "MohrCoulombLength",
- (ArithPointer) &particleExt->length - (ArithPointer) &materialPoint,
- Variable_DataType_Float );
-
- self->thickness = Swarm_NewScalarVariable(
- materialPointsSwarm,
- "MohrCoulombThickness",
- (ArithPointer) &particleExt->thickness - (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,
- "MohrCoulombTensileFailure",
- (ArithPointer) &particleExt->tensileFailure - (ArithPointer) &materialPoint,
- Variable_DataType_Char );
-}
-
-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 ){
- MohrCoulomb* self = (MohrCoulomb*)rheology;
- FeVariable* pressureField;
- MaterialPointsSwarm* materialPointsSwarm;
- FeVariable* velocityGradientsField;
- FiniteElementContext* context;
- Director* director;
-
- /* Construct Parent */
- _YieldRheology_Construct( self, cf );
-
- /* 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 );
-
- context = Stg_ComponentFactory_ConstructByName( cf,
- "context", FiniteElementContext, True );
- materialPointsSwarm = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "MaterialPointsSwarm", MaterialPointsSwarm, True );
- pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "PressureField", FeVariable, True );
- velocityGradientsField = Stg_ComponentFactory_ConstructByKey( cf, self->name,
- "VelocityGradientsField", FeVariable, True );
- director = Stg_ComponentFactory_ConstructByKey( cf, self->name, "Director", Director, True );
-
- _MohrCoulomb_Init(
- self,
- pressureField,
- velocityGradientsField,
- materialPointsSwarm,
- context,
- director,
- 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, "minimumYieldStress", 0.0 ),
- Stg_ComponentFactory_GetBool( cf, self->name, "ignoreOldOrientation", False ) );
-}
-
-void _MohrCoulomb_Build( void* rheology, void* data ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
-
- /* Build parent */
- _YieldRheology_Build( self, data );
-
- Stg_Component_Build( self->slipRate, data, False );
- Stg_Component_Build( self->slip, data, False );
- Stg_Component_Build( self->brightness, data, False );
- Stg_Component_Build( self->opacity, data, False );
- Stg_Component_Build( self->length, data, False );
- Stg_Component_Build( self->thickness, data, False );
- Stg_Component_Build( self->tensileFailure, data, False );
-
-}
-
-void _MohrCoulomb_Initialise( void* rheology, void* data ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- XYZ normal = { 0.0,1.0,0.0 };
- XYZ slip = { 1.0,0.0,0.0 };
- double* ptr;
- Particle_Index lParticle_I;
- Particle_Index particleLocalCount;
- double normalLength2;
- double invNormalLength;
- double initialDamageFraction;
- double* normalDirector;
- MaterialPoint* materialPoint;
- Index dof_I;
- Dimension_Index dim = self->materialPointsSwarm->dim;
- 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->director, data, False );
- Stg_Component_Initialise( self->strainWeakening, data, False );
-
- /* Initialise variables that I've created - (mainly just SwarmVariables)
- * This will run a Variable_Update for us */
- Stg_Component_Initialise( self->slipRate, data, False );
- Stg_Component_Initialise( self->slip, data, False );
- Stg_Component_Initialise( self->brightness, data, False );
- Stg_Component_Initialise( self->opacity, data, False );
- Stg_Component_Initialise( self->length, data, False );
- Stg_Component_Initialise( self->thickness, data, False );
- Stg_Component_Initialise( self->tensileFailure, 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 );
-
- Variable_SetValueDouble( self->slipRate->variable, lParticle_I, 0.0 );
-
- ptr = Variable_GetPtrDouble( self->slip->variable, lParticle_I );
-
- materialPoint = (MaterialPoint*)Swarm_ParticleAt( self->materialPointsSwarm, lParticle_I );
-
- normalDirector = Director_GetNormalPtr( self->director, materialPoint);
-
- initialDamageFraction = StrainWeakening_GetInitialDamageFraction( self->strainWeakening, materialPoint );
-
- if (drand48() < initialDamageFraction) {
- normalLength2 = 0.0;
-
- for( dof_I=0; dof_I < dim ; dof_I++) {
-
- normal[dof_I] = 1.0 - 2.0 * drand48();
- normalLength2 += normal[dof_I] * normal[dof_I];
- }
-
- invNormalLength = 1.0/sqrt(normalLength2);
-
- for( dof_I=0; dof_I < dim ; dof_I++){
- normal[dof_I] *= invNormalLength;
- }
-
- /*TODO : improve this initialisation (is it really needed ?)*/
- slip[0] = - normal[1];
- slip[1] = normal[0];
- slip[2] = 0.0;
-
- }
-
- memcpy( normalDirector, normal, sizeof(Coord) );
- memcpy( ptr, slip, sizeof(Coord) );
-
- Variable_SetValueDouble( self->slipRate->variable, lParticle_I, 0.0 );
-
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->length->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
- Variable_SetValueChar( self->tensileFailure->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;
- double postFailureWeakening;
- double yieldCriterion;
- double yieldIndicator; /* A materialPoint will yield if yieldCriterion < yieldIndicator */
-
- /* Don't want to yield on the first ever solve */
- if ( constitutiveMatrix->previousSolutionExists == False ) {
- return;
- }
-
- /* Calculate and store values of stress, pressure, velocity gradients, eigenvectors so they are only calculated once */
- _MohrCoulomb_StoreCurrentParameters( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint, xi );
-
- postFailureWeakening = StrainWeakening_GetPostFailureWeakening( self->strainWeakening, materialPoint );
-
- /* First part : treatment of the existing weakened directions */
- if ( !self->ignoreOldOrientation && postFailureWeakening > 0.0 ) {
-
- /* tryingOldOrientation is a flag used to know where we are. This is necessary because some parts are treated differently
- * if we are dealing with old orientation or with pristine material (see for example _MohrCoulomb_GetYieldIndicator) */
- self->tryingOldOrientation = True;
-
- if ( _MohrCoulomb_OldOrientationStillSoftening( self, materialPointsSwarm, materialPoint, constitutiveMatrix->dim ) ) {
- yieldCriterion = _MohrCoulomb_GetYieldCriterion( self, constitutiveMatrix, materialPointsSwarm,
- lElement_I, materialPoint, xi );
- yieldIndicator = _MohrCoulomb_GetYieldIndicator( self, constitutiveMatrix, materialPointsSwarm,
- lElement_I, materialPoint, xi );
-
- /* Set a bool to TRUE or FLAG depending on whether a materialPoint has failed or not */
- YieldRheology_SetParticleFlag( self, materialPointsSwarm, materialPoint, (yieldCriterion < yieldIndicator) );
-
- if ( yieldCriterion < yieldIndicator ) {
- StrainWeakening_AssignIncrement( self->strainWeakening, constitutiveMatrix,
- materialPointsSwarm, lElement_I, materialPoint, yieldCriterion, yieldIndicator );
- _MohrCoulomb_HasYielded( self, constitutiveMatrix, materialPointsSwarm, lElement_I, materialPoint,
- yieldCriterion, yieldIndicator );
- return;
- }
- }
- }
-
- /* Second part : treatment of pristine material with no direction favoured */
- self->tryingOldOrientation = False;
- _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;
- double* stress = self->currentStress;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double sigma_ns;
- Dimension_Index dim = constitutiveMatrix->dim;
- MohrCoulomb_Particle* particleExt;
- double stressMin;
- double stressMax;
- double theta;
- XYZ normal;
-
- Director_GetNormal( self->director, materialPoint, normal );
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
-
- if ( self->tryingOldOrientation ) {
- sigma_ns = SymmetricTensor_MultiplyByVectors( stress, normal, particleExt->slip, dim );
- }
- else {
- /* Failure criterion from stress field -
- * we can calculate whether the material has failed at the current stress from the principle stresses
- * (without having to compute the orientation) if it does fail,
- * we can then calculate the orientation for the orthotropy model and history */
- theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
-
- stressMin = eigenvectorList[0].eigenvalue;
- stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
-
- sigma_ns = 0.5 * sin( 2.0 * theta ) * ( stressMax - stressMin );
- }
-
- return fabs(sigma_ns);
-}
-
-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;
- double sigma_nn;
- MohrCoulomb_Particle* particleExt;
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
-
- /* Calculate frictional strength */
- effectiveFrictionCoefficient = _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint );
- effectiveCohesion = _MohrCoulomb_EffectiveCohesion( self, materialPoint );
- sigma_nn = _MohrCoulomb_Sigma_nn( self, materialPoint, constitutiveMatrix->dim );
-
- frictionalStrength = effectiveFrictionCoefficient * sigma_nn + effectiveCohesion;
-
- /* Check if it should break in tension (strictly : frictionalStrength < tensileStrength) */
- particleExt->tensileFailure = (frictionalStrength < 0.0);
-
- /* Make sure frictionalStrength is above the minimum */
- minimumYieldStress = self->minimumYieldStress;
-
- if ( frictionalStrength < minimumYieldStress)
- frictionalStrength = minimumYieldStress;
-
- return frictionalStrength;
-}
-
-void _MohrCoulomb_HasYielded(
- void* rheology,
- ConstitutiveMatrix* constitutiveMatrix,
- MaterialPointsSwarm* materialPointsSwarm,
- Element_LocalIndex lElement_I,
- MaterialPoint* materialPoint,
- double yieldCriterion,
- double yieldIndicator )
-{
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double beta;
- double* normal;
- double viscosity = ConstitutiveMatrix_GetIsotropicViscosity( constitutiveMatrix );
- MohrCoulomb_Particle* particleExt;
-
- normal = _MohrCoulomb_UpdateNormalDirection( self, materialPointsSwarm, materialPoint, constitutiveMatrix->dim );
-
- beta = 1.0 - yieldCriterion / yieldIndicator;
-
- ConstitutiveMatrix_SetSecondViscosity( constitutiveMatrix, viscosity * beta, normal );
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
- particleExt->slipRate = self->storedSlipRateValue;
-
- /* Set a flag to tell the director that particles which have failed don't need to have their normal updated */
- Director_SetDontUpdateParticleFlag( self->director, materialPoint, True );
-}
-
-Bool _MohrCoulomb_OldOrientationStillSoftening( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double* velocityGradient = self->currentVelocityGradient;
- double* stress = self->currentStress;
- MohrCoulomb_Particle* particleExt;
- double* n;
- double* s;
- double tau_nn;
- double traction[3];
- double dVparalleldXperpendicular;
- double dVparalleldXperpendicular1;
- XYZ normal;
-
- Director_GetNormal( self->director, materialPoint, normal );
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
- n = normal;
- s = particleExt->slip;
-
- /* If an existing weakened direction exists we test this direction for further failure
- * The slip direction should be set to the direction of maximum shear stress and stored
- * for plotting
- * (it does not need to be updated as a history variable like the director ) */
- if ( dim == 2 ) {
- traction[0] = stress[0] * n[0] + stress[2] * n[1];
- traction[1] = stress[2] * n[0] + stress[1] * n[1];
-
- tau_nn = n[0] * traction[0] + n[1] * traction[1];
-
- s[0] = traction[0] - tau_nn * n[0];
- s[1] = traction[1] - tau_nn * n[1];
- }
- else {
- traction[0] = stress[0] * n[0] + stress[3] * n[1] + stress[4] * n[2];
- traction[1] = stress[3] * n[0] + stress[1] * n[1] + stress[5] * n[2];
- traction[2] = stress[4] * n[0] + stress[5] * n[1] + stress[2] * n[2];
-
- tau_nn = n[0] * traction[0] + n[1] * traction[1] + n[2] * traction[2];
-
- s[0] = traction[0] - tau_nn * n[0];
- s[1] = traction[1] - tau_nn * n[1];
- s[2] = traction[2] - tau_nn * n[2];
- }
- StGermain_VectorNormalise( s, dim );
-
- /* Store tau_nn so that we can use it in calculating sigma_nn in _MohrCoulomb_Sigma_nn */
- self->tau_nn = tau_nn;
-
- /* Softening direction */
- dVparalleldXperpendicular = fabs(TensorArray_MultiplyByVectors( velocityGradient, s, n, dim )) ;
-
- /* Hardening Direction */
- dVparalleldXperpendicular1 = fabs(TensorArray_MultiplyByVectors( velocityGradient, n, s, dim )) ;
-
- /* Store this value on class in case we need to store it on the materialPoint if it fails */
- self->storedSlipRateValue = dVparalleldXperpendicular;
-
- /* We should only continue testing this plane if it is also in the
- * softening orientation with respect to the strain-rate gradient rather than hardening */
- return (dVparalleldXperpendicular1 < dVparalleldXperpendicular);
-}
-
-double* _MohrCoulomb_UpdateNormalDirection( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- MohrCoulomb_Particle* particleExt;
- double* velocityGradient = self->currentVelocityGradient;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double strainRate_ns[2];
- XYZ normal[2];
- XYZ slip[2];
- double theta;
- int favourablePlane;
- double* normalDirector;
- double tanPhi;
-
- normalDirector = Director_GetNormalPtr( self->director, materialPoint);
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
-
- /* This function is specific for pristine materials --
- * We don't need to calculate it for failure in old orientations */
- if ( self->tryingOldOrientation )
- return normalDirector;
-
- if ((fabs(tanPhi = _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint )))<0.000001)
- theta = M_PI / 4.0;
- else
- theta = 0.5 * atan( 1.0/ tanPhi );
-
- if (dim == 2){
-
- /* Identify potential failure directions */
- StGermain_RotateCoordinateAxis( slip[0], eigenvectorList[0].vector, K_AXIS, +theta);
- StGermain_RotateCoordinateAxis( slip[1], eigenvectorList[0].vector, K_AXIS, -theta);
- StGermain_RotateCoordinateAxis( normal[0], eigenvectorList[0].vector, K_AXIS, 0.5*M_PI + theta);
- StGermain_RotateCoordinateAxis( normal[1], eigenvectorList[0].vector, K_AXIS, 0.5*M_PI - theta);
- }
- else {
-
- /* Identify potential failure directions */
- StGermain_RotateVector( slip[0], eigenvectorList[0].vector, eigenvectorList[1].vector, + theta );
- StGermain_RotateVector( slip[1], eigenvectorList[0].vector, eigenvectorList[1].vector, - theta );
- StGermain_RotateVector( normal[0], eigenvectorList[0].vector, eigenvectorList[1].vector, 0.5*M_PI + theta);
- StGermain_RotateVector( normal[1], eigenvectorList[0].vector, eigenvectorList[1].vector, 0.5*M_PI - theta);
- }
-
- /* Resolve shear strain-rate for the potential failure planes */
- strainRate_ns[0] = fabs(TensorArray_MultiplyByVectors( velocityGradient, slip[0], normal[0], dim ));
- strainRate_ns[1] = fabs(TensorArray_MultiplyByVectors( velocityGradient, slip[1], normal[1], dim ));
-
- /* Choose the plane which is oriented favorably for continued slip */
- favourablePlane = strainRate_ns[0] > strainRate_ns[1] ? 0 : 1;
-
- memcpy( normalDirector, normal[favourablePlane], dim * sizeof(double) );
- memcpy( particleExt->slip, slip[favourablePlane], dim * sizeof(double) );
-
- self->storedSlipRateValue = strainRate_ns[ favourablePlane ];
-
- return normalDirector;
-}
-
-double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double effectiveCohesion;
-
- if ( self->tryingOldOrientation || self->ignoreOldOrientation ) {
- /* If this is the old orientation or old orientations are ignored
- * then the weakening should be considered here */
- double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- effectiveCohesion = self->cohesion * (1.0 - strainWeakeningRatio) +
- self->cohesionAfterSoftening * strainWeakeningRatio;
- }
- else {
- /* If old orientations are given priority and this is an unbroken direction
- * then consider material as completely pristine */
- effectiveCohesion = self->cohesion;
- }
-
- return effectiveCohesion;
-}
-
-double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- static double effectiveFrictionCoefficient = 0.0;
- static int prevTryingOldOrientationValue = -1;
- static void* prevParticle = NULL;
-
- /* See if we can return cached value */
- if ( self->tryingOldOrientation == (Bool) prevTryingOldOrientationValue && materialPoint == prevParticle )
- return effectiveFrictionCoefficient;
-
- /* This is a different situation to the cached case - therefore recalculate */
- if ( self->tryingOldOrientation || self->ignoreOldOrientation ) {
- /* If this is the old orientation or old orientations are ignored
- * then the weakening should be considered here */
- double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- effectiveFrictionCoefficient = self->frictionCoefficient * (1.0 - strainWeakeningRatio) +
- self->frictionCoefficientAfterSoftening * strainWeakeningRatio;
- }
- else {
- /* If old orientations are given priority and this is an unbroken direction
- * then consider material as completely pristine */
- effectiveFrictionCoefficient = self->frictionCoefficient;
- }
-
- /* Cache values */
- prevParticle = materialPoint;
- prevTryingOldOrientationValue = self->tryingOldOrientation;
-
- return effectiveFrictionCoefficient;
-}
-
-double _MohrCoulomb_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- double pressure = self->currentPressure;
- Eigenvector* eigenvectorList = self->currentEigenvectorList;
- double tau_nn;
- double sigma_nn;
-
- if ( self->tryingOldOrientation ) {
- /* tau_nn has already been calculated in _MohrCoulomb_OldOrientationStillSoftening */
- tau_nn = self->tau_nn;
- }
- else {
- double stressMin = eigenvectorList[0].eigenvalue;
- double stressMax = (dim == 2 ? eigenvectorList[1].eigenvalue : eigenvectorList[2].eigenvalue);
- double theta = 0.5 * atan( 1.0/ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint ) );
-
- tau_nn = 0.5 * (( stressMax + stressMin ) + cos( 2.0 * theta ) * ( stressMax - stressMin ));
- }
-
- sigma_nn = pressure - tau_nn;
-
- return sigma_nn;
-}
-
-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;
- MohrCoulomb_Particle* particleExt;
-
- particleExt = ExtensionManager_Get( materialPointsSwarm->particleExtensionMgr, materialPoint, self->particleExtHandle );
- particleExt->slipRate = 0.0;
-
- 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->currentStress );
-
- SymmetricTensor_CalcAllEigenvectors( self->currentStress, dim, self->currentEigenvectorList );
-
- Director_SetDontUpdateParticleFlag( self->director, materialPoint, False );
-}
-
-void _MohrCoulomb_UpdateDrawParameters( void* rheology ) {
- MohrCoulomb* self = (MohrCoulomb*) rheology;
- Particle_Index lParticle_I;
- Particle_Index particleLocalCount;
- StrainWeakening* strainWeakening = self->strainWeakening;
- StandardParticle* materialPoint;
- double slipRate;
-
- double length;
- double brightness;
- double opacity;
- double strainWeakeningRatio;
- double localMaxStrainIncrement;
- double localMaxSlipRate;
- double localMeanStrainIncrement;
- double localMeanSlipRate;
- Particle_Index localFailed;
-
- double globalMaxSlipRate;
- double globalMaxStrainIncrement;
- double globalMeanStrainIncrement;
- double globalMeanSlipRate;
- Particle_Index globalFailed;
-
- double averagedGlobalMaxSlipRate = 0.0;
- double averagedGlobalMaxStrainIncrement = 0.0;
-
- double oneOverGlobalMaxSlipRate;
- double oneOverGlobalMaxStrainIncrement;
- double postFailureWeakeningIncrement;
-
- /* This parameter is only needed for the alternative commented out set of parameters */
- /* double globalMaxStrainWeakeningRatio = 0.0; */
-
- /* 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 */
-
- localMaxStrainIncrement = 0.0;
- localMeanStrainIncrement = 0.0;
- localMaxSlipRate = 0.0;
- localMeanSlipRate = 0.0;
- localFailed = 0;
-
- /* Update all variables */
- Variable_Update( self->hasYieldedVariable->variable );
- Variable_Update( self->slipRate->variable );
- Variable_Update( self->brightness->variable );
- Variable_Update( self->opacity->variable );
- Variable_Update( self->length->variable );
- Variable_Update( self->thickness->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++;
- slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
-
- postFailureWeakeningIncrement =
- Variable_GetValueDouble( strainWeakening->postFailureWeakeningIncrement->variable, lParticle_I );
-
- localMeanSlipRate += slipRate;
- localMeanStrainIncrement += postFailureWeakeningIncrement;
-
- if(localMaxSlipRate < slipRate)
- localMaxSlipRate = slipRate;
-
- if(localMaxStrainIncrement < postFailureWeakeningIncrement)
- localMaxStrainIncrement = postFailureWeakeningIncrement;
- }
- }
-
- MPI_Allreduce( &localMaxSlipRate, &globalMaxSlipRate, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
- MPI_Allreduce( &localMaxStrainIncrement, &globalMaxStrainIncrement, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
- MPI_Allreduce( &localMeanSlipRate, &globalMeanSlipRate, 1, MPI_DOUBLE, MPI_SUM, 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;
- globalMeanSlipRate /= (double) globalFailed;
-
- averagedGlobalMaxStrainIncrement =
- 0.5 * averagedGlobalMaxStrainIncrement +
- 0.25 * globalMeanStrainIncrement +
- 0.25 * globalMaxStrainIncrement;
- averagedGlobalMaxSlipRate =
- 0.5 * averagedGlobalMaxSlipRate +
- 0.25 * globalMaxSlipRate +
- 0.25 * globalMeanSlipRate;
-
- /* Let's simply assume that twice the mean is a good place to truncate these values */
- oneOverGlobalMaxSlipRate = 1.0 / averagedGlobalMaxSlipRate;
- oneOverGlobalMaxStrainIncrement = 1.0 / averagedGlobalMaxStrainIncrement;
-
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
- 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->length->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
- continue;
- }
-
- slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
- strainWeakeningRatio = StrainWeakening_CalcRatio( strainWeakening, materialPoint );
-
-
- length = 0.005 + 0.005 * strainWeakeningRatio;
- length *= 10;
- brightness = strainWeakeningRatio * slipRate * oneOverGlobalMaxSlipRate;
- if( brightness > 1.0 )
- brightness = 1.0;
- opacity = 0.5 * pow(brightness,3.0);
-
- Variable_SetValueFloat( self->brightness->variable, lParticle_I, brightness );
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, opacity );
- Variable_SetValueFloat( self->length->variable, lParticle_I, (float) length );
- Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.5 + 10.0 * (float)length );
-
- /* This parameter is only needed for the alternative commented out set of parameters */
- /* if (strainWeakeningRatio>globalMaxStrainWeakeningRatio)
- globalMaxStrainWeakeningRatio = strainWeakeningRatio;
- */
- }
-
- /* This part is another set of parameters that has been used with shear problems.
- * it's kept here for now, but it shows that in principle we could define as many set of parameters
- * as we want (depending on the problem ?)*/
- /*
- for ( lParticle_I = 0 ; lParticle_I < particleLocalCount ; lParticle_I++ ) {
- materialPoint = Swarm_ParticleAt( strainWeakening->materialPointsSwarm, 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->length->variable, lParticle_I, 0.0 );
- Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.0 );
- continue;
- }
-
- slipRate = Variable_GetValueDouble( self->slipRate->variable, lParticle_I );
- strainWeakeningRatio = StrainWeakening_CalcRatio( strainWeakening, materialPoint );
-
- length = 0.03 + 0.03 * (strainWeakeningRatio/globalMaxStrainWeakeningRatio);
-
- brightness = 1.0-(strainWeakeningRatio/globalMaxStrainWeakeningRatio);
-
- if (brightness > 1.0)
- brightness = 1.0;
-
- opacity = (slipRate/globalMaxSlipRate);
-
- if (opacity > 0.90)
- opacity = 1.0;//this condition is to make sure we have enough planes that will be clearly seen.
-
- Variable_SetValueFloat( self->brightness->variable, lParticle_I, brightness );
- Variable_SetValueFloat( self->opacity->variable, lParticle_I, opacity );
- Variable_SetValueFloat( self->length->variable, lParticle_I, (float) length );
- Variable_SetValueFloat( self->thickness->variable, lParticle_I, 0.5 + 10.0 * (float)length );
- }
- */
-}
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2006-08-18 00:18:06 UTC (rev 4353)
@@ -1,178 +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;
-
- typedef struct {
- double slipRate;
- XYZ slip;
- float brightness;
- float opacity;
- float length;
- float thickness;
- Particle_Bool tensileFailure;
- } MohrCoulomb_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 __MohrCoulomb \
- /* Parent info */ \
- __YieldRheology \
- /* Virtual functions go here */ \
- /* Material Parameters */\
- SwarmVariable* slipRate; \
- SwarmVariable* slip; \
- SwarmVariable* brightness; \
- SwarmVariable* opacity; \
- SwarmVariable* length; \
- SwarmVariable* thickness; \
- SwarmVariable* tensileFailure; \
- 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 minimumYieldStress; \
- Bool ignoreOldOrientation; \
- /* Stored values that are calculated once for each particle */ \
- Eigenvector currentEigenvectorList[3]; \
- TensorArray currentVelocityGradient; \
- SymmetricTensor currentStress; \
- double currentPressure; \
- double tau_nn; \
- double storedSlipRateValue; \
- /* Flags set to tell functions what's going on */ \
- Bool tryingOldOrientation;
-
- 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 _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 );
-
- Bool _MohrCoulomb_OldOrientationStillSoftening( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) ;
- double* _MohrCoulomb_UpdateNormalDirection( void* rheology, MaterialPointsSwarm* materialPointsSwarm, void* materialPoint, Dimension_Index dim ) ;
-
- double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint ) ;
- double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint ) ;
- 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 ) ;
-
- void _MohrCoulomb_UpdateDrawParameters( void* rheology ) ;
-
-#endif
Deleted: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.meta 2006-08-18 00:18:06 UTC (rev 4353)
@@ -1,104 +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 a MohrCoulomb failure criterion</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>
- <struct>
- <param name="Name">ignoreOldOrientation</param>
- <param name="Type">Bool</param>
- <param name="Default">False</param>
- <param name="Description">Boolean to indicate if we want to take into account the oldOrientation</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">FiniteElementContext</param>
- <param name="Fallback Key">Context</param>
- <param name="Description">The context is essential because the MohrCoulomb_UpdateDrawParameters function is prepended to an existing Entry Point. </param>
- </struct>
-
- <struct>
- <param name="Essential">True</param>
- <param name="Type">Director</param>
- <param name="Description">This component allows to update the normal associated with a failure plane wrt the deformation. At the moment, it is only possible to update all normal at the same time, no matter if they correspond to active failure planes or not</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>
-
-
-</list>
-<!-- Add an example XML if possible -->
-<param name="Example">Underworld/InputFiles/ExtensionMC.xml</param>
-<param name="Example">Underworld/InputFiles/ExtensionMC3D.xml</param>
-
-</StGermainData>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Rheology.h 2006-08-18 00:18:06 UTC (rev 4353)
@@ -64,7 +64,7 @@
#include "VonMises.h"
#include "Byerlee.h"
#include "DruckerPrager.h"
- #include "MohrCoulomb.h"
+ #include "FaultingMoresiMuhlhaus2006.h"
#include "StrainWeakening.h"
#include "BuiterStrainWeakening.h"
#include "Director.h"
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h 2006-08-18 00:18:04 UTC (rev 4352)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/types.h 2006-08-18 00:18:06 UTC (rev 4353)
@@ -66,7 +66,7 @@
typedef struct Byerlee Byerlee;
typedef struct DruckerPrager DruckerPrager;
- typedef struct MohrCoulomb MohrCoulomb;
+ typedef struct FaultingMoresiMuhlhaus2006 FaultingMoresiMuhlhaus2006;
typedef struct StrainWeakening StrainWeakening;
typedef struct BuiterStrainWeakening BuiterStrainWeakening;
More information about the cig-commits
mailing list