[cig-commits] [commit] willic3/add-examples-grav2d: Simplify script to generate spatialdb file with initial stresses. (d55bda4)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/pylith

On branch  : willic3/add-examples-grav2d
Link       : https://github.com/geodynamics/pylith/compare/60696f9a25dbef5dc59363251b5bdb722f17da25...d55bda4e00b1c64bd6c2b0c984b928552092b7b9

>---------------------------------------------------------------

commit d55bda4e00b1c64bd6c2b0c984b928552092b7b9
Author: Brad Aagaard <baagaard at usgs.gov>
Date:   Fri Nov 21 08:14:27 2014 -0800

    Simplify script to generate spatialdb file with initial stresses.


>---------------------------------------------------------------

d55bda4e00b1c64bd6c2b0c984b928552092b7b9
 examples/2d/gravity/stress2spatialdb.py | 89 +++++++++++++--------------------
 1 file changed, 34 insertions(+), 55 deletions(-)

diff --git a/examples/2d/gravity/stress2spatialdb.py b/examples/2d/gravity/stress2spatialdb.py
index 6ec63e5..e7b0c1c 100755
--- a/examples/2d/gravity/stress2spatialdb.py
+++ b/examples/2d/gravity/stress2spatialdb.py
@@ -4,72 +4,51 @@ This script creates an initial stress spatial database from output stress
 results.
 """
 
+materials = ["elastic","viscoelastic"]
+
 # The code requires the numpy and h5py packages.
 import numpy
 import h5py
-import math
-# import pdb
-# pdb.set_trace()
 
-# Define input/output files.
-inFiles = ['output/grav_stress-elastic.h5',
-           'output/grav_stress-viscoelastic.h5']
-outFiles = ['grav_stress-elastic.spatialdb',
-            'grav_stress-viscoelastic.spatialdb']
-numFiles = len(inFiles)
 
-# ----------------------------------------------------------------------
-def convertStress(inFile, outFile):
-  """
-  Function to read HDF5 file and output stresses in a spatialdb.
-  """
+# Create coordinate system for spatial database
+from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+from spatialdata.geocoords.CSCart import CSCart
+cs = CSCart()
+cs._configure()
+cs.setSpaceDim(2)
+
+for material in materials:
+
+  filenameH5 = "output/grav_stress-%s.h5" % material
+  filenameDB = "grav_stress-%s.spatialdb" % material
 
   # Open HDF5 file and get coordinates, cells, and stress.
-  h5 = h5py.File(inFile, "r")
-  coords = numpy.array(h5['geometry/vertices'][:], dtype=numpy.float64)
+  h5 = h5py.File(filenameH5, "r")
+  vertices = h5['geometry/vertices'][:]
   cells = numpy.array(h5['topology/cells'][:], dtype=numpy.int)
-  stress = numpy.array(h5['cell_fields/stress'][0,:,:], dtype=numpy.float64)
-  numCells = cells.shape[0]
+  stress = h5['cell_fields/stress'][0,:,:]
   h5.close()
 
   # Get cell centers for output.
-  cellCoords = coords[cells,:]
+  cellCoords = vertices[cells,:]
   cellCenters = numpy.mean(cellCoords, axis=1)
 
-  # Stack coordinates and stresses for output.
-  outputArr = numpy.hstack((cellCenters, stress))
-
-  # Write spatial database.
-  dbHead1 = """// -*- C++ -*- (tell Emacs to use C++ mode)
-//
-// This spatial database gives the initial stresses for the model
-//
-#SPATIAL.ascii 1
-SimpleDB {
-  num-values = 3
-  value-names =  stress-xx stress-yy stress-xy
-  value-units =  Pa Pa Pa
-  num-locs = %d
-  """ % numCells
-
-  dbHead2 = """data-dim = 2
-  space-dim = 2
-  cs-data = cartesian {
-    to-meters = 1.0
-    space-dim = 2
-  }
-}
-"""
-  f = open(outFile, 'w')
-  f.write(dbHead1)
-  f.write(dbHead2)
-  numpy.savetxt(f, outputArr)
-  f.close()
-
-  return
-
-# ======================================================================
-# Loop over stress files.
-for fileNum in range(numFiles):
-  convertStress(inFiles[fileNum], outFiles[fileNum])
+  # Create writer for spatial database file
+  writer = SimpleIOAscii()
+  writer.inventory.filename = filenameDB
+  writer._configure()
+  writer.write({'points': cellCenters,
+              'coordsys': cs,
+              'data_dim': 2,
+              'values': [{'name': "stress-xx",
+                          'units': "Pa",
+                          'data': stress[:,0]},
+                         {'name': "stress-yy",
+                          'units': "Pa",
+                          'data': stress[:,1]},
+                         {'name': "stress-xy",
+                          'units': "Pa",
+                          'data': stress[:,2]},
+                       ]})