[cig-commits] r12582 - in long/3D/Gale/trunk: . src/StgFEM/plugins/StandardConditionFunctions

walter at geodynamics.org walter at geodynamics.org
Thu Aug 7 15:56:27 PDT 2008


Author: walter
Date: 2008-08-07 15:56:27 -0700 (Thu, 07 Aug 2008)
New Revision: 12582

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
   long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h
Log:
 r2296 at earth:  boo | 2008-08-07 13:33:12 -0700
 Add in a non-working PressureEquilibrium function.  It is not really needed, this is just for archiving purposes



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

Modified: long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c	2008-08-07 14:32:47 UTC (rev 12581)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.c	2008-08-07 22:56:27 UTC (rev 12582)
@@ -131,6 +131,9 @@
 	condFunc = ConditionFunction_New( StgFEM_StandardConditionFunctions_TemperatureProfile, "TemperatureProfile");
 	ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
 
+	condFunc = ConditionFunction_New( StgFEM_StandardConditionFunctions_PressureEquilibrium, "PressureEquilibrium");
+	ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
+
 	condFunc = ConditionFunction_New( StG_FEM_StandardConditionFunctions_Gaussian, "Gaussian");
 	ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
 
@@ -1109,7 +1112,100 @@
           }
 }
 
+/* Wall types */
+typedef enum
+  {
+    Wall_Back,
+    Wall_Left,
+    Wall_Bottom,
+    Wall_Right,
+    Wall_Top,
+    Wall_Front,
+    Wall_Size
+  } Wall;
 
+void StgFEM_StandardConditionFunctions_PressureEquilibrium( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
+	FiniteElementContext *	context            = (FiniteElementContext*)_context;
+	FeVariable             *pressure, *deviatoric_stress;
+	FeMesh*     mesh               = NULL;
+	Dictionary*             dictionary         = context->dictionary;
+	double*                 result             = (double*) _result;
+        double*                 coord;
+        IJK ijk;
+	unsigned	nDims, direction, boundary, basis[3];
+	Grid*		grid;
+        Bool include_boundary[3];
+        double normal_stress, tangential_stress[9], tangential_norm, n[3], p,
+          reference_pressure, rho, gravity;
+
+        pressure=(FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "NodalPressureField" );
+
+        FeVariable_GetValueAtNode(pressure,node_lI,&p);
+        /* If this is the first step and the nodal pressure has not
+           been set, then don't move the boundary. */
+        if(p==0 || context->dt==0.0)
+          {
+            *result=0;
+          }
+        else
+          {
+            deviatoric_stress=(FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "StressField" );
+            
+            mesh = pressure->feMesh;
+            nDims = Mesh_GetDimSize( mesh );
+            coord = Mesh_GetVertex( mesh, node_lI );
+            RegularMeshUtils_Node_1DTo3D
+              (mesh, Mesh_DomainToGlobal(mesh,MT_VERTEX,node_lI),ijk);
+            grid = *(Grid**)ExtensionManager_Get
+              ( mesh->info, mesh, 
+                ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+            
+            /* Hard code to the bottom of the box for now */
+            StaticFrictionVC_get_basis_vectors(Wall_Bottom,grid->sizes,
+                                               &direction,&boundary,basis);
+            
+            include_boundary[0]=True;
+            include_boundary[1]=True;
+            include_boundary[2]=True;
+            
+            StaticFrictionVC_get_interface_stresses
+              (deviatoric_stress, mesh, grid->sizes, node_lI, ijk, nDims, basis,
+               include_boundary,include_boundary,
+               &normal_stress,tangential_stress,&tangential_norm,n);
+            
+            reference_pressure=Dictionary_GetDouble_WithDefault
+              ( dictionary, "pressureEquilibriumPressure", 0.0 );
+            rho=Dictionary_GetDouble_WithDefault
+              ( dictionary, "pressureEquilibriumDensity", 0.0 );
+            gravity=Dictionary_GetDouble_WithDefault
+              ( dictionary, "pressureEquilibriumGravity", 0.0 );
+            if(reference_pressure==0.0 || rho==0.0 || gravity==0.0)
+              {
+                printf("pressureEquilibriumPressure, pressureEquilibriumDensity, and pressureEqulibriumGravity must be non-zero when\nusing the PressureEquilibirum function.\n");
+                abort();
+              }
+            else
+              {
+                /* We set the velocity so that with the old time step,
+                   distance traveled will be cancel out the pressure
+                   disequilibrium.  For stability, we multiply by 0.5, so
+                   that it does not overshoot the solution. */
+                *result=(reference_pressure-(p+normal_stress))*0.5
+                  /(context->dt*rho*gravity);
+
+                if(fabs(*result)>1.0)
+                  *result*=1/fabs(*result);
+              }
+
+            printf("pressure %d %lf %lf %lf %lf %lf %lf\n",
+                   node_lI,reference_pressure,
+                   p,normal_stress,context->dt,gravity,*result);
+          }
+
+        return;
+}
+
+
 Bool StgFEM_StandardConditionFunctions_Init( int* argc, char** argv[] ) {
   Stg_ComponentRegister* componentsRegister = Stg_ComponentRegister_Get_ComponentRegister();
   Stg_ComponentRegister_Add(componentsRegister,

Modified: long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h	2008-08-07 14:32:47 UTC (rev 12581)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/StandardConditionFunctions/StandardConditionFunctions.h	2008-08-07 22:56:27 UTC (rev 12582)
@@ -88,6 +88,7 @@
 void StG_FEM_StandardConditionFunctions_StepFunctionProduct3( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
 void StG_FEM_StandardConditionFunctions_StepFunctionProduct4( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
 void StgFEM_StandardConditionFunctions_TemperatureProfile( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
+void StgFEM_StandardConditionFunctions_PressureEquilibrium( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
 void StG_FEM_StandardConditionFunctions_Gaussian( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;
 
 void StgFEM_StandardConditionFunctions_ConvectionBenchmark( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) ;



More information about the cig-commits mailing list