[cig-commits] r14314 - in long/3D/Gale/trunk: . src/Underworld/plugins/Output/VTKOutput

[walter at geodynamics.org]

Author: walter
Date: 2009-03-13 13:57:23 -0700 (Fri, 13 Mar 2009)
New Revision: 14314

Modified:
   long/3D/Gale/trunk/
   long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
Log:
 r2552 at dante:  boo | 2009-03-13 13:56:45 -0700
 Make VTKOutput print out the complete pressure



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

Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c	2009-03-13 19:06:58 UTC (rev 14313)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/VTKOutput/VTKOutput.c	2009-03-13 20:57:23 UTC (rev 14314)
@@ -45,6 +45,7 @@
 #include <PICellerator/PICellerator.h>
 #include <Underworld/Underworld.h>
 #include "VTKOutput.h"
+#include <PICellerator/Utils/HydrostaticTerm.h>
 
 #include <stdlib.h>
 #include <string.h>
@@ -435,6 +436,7 @@
   int nDims, pn, i;
   int lower[3], upper[3], p_lower[3], p_upper[3];
 
+  HydrostaticTerm *hydrostaticTerm;
   Name field_filename;
   FILE *field_fp, *pfield_fp;
 
@@ -442,6 +444,9 @@
      the values. */
   Grid *elGrid, *vertGrid;
 
+  hydrostaticTerm = (HydrostaticTerm*)LiveComponentRegister_Get(
+                                                            self->CF->LCRegister,
+                                                            "hydrostaticTerm" );
   /* Open the file */
 
   Stg_asprintf( &field_filename, "%s/fields.%d.%05d.vts", self->outputPath,
@@ -460,6 +465,21 @@
       Memory_Free( field_filename );
     }
 
+  /* First, update the fields that are derived from particles.  This
+     is needed to get an accurate pressure, because it needs the
+     stresses. */
+
+  for ( var_I = 0; var_I < self->fieldVariable_Register->objects->count;
+        var_I++ ) {
+    FieldVariable* fieldVar;
+    fieldVar = FieldVariable_Register_GetByIndex( self->fieldVariable_Register,
+                                                  var_I );
+    if (Stg_Class_IsInstance( fieldVar, ParticleFeVariable_Type ))
+      {
+        ParticleFeVariable_Update( fieldVar );
+      }
+  }
+
   /* First, output the coordinates.  We have to do a huge song and
      dance just to get the extents of the mesh. */
 
@@ -478,11 +498,6 @@
 
       feVar=(FeVariable*)fieldVar;
 
-      if (Stg_Class_IsInstance( fieldVar, ParticleFeVariable_Type ))
-        {
-          ParticleFeVariable_Update( fieldVar );
-        }
-
       if(!header_printed)
         {
           Mesh *mesh;
@@ -619,9 +634,59 @@
                 case 3:
                 case 9:
                 case 27:
-                  for ( dof_I = 0; dof_I < dofAtEachNodeCount; dof_I++ ) {
-                    fprintf(field_fp, "%.15g ", variableValues[dof_I] );
-                  }
+                  /* Special case the pressure */
+                  if(!strcmp(feVar->name,"PressureField"))
+                    {
+                      double p;
+                      /* First add the trace of the stress */
+                      int stress_I;
+
+                      for(stress_I = 0;
+                          stress_I<self->fieldVariable_Register->objects->count;
+                          stress_I++ ) {
+                        FieldVariable* stressVar;
+                        stressVar =
+                          FieldVariable_Register_GetByIndex(self->fieldVariable_Register,
+                                                            stress_I );
+                        if(!strcmp(stressVar->name,"StressField"))
+                          {
+                            if (Stg_Class_IsInstance( stressVar,
+                                                      FeVariable_Type ))
+                              {
+                                FeVariable* sVar;
+                                double stressValues[MAX_FIELD_COMPONENTS];	
+                                sVar=(FeVariable*)stressVar;
+                                
+                                FeVariable_GetValueAtNode(sVar,local,stressValues);
+                                if(nDims==2)
+                                  {
+                                    p=(stressValues[0]+stressValues[1])/2;
+                                  }
+                                else
+                                  {
+                                    p=(stressValues[0]+stressValues[1]
+                                       +stressValues[2])/3;
+                                  }
+                              }
+                            break;
+                          }
+                      }
+                      /* Next add the hydrostatic term */
+                      if(hydrostaticTerm)
+                        {
+                          double *coord;
+                          coord=Mesh_GetVertex(feVar->feMesh,local);
+                          p+=HydrostaticTerm_Pressure(hydrostaticTerm,coord);
+                        }
+                      p+=variableValues[0];
+                      fprintf(field_fp, "%.15g ", p );
+                    }
+                  else
+                    {
+                      for ( dof_I = 0; dof_I < dofAtEachNodeCount; dof_I++ ) {
+                        fprintf(field_fp, "%.15g ", variableValues[dof_I] );
+                      }
+                    }
                   break;
                 case 4:
                   fprintf(field_fp, "%.15g %.15g 0", variableValues[0],