[cig-commits] r7516 - in cs/cigma/branches/cigma-0.9: . contrib

[luis at geodynamics.org]

Author: luis
Date: 2007-06-26 10:13:37 -0700 (Tue, 26 Jun 2007)
New Revision: 7516

Added:
   cs/cigma/branches/cigma-0.9/contrib/
   cs/cigma/branches/cigma-0.9/contrib/bmssnog_tet4.py
Log:
darcs patches:
  * Added contrib directory
  * Added script for importing benchmark meshes of various resolutions
  - Linear tetrahedral meshes for bmssnog benchmark (strike-slip no-gravity)



Added: cs/cigma/branches/cigma-0.9/contrib/bmssnog_tet4.py
===================================================================
--- cs/cigma/branches/cigma-0.9/contrib/bmssnog_tet4.py	2007-06-26 17:12:02 UTC (rev 7515)
+++ cs/cigma/branches/cigma-0.9/contrib/bmssnog_tet4.py	2007-06-26 17:13:37 UTC (rev 7516)
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+#
+# The data files for this example were taken from .tgz files in
+# http://geodynamics.org/cig/workinggroups/short/workarea/benchmarks/benchmark-strikeslip/pylith-0.8-input
+#
+
+import numpy
+import tables
+
+def readfromfile(fp, x, **kw):
+    (n,d) = x.shape
+    skip = kw.get('skip', 0)
+    factor = kw.get('factor', 1.0)
+    range = kw.get('range', slice(1,4))
+    dtype = kw.get('dtype', float)
+    offset = kw.get('offset', 1)
+    # skip the specified number of lines
+    for i in xrange(skip):
+        fp.readline()
+    # read n lines into array
+    for i in xrange(n):
+        line = fp.readline()
+        cols = line.split()
+        k = int(cols[0]) - offset
+        x[k] = map(dtype, cols[range])
+    # finally, apply the conversion factor
+    x *= factor
+    return
+
+def readfrom(filename, x, **kw):
+    print "Reading file", filename
+    fp = open(filename, 'r')
+    readfromfile(fp, x, **kw)
+    fp.close()
+    return
+
+###############################################################################
+
+if __name__ == '__main__':
+
+    filename = 'bmssnog_tet4.h5'
+    topdir = 'bmssnog_tet4_%s/'
+    coordinates_file = topdir + 'tet4_%s.coord'
+    connectivity_file = topdir + 'tet4_%s.connect'
+
+    resolutions = ['1000m', '0500m', '0250m']
+    nno_db = {'1000m': 15625, '0500m': 117649, '0250m': 912673}
+    nel_db = {'1000m': 82944, '0500m': 663552, '0250m': 5308416}
+
+
+    h5 = tables.openFile(filename, 'w')
+
+    for res in resolutions:
+        (nno,nel) = (nno_db[res], nel_db[res])
+        coordinates = numpy.zeros((nno,3), numpy.float32)
+        connectivity = numpy.zeros((nel,4), numpy.int32)
+        readfrom(coordinates_file % (res,res), coordinates,
+                 skip=1, factor=1000.0)
+        readfrom(connectivity_file % (res,res), connectivity,
+                 skip=0, range=slice(4,8), dtype=int)
+        connectivity -= 1
+        resgroup = h5.createGroup('/', 'mesh_'+res)
+        h5.createArray(resgroup, 'coordinates', coordinates)
+        h5.createArray(resgroup, 'connectivity', connectivity)
+
+    h5.close()