[cig-commits] r13853 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Mon Jan 12 16:03:36 PST 2009
Author: walter
Date: 2009-01-12 16:03:36 -0800 (Mon, 12 Jan 2009)
New Revision: 13853
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/Utils/src/DivergenceForce.c
long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c
long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h
long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c
long/3D/Gale/trunk/src/Gale/Utils/src/types.h
Log:
r2449 at dante: boo | 2009-01-12 10:50:49 -0800
Add HydrostaticTerm to StressBC
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2448
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2449
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/DivergenceForce.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/DivergenceForce.c 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/DivergenceForce.c 2009-01-13 00:03:36 UTC (rev 13853)
@@ -49,6 +49,7 @@
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
#include "types.h"
#include "DivergenceForce.h"
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StaticFrictionVC.c 2009-01-13 00:03:36 UTC (rev 13853)
@@ -26,6 +26,7 @@
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
#include "types.h"
#include "StaticFrictionVC.h"
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c 2009-01-13 00:03:36 UTC (rev 13853)
@@ -48,6 +48,7 @@
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
#include <PICellerator/Voronoi/Voronoi.h>
#include <PICellerator/PopulationControl/PopulationControl.h>
#include <PICellerator/Weights/Weights.h>
@@ -217,15 +218,15 @@
self->_wall = Wall_Size; /* invalid entry */
}
}
- _StressBC_GetValues(cf,self,"x");
- _StressBC_GetValues(cf,self,"y");
- _StressBC_GetValues(cf,self,"z");
+ _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)
+ char *direction, void *data)
{
StressBC* self = (StressBC*)stressBC;
char temp_str[20];
@@ -282,6 +283,15 @@
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 );
@@ -317,87 +327,121 @@
}
void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
- StressBC* self = (StressBC*) forceTerm;
- Element_NodeIndex elementNodeCount;
- 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;
+ StressBC* self = (StressBC*) forceTerm;
+ Element_NodeIndex elementNodeCount;
+ 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;
+
+ Node_DomainIndex* elementNodes = NULL;
- Node_DomainIndex* elementNodes = NULL;
+ 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" ) );
-
- area=StressBC_compute_face_area(self->_wall,mesh,lElement_I,
- &elementNodeCount,dim,&elementNodes);
- /* 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)
- {
- unsigned overcount=StressBC_get_overcount(dim,ijk,grid->sizes);
- 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;
- }
- /* 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. */
- elForceVec[ eNode_I * nodeDofCount
- + self->_entryTbl[entry_I].axis ] += stress*area/overcount;
- }
- }
- }
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ nodeDofCount = dim;
+
+ grid=*(Grid**)ExtensionManager_Get(mesh->info, mesh,
+ ExtensionManager_GetHandle(mesh->info,
+ "vertexGrid"));
+
+ area=StressBC_compute_face_area(self->_wall,mesh,lElement_I,
+ &elementNodeCount,dim,&elementNodes);
+ /* 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)
+ {
+ unsigned overcount=StressBC_get_overcount(dim,ijk,grid->sizes);
+ 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->context->nonLinearIteration_I)!=0) */
+/* { */
+/* double density=1; */
+/* double gravity=1; */
+/* double v[3]; */
+/* double damping; */
+/* FeVariable_GetValueAtNode(velocityField, */
+/* elementNodes[eNode_I],v); */
+
+/* printf("stress %d %d %d %g %g %g %g\n", */
+/* *(self->context->nonLinearIteration_I), */
+/* lElement_I, */
+/* elementNodes[eNode_I],stress, */
+/* density*gravity*self->context->dt*v[1], */
+/* /\* stress-density*gravity*1.09375*v[1],v[1]); *\/ */
+/* stress-density*gravity*self->context->dt*v[1],v[1]); */
+
+/* /\* damping=density*gravity*1.09375*v[1]; *\/ */
+/* damping=density*gravity*self->context->dt*v[1]; */
+/* /\* Make sure that we don't overcompensate. *\/ */
+/* /\* if(((damping>0 && stress>=0) || (damping<0 && stress<=0)) *\/ */
+/* /\* && fabs(damping)<fabs(stress)) *\/ */
+/* /\* stress+=damping; *\/ */
+/* } */
+ 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. */
+ elForceVec[eNode_I*nodeDofCount + self->_entryTbl[entry_I].axis]+=
+ stress*area/overcount;
+ }
+ }
+ }
}
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.h 2009-01-13 00:03:36 UTC (rev 13853)
@@ -105,7 +105,7 @@
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 _StressBC_GetValues(Stg_ComponentFactory* cf, void *stressBC, char *direction, void *data);
double StressBC_compute_face_area(Wall wall, FeMesh *mesh,
Index lElement_I,
Element_NodeIndex *elementNodeCount,
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/Wall.c 2009-01-13 00:03:36 UTC (rev 13853)
@@ -1,6 +1,7 @@
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
#include "types.h"
const char* WallEnumToStr[Wall_Size] = {
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2009-01-13 00:03:33 UTC (rev 13852)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2009-01-13 00:03:36 UTC (rev 13853)
@@ -56,7 +56,8 @@
typedef enum
{
StressBC_Double,
- StressBC_ConditionFunction
+ StressBC_ConditionFunction,
+ StressBC_HydrostaticTerm
} StressBC_Types;
typedef struct
@@ -65,6 +66,7 @@
double DoubleValue;
Index CFIndex;
Axis axis;
+ HydrostaticTerm *hydrostaticTerm;
} StressBC_Entry;
/* Wall types */
More information about the CIG-COMMITS
mailing list