[cig-commits] r14221 - in long/3D/Gale/trunk: . src/Underworld/Rheology/src
walter at geodynamics.org
walter at geodynamics.org
Wed Mar 4 09:09:32 PST 2009
Author: walter
Date: 2009-03-04 09:09:32 -0800 (Wed, 04 Mar 2009)
New Revision: 14221
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
Log:
r2523 at dante: boo | 2009-03-04 06:27:22 -0800
Make MohrCoulomb check if a particle is on the boundary by looking at which element it is in, rather than using a boundary shape
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2522
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2523
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-04 17:09:27 UTC (rev 14220)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.c 2009-03-04 17:09:32 UTC (rev 14221)
@@ -120,7 +120,12 @@
double boundaryCohesionAfterSoftening,
double boundaryFrictionCoefficient,
double boundaryFrictionCoefficientAfterSoftening,
- Stg_Shape *boundaryShape,
+ Bool boundaryBottom,
+ Bool boundaryTop,
+ Bool boundaryLeft,
+ Bool boundaryRight,
+ Bool boundaryFront,
+ Bool boundaryBack,
double minimumYieldStress,
HydrostaticTerm* hydrostaticTerm)
{
@@ -140,7 +145,12 @@
self->boundaryFrictionCoefficient = boundaryFrictionCoefficient;
self->boundaryCohesionAfterSoftening = boundaryCohesionAfterSoftening;
self->boundaryFrictionCoefficientAfterSoftening = boundaryFrictionCoefficientAfterSoftening;
- self->boundaryShape=boundaryShape;
+ self->boundaryBottom=boundaryBottom;
+ self->boundaryTop=boundaryTop;
+ self->boundaryLeft=boundaryLeft;
+ self->boundaryRight=boundaryRight;
+ self->boundaryFront=boundaryFront;
+ self->boundaryBack=boundaryBack;
self->minimumYieldStress = minimumYieldStress;
self->hydrostaticTerm=hydrostaticTerm;
@@ -203,7 +213,12 @@
Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryCohesionAfterSoftening", 0.0 ),
Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryFrictionCoefficient", 0.0 ),
Stg_ComponentFactory_GetDouble( cf, self->name, "boundaryFrictionCoefficientAfterSoftening", 0.0 ),
- Stg_ComponentFactory_ConstructByKey( cf, self->name, "boundaryShape", Stg_Shape, False, data ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryBottom", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryTop", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryLeft", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryRight", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryFront", False ),
+ Stg_ComponentFactory_GetBool( cf, self->name, "boundaryBack", False ),
Stg_ComponentFactory_GetDouble( cf, self->name, "minimumYieldStress", 0.0 ),
Stg_ComponentFactory_ConstructByKey( cf, self->name,
"HydrostaticTerm", HydrostaticTerm, False, data ) );
@@ -310,6 +325,10 @@
double factor;
Cell_Index cell_I;
Coord coord;
+ Element_GlobalIndex element_gI = 0;
+ unsigned inds[3];
+ Grid* elGrid;
+ Bool inside_boundary;
FeVariable_InterpolateWithinElement( pressureField, lElement_I, xi, &pressure );
@@ -324,25 +343,33 @@
/* Calculate frictional strength. We modify the friction and
cohesion because we have grouped terms from the normal
stresses and moved it to the yield indicator. */
+
+
+ /* Big song and dance to see if we are at a boundary that we care about */
+ elGrid = *(Grid**)ExtensionManager_Get(pressureField->feMesh->info,
+ pressureField->feMesh,
+ ExtensionManager_GetHandle(pressureField->feMesh->info, "elementGrid" ) );
+
+ element_gI = FeMesh_ElementDomainToGlobal( pressureField->feMesh, lElement_I );
+ RegularMeshUtils_Element_1DTo3D( pressureField->feMesh, element_gI, inds );
+
+ inside_boundary=(self->boundaryBottom && inds[1]==0)
+ || (self->boundaryTop && inds[1]==elGrid->sizes[1]-1)
+ || (self->boundaryLeft && inds[0]==0)
+ || (self->boundaryRight && inds[0]==elGrid->sizes[0]-1)
+ || (dim==3 && ((self->boundaryBack && inds[2]==0)
+ || (self->boundaryFront && inds[2]==elGrid->sizes[2]-1)));
+
effectiveFrictionCoefficient =
- _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint, coord );
- effectiveCohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint,coord);
+ _MohrCoulomb_EffectiveFrictionCoefficient( self, materialPoint,
+ inside_boundary );
+ effectiveCohesion = _MohrCoulomb_EffectiveCohesion(self,materialPoint,
+ inside_boundary);
- if(self->boundaryShape
- && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
- {
- factor=1;
-/* printf("strength %g %g %g %g %g %g %g\n",coord[0],coord[1], */
-/* frictionalStrength, */
-/* factor,effectiveFrictionCoefficient,pressure, */
-/* effectiveCohesion); */
- }
- else
- {
- double theta;
- theta = 0.5 * atan( 1.0/ effectiveFrictionCoefficient);
- factor=sin(2*theta);
- }
+ /* effectiveFrictionCoefficient=tan(phi).
+ If factor=sin(atan(1/tan(phi)))
+ => factor=cos(alpha)=1/sqrt(1+tan(phi)**2) */
+ factor=1/sqrt(1 + effectiveFrictionCoefficient*effectiveFrictionCoefficient);
if(dim==2)
{
@@ -397,14 +424,13 @@
}
double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint,
- double *coord ) {
+ Bool inside_boundary ) {
MohrCoulomb* self = (MohrCoulomb*) rheology;
double effectiveCohesion;
double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
-
- if(self->boundaryShape
- && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
+
+ if(inside_boundary)
{
effectiveCohesion=self->boundaryCohesion*(1.0-strainWeakeningRatio)+
self->boundaryCohesionAfterSoftening*strainWeakeningRatio;
@@ -418,14 +444,13 @@
return effectiveCohesion;
}
-double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, double *coord ) {
+double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, Bool inside_boundary ) {
MohrCoulomb* self = (MohrCoulomb*) rheology;
double effectiveFrictionCoefficient = 0.0;
double strainWeakeningRatio = StrainWeakening_CalcRatio( self->strainWeakening, materialPoint );
- if(self->boundaryShape
- && Stg_Shape_IsCoordInside(self->boundaryShape,coord))
+ if(inside_boundary)
{
effectiveFrictionCoefficient =
self->boundaryFrictionCoefficient * (1.0 - strainWeakeningRatio) +
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-04 17:09:27 UTC (rev 14220)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/MohrCoulomb.h 2009-03-04 17:09:32 UTC (rev 14221)
@@ -70,7 +70,12 @@
double boundaryCohesionAfterSoftening; \
double boundaryFrictionCoefficient; \
double boundaryFrictionCoefficientAfterSoftening; \
- Stg_Shape* boundaryShape; \
+ Bool boundaryTop; \
+ Bool boundaryBottom; \
+ Bool boundaryLeft; \
+ Bool boundaryRight; \
+ Bool boundaryFront; \
+ Bool boundaryBack; \
double minimumYieldStress; \
/* Stored values that are calculated once for each particle */ \
Eigenvector currentEigenvectorList[3]; \
@@ -145,8 +150,8 @@
double yieldCriterion,
double yieldIndicator );
- double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint, double *coord ) ;
- double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, double *coord ) ;
+ double _MohrCoulomb_EffectiveCohesion( void* rheology, void* materialPoint, Bool inside_boundary ) ;
+ double _MohrCoulomb_EffectiveFrictionCoefficient( void* rheology, void* materialPoint, Bool inside_boundary ) ;
double _MohrCoulomb_Sigma_nn( void* rheology, void* materialPoint, Dimension_Index dim ) ;
void _MohrCoulomb_StoreCurrentParameters(
More information about the CIG-COMMITS
mailing list