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

[walter at geodynamics.org]

Author: walter
Date: 2007-01-22 15:35:32 -0800 (Mon, 22 Jan 2007)
New Revision: 5859

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
Log:
 r1467 at earth:  boo | 2007-01-22 15:30:53 -0800
 Make stress boundaries work in 3D



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

Modified: long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c	2007-01-22 22:51:36 UTC (rev 5858)
+++ long/3D/Gale/trunk/src/PICellerator/Utils/src/StressBC.c	2007-01-22 23:35:32 UTC (rev 5859)
@@ -322,6 +322,25 @@
 	_ForceTerm_Destroy( forceTerm, data );
 }
 
+double cross(double *coord1, double *coord2, int normal)
+{
+  switch(normal)
+    {
+    case 0:
+      return coord1[1]*coord2[2]-coord1[2]*coord2[1];
+      break;
+    case 1:
+      return coord1[0]*coord2[2]-coord1[2]*coord2[0];
+      break;
+    case 2:
+      return coord1[0]*coord2[1]-coord1[1]*coord2[0];
+      break;
+    default:
+      printf("Bad call of cross() in StressBC\n");
+      abort();
+    }
+  return 0;
+}
 
 void _StressBC_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
 	StressBC*               self               = (StressBC*) forceTerm;
@@ -382,51 +401,79 @@
         else
           {
             double *coord1, *coord2, *coord3, *coord4;
-            int lower,upper,direction;
+            int first, second, third, fourth, normal;
+
+
+            Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
+                              MT_VERTEX,&elementNodeCount, &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+
+            */
+
             switch(self->_wall)
               {
               case StressBC_Wall_Left:
-                lower=0;
-                upper=3;
-                direction=1;
+                first=0;
+                second=2;
+                third=6;
+                fourth=4;
+                normal=0;
                 break;
               case StressBC_Wall_Right:
-                lower=1;
-                upper=2;
-                direction=1;
+                first=1;
+                second=3;
+                third=7;
+                fourth=5;
+                normal=0;
                 break;
               case StressBC_Wall_Bottom:
-                lower=0;
-                upper=1;
-                direction=0;
+                first=0;
+                second=1;
+                third=5;
+                fourth=4;
+                normal=1;
                 break;
               case StressBC_Wall_Top:
-                lower=3;
-                upper=2;
-                direction=0;
+                first=2;
+                second=3;
+                third=7;
+                fourth=6;
+                normal=1;
                 break;
               case StressBC_Wall_Front:
-                lower=0;
-                upper=3;
-                direction=1;
+                first=0;
+                second=1;
+                third=4;
+                fourth=3;
+                normal=2;
                 break;
               case StressBC_Wall_Back:
-                lower=0;
-                upper=3;
-                direction=1;
+                first=4;
+                second=5;
+                third=7;
+                fourth=6;
+                normal=2;
                 break;
               }
 
             Mesh_GetIncidence(mesh, Mesh_GetDimSize(mesh), lElement_I,
                               MT_VERTEX,&elementNodeCount, &elementNodes);
                               
-            coord1=Mesh_GetVertex(mesh,elementNodes[lower]);
-            coord2=Mesh_GetVertex(mesh,elementNodes[upper]);
-            area=coord2[direction]-coord1[direction];
-
-
-            printf("Need to implement area in StressBC.c for dim==3\n");
-            assert(0);
+            coord1=Mesh_GetVertex(mesh,elementNodes[first]);
+            coord2=Mesh_GetVertex(mesh,elementNodes[second]);
+            coord3=Mesh_GetVertex(mesh,elementNodes[third]);
+            coord4=Mesh_GetVertex(mesh,elementNodes[fourth]);
+            area=fabs(cross(coord1,coord2,normal) + cross(coord2,coord3,normal)
+                      + cross(coord3,coord4,normal)
+                      + cross(coord4,coord1,normal))/2;
           }
 
         /* Apply the stress */