[cig-commits] r19608 - in seismo/3D/SPECFEM3D/trunk: examples/Mount_StHelens src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Feb 10 06:17:27 PST 2012
Author: danielpeter
Date: 2012-02-10 06:17:26 -0800 (Fri, 10 Feb 2012)
New Revision: 19608
Added:
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount_stl.py
Modified:
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py
seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
Log:
fixes sign() function for double precision value; updates mesh_mount.py for Mount St.Helens example
Modified: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py 2012-02-09 19:25:48 UTC (rev 19607)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount.py 2012-02-10 14:17:26 UTC (rev 19608)
@@ -1,5 +1,13 @@
#!/usr/bin/env python
-
+#############################################################
+#
+# script uses ACIS surface formats
+#
+# note: this script seems to work with CUBIT version > 12.2
+# meshing takes about 15 minutes (without refinement)
+#
+#
+#############################################################
import cubit
import boundary_definition
import cubit2specfem3d
@@ -7,52 +15,46 @@
import os
import sys
import os.path
+import time
+# time stamp
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+
+# working directory
+cwd = os.getcwd()
+print "#current working directory: " + str(cwd)
+if cwd[len(cwd)-14:len(cwd)] != "Mount_StHelens":
+ print ""
+ print "#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/"
+ print ""
+
+cubit.cmd('version')
cubit.cmd('reset')
-os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
print "running meshing script..."
print ""
-print "note: this script uses topography surface in STL format"
-print " meshing will take around 2 min"
+print "note: this script uses topography surface in ACIS format"
+print " meshing will take around 15 min"
print ""
-# note: this is a workaround to use STL file formats rather than ACIS formats.
-# for our purpose to create a simple mesh, this STL formats are faster
+# uses developer commands
+cubit.cmd('set developer commands on')
+cubit.cmd('set import mesh tolerance 1')
#############################################################
#
-# 1. step: creates temporary brick volume in STL format
+# 0. step: loading topography surface
#
#############################################################
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#loading topo surface..."
-# creates temporary brick volume and exports in STL format
-# new brick volume (depth will become 1/2 * 20,000 = 10,000 m)
-cubit.cmd('brick x 15000 y 22000 z 20000')
-# moves volume to UTM coordinates of topography surface
-cubit.cmd('volume 1 move x 561738. y 5116370. z 0 ')
-# saves as STL body
-cubit.cmd('export stl ascii "topo_brick.stl" overwrite')
-# clear
-cubit.cmd('reset')
-#checks if new file available
-if not os.path.exists("topo_brick.stl"):
- print ""
- print "error creating new STL file topo_brick.stl, please check manually..."
- print ""
- cubit.cmd('pause')
-
-#############################################################
-#
-# 2. step: imports topography surface and brick volume in STL format
-#
-#############################################################
-
# topography surface
-if os.path.exists("topo.stl"):
+if os.path.exists("topo.cub"):
print "opening existing topography surface"
+ # topography surface
# previously run, just reopen the cubit file
- cubit.cmd('import stl "topo.stl" merge stitch')
+ cubit.cmd('open "topo.cub"')
else:
# topo surface doesn't exist yet, this creates it:
print "reading in topography surface"
@@ -61,73 +63,146 @@
# clear
cubit.cmd('reset')
# now reopen the cubit file
- cubit.cmd('import stl "topo.stl" merge stitch')
-# re-opens brick in STL format
-cubit.cmd('import stl "topo_brick.stl" merge stitch')
-# imprints topography surfaces into brick volume, creates connected surfaces
-# note: this is a develop feature
-cubit.cmd('set developer commands on')
+ cubit.cmd('open "topo.cub"')
+
+# healing the surface...
+cubit.cmd('Auto_clean volume 1 small_surfaces small_curve_size 10')
+cubit.cmd('regularize volume 1')
+
+#############################################################
+#
+# 1. step: creates temporary brick volume
+#
+#############################################################
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#creating brick..."
+
+# creates temporary brick volume
+# new brick volume (depth will become 1/2 * 20,000 = 10,000 m)
+cubit.cmd('brick x 15000 y 22000 z 20000')
+
+# moves volume to UTM coordinates of topography surface
+cubit.cmd('volume 2 move x 561738. y 5116370. z 0 ')
+
+# temporary backup
+cubit.cmd('save as "topo_1.cub" overwrite')
+
+#############################################################
+#
+# 2. step: creates volume with topography surface
+#
+#############################################################
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#creating volume with topography..."
+
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#imprinting volume, this will take around 1 min, please be patience..."
+cubit.cmd('merge all')
cubit.cmd('imprint all')
-# exports as STL only surfaces which will create new volume
-cubit.cmd('export stl ascii "topo_vol.stl" surface 9 5 6 8 11 13 overwrite')
-# clear
-cubit.cmd('reset')
-#checks if new file available
-if not os.path.exists("topo_vol.stl"):
- print ""
- print "error creating new STL file topo_vol.stl, please check manually..."
- print ""
- cubit.cmd('pause')
+# exports only surfaces which will create single volume
+cubit.cmd('export acis "topo_2.acis" surface 3 10 12 14 15 9 ascii overwrite')
+
+# backup
+cubit.cmd('save as "topo_2.cub" overwrite')
+
#############################################################
#
-# 3. step: manipulate STL file to create a single volume
+# 3. step: manipulate ACIS file to create a single volume
#
#############################################################
-
-os.system('awk \'BEGIN{print \"solid Body_1\";}{if($0 !~ /solid/) print $0;}END{print \"endsolid Body_1\";}\' topo_vol.stl > topo_vol2.stl')
-#checks if new file available
-if not os.path.exists("topo_vol2.stl"):
+# checks if new file available
+if not os.path.exists("topo_2.acis"):
print ""
- print "error creating new STL file topo_vol2.stl, please check manually..."
+ print "error creating new volume, please check manually..."
print ""
cubit.cmd('pause')
-
+# clears workspace
+cubit.cmd('reset')
+
+# imports surfaces and merges to single volume
+cubit.cmd('import acis "topo_2.acis" ascii merge_globally')
+
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#creating new volume, this will take another 2 min..."
+cubit.cmd('create volume surface all heal')
+
+# backup
+cubit.cmd('save as "topo_3.cub" overwrite')
+
#############################################################
#
-# 4. step: import STL volume and create mesh
+# 4. step: create mesh
#
#############################################################
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#initial meshing..."
+print "#(will take around 7 min)"
-# re-opens new volume in STL format
-cubit.cmd('import stl "topo_vol2.stl" merge stitch')
+# optional: refining mesh at surface
+#
+# note: refining using ACIS surface format takes a long time ... (up to 3 hours)
+DO_TOPO_REFINEMENT = False
+
# Meshing the volumes
-elementsize = 2000.0
-cubit.cmd('volume all size '+str(elementsize))
-# sets meshing type
-# uses a sweep algorithm in vertical (Z) direction
-cubit.cmd('volume all scheme sweep Vector 0 0 1')
-# initial coarse mesh
-cubit.cmd('mesh volume all')
-# draw/update mesh lines for visualization
-# this will draw also the tripling layer mesh lines
-cubit.cmd('draw volume all')
-# optional smoothing to improve mesh quality (takes up to 5 min)
-cubit.cmd('volume all smooth scheme condition number beta 1.0 cpu 5')
-cubit.cmd('smooth volume all')
-# refines global mesh
-cubit.cmd('refine volume 1 numsplit 1')
-cubit.cmd('draw volume all')
-# optional smoothing
-cubit.cmd('smooth volume all')
-# refines elements at topography surface
-cubit.cmd('refine surface 2 numsplit 1 bias 1.0 depth 1')
-cubit.cmd('draw volume all')
-# optional smoothing to improve mesh quality (takes up all 10 min if beta is 2.0)
-cubit.cmd('volume all smooth scheme condition number beta 2.3 cpu 10')
-cubit.cmd('smooth volume all')
-cubit.cmd('draw volume all')
+if DO_TOPO_REFINEMENT == False:
+ elementsize = 500.0
+ cubit.cmd('volume all size '+str(elementsize))
+ # note: we will mesh first the topography surface, then sweep down the mesh
+ # topography surface
+ #cubit.cmd('control skew surface 12')
+ cubit.cmd('surface 12 submap smooth off')
+ cubit.cmd('surface 12 scheme submap')
+ cubit.cmd('mesh surface 12')
+ # propagates mesh down for whole volume
+ cubit.cmd('volume 1 redistribute nodes off')
+ cubit.cmd('volume 1 scheme sweep source surface 12 target surface 7')
+ cubit.cmd('mesh volume all')
+ # draw/update mesh lines for visualization
+ # this will draw also the tripling layer mesh lines in case
+ cubit.cmd('draw volume all')
+else:
+ # optional surface refinement
+ # time stamp
+ print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+ print "#refining surface mesh..."
+ print "#(will take around 3 hours)"
+ # starts with a crude mesh
+ elementsize = 2000.0
+ cubit.cmd('volume all size '+str(elementsize))
+ # sets meshing type
+ # explicitly sets scheme for topography surface
+ cubit.cmd('surface 12 submap smooth off')
+ cubit.cmd('surface 12 scheme submap')
+ # uses a sweep algorithm in vertical (Z) direction
+ cubit.cmd('volume all scheme sweep Vector 0 0 1')
+ # initial coarse mesh
+ cubit.cmd('mesh volume all')
+ # optional smoothing to improve mesh quality (takes up to 5 min)
+ cubit.cmd('volume all smooth scheme condition number beta 1.0 cpu 5')
+ cubit.cmd('smooth volume all')
+ # refines global mesh
+ cubit.cmd('refine volume 1 numsplit 1')
+ cubit.cmd('draw volume all')
+ # optional smoothing
+ cubit.cmd('smooth volume all')
+ # refines elements at topography surface
+ cubit.cmd('refine surface 12 numsplit 1 bias 1.0 depth 1')
+ cubit.cmd('draw volume all')
+ # optional smoothing to improve mesh quality (takes up to 10 min)
+ cubit.cmd('volume all smooth scheme condition number beta 2.3 cpu 10')
+ cubit.cmd('smooth volume all')
+ # displays final mesh
+ cubit.cmd('draw volume all')
+
+# time stamp
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+print "#done meshing..."
+
+# backup
+cubit.cmd('save as "topo_4.cub" overwrite')
+
#### End of meshing
###### This is boundary_definition.py of GEOCUBIT
@@ -151,8 +226,6 @@
# optional saves backups
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
-# cleanup
-os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
@@ -160,3 +233,6 @@
cubit2specfem3d.export2SESAME('MESH')
# all files needed by SCOTCH are now in directory MESH
+
+# time stamp
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
Added: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount_stl.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount_stl.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount_stl.py 2012-02-10 14:17:26 UTC (rev 19608)
@@ -0,0 +1,198 @@
+#!/usr/bin/env python
+#############################################################
+#
+# script uses STL surface formats
+#
+# note: this script seems only to work with CUBIT version 12.2
+#
+# ( try out script mesh_mount.py when using a different CUBIT version)
+#
+#############################################################
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+import os.path
+
+# time stamp
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
+
+# working directory
+cwd = os.getcwd()
+print "#current working directory: " + str(cwd)
+if cwd[len(cwd)-14:len(cwd)] != "Mount_StHelens":
+ print ""
+ print "#please run this script from example directory: SPECFEM3D/example/Mount_StHelens/"
+ print ""
+
+cubit.cmd('version')
+cubit.cmd('reset')
+
+os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
+
+print "running meshing script..."
+print ""
+print "note: this script uses topography surface in STL format"
+print " meshing will take around 2 min"
+print ""
+
+# note: this is a workaround to use STL file formats rather than ACIS formats.
+# for our purpose to create a simple mesh, this STL formats are faster
+
+# uses developer commands
+cubit.cmd('set developer commands on')
+cubit.cmd('set import mesh tolerance 100')
+
+
+#############################################################
+#
+# 1. step: creates temporary brick volume in STL format
+#
+#############################################################
+
+# creates temporary brick volume and exports in STL format
+# new brick volume (depth will become 1/2 * 20,000 = 10,000 m)
+cubit.cmd('brick x 15000 y 22000 z 20000')
+# moves volume to UTM coordinates of topography surface
+cubit.cmd('volume 1 move x 561738. y 5116370. z 0 ')
+# saves as STL body
+cubit.cmd('export stl ascii "topo_brick.stl" overwrite')
+# clear
+cubit.cmd('reset')
+#checks if new file available
+if not os.path.exists("topo_brick.stl"):
+ print ""
+ print "error creating new STL file topo_brick.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 2. step: imports topography surface and brick volume in STL format
+#
+#############################################################
+
+# topography surface
+if os.path.exists("topo.stl"):
+ print "opening existing topography surface"
+ # previously run, just reopen the cubit file
+ cubit.cmd('import stl "topo.stl" merge stitch')
+else:
+ # topo surface doesn't exist yet, this creates it:
+ print "reading in topography surface"
+ # reads in topography points and creates sheet surface
+ execfile("./read_topo.py")
+ # clear
+ cubit.cmd('reset')
+ # now reopen the cubit file
+ cubit.cmd('import stl "topo.stl" merge stitch')
+
+# re-opens brick in STL format
+cubit.cmd('import stl "topo_brick.stl" merge stitch')
+# imprints topography surfaces into brick volume, creates connected surfaces
+# note: this is a develop feature
+cubit.cmd('imprint all')
+
+
+cubit.cmd('version')
+
+# cubit 12.2 specific
+# exports as STL only surfaces which will create new volume
+cubit.cmd('export stl ascii "topo_vol.stl" surface 9 5 6 8 11 13 overwrite')
+
+# clear
+cubit.cmd('reset')
+#checks if new file available
+if not os.path.exists("topo_vol.stl"):
+ print ""
+ print "error creating new STL file topo_vol.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 3. step: manipulate STL file to create a single volume
+#
+#############################################################
+
+os.system('awk \'BEGIN{print \"solid Body_1\";}{if($0 !~ /solid/) print $0;}END{print \"endsolid Body_1\";}\' topo_vol.stl > topo_vol2.stl')
+#checks if new file available
+if not os.path.exists("topo_vol2.stl"):
+ print ""
+ print "error creating new STL file topo_vol2.stl, please check manually..."
+ print ""
+ cubit.cmd('pause')
+
+#############################################################
+#
+# 4. step: import STL volume and create mesh
+#
+#############################################################
+
+# re-opens new volume in STL format
+cubit.cmd('import stl "topo_vol2.stl" merge stitch')
+
+# Meshing the volumes
+elementsize = 2000.0
+cubit.cmd('volume all size '+str(elementsize))
+# sets meshing type
+# uses a sweep algorithm in vertical (Z) direction
+cubit.cmd('volume all scheme sweep Vector 0 0 1')
+# initial coarse mesh
+cubit.cmd('mesh volume all')
+# draw/update mesh lines for visualization
+# this will draw also the tripling layer mesh lines
+cubit.cmd('draw volume all')
+# optional smoothing to improve mesh quality (takes up to 5 min)
+cubit.cmd('volume all smooth scheme condition number beta 1.0 cpu 5')
+cubit.cmd('smooth volume all')
+# refines global mesh
+cubit.cmd('refine volume 1 numsplit 1')
+cubit.cmd('draw volume all')
+# optional smoothing
+cubit.cmd('smooth volume all')
+# refines elements at topography surface
+cubit.cmd('refine surface 2 numsplit 1 bias 1.0 depth 1')
+cubit.cmd('draw volume all')
+# optional smoothing to improve mesh quality (takes up all 10 min if beta is 2.0)
+cubit.cmd('volume all smooth scheme condition number beta 2.3 cpu 10')
+cubit.cmd('smooth volume all')
+cubit.cmd('draw volume all')
+
+#### End of meshing
+
+###### This is boundary_definition.py of GEOCUBIT
+#..... which extracts the bounding faces and defines them into blocks
+boundary_definition.entities=['face']
+boundary_definition.define_bc(boundary_definition.entities,parallel=True)
+
+#### Define material properties for the 3 volumes ################
+cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
+
+
+cubit.cmd('block 1 name "elastic 1" ') # elastic material region
+cubit.cmd('block 1 attribute count 6')
+cubit.cmd('block 1 attribute index 1 1') # flag for material: 1 for 1. material
+cubit.cmd('block 1 attribute index 2 2800') # vp
+cubit.cmd('block 1 attribute index 3 1500') # vs
+cubit.cmd('block 1 attribute index 4 2300') # rho
+cubit.cmd('block 1 attribute index 5 150.0') # Qmu
+cubit.cmd('block 1 attribute index 6 0 ') # anisotropy_flag
+
+# optional saves backups
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+# cleanup
+os.system('rm -f topo_brick.stl topo_vol.stl topo_vol2.stl')
+
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+
+os.system('mkdir -p MESH')
+cubit2specfem3d.export2SESAME('MESH')
+
+# all files needed by SCOTCH are now in directory MESH
+
+# time stamp
+print "#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/mesh_mount_stl.py
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py 2012-02-09 19:25:48 UTC (rev 19607)
+++ seismo/3D/SPECFEM3D/trunk/examples/Mount_StHelens/read_topo.py 2012-02-10 14:17:26 UTC (rev 19608)
@@ -49,6 +49,7 @@
print '#done points: '+str(count)
print ''
+cubit.cmd('Display')
# creates smooth spline curves for surface
print '#creating curves...'
@@ -57,14 +58,21 @@
if i > 1 :
if i % nstep == 0 :
countcurves = countcurves + 1
- cubit.cmd('create curve spline vertex '+str(i-nstep+1) + ' to ' + str(i) )
+ cubit.cmd('create curve spline vertex '+str(i-nstep+1) + ' to ' + str(i) + ' delete' )
print '#done curves: '+str(countcurves)
print ''
+cubit.cmd('Display')
+cubit.cmd('pause')
+
# creates surface
print '#creating skin surface...'
cubit.cmd('create surface skin curve all')
+
+cubit.cmd('Display')
+cubit.cmd('pause')
+
# cleans up
cubit.cmd('merge all ')
cubit.cmd('delete vertex all')
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-02-09 19:25:48 UTC (rev 19607)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-02-10 14:17:26 UTC (rev 19608)
@@ -570,7 +570,7 @@
! determines factor +/-1 depending on sign of moment tensor
! (see e.g. Krebs et al., 2009. Fast full-wavefield seismic inversion using encoded sources,
! Geophysics, 74 (6), WCC177-WCC188.)
- pm1_source_encoding(isource) = sign(1.0,Mxx(isource))
+ pm1_source_encoding(isource) = sign(1.0d0,Mxx(isource))
! source array interpolated on all element gll points (only used for non point sources)
call compute_arrays_source_acoustic(xi_source(isource),eta_source(isource),gamma_source(isource),&
@@ -748,6 +748,7 @@
end subroutine setup_receivers_precompute_intp
+
!
!-------------------------------------------------------------------------------------------------
!
More information about the CIG-COMMITS
mailing list