[cig-commits] r5376 - in long/3D/Gale/trunk: . src/PICellerator/Utils/src

[walter at geodynamics.org]

Author: walter
Date: 2006-11-29 10:08:18 -0800 (Wed, 29 Nov 2006)
New Revision: 5376

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
Log:
 r1121 at earth:  boo | 2006-11-29 01:29:43 -0800
 Pull edge detection out of StressBC_CalcForce



Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
   - 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1120
   + 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1121

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:08:10 UTC (rev 5375)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c	2006-11-29 18:08:18 UTC (rev 5376)
@@ -296,10 +296,6 @@
 
 void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
 	StressBC*               self               = (StressBC*) forceTerm;
-	IntegrationPoint*                particle;
-	StressBC_MaterialExt*   materialExt;
-	Particle_InCellIndex             cParticle_I;
-	Particle_InCellIndex             cellParticleCount;
 	Element_NodeIndex                elementNodeCount;
 	Dimension_Index                  dim                = forceVector->dim;
 	FiniteElement_Mesh*              mesh               = forceVector->feVariable->feMesh;
@@ -308,17 +304,23 @@
 	Dof_Index                        nodeDofCount;
 	double                           stress, area;
 
+	MeshLayout*		meshLayout = mesh->layout;
+	HexaMD*			decomp = (HexaMD*)meshLayout->decomp;
+
 	elementType       = FeMesh_ElementTypeAt( mesh, lElement_I );
 	elementNodeCount  = elementType->nodeCount;
 	nodeDofCount      = dim;
+        
+        IJK			ijk;
 
-        Coord *coord1, *coord2;
-
-        coord1=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][0]));
-        coord2=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][1]));
-
         if(dim==2)
-          area=(*coord2)[0]-(*coord1)[0];
+          {
+            Coord *coord1, *coord2;
+            
+            coord1=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][0]));
+            coord2=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][1]));
+            area=(*coord2)[0]-(*coord1)[0];
+          }
         else
           {
             printf("Need to implement area in StressBC.c for dim==3\n");
@@ -327,9 +329,13 @@
 
         /* Apply force in verticle direction */
         for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) { 		
-          stress = StressBC_CalcForce(self,mesh,lElement_I,eNode_I);
+          /* Make sure that we are on the boundary */
+          RegularMeshUtils_Node_1DTo3D( decomp, mesh->elementNodeTbl[lElement_I][eNode_I],
+                                        &ijk[0], &ijk[1], &ijk[2] );      
+          if( ijk[1] == 0 ) {
+            stress = StressBC_CalcForce(self,mesh,lElement_I,eNode_I);
 
-          elForceVec[ eNode_I * nodeDofCount + I_AXIS ] += stress*area;
+            elForceVec[ eNode_I * nodeDofCount + I_AXIS ] += stress*area;
 /*           elForceVec[ eNode_I * nodeDofCount + J_AXIS ] += force; */
 
 
@@ -338,24 +344,16 @@
 /*                  Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][eNode_I])[1], */
 /*                  stress,area,elForceVec[ eNode_I * nodeDofCount + I_AXIS ]); */
 
+          }
         }
 }
 
 
 double _StressBC_CalcForce( void* forceTerm, Mesh* mesh, Element_DomainIndex lElement_I, Index node_I ) {
 	StressBC*               self               = (StressBC*) forceTerm;
-	MeshLayout*		meshLayout = mesh->layout;
-	HexaMD*			decomp = (HexaMD*)meshLayout->decomp;
-        
-        IJK			ijk;
 		
-        RegularMeshUtils_Node_1DTo3D( decomp, mesh->elementNodeTbl[lElement_I][node_I],
-                                      &ijk[0], &ijk[1], &ijk[2] );      
-        if( ijk[1] == 0 ) {
-          Coord *coord;
 /*             return self->force; */
-          coord=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][node_I]));
-          return exp(-(0.5-(*coord)[0])*(0.5-(*coord)[0])/0.01)*self->force;
-        }
-	return 0.0;
+        Coord *coord;
+        coord=&(Mesh_CoordAt(mesh,mesh->elementNodeTbl[lElement_I][node_I]));
+        return exp(-(0.5-(*coord)[0])*(0.5-(*coord)[0])/0.01)*self->force;
 }