[cig-commits] [commit] master: Started work on tpv28. (ba4ce6e)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Feb 18 17:31:52 PST 2014
Repository : ssh://geoshell/pylith_benchmarks
On branch : master
Link : https://github.com/geodynamics/pylith_benchmarks/compare/ad0cd8a9f840a53e433cb88f0da02308038cc20b...ba4ce6e681eb615c048e28577878b57d81946529
>---------------------------------------------------------------
commit ba4ce6e681eb615c048e28577878b57d81946529
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Tue Feb 18 17:31:51 2014 -0800
Started work on tpv28.
>---------------------------------------------------------------
ba4ce6e681eb615c048e28577878b57d81946529
dynamic/scecdynrup/tpv28/create_faultsurf.py | 62 ++++++++++++++
dynamic/scecdynrup/tpv28/cubit_io.py | 117 +++++++++++++++++++++++++++
2 files changed, 179 insertions(+)
diff --git a/dynamic/scecdynrup/tpv28/create_faultsurf.py b/dynamic/scecdynrup/tpv28/create_faultsurf.py
new file mode 100755
index 0000000..bfdbfb5
--- /dev/null
+++ b/dynamic/scecdynrup/tpv28/create_faultsurf.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+# Python script to create Exodus II mapped mesh of fault geometry for
+# import into CUBIT.
+
+# Filename for ExodusII file (output)
+filenameEXO = "fault.exo"
+
+# Geometry of simulation domain
+domainLength = 64.0e+3
+domainHeight = 32.0e+3
+buffer = 2.0e+3 # Extend fault plane beyond domain
+
+# Geometry of bump
+bumpRadius = 3.0e+3
+bumpAmplitude = 6.0e+2
+bumpY = 10.5e+3
+bumpZ = -7.5e+3
+bumpDx = 1.0e+3
+
+# ======================================================================
+import os
+import numpy
+
+from cubit_io import write_exodus_file
+from netCDF4 import Dataset
+
+# ----------------------------------------------------------------------
+def bump(yz):
+ y = yz[:,0]
+ z = yz[:,1]
+ r1 = ((y+bumpY)**2 + (z-bumpZ)**2)**0.5
+ r2 = ((y-bumpY)**2 + (z-bumpZ)**2)**0.5
+ mask1 = r1 <= 3.0e+3
+ mask2 = r2 <= 3.0e+3
+ npts = yz.shape[0]
+ x = mask1*-0.5*bumpAmplitude*(1.0+numpy.cos(numpy.pi*r1/3.0e+3)) + \
+ mask2*-0.8*bumpAmplitude*(1.0+numpy.cos(numpy.pi*r2/3.0e+3))
+ return x
+
+
+# ----------------------------------------------------------------------
+yA = numpy.linspace(-0.5*domainLength-buffer, -bumpY-bumpRadius, 2)
+yB = numpy.arange(-bumpY-bumpRadius+bumpDx, -bumpY+bumpRadius, bumpDx)
+yC = numpy.linspace(-bumpY+bumpRadius, +bumpY-bumpRadius, 2)
+yD = numpy.arange(+bumpY-bumpRadius+bumpDx, +bumpY+bumpRadius, bumpDx)
+yE = numpy.linspace(+bumpY+bumpRadius, +0.5*domainLength+buffer, 2)
+y = numpy.hstack((yA, yB, yC, yD, yE))
+numY = y.shape[0]
+
+zA = numpy.linspace(-domainHeight-buffer, bumpZ-bumpRadius, 2)
+zB = numpy.arange(bumpZ-bumpRadius+bumpDx, bumpZ+bumpRadius, bumpDx)
+zC = numpy.linspace(bumpZ+bumpRadius, 0.0, 2)
+z = numpy.hstack((zA, zB, zC))
+numZ = z.shape[0]
+
+yz = numpy.zeros((numY*numZ, 2), dtype=numpy.float64)
+
+#write_exodus_file(filenameEXO, cells, vertices)
+
+
+# End of file
diff --git a/dynamic/scecdynrup/tpv28/cubit_io.py b/dynamic/scecdynrup/tpv28/cubit_io.py
new file mode 100644
index 0000000..ae8ed45
--- /dev/null
+++ b/dynamic/scecdynrup/tpv28/cubit_io.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+
+# ----------------------------------------------------------------------
+# write_exodus_file
+# ----------------------------------------------------------------------
+def write_exodus_file(filename, cells, vertices, shape="SHELL4"):
+ """
+ Write Exodus-II file compatible with CUBIT.
+
+ cells is a 0-based array (ncells, ncorners).
+
+ vertices is (nvertices, dim).
+
+ All cells are placed in a single block.
+
+ Requires netCDF4 module.
+ """
+ import numpy
+ from netCDF4 import Dataset
+
+ len_string = 33
+
+ root = Dataset(filename, 'w', format='NETCDF3_CLASSIC')
+
+ # Set global attributes
+ root.api_version = 4.98
+ root.version = 4.98
+ root.floating_point_word_size = 8
+ root.file_size = 0
+ root.title = "cubit"
+
+ # Setup dimensions
+
+ # Generic information
+ root.createDimension('len_string', len_string)
+ root.createDimension('len_line', 81)
+ root.createDimension('four', 4)
+ root.createDimension('num_qa_rec', 1)
+ root.createDimension('time_step', None)
+
+ # Mesh specific information
+ (ncells, ncorners) = cells.shape
+ (nvertices, dim) = vertices.shape
+ root.createDimension('num_dim', dim)
+ root.createDimension('num_el_blk', 1)
+ root.createDimension('num_nod_per_el1', ncorners)
+ root.createDimension('num_att_in_blk1', 1)
+
+ root.createDimension('num_nodes', nvertices)
+ root.createDimension('num_elem', ncells)
+ root.createDimension('num_el_in_blk1', ncells)
+
+ # Setup variables
+ connect1 = root.createVariable('connect1', numpy.int32,
+ ('num_el_in_blk1', 'num_nod_per_el1',))
+
+ coord = root.createVariable('coord', numpy.float64,
+ ('num_dim', 'num_nodes',))
+
+ time_whole = root.createVariable('time_whole', numpy.float64,
+ ('time_step',))
+
+ coor_names = root.createVariable('coor_names', 'S1',
+ ('num_dim', 'len_string',))
+
+ qa_records = root.createVariable('qa_records', 'S1',
+ ('num_qa_rec', 'four', 'len_string',))
+
+ eb_names = root.createVariable('eb_names', 'S1',
+ ('num_el_blk', 'len_string',))
+
+ elem_map = root.createVariable('elem_map', numpy.int32,
+ ('num_elem',))
+
+ eb_status = root.createVariable('eb_status', numpy.int32,
+ ('num_el_blk',))
+
+ eb_prop1 = root.createVariable('eb_prop1', numpy.int32,
+ ('num_el_blk',))
+
+ attrib1 = root.createVariable('attrib1', numpy.float64,
+ ('num_el_in_blk1', 'num_att_in_blk1',))
+
+ # Set variable values
+ connect1[:] = 1+cells[:]
+ connect1.elem_type = shape
+
+ coord[:] = vertices.transpose()[:]
+
+ from netCDF4 import stringtoarr
+ if dim == 2:
+ coor_names[0,:] = stringtoarr("x", len_string)
+ coor_names[1,:] = stringtoarr("y", len_string)
+ elif dim == 3:
+ coor_names[0,:] = stringtoarr("x", len_string)
+ coor_names[1,:] = stringtoarr("y", len_string)
+ coor_names[2,:] = stringtoarr("z", len_string)
+
+
+ qa_records[0,0,:] = stringtoarr("CUBIT", len_string)
+ qa_records[0,1,:] = stringtoarr("11.0", len_string)
+ qa_records[0,2,:] = stringtoarr("01/01/2000", len_string)
+ qa_records[0,3,:] = stringtoarr("12:00:00", len_string)
+
+ elem_map[:] = numpy.arange(1, ncells+1, dtype=numpy.int32)[:]
+
+ eb_status[:] = numpy.ones( (1,), dtype=numpy.int32)[:]
+
+ eb_prop1[:] = numpy.ones( (1,), dtype=numpy.int32)[:]
+ eb_prop1.name = "ID"
+
+ attrib1[:] = numpy.ones( (1, ncells), dtype=numpy.int32)[:]
+
+ root.close()
+
+
+# End of file
More information about the CIG-COMMITS
mailing list