[cig-commits] r7644 - in long/3D/Gale/trunk: . input input/cookbook
src/StGermain/Base/Automation/src
src/StGermain/Discretisation/Utils/src
src/StgFEM/Discretisation/src src/StgFEM/SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Wed Jul 11 14:02:09 PDT 2007
Author: walter
Date: 2007-07-11 14:02:00 -0700 (Wed, 11 Jul 2007)
New Revision: 7644
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/input/cookbook/viscous_extension.xml
long/3D/Gale/trunk/input/extension.xml
long/3D/Gale/trunk/src/StGermain/Base/Automation/src/VariableCondition.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.h
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
Log:
r1860 at earth: boo | 2007-07-08 02:19:36 -0700
Rough version of static friction. It only works on the bottom and the friction coefficient is hard coded
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1858
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1860
Modified: long/3D/Gale/trunk/input/cookbook/viscous_extension.xml
===================================================================
--- long/3D/Gale/trunk/input/cookbook/viscous_extension.xml 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/input/cookbook/viscous_extension.xml 2007-07-11 21:02:00 UTC (rev 7644)
@@ -125,6 +125,17 @@
<param name="DofLayout">pressureDofLayout</param>
<param name="LinkedDofInfo">pressureLinkedDofs</param>
</struct>
+
+ <struct name="StressField">
+ <param name="Type">StressField</param>
+ <param name="StrainRateField">StrainRateField</param>
+ <param name="Context">context</param>
+ <param name="ConstitutiveMatrix">constitutiveMatrix</param>
+ <param name="Swarm">picIntegrationPoints</param>
+ <param name="Mesh">mesh-linear</param>
+ <param name="IC">stressICs</param>
+ </struct>
+
<struct name="cellLayout">
<param name="Type">SingleCellLayout</param>
</struct>
@@ -474,5 +485,23 @@
</struct>
</list>
</struct>
+
+
+ <struct name="stressICs">
+ <param name="type">CompositeVC</param>
+ <list name="vcList">
+ <struct>
+ <param name="type">AllNodesVC</param>
+ <list name="variables">
+ <struct>
+ <param name="name">stress</param>
+ <param name="type">double</param>
+ <param name="value">0.0</param>
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
<param name="checkpointEvery">1</param>
</StGermainData>
Modified: long/3D/Gale/trunk/input/extension.xml
===================================================================
--- long/3D/Gale/trunk/input/extension.xml 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/input/extension.xml 2007-07-11 21:02:00 UTC (rev 7644)
@@ -122,6 +122,17 @@
<param name="DofLayout">pressureDofLayout</param>
<param name="LinkedDofInfo">pressureLinkedDofs</param>
</struct>
+
+ <struct name="StressField">
+ <param name="Type">StressField</param>
+ <param name="StrainRateField">StrainRateField</param>
+ <param name="Context">context</param>
+ <param name="ConstitutiveMatrix">constitutiveMatrix</param>
+ <param name="Swarm">picIntegrationPoints</param>
+ <param name="Mesh">mesh-linear</param>
+ <param name="IC">stressICs</param>
+ </struct>
+
<struct name="cellLayout">
<param name="Type">SingleCellLayout</param>
</struct>
@@ -288,6 +299,12 @@
<param name="nonLinearTolerance">nonLinearTolerance</param>
<param name="makeConvergenceFile">false</param>
</struct>
+ <struct name="buoyancyForceTerm">
+ <param name="Type">BuoyancyForceTerm</param>
+ <param name="ForceVector">mom_force</param>
+ <param name="Swarm">picIntegrationPoints</param>
+ <param name="gravity">1.0</param>
+ </struct>
<struct name="backgroundShape">
<param name="Type">Everywhere</param>
</struct>
@@ -349,6 +366,7 @@
<struct name="crust">
<param name="Type">RheologyMaterial</param>
<param name="Shape">crustShape</param>
+ <param name="density">100.0</param>
<list name="Rheology">
<param>crustViscosity</param>
<param>yielding</param>
@@ -359,6 +377,7 @@
<struct name="pdms">
<param name="Type">RheologyMaterial</param>
<param name="Shape">pdmsShape</param>
+ <param name="density">100.0</param>
<list name="Rheology">
<param>pdmsViscosity</param>
<param>storeViscosity</param>
@@ -517,6 +536,7 @@
</struct>
</list>
</struct>
+
<struct>
<param name="type">WallVC</param>
<param name="wall">bottom</param>
@@ -526,6 +546,13 @@
<param name="type">double</param>
<param name="value">0</param>
</struct>
+ </list>
+ </struct>
+
+ <struct>
+ <param name="type">FrictionVC</param>
+ <param name="wall">bottom</param>
+ <list name="variables">
<struct>
<param name="name">vx</param>
<param name="type">func</param>
@@ -535,8 +562,25 @@
</struct>
</list>
</struct>
- <param name="StepFunctionLowerOffset">0.8</param>
- <param name="StepFunctionUpperOffset">1.2</param>
+
+ <struct name="stressICs">
+ <param name="type">CompositeVC</param>
+ <list name="vcList">
+ <struct>
+ <param name="type">AllNodesVC</param>
+ <list name="variables">
+ <struct>
+ <param name="name">stress</param>
+ <param name="type">double</param>
+ <param name="value">0.0</param>
+ </struct>
+ </list>
+ </struct>
+ </list>
+ </struct>
+
+ <param name="StepFunctionLowerOffset">0.99</param>
+ <param name="StepFunctionUpperOffset">1.01</param>
<param name="StepFunctionValue">1.0</param>
<param name="StepFunctionDim">0</param>
<param name="StepFunctionLessThan">False</param>
Modified: long/3D/Gale/trunk/src/StGermain/Base/Automation/src/VariableCondition.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Base/Automation/src/VariableCondition.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StGermain/Base/Automation/src/VariableCondition.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -387,7 +387,8 @@
/* Force the building of the variable (to be safe) */
var = self->variable_Register->_variable[self->vcTbl[i][vcVar_I].varIndex];
- Build( var, data, False );
+/* Build( var, data, False ); */
+ Build( var, data, True );
}
}
}
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -32,6 +32,8 @@
#include "types.h"
#include "FrictionVC.h"
#include "RegularMeshUtils.h"
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
#include <string.h>
#include <assert.h>
@@ -262,6 +264,7 @@
self->_wall = FrictionVC_Wall_Size;
self->_entryTbl = 0;
self->_entryCount = 0;
+ self->context=NULL;
}
@@ -309,7 +312,6 @@
}
/* Obtain the variable entries */
- self->_entryCount = 0;
self->_entryCount = Dictionary_Entry_Value_GetCount(Dictionary_Entry_Value_GetMember(vcDictVal, "variables"));
self->_entryTbl = Memory_Alloc_Array( FrictionVC_Entry, self->_entryCount, "FrictionVC->_entryTbl" );
varsVal = Dictionary_Entry_Value_GetMember(vcDictVal, "variables");
@@ -403,6 +405,7 @@
self->_wall = FrictionVC_Wall_Size;
self->_entryCount = 0;
self->_entryTbl = NULL;
+ self->context=NULL;
}
}
@@ -518,6 +521,7 @@
newFrictionVC->_dictionaryEntryName = self->_dictionaryEntryName;
newFrictionVC->_wall = self->_wall;
newFrictionVC->_entryCount = self->_entryCount;
+ newFrictionVC->context = self->context;
if( deep ) {
newFrictionVC->_mesh = (Mesh*)Stg_Class_Copy( self->_mesh, NULL, deep, nameExt, map );
@@ -566,7 +570,8 @@
void _FrictionVC_BuildSelf( void* wallVC, void* data ) {
FrictionVC* self = (FrictionVC*)wallVC;
-
+
+ self->context=(FiniteElementContext *)data;
if( self->_mesh ) {
Build( self->_mesh, data, False );
}
@@ -579,44 +584,41 @@
IndexSet *set = NULL;
Stream* warningStr = Journal_Register( Error_Type, self->type );
unsigned nDims;
- unsigned* gSize;
+ Grid* vertGrid;
+ FiniteElementContext * context = self->context;
+ FeVariable* deviatoric_stress = NULL;
+ FeVariable* pressure = NULL;
+ deviatoric_stress = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "StressField" );
+ pressure = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "PressureField" );
+
nDims = Mesh_GetDimSize( self->_mesh );
- gSize = (unsigned*)ExtensionManager_Get( self->_mesh->info, self->_mesh,
- ExtensionManager_GetHandle( self->_mesh->info,
- "cartesianGlobalSize" ) );
+ vertGrid = *(Grid**)ExtensionManager_Get( self->_mesh->info, self->_mesh,
+ ExtensionManager_GetHandle( self->_mesh->info,
+ "vertexGrid" ) );
switch (self->_wall) {
case FrictionVC_Wall_Front:
- if ( nDims < 3 || !gSize[2] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "K" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if ( nDims < 3 || !vertGrid->sizes[2] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
- set = RegularMeshUtils_CreateGlobalFrontSet(self->_mesh);
+ set = RegularMeshUtils_CreateGlobalFrontSet( self->_mesh );
}
break;
case FrictionVC_Wall_Back:
- if ( nDims < 3 || !gSize[2] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "K" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if ( nDims < 3 || !vertGrid->sizes[2] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
- set = RegularMeshUtils_CreateGlobalBackSet(self->_mesh);
+ set = RegularMeshUtils_CreateGlobalBackSet( self->_mesh );
}
break;
case FrictionVC_Wall_Top:
- if ( nDims < 2 || !gSize[1] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "J" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if ( nDims < 2 || !vertGrid->sizes[1] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
set = RegularMeshUtils_CreateGlobalTopSet(self->_mesh);
@@ -624,23 +626,44 @@
break;
case FrictionVC_Wall_Bottom:
- if ( nDims < 2 || !gSize[1] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "J" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if ( nDims < 2 || !vertGrid->sizes[1] ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
- set = RegularMeshUtils_CreateGlobalBottomSet(self->_mesh);
- }
+ if(deviatoric_stress!=NULL && deviatoric_stress->dofLayout!=NULL)
+ {
+ Node_LocalIndex n_i;
+ unsigned nNodes;
+ IJK ijk;
+
+ nNodes = Mesh_GetDomainSize( self->_mesh, MT_VERTEX );
+ set = IndexSet_New( nNodes );
+ for( n_i = 0; n_i < nNodes; n_i++ ) {
+ RegularMeshUtils_Node_1DTo3D( self->_mesh, Mesh_DomainToGlobal( self->_mesh, MT_VERTEX, n_i ), ijk );
+ if(ijk[1]==0)
+ {
+ double str[6], p;
+ double friction;
+ friction=10.0;
+ FeVariable_GetValueAtNode(deviatoric_stress,n_i,str);
+ FeVariable_GetValueAtNode(pressure,n_i,&p);
+
+ /* For now, assume that the bottom is flat and not inclined */
+ if(str[1]==0 || p-str[2]>friction*fabs(str[1]))
+ IndexSet_Add(set,n_i);
+ }
+ }
+ }
+ else
+ {
+ set = RegularMeshUtils_CreateGlobalBottomSet(self->_mesh);
+ }
+ }
break;
case FrictionVC_Wall_Left:
- if ( !gSize[0] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "I" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if ( nDims < 1 ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
set = RegularMeshUtils_CreateGlobalLeftSet(self->_mesh);
@@ -648,11 +671,8 @@
break;
case FrictionVC_Wall_Right:
- if ( !gSize[0] ) {
- Journal_Printf( warningStr, "Warning - in %s: Can't build a %s wall VC "
- "when mesh has no elements in the %s axis. Returning an empty set.\n", __func__,
- FrictionVC_WallEnumToStr[self->_wall], "I" );
- set = IndexSet_New(Mesh_GetDomainSize( self->_mesh, MT_VERTEX ));
+ if( nDims < 1 ) {
+ set = IndexSet_New( Mesh_GetDomainSize( self->_mesh, MT_VERTEX ) );
}
else {
set = RegularMeshUtils_CreateGlobalRightSet(self->_mesh);
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.h 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Utils/src/FrictionVC.h 2007-07-11 21:02:00 UTC (rev 7644)
@@ -56,7 +56,11 @@
FrictionVC_Wall _wall; \
FrictionVC_Entry_Index _entryCount; \
FrictionVC_Entry* _entryTbl; \
- Mesh* _mesh;
+ Mesh* _mesh; \
+ /* I would like to make this a FiniteElementContext*, but */ \
+ /* then there are problems compiling because this is in */ \
+ /* StGermain, and we do not have access to StgFEM yet. */ \
+ void* context;
struct _FrictionVC { __FrictionVC };
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -484,10 +484,10 @@
void FeEquationNumber_Build( void* feEquationNumber ) {
FeEquationNumber* self = (FeEquationNumber*)feEquationNumber;
- if ( False == self->isBuilt ) {
+/* if ( False == self->isBuilt ) { */
self->_build( self, NULL );
self->isBuilt = True;
- }
+/* } */
}
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -513,17 +513,16 @@
FeVariable* self = (FeVariable*)variable;
DiscretisationContext* context = (DiscretisationContext*)data;
- if ( False == self->isBuilt ) {
+/* if ( False == self->isBuilt ) { */
self->isBuilt = True;
Journal_DPrintf( self->debug, "In %s- for %s:\n", __func__, self->name );
Stream_IndentBranch( StgFEM_Debug );
- /* build the BCs */
Build( self->feMesh, data, False );
Build( self->dofLayout, data, False );
if ( self->bcs ){
- Build( self->bcs, data, False );
+ Build( self->bcs, data, True );
}
/* only bother building the ics specified via XML/construct if we are not in restart mode
- otherwise, we will use the checkpointed values anyway */
@@ -541,7 +540,7 @@
FeEquationNumber_Build( self->eqNum );
Stream_UnIndentBranch( StgFEM_Debug );
- }
+/* } */
}
void _FeVariable_Construct( void* variable, Stg_ComponentFactory* cf, void* data )
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -516,7 +516,8 @@
/* ensure variables are built */
if( self->rowVariable )
- Build( self->rowVariable, data, False );
+ Build( self->rowVariable, data, True );
+/* Build( self->rowVariable, data, False ); */
/* If we don't have a communicator, grab one off the mesh. */
if( !self->comm ) {
@@ -526,7 +527,8 @@
}
if( self->columnVariable )
- Build( self->columnVariable, data, False );
+ Build( self->columnVariable, data, True );
+/* Build( self->columnVariable, data, False ); */
#if DEBUG
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2007-07-11 20:17:32 UTC (rev 7643)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2007-07-11 21:02:00 UTC (rev 7644)
@@ -403,23 +403,27 @@
/* build the matrices */
for ( index = 0; index < self->stiffnessMatrices->count; index++ ) {
/* Update rowSize and colSize if boundary conditions have been applied */
- Build( self->stiffnessMatrices->data[index], _context, False );
+/* Build( self->stiffnessMatrices->data[index], _context, False ); */
+ Build( self->stiffnessMatrices->data[index], _context, True );
}
/* and the vectors */
for ( index = 0; index < self->forceVectors->count; index++ ) {
/* Build the force vectors - includes updateing matrix size based on Dofs */
- Build( self->forceVectors->data[index], _context, False );
+/* Build( self->forceVectors->data[index], _context, False ); */
+ Build( self->forceVectors->data[index], _context, True );
}
/* and the solutions */
for ( index = 0; index < self->solutionVectors->count; index++ ) {
/* Build the force vectors - includes updateing matrix size based on Dofs */
- Build( self->solutionVectors->data[index], _context, False );
+/* Build( self->solutionVectors->data[index], _context, False ); */
+ Build( self->solutionVectors->data[index], _context, True );
}
/* lastly, the solver */
- Build( self->solver, self, False );
+/* Build( self->solver, self, False ); */
+ Build( self->solver, self, True );
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -652,6 +656,7 @@
Stream* errorStream = Journal_Register( Error_Type, self->type );
double wallTime;
Iteration_Index minIterations = self->nonLinearMinIterations;
+ FeVariable* stress = NULL;
Journal_Printf( self->info, "In %s\n", __func__ );
Stream_IndentBranch( StgFEM_Debug );
@@ -662,6 +667,11 @@
self->nonLinearIteration_I = 0;
Journal_Printf(self->info,"\nNon linear solver - iteration %d\n", self->nonLinearIteration_I);
+ /* Add this to create a new set of bcs dependent on the stress used in friction bc's. */
+ stress = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "StressField" );
+ ParticleFeVariable_Update( stress );
+ _SystemLinearEquations_Build(sle,_context);
+
self->linearExecute( self, _context );
self->hasExecuted = True;
@@ -671,6 +681,12 @@
Vector_SetLocalSize( previousVector, Vector_GetLocalSize( currentVector ) );
for ( self->nonLinearIteration_I = 1 ; self->nonLinearIteration_I < maxIterations ; self->nonLinearIteration_I++ ) {
+
+ /* Add this to create a new set of bcs dependent on the stress used in friction bc's. */
+ ParticleFeVariable_Update( stress );
+ _SystemLinearEquations_Build(sle,_context);
+
+
Vector_CopyEntries( currentVector, previousVector );
Journal_Printf(self->info,"Non linear solver - iteration %d\n", self->nonLinearIteration_I);
More information about the cig-commits
mailing list