[cig-commits] r5374 - in long/3D/Gale/trunk: .
src/PICellerator/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Wed Nov 29 10:08:07 PST 2006
Author: walter
Date: 2006-11-29 10:08:07 -0800 (Wed, 29 Nov 2006)
New Revision: 5374
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/fountain.xml
long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h
Log:
r1119 at earth: boo | 2006-11-28 16:36:18 -0800
Very rough and hard coded stress BC's. It only seems to work for lateral stress. Normal stress causes bad convergence
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1115
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1119
Modified: long/3D/Gale/trunk/fountain.xml
===================================================================
--- long/3D/Gale/trunk/fountain.xml 2006-11-29 18:03:27 UTC (rev 5373)
+++ long/3D/Gale/trunk/fountain.xml 2006-11-29 18:08:07 UTC (rev 5374)
@@ -358,17 +358,17 @@
<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">gravity</param>
- </struct>
+<!-- <struct name="buoyancyForceTerm"> -->
+<!-- <param name="Type">BuoyancyForceTerm</param> -->
+<!-- <param name="ForceVector">mom_force</param> -->
+<!-- <param name="Swarm">picIntegrationPoints</param> -->
+<!-- <param name="gravity">gravity</param> -->
+<!-- </struct> -->
<struct name="stressBC">
<param name="Type">StressBC</param>
<param name="ForceVector">mom_force</param>
<param name="Swarm">picIntegrationPoints</param>
- <param name="force">10.0</param>
+ <param name="force">-1.0</param>
</struct>
<struct name="background">
<param name="Type">Everywhere</param>
@@ -470,6 +470,7 @@
<param name="remesher">velocityRemesher</param>
<param name="velocityField">VelocityField</param>
<param name="wrapTop">True</param>
+ <param name="staticSides">True</param>
<list name="fields">
<struct>
<param name="field">VelocityField</param>
@@ -541,4 +542,5 @@
</struct>
<param name="checkpointEvery">1</param>
<param name="gravity">1.</param>
+<!-- <param name="journal.info">True</param> -->
</StGermainData>
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c 2006-11-29 18:03:27 UTC (rev 5373)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c 2006-11-29 18:08:07 UTC (rev 5374)
@@ -302,80 +302,56 @@
Particle_InCellIndex cellParticleCount;
Element_NodeIndex elementNodeCount;
Dimension_Index dim = forceVector->dim;
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
FiniteElement_Mesh* mesh = forceVector->feVariable->feMesh;
Node_ElementLocalIndex eNode_I;
- Cell_Index cell_I;
ElementType* elementType;
Dof_Index nodeDofCount;
- double force;
- double detJac = 0.0;
- double factor;
- double Ni[8];
- double* xi;
- Material* material;
- MaterialPoint* materialparticle;
- MaterialPointsSwarm* materialSwarm;
+ double stress, area;
- double totalWeight = 0.0;
- double adjustFactor = 0.0;
-
elementType = FeMesh_ElementTypeAt( mesh, lElement_I );
elementNodeCount = elementType->nodeCount;
nodeDofCount = dim;
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[cell_I];
- /* adjust & adjustFactor -- 20060411 Alan
- *
- * The adjust decides whether an adjustFactor should be applied to the resulting factor.
- * If on, the total weight of the particles in the cell are scaled to the cell local volume.
- *
- * This is designed to be used when integrating with swarms which do not cover the whole domain
- * (ie - I used it to do dave.m's test of 1 swarm for blob, 1 swarm for background)
- */
- if ( self->adjust ) {
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- totalWeight += particle->weight;
- }
- adjustFactor = swarm->weights->cellLocalVolume / totalWeight;
- }
- else {
- adjustFactor = 1.0;
- }
-
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- xi = particle->xi;
+ Coord *coord1, *coord2;
- detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
- ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+ coord1=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][0]));
+ coord2=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][1]));
- /* Get parameters */
- material = IntegrationPointsSwarm_GetMaterialOn( (IntegrationPointsSwarm*) swarm, particle );
- materialparticle = OneToOneMapper_GetMaterialPoint( swarm->mapper, particle, &materialSwarm );
+ if(dim==2)
+ area=(*coord2)[0]-(*coord1)[0];
+ else
+ {
+ printf("Need to implement area in StressBC.c for dim==3\n");
+ assert(0);
+ }
- /* Calculate Force */
- force = StressBC_CalcForce( self, (Swarm*)swarm, lElement_I, materialparticle );
- factor = detJac * particle->weight * adjustFactor * force;
+ /* Apply force in verticle direction */
+ for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
+ stress = StressBC_CalcForce(self,mesh,lElement_I,eNode_I);
- /* Apply force in verticle direction */
- for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
- elForceVec[ eNode_I * nodeDofCount + J_AXIS ] += factor * Ni[ eNode_I ] ;
- }
- }
-
+ elForceVec[ eNode_I * nodeDofCount + I_AXIS ] += stress*area;
+/* elForceVec[ eNode_I * nodeDofCount + J_AXIS ] += force; */
+
+
+/* printf("force %d %d %lf %lf %lf %lf %lf\n",eNode_I,lElement_I, */
+/* Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][eNode_I])[0], */
+/* Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][eNode_I])[1], */
+/* stress,area,elForceVec[ eNode_I * nodeDofCount + I_AXIS ]); */
+
+ }
}
-/* The default implementation is for the force to be constant. */
-double _StressBC_CalcForce( void* forceTerm, Swarm* swarm, Element_DomainIndex dElement_I, void* particle ) {
+
+double _StressBC_CalcForce( void* forceTerm, Mesh* mesh, Element_DomainIndex lElement_I, Index node_I ) {
StressBC* self = (StressBC*) forceTerm;
- double *coord;
+ Coord *coord;
- coord=((MaterialPoint*)particle)->coord;
+ coord=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][node_I]));
- if(coord[0]>0.4 && coord[0]<0.6 && coord[1]<0.1)
- return self->force;
+ if((*coord)[1]<0.001)
+ {
+/* return self->force; */
+ return exp(-(0.5-(*coord)[0])*(0.5-(*coord)[0])/0.01)*self->force;
+ }
return 0.0;
}
Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h 2006-11-29 18:03:27 UTC (rev 5373)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.h 2006-11-29 18:08:07 UTC (rev 5374)
@@ -48,9 +48,9 @@
typedef double (StressBC_CalcForceFunction) (
void* forceTerm,
- Swarm* swarm,
+ Mesh* mesh,
Element_LocalIndex lElement_I,
- void* particle );
+ Index node_I );
/** Textual name of this class */
extern const Type StressBC_Type;
@@ -122,9 +122,9 @@
void _StressBC_Destroy( void* forceTerm, void* data ) ;
void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
- double _StressBC_CalcForce( void* forceTerm, Swarm* swarm, Element_LocalIndex lElement_I, void* particle ) ;
+ double _StressBC_CalcForce( void* forceTerm, Mesh* mesh, Element_DomainIndex dElement_I, Index node_I ) ;
- #define StressBC_CalcForce( forceTerm, swarm, lElement_I, particle )\
- (( (StressBC*) forceTerm )->_calcForce( forceTerm, swarm, lElement_I, particle ) )
+ #define StressBC_CalcForce( forceTerm, mesh, lElement_I, node_I )\
+ (( (StressBC*) forceTerm )->_calcForce( forceTerm, mesh, lElement_I, node_I ) )
#endif
More information about the cig-commits
mailing list