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

[walter at geodynamics.org]

Author: walter
Date: 2009-04-17 15:20:31 -0700 (Fri, 17 Apr 2009)
New Revision: 14747

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
Log:
 r2658 at dante:  boo | 2009-04-17 15:20:09 -0700
 Fix hydrostatic StressBC's



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

Modified: long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c	2009-04-17 22:20:02 UTC (rev 14746)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/StressBC.c	2009-04-17 22:20:31 UTC (rev 14747)
@@ -414,15 +414,81 @@
               case StressBC_HydrostaticTerm:
                 stress=
                   -HydrostaticTerm_Pressure(self->_entryTbl[entry_I].hydrostaticTerm,
-                                           Mesh_GetVertex(mesh,elementNodes[eNode_I]));
+                                            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. */
-            elForceVec[eNode_I*nodeDofCount + self->_entryTbl[entry_I].axis]+=
-              stress*area/overcount;
+            if(self->_entryTbl[entry_I].type!=StressBC_HydrostaticTerm)
+              elForceVec[eNode_I*nodeDofCount + self->_entryTbl[entry_I].axis]+=
+                stress*area/overcount;
           }
       }
   }