[cig-commits] commit: Add StressBC
Mercurial
hg at geodynamics.org
Mon Dec 21 13:28:42 PST 2009
changeset: 689:d058abdcb8f2
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Mon Dec 21 13:27:10 2009 -0800
files: Utils/src/Init.c Utils/src/RadiogenicHeatingTerm.c Utils/src/StressBC.c Utils/src/StressBC.h Utils/src/StressBC.meta Utils/src/Utils.h Utils/src/Wall.c Utils/src/types.h
description:
Add StressBC
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/Init.c
--- a/Utils/src/Init.c Mon Dec 07 13:47:23 2009 -0800
+++ b/Utils/src/Init.c Mon Dec 21 13:27:10 2009 -0800
@@ -74,6 +74,8 @@ Bool Underworld_Utils_Init( int* argc, c
Stg_ComponentRegister_Add( componentRegister, DensityField_Type , "0", _DensityField_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, MixedStabiliserTerm_Type,
"0", _MixedStabiliserTerm_DefaultNew );
+ Stg_ComponentRegister_Add( componentRegister, StressBC_Type,
+ "0", _StressBC_DefaultNew );
RegisterParent( UnderworldContext_Type, PICelleratorContext_Type );
RegisterParent( PressureTemperatureOutput_Type, SwarmOutput_Type );
@@ -85,6 +87,7 @@ Bool Underworld_Utils_Init( int* argc, c
RegisterParent( ViscosityField_Type, ParticleFeVariable_Type );
RegisterParent( DensityField_Type, ParticleFeVariable_Type );
RegisterParent( MixedStabiliserTerm_Type, StiffnessMatrixTerm_Type );
+ RegisterParent( StressBC_Type, ForceTerm_Type );
return True;
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/RadiogenicHeatingTerm.c
--- a/Utils/src/RadiogenicHeatingTerm.c Mon Dec 07 13:47:23 2009 -0800
+++ b/Utils/src/RadiogenicHeatingTerm.c Mon Dec 21 13:27:10 2009 -0800
@@ -7,6 +7,7 @@
#include <PICellerator/PopulationControl/PopulationControl.h>
#include <PICellerator/Weights/Weights.h>
#include <PICellerator/MaterialPoints/MaterialPoints.h>
+#include <PICellerator/Utils/Utils.h>
#include "types.h"
#include "RadiogenicHeatingTerm.h"
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/StressBC.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/StressBC.c Mon Dec 21 13:27:10 2009 -0800
@@ -0,0 +1,615 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Modified 2006 by Walter Landry (California Institute of Technology)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id$
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgDomain/StgDomain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
+#include <PICellerator/PopulationControl/PopulationControl.h>
+#include <PICellerator/Weights/Weights.h>
+#include <PICellerator/MaterialPoints/MaterialPoints.h>
+
+#include "PICellerator/Utils/types.h"
+#include "PICellerator/Utils/MaterialSwarmVariable.h"
+
+#include "types.h"
+#include "StressBC.h"
+
+#include <assert.h>
+#include <string.h>
+
+/* Textual name of this class */
+const Type StressBC_Type = "StressBC";
+
+extern const char* WallEnumToStr[Wall_Size];
+
+StressBC* StressBC_New(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ StressBC* self = (StressBC*) _StressBC_DefaultNew( name );
+
+ StressBC_InitAll(
+ self,
+ forceVector,
+ integrationSwarm,
+ conFunc_Register,
+ context);
+
+ return self;
+}
+
+/* Creation implementation / Virtual constructor */
+StressBC* _StressBC_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,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context)
+{
+ StressBC* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(StressBC) );
+ self = (StressBC*) _ForceTerm_New(
+ sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _assembleElement,
+ name );
+
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+
+ return self;
+}
+
+void _StressBC_Init(
+ StressBC* self,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ self->numEntries = 0;
+ self->conFunc_Register=conFunc_Register;
+ self->context=context;
+}
+
+void StressBC_InitAll(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context)
+{
+ StressBC* self = (StressBC*) forceTerm;
+
+ ForceTerm_InitAll( self, forceVector, integrationSwarm, NULL );
+ _StressBC_Init( self, conFunc_Register, context );
+}
+
+void _StressBC_Delete( void* forceTerm ) {
+ StressBC* self = (StressBC*)forceTerm;
+
+ _ForceTerm_Delete( self );
+}
+
+void _StressBC_Print( void* forceTerm, Stream* stream ) {
+ StressBC* self = (StressBC*)forceTerm;
+
+ _ForceTerm_Print( self, stream );
+}
+
+void* _StressBC_DefaultNew( Name name ) {
+ return (void*)_StressBC_New(
+ sizeof(StressBC),
+ StressBC_Type,
+ _StressBC_Delete,
+ _StressBC_Print,
+ NULL,
+ _StressBC_DefaultNew,
+ _StressBC_Construct,
+ _StressBC_Build,
+ _StressBC_Initialise,
+ _StressBC_Execute,
+ _StressBC_Destroy,
+ _StressBC_AssembleElement,
+ NULL,
+ name,
+ NULL);
+}
+
+void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ FiniteElementContext* context;
+
+ /* Construct Parent */
+ _ForceTerm_Construct( self, cf, data );
+
+ context = (FiniteElementContext*)Stg_ComponentFactory_ConstructByName
+ ( cf, "context", FiniteElementContext, True, data ) ;
+
+ _StressBC_Init( self, context->condFunc_Register, context );
+ {
+ char* wallStr;
+
+ /* Obtain which wall */
+ wallStr = Stg_ComponentFactory_GetString( cf, self->name,
+ "wall", "");
+ if (!strcasecmp(wallStr, "back"))
+ self->_wall = Wall_Back;
+ else if (!strcasecmp(wallStr, "left"))
+ self->_wall = Wall_Left;
+ else if (!strcasecmp(wallStr, "bottom"))
+ self->_wall = Wall_Bottom;
+ else if (!strcasecmp(wallStr, "right"))
+ self->_wall = Wall_Right;
+ else if (!strcasecmp(wallStr, "top"))
+ self->_wall = Wall_Top;
+ else if (!strcasecmp(wallStr, "front"))
+ self->_wall = Wall_Front;
+ else {
+ assert( 0 );
+ self->_wall = Wall_Size; /* invalid entry */
+ }
+ }
+ _StressBC_GetValues(cf,self,"x",data);
+ _StressBC_GetValues(cf,self,"y",data);
+ _StressBC_GetValues(cf,self,"z",data);
+
+}
+
+/* Gets the actual values used by the StressBC (e.g. a float or a function). */
+void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC,
+ char *direction, void *data)
+{
+ StressBC* self = (StressBC*)stressBC;
+ char temp_str[20];
+ char *type;
+
+ switch(*direction)
+ {
+ case 'x':
+ self->_entryTbl[self->numEntries].axis=I_AXIS;
+ break;
+ case 'y':
+ self->_entryTbl[self->numEntries].axis=J_AXIS;
+ break;
+ case 'z':
+ self->_entryTbl[self->numEntries].axis=K_AXIS;
+ break;
+ }
+
+ strcat(strcpy(temp_str,direction),"_type");
+ type=Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ strcat(strcpy(temp_str,direction),"_value");
+
+ if(!strcasecmp(type,"double") || !strcasecmp(type,"float"))
+ {
+ self->_entryTbl[self->numEntries].type = StressBC_Double;
+ self->_entryTbl[self->numEntries].DoubleValue =
+ Stg_ComponentFactory_GetDouble( cf, self->name, temp_str, 0.0);
+ (self->numEntries)++;
+ }
+ else if(!strcasecmp(type,"func"))
+ {
+ char *funcName =
+ Stg_ComponentFactory_GetString( cf, self->name, temp_str, "");
+ Index cfIndex;
+
+ cfIndex = ConditionFunction_Register_GetIndex
+ ( self->conFunc_Register, funcName);
+ self->_entryTbl[self->numEntries].type =
+ StressBC_ConditionFunction;
+ if ( cfIndex == (unsigned)-1 ) {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of StressBC (applies to wall \"%s\"), the cond. func. applied to "
+ "direction \"%s\" - \"%s\" - wasn't found in the c.f. register.\n",
+ __func__, WallEnumToStr[self->_wall],
+ direction, funcName );
+ Journal_Printf( errorStr, "(Available functions in the C.F. register are: ");
+ ConditionFunction_Register_PrintNameOfEachFunc
+ ( self->conFunc_Register, errorStr );
+ Journal_Printf( errorStr, ")\n");
+ assert(0);
+ }
+ self->_entryTbl[self->numEntries].CFIndex = cfIndex;
+ (self->numEntries)++;
+ }
+ else if(!strcasecmp(type,"HydrostaticTerm"))
+ {
+ self->_entryTbl[self->numEntries].type = StressBC_HydrostaticTerm;
+ self->_entryTbl[self->numEntries].hydrostaticTerm =
+ Stg_ComponentFactory_ConstructByKey( cf, self->name, temp_str,
+ HydrostaticTerm, True,
+ data ) ;
+ (self->numEntries)++;
+ }
+ else if(strlen(type)!=0)
+ {
+ Stream* errorStr = Journal_Register( Error_Type, self->type );
+ Journal_Printf( errorStr, "Error- in %s: While parsing "
+ "definition of StressBC (applies to wall \"%s\"), the type of condition \"%s\"\nis not supported. Supported types are \"double\" and \"function\".",
+ __func__, WallEnumToStr[self->_wall],
+ type );
+ assert(0);
+ }
+}
+
+void _StressBC_Build( void* forceTerm, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ Name name;
+
+ _ForceTerm_Build( self, data );
+}
+
+void _StressBC_Initialise( void* forceTerm, void* data ) {
+ StressBC* self = (StressBC*)forceTerm;
+ Index i;
+
+ _ForceTerm_Initialise( self, data );
+
+}
+
+void _StressBC_Execute( void* forceTerm, void* data ) {
+ _ForceTerm_Execute( forceTerm, data );
+}
+
+void _StressBC_Destroy( void* forceTerm, void* data ) {
+ _ForceTerm_Destroy( forceTerm, data );
+}
+
+void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
+ StressBC* self = (StressBC*) forceTerm;
+ Dimension_Index dim = forceVector->dim;
+ FeMesh* mesh = forceVector->feVariable->feMesh;
+ Grid* grid;
+ Node_ElementLocalIndex eNode_I;
+ ElementType* elementType;
+ Dof_Index nodeDofCount;
+ double stress, area;
+ IJK ijk;
+ int overcount;
+ IArray *incidence;
+ int elementNodeCount;
+ int *elementNodes;
+
+ FeVariable* velocityField = NULL;
+ velocityField = (FeVariable*)FieldVariable_Register_GetByName
+ ( self->context->fieldVariable_Register, "VelocityField" );
+
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ nodeDofCount = dim;
+
+ grid=*(Grid**)ExtensionManager_Get(mesh->info, mesh,
+ ExtensionManager_GetHandle(mesh->info,
+ "vertexGrid"));
+
+ incidence=IArray_New();
+ area=StressBC_compute_face_area(self->_wall,mesh,lElement_I,dim,incidence);
+ elementNodeCount=IArray_GetSize(incidence);
+ elementNodes=IArray_GetPtr(incidence);
+
+ if(dim==2)
+ {
+ overcount=2;
+ }
+ else
+ {
+ overcount=4;
+ }
+ /* Apply the stress */
+ for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
+ /* Make sure that we are on the boundary */
+ int condition, entry_I;
+ ConditionFunction* cf;
+ RegularMeshUtils_Node_1DTo3D
+ (mesh,Mesh_DomainToGlobal(mesh,MT_VERTEX,elementNodes[eNode_I]),ijk);
+
+ switch(self->_wall)
+ {
+ case Wall_Left:
+ condition=(ijk[0] == 0);
+ break;
+ case Wall_Right:
+ condition=(ijk[0] == ( grid->sizes[0] - 1 ));
+ break;
+ case Wall_Bottom:
+ condition=(ijk[1] == 0);
+ break;
+ case Wall_Top:
+ condition=(ijk[1] == ( grid->sizes[1] - 1 ));
+ break;
+ case Wall_Front:
+ condition=(ijk[2] == 0);
+ break;
+ case Wall_Back:
+ condition=(ijk[2] == ( grid->sizes[2] - 1 ));
+ break;
+ }
+
+ if(condition)
+ {
+ for(entry_I=0; entry_I<self->numEntries; ++entry_I)
+ {
+ switch(self->_entryTbl[entry_I].type)
+ {
+ case StressBC_Double:
+ stress=self->_entryTbl[entry_I].DoubleValue;
+ break;
+ case StressBC_ConditionFunction:
+ cf=self->conFunc_Register->
+ _cf[self->_entryTbl[entry_I].CFIndex];
+
+ /* We use a variable number of zero "0", because
+ we don't use the variable number and that one
+ is always going to exist. */
+ ConditionFunction_Apply(cf,elementNodes[eNode_I],
+ 0,self->context,&stress);
+ break;
+ case StressBC_HydrostaticTerm:
+ stress=
+ -HydrostaticTerm_Pressure(self->_entryTbl[entry_I].hydrostaticTerm,
+ Mesh_GetVertex(mesh,elementNodes[eNode_I]));
+ if(self->_wall!=Wall_Top)
+ {
+ Stream* errorStr=Journal_Register( Error_Type, self->type );
+ Journal_Firewall(0,errorStr,"You can only apply a HydrostaticTerm StressBC to the top wall.\nYou applied it to the %s wall",WallVC_WallEnumToStr[self->_wall]);
+ }
+
+ if(dim==2)
+ {
+ double dx, dy, *coord0, *coord1;
+ /* StGermain uses the ordering
+
+ 0:x,y
+ 1:x+,y
+ 2:x+,y+
+ 3:x,y+
+
+ So the top two vertices are 2 and 3 */
+
+ coord0=Mesh_GetVertex(mesh,elementNodes[3]);
+ coord1=Mesh_GetVertex(mesh,elementNodes[2]);
+
+ dx=coord1[0]-coord0[0];
+ dy=coord1[1]-coord0[1];
+
+ elForceVec[eNode_I*nodeDofCount + I_AXIS]+=
+ stress*dy/overcount;
+ elForceVec[eNode_I*nodeDofCount + J_AXIS]+=
+ stress*dx/overcount;
+ }
+ else
+ {
+ double *coord0,*coord1,*coord2,*coord3,
+ vector0[3],vector1[3],normal[3];
+
+ /* Decompose the quadrilateral into two triangles.
+ For each triangle, take the cross product and
+ apply the force in that direction. */
+
+ coord0=Mesh_GetVertex(mesh,elementNodes[2]);
+ coord1=Mesh_GetVertex(mesh,elementNodes[3]);
+ coord2=Mesh_GetVertex(mesh,elementNodes[7]);
+ coord3=Mesh_GetVertex(mesh,elementNodes[6]);
+
+ StGermain_VectorSubtraction( vector0, coord1, coord0, dim) ;
+ StGermain_VectorSubtraction( vector1, coord2, coord1, dim) ;
+ StGermain_VectorCrossProduct( normal, vector0, vector1 );
+
+ elForceVec[eNode_I*nodeDofCount + I_AXIS]+=
+ stress*normal[0]/(2*overcount);
+ elForceVec[eNode_I*nodeDofCount + J_AXIS]+=
+ stress*normal[1]/(2*overcount);
+ elForceVec[eNode_I*nodeDofCount + K_AXIS]+=
+ stress*normal[2]/(2*overcount);
+
+ StGermain_VectorSubtraction( vector0, coord2, coord1, dim) ;
+ StGermain_VectorSubtraction( vector1, coord3, coord2, dim) ;
+ StGermain_VectorCrossProduct( normal, vector0, vector1 );
+
+ elForceVec[eNode_I*nodeDofCount + I_AXIS]+=
+ stress*normal[0]/(2*overcount);;
+ elForceVec[eNode_I*nodeDofCount + J_AXIS]+=
+ stress*normal[1]/(2*overcount);;
+ elForceVec[eNode_I*nodeDofCount + K_AXIS]+=
+ stress*normal[2]/(2*overcount);;
+ }
+ break;
+ }
+ /* We have to divide by an overcount_factor, because
+ otherwise different elements will count the same node
+ more than once. In 2D, nodes are counted twice, in 3D,
+ nodes are counted four times. */
+ if(self->_entryTbl[entry_I].type!=StressBC_HydrostaticTerm)
+ elForceVec[eNode_I*nodeDofCount + self->_entryTbl[entry_I].axis]+=
+ stress*area/overcount;
+ }
+ }
+ }
+ NewClass_Delete(incidence);
+}
+
+
+double StressBC_compute_face_area(Wall wall, FeMesh *mesh, Index lElement_I,
+ Dimension_Index dim,
+ IArray *incidence)
+{
+ /* Compute the area of the face. */
+ if(dim==2)
+ {
+ double *coord1, *coord2;
+ int lower,upper,direction;
+ int *elementNodes;
+ switch(wall)
+ {
+ case Wall_Left:
+ lower=0;
+ upper=3;
+ direction=1;
+ break;
+ case Wall_Right:
+ lower=1;
+ upper=2;
+ direction=1;
+ break;
+ case Wall_Bottom:
+ lower=0;
+ upper=1;
+ direction=0;
+ break;
+ case Wall_Top:
+ lower=3;
+ upper=2;
+ direction=0;
+ break;
+ }
+
+ Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
+ MT_VERTEX,incidence);
+ elementNodes=IArray_GetPtr(incidence);
+
+ coord1=Mesh_GetVertex(mesh,elementNodes[lower]);
+ coord2=Mesh_GetVertex(mesh,elementNodes[upper]);
+ return coord2[direction]-coord1[direction];
+ }
+ else
+ {
+ double *coord1, *coord2, *coord3, *coord4;
+ int first, second, third, fourth;
+ int *elementNodes;
+
+ /* StGermain uses the ordering
+ 0: x,y,z
+ 1: x+,y,z
+ 2: x,y+,z
+ 3: x+,y+,z
+ 4: x,y,z+
+ 5: x+,y,z+
+ 6: x,y+,z+
+ 7: x+,y+,z+
+ */
+
+ /* Get the indices for which wall we want to get the area of. */
+ switch(wall)
+ {
+ case Wall_Left:
+ first=0;
+ second=2;
+ third=6;
+ fourth=4;
+ break;
+ case Wall_Right:
+ first=1;
+ second=3;
+ third=7;
+ fourth=5;
+ break;
+ case Wall_Bottom:
+ first=0;
+ second=1;
+ third=5;
+ fourth=4;
+ break;
+ case Wall_Top:
+ first=2;
+ second=3;
+ third=7;
+ fourth=6;
+ break;
+ case Wall_Front:
+ first=0;
+ second=1;
+ third=4;
+ fourth=3;
+ break;
+ case Wall_Back:
+ first=4;
+ second=5;
+ third=7;
+ fourth=6;
+ break;
+ }
+
+ Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
+ MT_VERTEX,incidence);
+ elementNodes=IArray_GetPtr(incidence);
+
+ coord1=Mesh_GetVertex(mesh,elementNodes[first]);
+ coord2=Mesh_GetVertex(mesh,elementNodes[second]);
+ coord3=Mesh_GetVertex(mesh,elementNodes[third]);
+ coord4=Mesh_GetVertex(mesh,elementNodes[fourth]);
+
+ return StGermain_ConvexQuadrilateralArea(coord1,coord2,coord3,coord4,
+ dim);
+ }
+}
+
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/StressBC.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/StressBC.h Mon Dec 21 13:27:10 2009 -0800
@@ -0,0 +1,115 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+** Copyright (c) 2005-2006, Monash Cluster Computing, Building 28, Monash University Clayton Campus,
+** Victoria, 3800, Australia
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+**
+** Contributors:
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Patrick D. Sunter, Software Engineer, VPAC. (patrick at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** David Stegman, Postdoctoral Fellow, Monash University. (david.stegman at sci.monash.edu.au)
+** Wendy Sharples, PhD Student, Monash University (wendy.sharples at sci.monash.edu.au)
+**
+** Modified 2006 by Walter Landry (California Institute of Technology)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+*/
+
+
+#ifndef __Gale_Utils_StressBC_h__
+#define __Gale_Utils_StressBC_h__
+
+ /** Textual name of this class */
+ extern const Type StressBC_Type;
+
+ /** StressBC class contents */
+ #define __StressBC \
+ /* General info */ \
+ __ForceTerm \
+ \
+ /* StressBC info */ \
+ Wall _wall; \
+ StressBC_Entry _entryTbl[3]; \
+ int numEntries; \
+ ConditionFunction_Register* conFunc_Register; \
+ FiniteElementContext* context; \
+
+ struct StressBC { __StressBC };
+
+ StressBC* StressBC_New(
+ Name name,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+ StressBC* _StressBC_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,
+ ForceTerm_AssembleElementFunction* _assembleElement,
+ ConditionFunction_Register* conFunc_Register,
+ Name name,
+ FiniteElementContext* context);
+
+ void StressBC_InitAll(
+ void* forceTerm,
+ ForceVector* forceVector,
+ Swarm* integrationSwarm,
+ ConditionFunction_Register* conFunc_Register,
+ FiniteElementContext* context);
+
+ void _StressBC_Delete( void* forceTerm );
+ void _StressBC_Print( void* forceTerm, Stream* stream );
+
+ void* _StressBC_DefaultNew( Name name ) ;
+ void _StressBC_Construct( void* forceTerm, Stg_ComponentFactory* cf, void* data ) ;
+ void _StressBC_Build( void* forceTerm, void* data ) ;
+ void _StressBC_Initialise( void* forceTerm, void* data ) ;
+ void _StressBC_Execute( void* forceTerm, void* data ) ;
+ void _StressBC_Destroy( void* forceTerm, void* data ) ;
+
+ void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
+ void _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC, char *direction, void *data);
+ double StressBC_compute_face_area(Wall wall, FeMesh *mesh,
+ Index lElement_I,
+ Dimension_Index dim,
+ IArray *incidence);
+ unsigned StressBC_get_overcount(Dimension_Index dim, IJK ijk,
+ unsigned sizes[]);
+
+#endif
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/StressBC.meta
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/StressBC.meta Mon Dec 21 13:27:10 2009 -0800
@@ -0,0 +1,36 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">StressBC</param>
+<param name="Organisation">VPAC and MCC</param>
+<param name="Project">Gale</param>
+<param name="Location">./Gale/Utils/src/</param>
+<param name="Project Web">https://csd.vpac.org/twiki/bin/view/PICellerator/WebHome</param>
+<param name="Copyright">Copyright (C) 2005 VPAC and Monash Cluster Computing.</param>
+<param name="License">https://csd.vpac.org/twiki/bin/view/Stgermain/SoftwareLicense http://www.opensource.org/licenses/bsd-license.php</param>
+<param name="Parent">ForceTerm</param>
+<param name="Description">This adds a stress boundary condition.</param>
+
+<!--Now the interesting stuff-->
+
+
+<list name="Params">
+ <struct>
+ </struct>
+
+</list>
+
+<!-- Add an exmaple XML if possible -->
+<!-- <param name="Example"> -->
+<!-- <struct name="stressBC"> -->
+<!-- <param name="Type">StressBC</param> -->
+<!-- <param name="ForceVector">mom_force</param> -->
+<!-- <param name="Swarm">materialPoints</param> -->
+<!-- <param name="wall">bottom</param> -->
+<!-- <param name="x_type">double</param> -->
+<!-- <param name="x_value">1.0</param> -->
+<!-- </struct> -->
+<!-- </param> -->
+
+</StGermainData>
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/Utils.h
--- a/Utils/src/Utils.h Mon Dec 07 13:47:23 2009 -0800
+++ b/Utils/src/Utils.h Mon Dec 21 13:27:10 2009 -0800
@@ -58,6 +58,7 @@
#include "Underworld_SwarmOutput.h"
#include "RadiogenicHeatingTerm.h"
#include "XDMFGenerator.h"
+ #include "StressBC.h"
#include "Init.h"
#include "Finalise.h"
#endif
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/Wall.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/Wall.c Mon Dec 21 13:27:10 2009 -0800
@@ -0,0 +1,14 @@
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgDomain/StgDomain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
+#include "types.h"
+
+const char* WallEnumToStr[Wall_Size] = {
+ "back",
+ "left",
+ "bottom",
+ "right",
+ "top",
+ "front" };
diff -r 1420849a455e -r d058abdcb8f2 Utils/src/types.h
--- a/Utils/src/types.h Mon Dec 07 13:47:23 2009 -0800
+++ b/Utils/src/types.h Mon Dec 21 13:27:10 2009 -0800
@@ -56,5 +56,34 @@
typedef struct ViscosityField ViscosityField;
typedef struct DensityField DensityField;
typedef struct MixedStabiliserTerm MixedStabiliserTerm;
+ typedef struct StressBC StressBC;
+
+typedef enum
+ {
+ StressBC_Double,
+ StressBC_ConditionFunction,
+ StressBC_HydrostaticTerm
+ } StressBC_Types;
+
+typedef struct
+{
+ StressBC_Types type;
+ double DoubleValue;
+ Index CFIndex;
+ Axis axis;
+ HydrostaticTerm *hydrostaticTerm;
+} StressBC_Entry;
+
+/* Wall types */
+typedef enum
+ {
+ Wall_Back,
+ Wall_Left,
+ Wall_Bottom,
+ Wall_Right,
+ Wall_Top,
+ Wall_Front,
+ Wall_Size
+ } Wall;
#endif
More information about the CIG-COMMITS
mailing list