[cig-commits] r12000 - in long/3D/Gale/trunk: . src/StgFEM/SLE/SystemSetup/src

[walter at geodynamics.org]

Author: walter
Date: 2008-05-22 11:20:25 -0700 (Thu, 22 May 2008)
New Revision: 12000

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
Log:
 r2175 at earth:  boo | 2008-05-22 11:20:04 -0700
 Update nodal pressures, if defined



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

Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c	2008-05-22 18:16:23 UTC (rev 11999)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c	2008-05-22 18:20:25 UTC (rev 12000)
@@ -644,7 +644,59 @@
 	}
 }
 
+void update_nodal_pressure(void *_context)
+{
+  FeVariable *nodal_pressure, *pressure;
+  nodal_pressure = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "NodalPressureField" );
+  pressure = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "PressureField" );
+  if(nodal_pressure!=NULL)
+    {
+      Node_LocalIndex n_i;
+      unsigned	N, nDims;
+            
+      nDims = Mesh_GetDimSize( nodal_pressure->feMesh );
 
+        /* Average the pressure over neighboring cells to get the
+           pressure at the node. */
+      N = FeMesh_GetNodeLocalSize( nodal_pressure->feMesh);
+      for( n_i = 0; n_i < N; n_i++ ) {
+        unsigned nInc;
+        unsigned *inc;
+        unsigned nn;
+        double p,temp;
+        
+        Mesh_GetIncidence(nodal_pressure->feMesh, MT_VERTEX, n_i, nDims, &nInc, &inc);
+        p=0;
+        for(nn=0; nn<nInc; ++nn)
+          {
+            FeVariable_GetValueAtNode(pressure,inc[nn],&temp);
+            p+=temp;
+          }
+        p/=nInc;
+        FeVariable_SetValueAtNode(nodal_pressure,n_i,&p);
+      }
+
+      /* Average the nodal pressures to get a smoothed element pressure */
+      N = FeMesh_GetElementLocalSize( pressure->feMesh);
+      for( n_i = 0; n_i < N; n_i++ ) {
+        unsigned nInc;
+        unsigned *inc;
+        unsigned nn;
+        double p,temp;
+        
+/*         Mesh_GetIncidence(pressure->feMesh, nDims, n_i, MT_VERTEX, &nInc, &inc); */
+        p=0;
+/*         for(nn=0; nn<nInc; ++nn) */
+/*           { */
+/*             FeVariable_GetValueAtNode(nodal_pressure,inc[nn],&temp); */
+/*             p+=temp; */
+/*           } */
+/*         p/=nInc; */
+        FeVariable_SetValueAtNode(pressure,n_i,&p);
+      }
+    }
+}
+
 void SystemLinearEquations_NonLinearExecute( void* sle, void* _context ) {
 	SystemLinearEquations*	self            = (SystemLinearEquations*) sle;
 	Vector*                 previousVector  = NULL;
@@ -665,6 +717,7 @@
 	wallTime = MPI_Wtime();		
         
         /* Set up the current BC's (only really useful for friction BC's). */
+
         stress     = (FeVariable*)FieldVariable_Register_GetByName( ((FiniteElementContext *)_context)->fieldVariable_Register, "StressField" );
         if(stress!=NULL)
           {
@@ -677,6 +730,8 @@
         currentVector =
           SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
         
+        update_nodal_pressure(_context);
+
         /* Non linear iterations */
         converged=0;
         for ( self->nonLinearIteration_I = 1 ;
@@ -697,6 +752,8 @@
           self->linearExecute( self, _context );
           self->hasExecuted = True;
           
+          update_nodal_pressure(_context);
+
           /* TODO - Give option which solution vector to test */
           currentVector =
             SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;