[cig-commits] [commit] master: Created scripts to generate initial mesh. (385c7d3)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Feb 19 16:29:01 PST 2014
Repository : ssh://geoshell/pylith_benchmarks
On branch : master
Link : https://github.com/geodynamics/pylith_benchmarks/compare/997d255b545bda4a5ee4abf4d2d9ce5089eaacdb...385c7d351cee0569252f2d7ac9caaacc1d6b3888
>---------------------------------------------------------------
commit 385c7d351cee0569252f2d7ac9caaacc1d6b3888
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Wed Feb 19 16:29:04 2014 -0800
Created scripts to generate initial mesh.
>---------------------------------------------------------------
385c7d351cee0569252f2d7ac9caaacc1d6b3888
dynamic/scecdynrup/tpv28/adjust_faultgeom.py | 46 ++++++++++++++++++++++
dynamic/scecdynrup/tpv28/create_faultsurf.jou | 25 ++++++++++++
dynamic/scecdynrup/tpv28/create_faultsurf.py | 32 +++++++++------
dynamic/scecdynrup/{tpv102 => tpv28}/createbc.jou | 36 ++++++-----------
dynamic/scecdynrup/{tpv102 => tpv28}/geometry.jou | 35 ++++++----------
dynamic/scecdynrup/tpv28/gradient.jou | 24 +++++++++++
dynamic/scecdynrup/{tpv102 => tpv28}/tet4_200m.jou | 2 +-
7 files changed, 142 insertions(+), 58 deletions(-)
diff --git a/dynamic/scecdynrup/tpv28/adjust_faultgeom.py b/dynamic/scecdynrup/tpv28/adjust_faultgeom.py
new file mode 100755
index 0000000..fef7987
--- /dev/null
+++ b/dynamic/scecdynrup/tpv28/adjust_faultgeom.py
@@ -0,0 +1,46 @@
+#!/usr/bin/env python
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2014 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+# PREREQUISITES: numpy, netCDF4
+
+filenameEXO = "tet4_200m.exo"
+
+# ======================================================================
+import numpy
+import netCDF4
+
+from create_faultsurf import bump
+
+exodus = netCDF4.Dataset(filenameEXO, 'a')
+coordx = exodus.variables['coordx'][:]
+coordy = exodus.variables['coordy'][:]
+coordz = exodus.variables['coordz'][:]
+
+nodesetNames = exodus.variables['ns_names']
+for i in range(nodesetNames.shape[0]):
+ if netCDF4.chartostring(nodesetNames[i,:]) == "fault":
+ indices = exodus.variables['node_ns%d' % i][:]
+
+faulty = coordy[indices]
+faultz = coordz[indices]
+faultx = bump(faulty, faultz)
+coordx[indices] = faultx
+exodus.variables['coordx'][:] = coordx
+
+exodus.close()
+
+
+# End of file
diff --git a/dynamic/scecdynrup/tpv28/create_faultsurf.jou b/dynamic/scecdynrup/tpv28/create_faultsurf.jou
new file mode 100644
index 0000000..714f2b6
--- /dev/null
+++ b/dynamic/scecdynrup/tpv28/create_faultsurf.jou
@@ -0,0 +1,25 @@
+# -*- Python -*- (syntax highlighting)
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+#
+# ----------------------------------------------------------------------
+#
+# CUBIT journal file with geometry for TPV210.
+#
+# ----------------------------------------------------------------------
+# Set units to SI.
+# ----------------------------------------------------------------------
+#{Units('si')}
+#
+# ----------------------------------------------------------------------
+# Reset geometry.
+# ----------------------------------------------------------------------
+reset
+
+import mesh geometry "fault.exo" all
+create surface net from mapped surface 1 heal
+export acis "fault.sat" body 2 overwrite
+
+# End of file
+
diff --git a/dynamic/scecdynrup/tpv28/create_faultsurf.py b/dynamic/scecdynrup/tpv28/create_faultsurf.py
index 5585ec3..aa5f65b 100755
--- a/dynamic/scecdynrup/tpv28/create_faultsurf.py
+++ b/dynamic/scecdynrup/tpv28/create_faultsurf.py
@@ -9,7 +9,7 @@ filenameEXO = "fault.exo"
# Geometry of simulation domain
domainLength = 64.0e+3
domainHeight = 32.0e+3
-buffer = 2.0e+3 # Extend fault plane beyond domain
+buffer = 1.5e+3 # Extend fault plane beyond domain
# Geometry of bump
bumpRadius = 3.0e+3
@@ -26,7 +26,7 @@ from cubit_io import write_exodus_file
from netCDF4 import Dataset
# ----------------------------------------------------------------------
-def bump(yz):
+def bump(y, z):
y = yz[:,0]
z = yz[:,1]
r1 = ((y+bumpY)**2 + (z-bumpZ)**2)**0.5
@@ -40,17 +40,27 @@ def bump(yz):
# ----------------------------------------------------------------------
-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)
+y0 = -0.5*domainLength-buffer; y1 = -bumpY-bumpRadius; ny = int(1+(y1-y0)/(8.0*bumpDx))
+yA = numpy.linspace(y0, y1, ny)
+y0 = y1 + bumpDx; y1 = -bumpY+bumpRadius
+yB = numpy.arange(y0, y1, bumpDx)
+
+y0 = y1+bumpDx; y1 = +bumpY-bumpRadius; ny = int(1+(y1-y0)/(8.0*bumpDx))
+yC = numpy.linspace(y0, y1, ny)
+y0 = y1 + bumpDx; y1 = +bumpY+bumpRadius
+yD = numpy.arange(y0, y1, bumpDx)
+
+y0 = y1+bumpDx; y1 = +0.5*domainLength+buffer; ny = int(1+(y1-y0)/(8.0*bumpDx))
+yE = numpy.linspace(y0, y1, ny)
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)
+z0 = -domainHeight-buffer; z1 = bumpZ-bumpRadius; nz = (1+(z1-z0)/(8.0*bumpDx))
+zA = numpy.linspace(z0, z1, nz)
+zB = numpy.arange(z1+bumpDx, bumpZ+bumpRadius, bumpDx)
+
+z0 = bumpZ+bumpRadius+bumpDx; z1 = +buffer; nz = int(1+(z1-z0)/(8.0*bumpDx))
+zC = numpy.linspace(z0, z1, nz)
z = numpy.hstack((zA, zB, zC))
numZ = z.shape[0]
@@ -60,7 +70,7 @@ for iY in xrange(numY):
yz[iY*numZ:(iY+1)*numZ,1] = z
yz[iY*numZ:(iY+1)*numZ,0] = y[iY]
-x = bump(yz)
+x = bump(yz[:,0],yz[:,1])
vertices = numpy.zeros((numY*numZ,3), dtype=numpy.float64)
vertices[:,0] = x
vertices[:,1] = yz[:,0]
diff --git a/dynamic/scecdynrup/tpv102/createbc.jou b/dynamic/scecdynrup/tpv28/createbc.jou
similarity index 66%
copy from dynamic/scecdynrup/tpv102/createbc.jou
copy to dynamic/scecdynrup/tpv28/createbc.jou
index c16ed55..4dfec67 100644
--- a/dynamic/scecdynrup/tpv102/createbc.jou
+++ b/dynamic/scecdynrup/tpv28/createbc.jou
@@ -1,7 +1,7 @@
# ----------------------------------------------------------------------
# Create blocks for materials
# ----------------------------------------------------------------------
-block 1 volume 1 3
+block 1 volume domain domain at A
block 1 name "elastic"
# ----------------------------------------------------------------------
@@ -28,16 +28,16 @@ nodeset 11 name "face_xneg"
# ----------------------------------------------------------------------
# Create nodeset for +y face
# ----------------------------------------------------------------------
-group "face_ypos" add node in surface 12
-group "face_ypos" add node in surface 14
+group "face_ypos" add node in surface 10
+group "face_ypos" add node in surface 17
nodeset 12 group face_ypos
nodeset 12 name "face_ypos"
# ----------------------------------------------------------------------
# Create nodeset for -y face
# ----------------------------------------------------------------------
-group "face_yneg" add node in surface 10
-group "face_yneg" add node in surface 17
+group "face_yneg" add node in surface 12
+group "face_yneg" add node in surface 15
nodeset 13 group face_yneg
nodeset 13 name "face_yneg"
@@ -45,28 +45,16 @@ nodeset 13 name "face_yneg"
# Create nodeset for -z face
# ----------------------------------------------------------------------
group "face_zneg" add node in surface 11
-group "face_zneg" add node in surface 15
-nodeset 15 group face_zneg
-nodeset 15 name "face_zneg"
+group "face_zneg" add node in surface 16
+nodeset 14 group face_zneg
+nodeset 14 name "face_zneg"
# ----------------------------------------------------------------------
# Create nodeset for +z face
# ----------------------------------------------------------------------
group "face_zpos" add node in surface 9
-group "face_zpos" add node in surface 16
-nodeset 16 group face_zpos
-nodeset 16 name "face_zpos"
-
-# ----------------------------------------------------------------------
-# Create nodeset for subset of +z face
-# ----------------------------------------------------------------------
-group "face_zpos_subset" add node in group face_zpos
-group "face_zpos_subset" remove node with y_coord > +15.001e+3
-group "face_zpos_subset" remove node with y_coord < -15.001e+3
-group "face_zpos_subset" remove node with x_coord > +10.001e+3
-group "face_zpos_subset" remove node with x_coord < -10.001e+3
-nodeset 17 group face_zpos_subset
-nodeset 17 name "face zpos_subset"
-
-
+group "face_zpos" add node in surface 14
+nodeset 15 group face_zpos
+nodeset 15 name "face_zpos"
+# End of file
diff --git a/dynamic/scecdynrup/tpv102/geometry.jou b/dynamic/scecdynrup/tpv28/geometry.jou
similarity index 70%
copy from dynamic/scecdynrup/tpv102/geometry.jou
copy to dynamic/scecdynrup/tpv28/geometry.jou
index 452b4cd..26cf0e9 100644
--- a/dynamic/scecdynrup/tpv102/geometry.jou
+++ b/dynamic/scecdynrup/tpv28/geometry.jou
@@ -30,51 +30,42 @@ reset
#{blockHeight=32.0*km}
#
#{faultLength=36.0*km}
-#{faultWidth=18.0*km}
+#{faultWidth=15.0*km}
+#{buffer=1.5*km}
brick x {blockWidth} y {blockLength} z {blockHeight}
# Translate block so the top is at z=0
volume 1 move x 0 y 0 z {-0.5*blockHeight}
+volume {Id("volume")} name "domain"
# ----------------------------------------------------------------------
# Create interface surfaces
# ----------------------------------------------------------------------
-create planar surface with plane xplane offset 0
-#{sI=Id("surface")}
+import acis "fault.sat"
+#{vFault=Id("body")}
# ----------------------------------------------------------------------
# Divide volumes using interface surfaces
# ----------------------------------------------------------------------
-webcut volume 1 with plane surface {sI}
+webcut volume domain tool volume {vFault}
# ----------------------------------------------------------------------
# Split fault surface
# ----------------------------------------------------------------------
-create vertex 0 {-0.5*faultLength} 0
-#{vAn=Id("vertex")}
-create vertex 0 {+0.5*faultLength} 0
-#{vAp=Id("vertex")}
-create vertex 0 {-0.5*faultLength} {-blockHeight}
-#{vBn=Id("vertex")}
-create vertex 0 {+0.5*faultLength} {-blockHeight}
-#{vBp=Id("vertex")}
-create vertex 0 {-0.5*faultLength} {-faultWidth}
-#{vCn=Id("vertex")}
-create vertex 0 {+0.5*faultLength} {-faultWidth}
-#{vCp=Id("vertex")}
-split surface 7 across location vertex {vAn} location vertex {vBn}
-split surface 18 across location vertex {vAp} location vertex {vBp}
-split surface 21 across location vertex {vCn} location vertex {vCp}
-surface 22 name "fault_surface"
-delete vertex all
+brick x {3.0*km} y {faultLength} z {faultWidth+2*buffer}
+{vChopper=Id("volume")}
+volume {vChopper} move x 0 y 0 z {-0.5*faultWidth+buffer}
+chop volume {vFault} with volume {vChopper} keep
+delete volume {vChopper}
# ----------------------------------------------------------------------
# Imprint all volumes, merging surfaces
# ----------------------------------------------------------------------
imprint all
merge all
-delete body 2
+delete body {vFault} 5 6
+surface 26 name "fault_surface"
diff --git a/dynamic/scecdynrup/tpv28/gradient.jou b/dynamic/scecdynrup/tpv28/gradient.jou
new file mode 100644
index 0000000..6a92fd6
--- /dev/null
+++ b/dynamic/scecdynrup/tpv28/gradient.jou
@@ -0,0 +1,24 @@
+# ----------------------------------------------------------------------
+# Set vertex spacing with increasing spacing away from fault
+# ----------------------------------------------------------------------
+#{bias_factor=1.02}
+
+# ----------------------------------------------------------------------
+# Set sizes
+# ----------------------------------------------------------------------
+curve all scheme default
+surface all sizing function none
+curve all sizing function none
+surface fault_surface size {dx}
+
+surface 27 sizing function type bias start curve 57 46 59 factor {bias_factor}
+surface 9 14 sizing function type bias start curve 60 factor {bias_factor}
+surface 12 15 sizing function type bias start curve 20 factor {bias_factor}
+surface 10 17 sizing function type bias start curve 18 factor {bias_factor}
+surface 11 16 sizing function type bias start curve 17 factor {bias_factor}
+
+# End of file
+
+
+
+
diff --git a/dynamic/scecdynrup/tpv102/tet4_200m.jou b/dynamic/scecdynrup/tpv28/tet4_200m.jou
similarity index 99%
copy from dynamic/scecdynrup/tpv102/tet4_200m.jou
copy to dynamic/scecdynrup/tpv28/tet4_200m.jou
index c60434d..7a095c2 100644
--- a/dynamic/scecdynrup/tpv102/tet4_200m.jou
+++ b/dynamic/scecdynrup/tpv28/tet4_200m.jou
@@ -45,6 +45,6 @@ playback 'createbc.jou'
# ----------------------------------------------------------------------
export mesh "tet4_200m.exo" dimension 3 overwrite
-
+q
# End of file
More information about the CIG-COMMITS
mailing list