[cig-commits] r16024 - in seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES: . waterlayered_halfspace
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Nov 20 18:26:52 PST 2009
Author: danielpeter
Date: 2009-11-20 18:26:51 -0800 (Fri, 20 Nov 2009)
New Revision: 16024
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/CMTSOLUTION
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/Par_file
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/STATIONS
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/boundary_definition.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8-nodoubling.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py
Log:
added acoustic examples in folder EXAMPLES/waterlayered_halfspace
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/CMTSOLUTION
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/CMTSOLUTION (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/CMTSOLUTION 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,13 @@
+PDE 1999 01 01 00 00 00.00 67000 67000 -25000 4.2 4.2 FIG8
+event name: FIG8_vertical
+time shift: 0.0000
+half duration: 2.0
+latitude: 67000.0
+longitude: 67000.0
+depth: -25000.0
+Mrr: 1.000000e+23
+Mtt: 1.000000e+23
+Mpp: 1.000000e+23
+Mrt: 0.000000
+Mrp: 0.000000
+Mtp: 0.000000
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/Par_file (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/Par_file 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,52 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
+SAVE_FORWARD = .false.
+
+# UTM projection parameters
+UTM_PROJECTION_ZONE = 11
+SUPPRESS_UTM_PROJECTION = .true.
+
+# number of MPI processors along xi and eta (can be different)
+NPROC_XI = 2
+NPROC_ETA = 2
+
+# time step parameters
+NSTEP = 4500
+DT = 0.005d0
+
+# parameters describing the model
+OCEANS = .false.
+TOPOGRAPHY = .false.
+ATTENUATION = .false.
+USE_OLSEN_ATTENUATION = .false.
+ANISOTROPY = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .false.
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .false.
+MOVIE_VOLUME = .false.
+NTSTEP_BETWEEN_FRAMES = 200
+CREATE_SHAKEMAP = .false.
+SAVE_DISPLACEMENT = .false.
+USE_HIGHRES_FOR_MOVIES = .false.
+HDUR_MOVIE = 0.0
+
+# save AVS or OpenDX mesh files to check the mesh
+SAVE_MESH_FILES = .false.
+
+# path to store the local database file on each node
+LOCAL_PATH = DATABASES_MPI
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 500
+
+# interval in time steps for writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 10000
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/README 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,59 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+
+step-by-step tutorial:
+
+1. create mesh:
+
+ - run cubit GUI:
+ > claro
+
+ then run meshing script:
+ claro -> Menu "Tools" -> "Play Journal File" ... and select file: "waterlayer_mesh_boundary_fig8.py"
+
+ (or select: "waterlayer_mesh_boundary_fig8-nodoubling.py" for regular mesh without doubling layer)
+
+ this creates all the mesh files in subdirectory MESH/
+
+
+2. decompose mesh files:
+
+ - run decomposer in directory decompose_mesh_SCOTCH/:
+ (example assumes 4 partitions with mesh files in OUTPUT_FILES/)
+
+ > make
+ > ./xdecompose_mesh_SCOTCH 4 ../EXAMPLES/waterlayered_halfspace/MESH/ ../DATABASES_MPI/
+
+ which creates mesh partitions "proc0000***_Database" in directory "DATABASES_MPI".
+ you can then specify "DATABASES_MPI" in "Par_file" for your "LOCAL_PATH"
+
+
+3. generate databases:
+
+ - copy Par_file, CMTSOLUTION & STATIONS files provided in this examples directory to DATA/ directory in
+ the SPECFEM3D_SESAME root directory
+
+ - compile generate_databases:
+ > make xgenerate_databases
+
+ - submit job script:
+ > qsub go_generate_database_pbs.sesame.bash
+
+ (note: if execution fails due to memory shortage - most likely fails when calling routine to
+ create regional mesh - then try to increase process memory stack size: ulimit -s 2000000 (2GB) )
+
+
+4. run simulation:
+
+ - compile specfem3D:
+ > make xspecfem3D
+
+ - submit job script:
+ > qsub go_solver_pbs.sesame.bash
+
+
+
+
+
\ No newline at end of file
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/STATIONS
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/STATIONS (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/STATIONS 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,230 @@
+Z1 DB 67000.00 67000.00 0.0 -1000
+Z2 DB 67000.00 67000.00 0.0 -3000.0
+Z3 DB 67000.00 67000.00 0.0 -6000.0
+Z4 DB 67000.00 67000.00 0.0 -54000.0
+Z5 DB 67000.00 67000.00 0.0 -57000.0
+Z6 DB 67000.00 67000.00 0.0 -60000.0
+X1 DB 67000.00 0.000000 0.0 -1000
+X2 DB 67000.00 1196.429 0.0 -1000
+X3 DB 67000.00 2392.857 0.0 -1000
+X4 DB 67000.00 3589.286 0.0 -1000
+X5 DB 67000.00 4785.714 0.0 -1000
+X6 DB 67000.00 5982.143 0.0 -1000
+X7 DB 67000.00 7178.571 0.0 -1000
+X8 DB 67000.00 8375.000 0.0 -1000
+X9 DB 67000.00 9571.429 0.0 -1000
+X10 DB 67000.00 10767.86 0.0 -1000
+X11 DB 67000.00 11964.29 0.0 -1000
+X12 DB 67000.00 13160.71 0.0 -1000
+X13 DB 67000.00 14357.14 0.0 -1000
+X14 DB 67000.00 15553.57 0.0 -1000
+X15 DB 67000.00 16750.00 0.0 -1000
+X16 DB 67000.00 17946.43 0.0 -1000
+X17 DB 67000.00 19142.86 0.0 -1000
+X18 DB 67000.00 20339.29 0.0 -1000
+X19 DB 67000.00 21535.71 0.0 -1000
+X20 DB 67000.00 22732.14 0.0 -1000
+X21 DB 67000.00 23928.57 0.0 -1000
+X22 DB 67000.00 25125.00 0.0 -1000
+X23 DB 67000.00 26321.43 0.0 -1000
+X24 DB 67000.00 27517.86 0.0 -1000
+X25 DB 67000.00 28714.29 0.0 -1000
+X26 DB 67000.00 29910.71 0.0 -1000
+X27 DB 67000.00 31107.14 0.0 -1000
+X28 DB 67000.00 32303.57 0.0 -1000
+X29 DB 67000.00 33500.00 0.0 -1000
+X30 DB 67000.00 34696.43 0.0 -1000
+X31 DB 67000.00 35892.86 0.0 -1000
+X32 DB 67000.00 37089.29 0.0 -1000
+X33 DB 67000.00 38285.71 0.0 -1000
+X34 DB 67000.00 39482.14 0.0 -1000
+X35 DB 67000.00 40678.57 0.0 -1000
+X36 DB 67000.00 41875.00 0.0 -1000
+X37 DB 67000.00 43071.43 0.0 -1000
+X38 DB 67000.00 44267.86 0.0 -1000
+X39 DB 67000.00 45464.29 0.0 -1000
+X40 DB 67000.00 46660.71 0.0 -1000
+X41 DB 67000.00 47857.14 0.0 -1000
+X42 DB 67000.00 49053.57 0.0 -1000
+X43 DB 67000.00 50250.00 0.0 -1000
+X44 DB 67000.00 51446.43 0.0 -1000
+X45 DB 67000.00 52642.86 0.0 -1000
+X46 DB 67000.00 53839.29 0.0 -1000
+X47 DB 67000.00 55035.71 0.0 -1000
+X48 DB 67000.00 56232.14 0.0 -1000
+X49 DB 67000.00 57428.57 0.0 -1000
+X50 DB 67000.00 58625.00 0.0 -1000
+X51 DB 67000.00 59821.43 0.0 -1000
+X52 DB 67000.00 61017.86 0.0 -1000
+X53 DB 67000.00 62214.29 0.0 -1000
+X54 DB 67000.00 63410.71 0.0 -1000
+X55 DB 67000.00 64607.14 0.0 -1000
+X56 DB 67000.00 65803.57 0.0 -1000
+X57 DB 67000.00 67000.00 0.0 -1000
+X58 DB 67000.00 68196.43 0.0 -1000
+X59 DB 67000.00 69392.86 0.0 -1000
+X60 DB 67000.00 70589.29 0.0 -1000
+X61 DB 67000.00 71785.71 0.0 -1000
+X62 DB 67000.00 72982.14 0.0 -1000
+X63 DB 67000.00 74178.57 0.0 -1000
+X64 DB 67000.00 75375.00 0.0 -1000
+X65 DB 67000.00 76571.43 0.0 -1000
+X66 DB 67000.00 77767.86 0.0 -1000
+X67 DB 67000.00 78964.29 0.0 -1000
+X68 DB 67000.00 80160.71 0.0 -1000
+X69 DB 67000.00 81357.14 0.0 -1000
+X70 DB 67000.00 82553.57 0.0 -1000
+X71 DB 67000.00 83750.00 0.0 -1000
+X72 DB 67000.00 84946.43 0.0 -1000
+X73 DB 67000.00 86142.86 0.0 -1000
+X74 DB 67000.00 87339.29 0.0 -1000
+X75 DB 67000.00 88535.71 0.0 -1000
+X76 DB 67000.00 89732.14 0.0 -1000
+X77 DB 67000.00 90928.57 0.0 -1000
+X78 DB 67000.00 92125.00 0.0 -1000
+X79 DB 67000.00 93321.43 0.0 -1000
+X80 DB 67000.00 94517.86 0.0 -1000
+X81 DB 67000.00 95714.29 0.0 -1000
+X82 DB 67000.00 96910.71 0.0 -1000
+X83 DB 67000.00 98107.14 0.0 -1000
+X84 DB 67000.00 99303.57 0.0 -1000
+X85 DB 67000.00 100500.0 0.0 -1000
+X86 DB 67000.00 101696.4 0.0 -1000
+X87 DB 67000.00 102892.9 0.0 -1000
+X88 DB 67000.00 104089.3 0.0 -1000
+X89 DB 67000.00 105285.7 0.0 -1000
+X90 DB 67000.00 106482.1 0.0 -1000
+X91 DB 67000.00 107678.6 0.0 -1000
+X92 DB 67000.00 108875.0 0.0 -1000
+X93 DB 67000.00 110071.4 0.0 -1000
+X94 DB 67000.00 111267.9 0.0 -1000
+X95 DB 67000.00 112464.3 0.0 -1000
+X96 DB 67000.00 113660.7 0.0 -1000
+X97 DB 67000.00 114857.1 0.0 -1000
+X98 DB 67000.00 116053.6 0.0 -1000
+X99 DB 67000.00 117250.0 0.0 -1000
+X100 DB 67000.00 118446.4 0.0 -1000
+X101 DB 67000.00 119642.9 0.0 -1000
+X102 DB 67000.00 120839.3 0.0 -1000
+X103 DB 67000.00 122035.7 0.0 -1000
+X104 DB 67000.00 123232.1 0.0 -1000
+X105 DB 67000.00 124428.6 0.0 -1000
+X106 DB 67000.00 125625.0 0.0 -1000
+X107 DB 67000.00 126821.4 0.0 -1000
+X108 DB 67000.00 128017.9 0.0 -1000
+X109 DB 67000.00 129214.3 0.0 -1000
+X110 DB 67000.00 130410.7 0.0 -1000
+X111 DB 67000.00 131607.1 0.0 -1000
+X112 DB 67000.00 132803.6 0.0 -1000
+Y1 DB 0.000000 67000.00 0.0 -1000
+Y2 DB 1196.429 67000.00 0.0 -1000
+Y3 DB 2392.857 67000.00 0.0 -1000
+Y4 DB 3589.286 67000.00 0.0 -1000
+Y5 DB 4785.714 67000.00 0.0 -1000
+Y6 DB 5982.143 67000.00 0.0 -1000
+Y7 DB 7178.571 67000.00 0.0 -1000
+Y8 DB 8375.000 67000.00 0.0 -1000
+Y9 DB 9571.429 67000.00 0.0 -1000
+Y10 DB 10767.86 67000.00 0.0 -1000
+Y11 DB 11964.29 67000.00 0.0 -1000
+Y12 DB 13160.71 67000.00 0.0 -1000
+Y13 DB 14357.14 67000.00 0.0 -1000
+Y14 DB 15553.57 67000.00 0.0 -1000
+Y15 DB 16750.00 67000.00 0.0 -1000
+Y16 DB 17946.43 67000.00 0.0 -1000
+Y17 DB 19142.86 67000.00 0.0 -1000
+Y18 DB 20339.29 67000.00 0.0 -1000
+Y19 DB 21535.71 67000.00 0.0 -1000
+Y20 DB 22732.14 67000.00 0.0 -1000
+Y21 DB 23928.57 67000.00 0.0 -1000
+Y22 DB 25125.00 67000.00 0.0 -1000
+Y23 DB 26321.43 67000.00 0.0 -1000
+Y24 DB 27517.86 67000.00 0.0 -1000
+Y25 DB 28714.29 67000.00 0.0 -1000
+Y26 DB 29910.71 67000.00 0.0 -1000
+Y27 DB 31107.14 67000.00 0.0 -1000
+Y28 DB 32303.57 67000.00 0.0 -1000
+Y29 DB 33500.00 67000.00 0.0 -1000
+Y30 DB 34696.43 67000.00 0.0 -1000
+Y31 DB 35892.86 67000.00 0.0 -1000
+Y32 DB 37089.29 67000.00 0.0 -1000
+Y33 DB 38285.71 67000.00 0.0 -1000
+Y34 DB 39482.14 67000.00 0.0 -1000
+Y35 DB 40678.57 67000.00 0.0 -1000
+Y36 DB 41875.00 67000.00 0.0 -1000
+Y37 DB 43071.43 67000.00 0.0 -1000
+Y38 DB 44267.86 67000.00 0.0 -1000
+Y39 DB 45464.29 67000.00 0.0 -1000
+Y40 DB 46660.71 67000.00 0.0 -1000
+Y41 DB 47857.14 67000.00 0.0 -1000
+Y42 DB 49053.57 67000.00 0.0 -1000
+Y43 DB 50250.00 67000.00 0.0 -1000
+Y44 DB 51446.43 67000.00 0.0 -1000
+Y45 DB 52642.86 67000.00 0.0 -1000
+Y46 DB 53839.29 67000.00 0.0 -1000
+Y47 DB 55035.71 67000.00 0.0 -1000
+Y48 DB 56232.14 67000.00 0.0 -1000
+Y49 DB 57428.57 67000.00 0.0 -1000
+Y50 DB 58625.00 67000.00 0.0 -1000
+Y51 DB 59821.43 67000.00 0.0 -1000
+Y52 DB 61017.86 67000.00 0.0 -1000
+Y53 DB 62214.29 67000.00 0.0 -1000
+Y54 DB 63410.71 67000.00 0.0 -1000
+Y55 DB 64607.14 67000.00 0.0 -1000
+Y56 DB 65803.57 67000.00 0.0 -1000
+Y57 DB 67000.00 67000.00 0.0 -1000
+Y58 DB 68196.43 67000.00 0.0 -1000
+Y59 DB 69392.86 67000.00 0.0 -1000
+Y60 DB 70589.29 67000.00 0.0 -1000
+Y61 DB 71785.71 67000.00 0.0 -1000
+Y62 DB 72982.14 67000.00 0.0 -1000
+Y63 DB 74178.57 67000.00 0.0 -1000
+Y64 DB 75375.00 67000.00 0.0 -1000
+Y65 DB 76571.43 67000.00 0.0 -1000
+Y66 DB 77767.86 67000.00 0.0 -1000
+Y67 DB 78964.29 67000.00 0.0 -1000
+Y68 DB 80160.71 67000.00 0.0 -1000
+Y69 DB 81357.14 67000.00 0.0 -1000
+Y70 DB 82553.57 67000.00 0.0 -1000
+Y71 DB 83750.00 67000.00 0.0 -1000
+Y72 DB 84946.43 67000.00 0.0 -1000
+Y73 DB 86142.86 67000.00 0.0 -1000
+Y74 DB 87339.29 67000.00 0.0 -1000
+Y75 DB 88535.71 67000.00 0.0 -1000
+Y76 DB 89732.14 67000.00 0.0 -1000
+Y77 DB 90928.57 67000.00 0.0 -1000
+Y78 DB 92125.00 67000.00 0.0 -1000
+Y79 DB 93321.43 67000.00 0.0 -1000
+Y80 DB 94517.86 67000.00 0.0 -1000
+Y81 DB 95714.29 67000.00 0.0 -1000
+Y82 DB 96910.71 67000.00 0.0 -1000
+Y83 DB 98107.14 67000.00 0.0 -1000
+Y84 DB 99303.57 67000.00 0.0 -1000
+Y85 DB 100500.0 67000.00 0.0 -1000
+Y86 DB 101696.4 67000.00 0.0 -1000
+Y87 DB 102892.9 67000.00 0.0 -1000
+Y88 DB 104089.3 67000.00 0.0 -1000
+Y89 DB 105285.7 67000.00 0.0 -1000
+Y90 DB 106482.1 67000.00 0.0 -1000
+Y91 DB 107678.6 67000.00 0.0 -1000
+Y92 DB 108875.0 67000.00 0.0 -1000
+Y93 DB 110071.4 67000.00 0.0 -1000
+Y94 DB 111267.9 67000.00 0.0 -1000
+Y95 DB 112464.3 67000.00 0.0 -1000
+Y96 DB 113660.7 67000.00 0.0 -1000
+Y97 DB 114857.1 67000.00 0.0 -1000
+Y98 DB 116053.6 67000.00 0.0 -1000
+Y99 DB 117250.0 67000.00 0.0 -1000
+Y100 DB 118446.4 67000.00 0.0 -1000
+Y101 DB 119642.9 67000.00 0.0 -1000
+Y102 DB 120839.3 67000.00 0.0 -1000
+Y103 DB 122035.7 67000.00 0.0 -1000
+Y104 DB 123232.1 67000.00 0.0 -1000
+Y105 DB 124428.6 67000.00 0.0 -1000
+Y106 DB 125625.0 67000.00 0.0 -1000
+Y107 DB 126821.4 67000.00 0.0 -1000
+Y108 DB 128017.9 67000.00 0.0 -1000
+Y109 DB 129214.3 67000.00 0.0 -1000
+Y110 DB 130410.7 67000.00 0.0 -1000
+Y111 DB 131607.1 67000.00 0.0 -1000
+Y112 DB 132803.6 67000.00 0.0 -1000
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/boundary_definition.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/boundary_definition.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,378 @@
+#############################################################################
+# boundary_definition.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+def define_absorbing_surf():
+ """
+ define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
+ it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
+ ...
+ absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+
+
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+# for k in list_surf:
+# center_point = cubit.get_center_point("surface", k)
+# if abs((center_point[0] - xmin_box)/xmin_box) <= 0.005:
+# absorbing_surf_xmin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[0] - xmax_box)/xmax_box) <= 0.005:
+# absorbing_surf_xmax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymin_box)/ymin_box) <= 0.005:
+# absorbing_surf_ymin.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[1] - ymax_box)/ymax_box) <= 0.005:
+# absorbing_surf_ymax.append(k)
+# absorbing_surf.append(k)
+# elif abs((center_point[2] - zmin_box)/zmin_box) <= 0.005:
+# absorbing_surf_bottom.append(k)
+# absorbing_surf.append(k)
+# else:
+# sbox=cubit.get_bounding_box('surface',k)
+# dz=abs((sbox[7] - zmax_box)/zmax_box)
+# normal=cubit.get_surface_normal(k)
+# zn=normal[2]
+# dn=abs(zn-1)
+# if dz <= 0.001 and dn < 0.2:
+# top_surf.append(k)
+
+ #box lengths
+ x_len = abs( xmax_box - xmin_box)
+ y_len = abs( ymax_box - ymin_box)
+ z_len = abs( zmax_box - zmin_box)
+
+ print '##boundary box: '
+ print '## x length: ' + str(x_len)
+ print '## y length: ' + str(y_len)
+ print '## z length: ' + str(z_len)
+
+ # tolerance parameters
+ absorbing_surface_distance_tolerance=0.005
+ topographic_surface_distance_tolerance=0.001
+ topographic_surface_normal_tolerance=0.2
+
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ if abs((center_point[0] - xmin_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[0] - xmax_box)/x_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_xmax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymin_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymin.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[1] - ymax_box)/y_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_ymax.append(k)
+ absorbing_surf.append(k)
+ elif abs((center_point[2] - zmin_box)/z_len) <= absorbing_surface_distance_tolerance:
+ absorbing_surf_bottom.append(k)
+ absorbing_surf.append(k)
+ else:
+ sbox=cubit.get_bounding_box('surface',k)
+ dz=abs((sbox[7] - zmax_box)/z_len)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ dn=abs(zn-1)
+ if dz <= topographic_surface_distance_tolerance and dn < topographic_surface_normal_tolerance:
+ top_surf.append(k)
+
+ return absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,top_surf
+
+def define_absorbing_surf_nopar():
+ """
+ define the absorbing surfaces for a layered topological box where boundary surfaces are not parallel to the axis.
+ it returns absorbing_surf,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ """
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ from sets import Set
+ def product(*args, **kwds):
+ # product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
+ # product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
+ pools = map(tuple, args) * kwds.get('repeat', 1)
+ result = [[]]
+ for pool in pools:
+ result = [x+[y] for x in result for y in pool]
+ return result
+ absorbing_surf=[]
+ absorbing_surf_xmin=[]
+ absorbing_surf_xmax=[]
+ absorbing_surf_ymin=[]
+ absorbing_surf_ymax=[]
+ absorbing_surf_bottom=[]
+ top_surf=[]
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ zmax_box=cubit.get_total_bounding_box("volume",list_vol)[7]
+ zmin_box=cubit.get_total_bounding_box("volume",list_vol)[6] #it is the z_min of the box ... box= xmin,xmax,d,ymin,ymax,d,zmin...
+ xmin_box=cubit.get_total_bounding_box("volume",list_vol)[0]
+ xmax_box=cubit.get_total_bounding_box("volume",list_vol)[1]
+ ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
+ ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ lv=[]
+ for k in list_surf:
+ sbox=cubit.get_bounding_box('surface',k)
+ dzmax=abs((sbox[7] - zmax_box)/zmax_box)
+ dzmin=abs((sbox[6] - zmin_box)/zmin_box)
+ normal=cubit.get_surface_normal(k)
+ zn=normal[2]
+ if dzmax <= 0.001 and zn > 0.7:
+ top_surf.append(k)
+ list_vertex=cubit.get_relatives('surface',k,'vertex')
+ for v in list_vertex:
+ valence=cubit.get_valence(v)
+ if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
+ lv.append(v)
+ elif dzmin <= 0.001 and zn < -0.7:
+ absorbing_surf.append(k)
+ lp=[]
+ combs=product(lv,lv)
+ for comb in combs:
+ v1=comb[0]
+ v2=comb[1]
+ c=Set(cubit.get_relatives("vertex",v1,"curve")) & Set(cubit.get_relatives("vertex",v2,"curve"))
+ if len(c) == 1:
+ p=cubit.get_center_point("curve",list(c)[0])
+ lp.append(p)
+ for k in list_surf:
+ center_point = cubit.get_center_point("surface", k)
+ for p in lp:
+ if abs((center_point[0] - p[0])/p[0]) <= 0.005 and abs((center_point[1] - p[1])/p[1]) <= 0.005:
+ absorbing_surf.append(k)
+ break
+ return absorbing_surf,top_surf
+
+def define_absorbing_surf_sphere():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ surf=[]
+ list_surf=cubit.parse_cubit_list("surface","all")
+ for s in list_surf:
+ v=cubit.get_relatives('surface',s,'volume')
+ if len(v) == 1:
+ surf.append(s)
+ return surf
+
+def define_block():
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ list_vol=cubit.parse_cubit_list("volume","all")
+ init_n_vol=len(list_vol)
+ list_name=map(lambda x: 'vol'+x,map(str,list_vol))
+ return list_vol,list_name
+
+def build_block(vol_list,name):
+ from sets import Set
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ block_list=cubit.get_block_id_list()
+ if len(block_list) > 0:
+ id_block=max(block_list)
+ else:
+ id_block=0
+ for v,n in zip(vol_list,name):
+ id_block+=1
+ v_other=Set(vol_list)-Set([v])
+ #command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ command= 'block '+str(id_block)+' hex in vol '+str(v)
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
+
+def build_block_side(surf_list,name,obj='surface'):
+ try:
+ cubit.cmd('comment')
+ except:
+ try:
+ import cubit
+ cubit.init([""])
+ except:
+ print 'error importing cubit'
+ import sys
+ sys.exit()
+ id_nodeset=cubit.get_next_nodeset_id()
+ id_block=cubit.get_next_block_id()
+
+
+ if obj == 'hex':
+ txt='hex in node in surface'
+ txt1='block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt2="block "+str(id_block)+" name '"+name+"'"
+ txt1=txt1.replace("["," ").replace("]"," ")
+ elif obj == 'node':
+ txt=obj+' in surface'
+ txt1= 'nodeset '+str(id_nodeset)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "nodeset "+str(id_nodeset)+" name '"+name+"'"
+ elif obj == 'face' or obj == 'edge':
+ txt=obj+' in surface'
+ txt1= 'block '+str(id_block)+ ' '+ txt +' '+str(list(surf_list))
+ txt1 = txt1.replace("["," ").replace("]"," ")
+ txt2 = "block "+str(id_block)+" name '"+name+"'"
+ else:
+ txt1=''
+ # do not execute: block id might be wrong
+ print "##block "+str(id_block)+" name '"+name+"_notsupported (only hex,face,edge,node)'"
+ txt2=''
+
+
+ cubit.cmd(txt1)
+ cubit.cmd(txt2)
+
+def define_bc(*args,**keys):
+ parallel=keys.get('parallel',True)
+ closed=keys.get('closed',False)
+ if not closed:
+ print "##open region"
+
+ if parallel:
+ surf,xmin,xmax,ymin,ymax,bottom,topo=define_absorbing_surf()
+ else:
+ surf,topo=define_absorbing_surf_nopar()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ entities=args[0]
+ print entities
+ for entity in entities:
+ print "##entity: "+str(entity)
+
+ build_block_side(topo,entity+'_topo',obj=entity)
+ build_block_side(surf,entity+'_abs',obj=entity)
+ if parallel:
+ build_block_side(xmin,entity+'_abs_xmin',obj=entity)
+ build_block_side(xmax,entity+'_abs_xmax',obj=entity)
+ build_block_side(ymin,entity+'_abs_ymin',obj=entity)
+ build_block_side(ymax,entity+'_abs_ymax',obj=entity)
+ build_block_side(bottom,entity+'_abs_bottom',obj=entity)
+ else:
+ print "##closed region"
+
+ surf=define_absorbing_surf_sphere()
+ v_list,name_list=define_block()
+ build_block(v_list,name_list)
+ entities=args[0]
+ for entity in entities:
+ build_block_side(surf,entity+'_closedvol',obj=entity)
+
+
+
+#entities=['surface','face','edge','node','hex']
+#define_bc(entities,parallel=True)
+#define_bc(entities,parallel=False)
+#define_bc(entities,parallel=False,closed=True)
+
+
+# call
+#entities=['surface','face']
+#define_bc(entities,parallel=True)
+
+#block 1 attribute count 5
+#block 2 attribute count 0
+#block 2 name '/prova/interface1'
+#block 2 attribute count 3
+#block 3 attribute count 5
+#block 1 attribute index 2 1500
+#block 1 attribute count 2
+#block 3 attribute index 1 2
+#block 3 attribute index 2 5800
+#block 3 attribute index 3 3900
+#block 3 attribute index 4 1500
+#block 3 attribute index 5 3.5
+#block 1 name 'top'
+#block 3 name 'bottom'
+#block 2 attribute index 1 -1
+#block 2 attribute index 2 2
+#block 2 attribute count 4
+#block 2 attribute index 2 1
+#block 2 attribute index 3 2
+
+
+
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/cubit2specfem3d.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,698 @@
+#!python
+#############################################################################
+# cubit2specfem3d.py #
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+#for a complete definition of the format of the mesh in SPECFEM3D_SESAME check the manual (http:/):
+#
+#USAGE
+#
+#############################################################################
+#PREREQUISITE
+#The mesh must be prepared
+# automatically using the module boundary_definition (see boundary_definition.py for more information)
+#or
+# manually following the convention:
+# - each material should have a block defined by:
+# material domain_flag (acoustic/elastic/poroelastic)name,flag of the material (integer),p velocity
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
+# - each mesh should have the block definition for the face on the free_surface (topography),
+# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
+# - each mesh should have the block definition for the faces on the absorbing boundaries,
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+#
+#############################################################################
+#RUN
+#In a python script or in the cubit python tab call:
+#
+# export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME)
+#
+#the modele create a python class for the mesh: ex. profile=mesh()
+#and it export the files of the mesh needed by the partitioner of SESAME
+#
+#############################################################################
+#OUTPUT
+#The default output are 11 ASCII files:
+#__________________________________________________________________________________________
+#mesh_name='mesh_file' -> the file that contains the connectity of the all mesh
+# format:
+# number of elements
+# id_elements id_node1 id_node2 id_node3 id_node4 id_node5 id_node6 id_node7 id_node8
+# .....
+#
+#__________________________________________________________________________________________
+##nodecoord_name='nodes_coords_file' -> the file that contains the coordinates of the nodes of the all mesh
+# format:
+# number of nodes
+# id_node x_coordinate y_coordinate z_coordinate
+# .....
+#
+#__________________________________________________________________________________________
+##material_name='materials_file' -> the file that contains the material flag of the elements
+# format:
+# id_element flag
+# .....
+#
+#__________________________________________________________________________________________
+##nummaterial_name='nummaterial_velocity_file' -> table of the material properties
+# format:
+# flag rho vp vs 0 0 #full definition of the properties, flag > 0
+# .....
+# flag 'tomography' file_name #for interpolation with tomography
+# .....
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
+#__________________________________________________________________________________________
+##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+# ....
+#
+#__________________________________________________________________________________________
+##freename='free_surface_file' -> file with the hex on the free surface (usually the topography)
+# format:
+# number of faces
+# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
+#
+#__________________________________________________________________________________________
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
+#
+#############################################################################
+
+
+import cubit
+
+
+class mtools(object):
+ """docstring for ciao"""
+ def __init__(self,frequency,list_surf,list_vp):
+ super(mtools, self).__init__()
+ self.frequency = frequency
+ self.list_surf = list_surf
+ self.list_vp = list_vp
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ def __repr__(self):
+ txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ return txt
+ def freq2meshsize(self,vp):
+ velocity=vp*.5
+ self.size=(1/2.5)*velocity/self.frequency*(self.ngll-1)/self.point_wavelength
+ self.dt=.4*self.size/vp*self.percent_gll
+ return self.size,self.dt
+ def mesh_it(self):
+ for surf,vp in zip(self.list_surf,self.list_vp):
+ command = "surface "+str(surf)+" size "+str(self.freq2meshsize(vp)[0])
+ cubit.cmd(command)
+ command = "surface "+str(surf)+ 'scheme pave'
+ cubit.cmd(command)
+ command = "mesh surf "+str(surf)
+ cubit.cmd(command)
+
+class block_tools:
+ def __int__(self):
+ pass
+ def create_blocks(self,mesh_entity,list_entity=None,):
+ if mesh_entity =='surface':
+ txt=' face in surface '
+ elif mesh_entity == 'curve':
+ txt=' edge in curve '
+ elif mesh_entity == 'group':
+ txt=' face in group '
+ if list_entity:
+ if not isinstance(list_entity,list):
+ list_entity=[list_entity]
+ for entity in list_entity:
+ iblock=cubit.get_next_block_id()
+ command = "block "+str(iblock)+ txt +str(entity)
+ cubit.cmd(command)
+ def material_file(self,filename):
+ matfile=open(filename,'w')
+ material=[]
+ for record in matfile:
+ mat_name,vp_str=record.split()
+ vp=float(vp_str)
+ material.append([mat_name,vp])
+ self.material=dict(material)
+ def assign_block_material(self,id_block,mat_name,vp=None):
+ try:
+ material=self.material
+ except:
+ material=None
+ cubit.cmd('block '+str(id_block)+' attribute count 2')
+ cubit.cmd('block '+str(id_block)+' attribute index 1 '+str(id_block))
+ if material:
+ if material.has_key(mat_name):
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(material[mat_name]))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(material[mat_name])+' from database'
+ elif vp:
+ cubit.cmd('block '+str(id_block)+' attribute index 2 '+str(vp))
+ print 'block '+str(id_block)+' - material '+mat_name+' - vp '+str(vp)
+ else:
+ print 'assignment impossible: check if '+mat_name+' is in the database or specify vp'
+
+class mesh_tools(block_tools):
+ """Tools for the mesh
+ #########
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ #########
+ length=edge_length(edge)
+ return the length of a edge
+ #########
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ #########
+ """
+ def __int__(self):
+ pass
+ def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
+ """
+ dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
+ """
+ ratiostore=1e10
+ dtstore=1e10
+ edgedtstore=-1
+ edgeratiostore=-1
+ for edge in edges:
+ d=self.edge_length(edge)
+ ratio=(1/2.5)*velocity/d*(self.ngll-1)/self.point_wavelength
+ dt=.4*d/velocity*self.percent_gll
+ if dt<dtstore:
+ dtstore=dt
+ edgedtstore=edge
+ if ratio < ratiostore:
+ ratiostore=ratio
+ edgeratiostore=edge
+ try:
+ for bin_d,bin_u,side in zip(bins_d,bins_u,sidelist):
+ if ratio >= bin_d and ratio < bin_u:
+ command = "sideset "+str(side)+" edge "+str(edge)
+ cubit.cmd(command)
+ #print command
+ break
+ except:
+ pass
+ return dtstore,edgedtstore,ratiostore,edgeratiostore
+ def edge_length(self,edge):
+ """
+ length=edge_length(edge)
+ return the length of a edge
+ """
+ from math import sqrt
+ nodes=cubit.get_connectivity('Edge',edge)
+ x0,y0,z0=cubit.get_nodal_coordinates(nodes[0])
+ x1,y1,z1=cubit.get_nodal_coordinates(nodes[1])
+ d=sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
+ return d
+ def edge_min_length(self,surface):
+ """
+ edge_min,length=edge_min_length(surface)
+ given the cubit id of a surface, it return the edge with minimun length
+ """
+ from math import sqrt
+ self.dmin=99999
+ edge_store=0
+ command = "group 'list_edge' add edge in surf "+str(surface)
+ command = command.replace("["," ").replace("]"," ")
+ #print command
+ cubit.cmd(command)
+ group=cubit.get_id_from_name("list_edge")
+ edges=cubit.get_group_edges(group)
+ command = "delete group "+ str(group)
+ cubit.cmd(command)
+ for edge in edges:
+ d=self.edge_length(edge)
+ if d<dmin:
+ self.dmin=d
+ edge_store=edge
+ self.edgemin=edge_store
+ return self.edgemin,self.dmin
+ def normal_check(self,nodes,normal):
+ tres=.2
+ p0=cubit.get_nodal_coordinates(nodes[0])
+ p1=cubit.get_nodal_coordinates(nodes[1])
+ p2=cubit.get_nodal_coordinates(nodes[2])
+ a=[p1[0]-p0[0],p1[1]-p0[1],p1[2]-p0[2]]
+ b=[p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]
+ axb=[a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1] - a[1]*b[0]]
+ dot=0.0
+ for i in (0,1,2):
+ dot=dot+axb[i]*normal[i]
+ if dot > 0:
+ return nodes
+ elif dot < 0:
+ return nodes[0],nodes[3],nodes[2],nodes[1]
+ else:
+ print 'error, dot=0', axb,normal,dot,p0,p1,p2
+ def mesh_analysis(self,frequency):
+ from sets import Set
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ bins_d=[0.0001]+range(0,int(frequency)+1)+[1000]
+ bins_u=bins_d[1:]
+ dt=[]
+ ed_dt=[]
+ r=[]
+ ed_r=[]
+ nstart=cubit.get_next_sideset_id()
+ command = "del sideset all"
+ cubit.cmd(command)
+ for bin_d,bin_u in zip(bins_d,bins_u):
+ nsideset=cubit.get_next_sideset_id()
+ command='create sideset '+str(nsideset)
+ cubit.cmd(command)
+ command = "sideset "+str(nsideset)+ " name "+ "'ratio-["+str(bin_d)+"_"+str(bin_u)+"['"
+ cubit.cmd(command)
+ nend=cubit.get_next_sideset_id()
+ sidelist=range(nstart,nend)
+ for block in self.block_mat:
+ name=cubit.get_exodus_entity_name('block',block)
+ velocity=self.material[name][1]
+ if velocity > 0:
+ faces=cubit.get_block_faces(block)
+ edges=[]
+ for face in faces:
+ es=cubit.get_sub_elements("face", face, 1)
+ edges=edges+list(es)
+ edges=Set(edges)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
+ dt.append(dtstore)
+ ed_dt.append(edgedtstore)
+ r.append(ratiostore)
+ ed_r.append(edgeratiostore)
+ self.ddt=zip(ed_dt,dt)
+ self.dr=zip(ed_r,r)
+ def sorter(x, y):
+ return cmp(x[1],y[1])
+ self.ddt.sort(sorter)
+ self.dr.sort(sorter)
+ print self.ddt,self.dr
+ print 'Deltat minimum => edge:'+str(self.ddt[0][0])+' dt: '+str(self.ddt[0][1])
+ print 'minimum frequency resolved => edge:'+str(self.dr[0][0])+' frequency: '+str(self.dr[0][1])
+ return self.ddt[0],self.dr[0]
+
+class mesh(object,mesh_tools):
+ def __init__(self):
+ super(mesh, self).__init__()
+ self.mesh_name='mesh_file'
+ self.nodecoord_name='nodes_coords_file'
+ self.material_name='materials_file'
+ self.nummaterial_name='nummaterial_velocity_file'
+ self.absname='absorbing_surface_file'
+ self.freename='free_surface_file'
+ self.recname='STATIONS'
+ self.face='QUAD4'
+ self.hex='HEX8'
+ self.edge='BAR2'
+ self.topo='face_topo'
+ self.rec='receivers'
+ self.block_definition()
+ self.ngll=5
+ self.percent_gll=0.172
+ self.point_wavelength=5
+ cubit.cmd('compress')
+ def __repr__(self):
+ pass
+ def block_definition(self):
+ block_flag=[]
+ block_mat=[]
+ block_bc=[]
+ block_bc_flag=[]
+ material={}
+ bc={}
+ blocks=cubit.get_block_id_list()
+ for block in blocks:
+ name=cubit.get_exodus_entity_name('block',block)
+ ty=cubit.get_block_element_type(block)
+ print block,name,blocks,ty,self.hex,self.face
+ if ty == self.hex:
+ flag=None
+ vel=None
+ vs=None
+ rho=None
+ q=0
+ ani=0
+ # material domain id
+ if name == "acoustic" :
+ imaterial = 1
+ elif name == "elastic" :
+ imaterial = 2
+ elif name == "poroelastic" :
+ imaterial = 3
+ else :
+ imaterial = 0
+
+ nattrib=cubit.get_block_attribute_count(block)
+ if nattrib != 0:
+ # material flag:
+ # positive => material properties,
+ # negative => interface/tomography domain
+ flag=int(cubit.get_block_attribute_value(block,0))
+ if flag > 0 and nattrib >= 2:
+ # vp
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ # vs
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
+ elif flag < 0:
+ # velocity model
+ vel=name
+ attrib=cubit.get_block_attribute_value(block,1)
+ if attrib == 1:
+ kind='interface'
+ flag_down=cubit.get_block_attribute_value(block,2)
+ flag_up=cubit.get_block_attribute_value(block,3)
+ elif attrib == 2:
+ kind='tomography'
+ else:
+ flag=block
+ vel,vs,rho,q,ani=(name,0,0,0,0)
+ block_flag.append(int(flag))
+ block_mat.append(block)
+ if flag > 0:
+ par=tuple([imaterial,flag,vel,vs,rho,q,ani])
+ elif flag < 0:
+ if kind=='interface':
+ par=tuple([imaterial,flag,kind,name,flag_down,flag_up])
+ elif kind=='tomography':
+ par=tuple([imaterial,flag,kind,name])
+ elif flag==0:
+ par=tuple([imaterial,flag,name])
+ material[block]=par
+ elif ty == self.face: #Stacey condition, we need hex here for pml
+ block_bc_flag.append(4)
+ block_bc.append(block)
+ bc[block]=4 #face has connectivity = 4
+ if name == self.topo: topography_face=block
+ else:
+ print 'blocks no properly defined',ty
+ return None, None,None,None,None,None,None,None
+ nsets=cubit.get_nodeset_id_list()
+ if len(nsets) == 0: self.receivers=None
+ for nset in nsets:
+ name=cubit.get_exodus_entity_name('nodeset',nset)
+ if name == self.rec:
+ self.receivers=nset
+ else:
+ print 'nodeset '+name+' not defined'
+ self.receivers=None
+ try:
+ self.block_mat=block_mat
+ self.block_flag=block_flag
+ self.block_bc=block_bc
+ self.block_bc_flag=block_bc_flag
+ self.material=material
+ self.bc=bc
+ self.topography=topography_face
+ except:
+ print 'blocks no properly defined'
+ print block_mat
+ print block_flag
+ print block_bc
+ print block_bc_flag
+ print material
+ print bc
+ print topography
+ def mat_parameter(self,properties):
+ #TODO: material property acoustic/elastic/poroelastic ? .... where?
+ print "#material properties:"
+ print properties
+ imaterial=properties[0]
+ flag=properties[1]
+ if flag > 0:
+ vel=properties[2]
+ if properties[3] is None and type(vel) != str:
+ # velocity model scales with given vp value
+ if vel >= 30:
+ m2km=1000.
+ else:
+ m2km=1.
+ vp=vel/m2km
+ rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
+ txt='%1i %3i %20f %20f %20f %1i %1i\n' % (properties[0],properties[1],rho,vel,vel/(3**.5),0,0)
+ elif type(vel) != str:
+ # velocity model given as vp,vs,rho,..
+ #format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%1i %3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[1],properties[4], \
+ properties[2],properties[3],properties[5],properties[6])
+ else:
+ txt='%1i %3i %s \n' % (properties[0],properties[1],properties[2])
+ elif flag < 0:
+ if properties[2] == 'tomography':
+ txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],properties[3])
+ elif properties[2] == 'interface':
+ txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],\
+ properties[4],properties[5])
+ return txt
+ def nummaterial_write(self,nummaterial_name):
+ print 'Writing '+nummaterial_name+'.....'
+ nummaterial=open(nummaterial_name,'w')
+ for block in self.block_mat:
+ #name=cubit.get_exodus_entity_name('block',block)
+ nummaterial.write(self.mat_parameter(self.material[block]))
+ nummaterial.close()
+ print 'Ok'
+ def mesh_write(self,mesh_name):
+ meshfile=open(mesh_name,'w')
+ print 'Writing '+mesh_name+'.....'
+ num_elems=cubit.get_hex_count()
+ print num_elems
+ meshfile.write(str(num_elems)+'\n')
+ num_write=0
+ for block,flag in zip(self.block_mat,self.block_flag):
+ #print block,flag
+ hexes=cubit.get_block_hexes(block)
+ #print len(hexes)
+ for hexa in hexes:
+ #print hexa
+ nodes=cubit.get_connectivity('Hex',hexa)
+ #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
+ txt=('%10i ')% hexa
+ txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ meshfile.write(txt)
+ meshfile.close()
+ print 'Ok',str(num_elems)
+ def material_write(self,mat_name):
+ mat=open(mat_name,'w')
+ print 'Writing '+mat_name+'.....'
+ for block,flag in zip(self.block_mat,self.block_flag):
+ hexes=cubit.get_block_hexes(block)
+ for hexa in hexes:
+ mat.write(('%10i %10i\n') % (hexa,flag))
+ mat.close()
+ print 'Ok'
+ def nodescoord_write(self,nodecoord_name):
+ nodecoord=open(nodecoord_name,'w')
+ print 'Writing '+nodecoord_name+'.....'
+ node_list=cubit.parse_cubit_list('node','all')
+ num_nodes=len(node_list)
+ nodecoord.write('%10i\n' % num_nodes)
+ #
+ for node in node_list:
+ x,y,z=cubit.get_nodal_coordinates(node)
+ txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
+ nodecoord.write(txt)
+ nodecoord.close()
+ print 'Ok'
+ def free_write(self,freename=None):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ normal=(0,0,1)
+ if not freename: freename=self.freename
+ freehex=open(freename,'w')
+ print 'Writing '+freename+'.....'
+ #
+ #
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block == self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print name,block
+ quads_all=cubit.get_block_faces(block)
+ print 'face = ',len(quads_all)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ freehex.write('%10i\n' % len(quads_all))
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ #print f
+ nodes=cubit.get_connectivity('Face',f)
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ freehex.write(txt)
+ freehex.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def abs_write(self,absname=None):
+ import re
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ from sets import Set
+ if not absname: absname=self.absname
+ #
+ #
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print name,block
+ absflag=False
+ if re.search('xmin',name):
+ print 'xmin'
+ abshex_local=open(absname+'_xmin','w')
+ normal=(-1,0,0)
+ elif re.search('xmax',name):
+ print "xmax"
+ abshex_local=open(absname+'_xmax','w')
+ normal=(1,0,0)
+ elif re.search('ymin',name):
+ print "ymin"
+ abshex_local=open(absname+'_ymin','w')
+ normal=(0,-1,0)
+ elif re.search('ymax',name):
+ print "ymax"
+ abshex_local=open(absname+'_ymax','w')
+ normal=(0,1,0)
+ elif re.search('bottom',name):
+ print "bottom"
+ abshex_local=open(absname+'_bottom','w')
+ normal=(0,0,-1)
+ elif re.search('abs',name):
+ print "abs all - no implemented yet"
+ absflag=True
+ abshex_local=open(absname,'w')
+ #
+ #
+ quads_all=cubit.get_block_faces(block)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ abshex_local.write('%10i\n' % len(quads_all))
+ #command = "group 'list_hex' add hex in face "+str(quads_all)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ #group=cubit.get_id_from_name("list_hex")
+ #list_hex=cubit.get_group_hexes(group)
+ #command = "delete group "+ str(group)
+ #cubit.cmd(command)
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ nodes=cubit.get_connectivity('Face',f)
+ if not absflag:
+ nodes_ok=self.normal_check(nodes,normal)
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ else:
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
+ abshex_local.write(txt)
+ abshex_local.close()
+ print 'Ok'
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ def rec_write(self,recname):
+ print 'Writing '+self.recname+'.....'
+ recfile=open(self.recname,'w')
+ nodes=cubit.get_nodeset_nodes(self.receivers)
+ for i,n in enumerate(nodes):
+ x,y,z=cubit.get_nodal_coordinates(n)
+ recfile.write('ST%i XX %20f %20f 0.0 0.0 \n' % (i,x,z))
+ recfile.close()
+ print 'Ok'
+ def write(self,path=''):
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ if len(path) != 0:
+ if path[-1] != '/': path=path+'/'
+ self.mesh_write(path+self.mesh_name)
+ self.material_write(path+self.material_name)
+ self.nodescoord_write(path+self.nodecoord_name)
+ self.free_write(path+self.freename)
+ self.abs_write(path+self.absname)
+ self.nummaterial_write(path+self.nummaterial_name)
+ if self.receivers: self.rec_write(path+self.recname)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+
+def export2SESAME(path_exporting_mesh_SPECFEM3D_SESAME):
+ sem_mesh=mesh()
+ sem_mesh.write(path=path_exporting_mesh_SPECFEM3D_SESAME)
+
+
+
+if __name__ == '__main__':
+ path='/Users/emanuele/Desktop/'
+ export2SESAME(path)
+
+#TODO: change the algorithm for the abs detection of the hex
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8-nodoubling.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8-nodoubling.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8-nodoubling.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,118 @@
+#!/usr/bin/env python
+
+###########################################################################
+#### TNM: This is the mesh generation, adapted from a journal file
+#### specific to the settings of Komatitsch and Tromp 1999, Fig.8
+#### Aug 2009
+###########################################################################
+#
+# modified version: top layer is water
+#
+###########################################################################
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+cubit.cmd('reset')
+cubit.cmd('brick x 134000 y 134000 z 60000')
+
+# This seems to conflict with boundary_definition.py
+# ....which leaves the model space at e.g. x=[-67,67] km
+cubit.cmd('volume 1 move x 67000 y 67000 z -30000')
+
+# create vertices for discontinuity
+cubit.cmd('split curve 9 distance 3000')
+cubit.cmd('split curve 10 distance 3000')
+cubit.cmd('split curve 11 distance 3000')
+cubit.cmd('split curve 12 distance 3000')
+
+# create surface for interface
+cubit.cmd('create surface vertex 9 10 12 11')
+
+cubit.cmd('section volume 1 with surface 7 keep normal')
+cubit.cmd('section volume 1 with surface 7 reverse')
+
+# create vertices for auxiliary interface to allow for refinement
+cubit.cmd('split curve 29 distance 9000')
+cubit.cmd('split curve 31 distance 9000')
+cubit.cmd('split curve 32 distance 9000')
+cubit.cmd('split curve 36 distance 9000')
+
+# create surface for buffer interface to refine BELOW the discontinuity
+cubit.cmd('create surface vertex 25 26 28 27')
+
+cubit.cmd('section volume 3 with surface 19 keep normal')
+cubit.cmd('section volume 3 with surface 19 reverse')
+
+cubit.cmd('delete volume 2 4')
+
+cubit.cmd('merge all')
+cubit.cmd('imprint all')
+
+# Meshing the volumes
+#elementsize = 1196.4 #hi-resolution
+#elementsize = 1500.0 # mid-resolution
+elementsize = 3000.0 # low-resolution
+
+cubit.cmd('volume 3 size '+str(elementsize))
+cubit.cmd('mesh volume 3')
+
+#cubit.cmd('refine surface 8 numsplit 1 bias 1.0 depth 1')
+cubit.cmd('volume 1 size '+str(elementsize))
+cubit.cmd('mesh volume 1')
+
+cubit.cmd('volume 5 size '+str(elementsize))
+cubit.cmd('mesh volume 5')
+
+#### 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 "acoustic" ') # acoustic material region
+cubit.cmd('block 1 attribute count 4')
+cubit.cmd('block 1 attribute index 1 1 ') # material 1
+cubit.cmd('block 1 attribute index 2 1480 ') # vp
+cubit.cmd('block 1 attribute index 3 0 ') # vs
+cubit.cmd('block 1 attribute index 4 1028 ') # rho (ocean salt water density:
+ # http://www.windows.ucar.edu/tour/link=/earth/Water/density.html
+
+cubit.cmd('block 2 name "elastic" ') # elastic material region
+cubit.cmd('block 2 attribute count 6')
+cubit.cmd('block 2 attribute index 1 2 ') # material 2
+cubit.cmd('block 2 attribute index 2 7500 ') # vp
+cubit.cmd('block 2 attribute index 3 4300 ') # vs
+cubit.cmd('block 2 attribute index 4 3200 ') # rho
+cubit.cmd('block 2 attribute index 5 6') # Q_flag
+cubit.cmd('block 2 attribute index 6 0 ') # anisotropy_flag
+
+cubit.cmd('block 3 name "elastic" ') # elastic material region
+cubit.cmd('block 3 attribute count 6')
+cubit.cmd('block 3 attribute index 1 3 ') # same properties as material 2
+cubit.cmd('block 3 attribute index 2 7500 ') # vp
+cubit.cmd('block 3 attribute index 3 4300 ') # vs
+cubit.cmd('block 3 attribute index 4 3200 ') # rho
+cubit.cmd('block 3 attribute index 5 6') # Q_flag
+cubit.cmd('block 3 attribute index 6 0 ') # anisotropy_flag
+
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+
+#### 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
+
+
+
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8-nodoubling.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,109 @@
+#!/usr/bin/env python
+
+###########################################################################
+#### TNM: This is the mesh generation, adapted from a journal file
+#### specific to the settings of Komatitsch and Tromp 1999, Fig.8
+#### Aug 2009
+###########################################################################
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+cubit.cmd('reset')
+cubit.cmd('brick x 134000 y 134000 z 60000')
+
+# This seems to conflict with boundary_definition.py
+# ....which leaves the model space at e.g. x=[-67,67] km
+cubit.cmd('volume 1 move x 67000 y 67000 z -30000')
+
+# create vertices for discontinuity
+cubit.cmd('split curve 9 distance 3000')
+cubit.cmd('split curve 10 distance 3000')
+cubit.cmd('split curve 11 distance 3000')
+cubit.cmd('split curve 12 distance 3000')
+
+# create surface for interface
+cubit.cmd('create surface vertex 9 10 12 11')
+
+cubit.cmd('section volume 1 with surface 7 keep normal')
+cubit.cmd('section volume 1 with surface 7 reverse')
+
+# create vertices for auxiliary interface to allow for refinement
+cubit.cmd('split curve 29 distance 9000')
+cubit.cmd('split curve 31 distance 9000')
+cubit.cmd('split curve 32 distance 9000')
+cubit.cmd('split curve 36 distance 9000')
+
+# create surface for buffer interface to refine BELOW the discontinuity
+cubit.cmd('create surface vertex 25 26 28 27')
+
+cubit.cmd('section volume 3 with surface 19 keep normal')
+cubit.cmd('section volume 3 with surface 19 reverse')
+
+cubit.cmd('delete volume 2 4')
+
+cubit.cmd('merge all')
+cubit.cmd('imprint all')
+
+# Meshing the volumes
+cubit.cmd('volume 3 size 3589.2')
+cubit.cmd('mesh volume 3')
+
+cubit.cmd('refine surface 8 numsplit 1 bias 1.0 depth 1')
+
+cubit.cmd('volume 1 size 1196.4')
+cubit.cmd('mesh volume 1')
+
+cubit.cmd('volume 5 size 4785.71')
+cubit.cmd('mesh volume 5')
+
+#### 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 attribute count 6')
+cubit.cmd('block 1 attribute index 1 1 ') # volume 1
+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 6 ') # Q_flag
+cubit.cmd('block 1 attribute index 6 0 ') # anisotropy_flag
+
+cubit.cmd('block 2 attribute count 6')
+cubit.cmd('block 2 attribute index 1 2 ') # volume 2
+cubit.cmd('block 2 attribute index 2 7500 ')
+cubit.cmd('block 2 attribute index 3 4300 ')
+cubit.cmd('block 2 attribute index 4 3200 ')
+cubit.cmd('block 2 attribute index 5 6 ')
+cubit.cmd('block 2 attribute index 6 0 ') # anisotropy_flag
+
+cubit.cmd('block 3 attribute count 6')
+cubit.cmd('block 3 attribute index 1 2 ') # same properties as for volume 2
+cubit.cmd('block 3 attribute index 2 7500 ')
+cubit.cmd('block 3 attribute index 3 4300 ')
+cubit.cmd('block 3 attribute index 4 3200 ')
+cubit.cmd('block 3 attribute index 5 6 ')
+cubit.cmd('block 3 attribute index 6 0 ') # anisotropy_flag
+
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+
+#### 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
+
+
+
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,60 @@
+#!/usr/bin/env python
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+# two volumes separating 134000x134000x60000 block in middle vertically
+cubit.cmd('reset')
+cubit.cmd('brick x 67000 y 134000 z 60000')
+cubit.cmd('volume 1 move x 33500 y 67000 z -30000')
+cubit.cmd('brick x 67000 y 134000 z 60000')
+cubit.cmd('volume 2 move x 100500 y 67000 z -30000')
+cubit.cmd('merge all')
+
+# Meshing the volumes
+elementsize = 3000.0
+
+cubit.cmd('volume 1 size '+str(elementsize))
+cubit.cmd('volume 2 size '+str(elementsize))
+cubit.cmd('mesh volume 1 2')
+
+
+#### 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 "acoustic" ') # elastic material region
+cubit.cmd('block 1 attribute count 4')
+cubit.cmd('block 1 attribute index 1 1') # flag for material: 1 for 1. material
+cubit.cmd('block 1 attribute index 2 3000') # vp
+cubit.cmd('block 1 attribute index 3 0') # vs
+cubit.cmd('block 1 attribute index 4 2300') # rho
+
+cubit.cmd('block 2 name "elastic" ') # elastic material region
+cubit.cmd('block 2 attribute count 6')
+cubit.cmd('block 2 attribute index 1 2') # flag for material: 2 for 2. material
+cubit.cmd('block 2 attribute index 2 3000') # vp
+cubit.cmd('block 2 attribute index 3 0.1') # vs
+cubit.cmd('block 2 attribute index 4 2300') # rho
+cubit.cmd('block 2 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+cubit.cmd('block 2 attribute index 6 0 ') # anisotropy_flag
+
+
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+
+#### 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
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_mesh_boundary_vertical.py
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py 2009-11-21 02:26:51 UTC (rev 16024)
@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+
+###########################################################################
+#### TNM: This is the mesh generation, adapted from a journal file
+#### specific to the settings of Komatitsch and Tromp 1999, Fig.8
+#### Aug 2009
+###########################################################################
+#
+# modified version: top layer is water
+#
+###########################################################################
+
+import cubit
+import boundary_definition
+import cubit2specfem3d
+
+import os
+import sys
+
+cubit.cmd('reset')
+cubit.cmd('brick x 134000 y 134000 z 60000')
+
+# This seems to conflict with boundary_definition.py
+# ....which leaves the model space at e.g. x=[-67,67] km
+cubit.cmd('volume 1 move x 67000 y 67000 z -30000')
+
+# create vertices for discontinuity
+cubit.cmd('split curve 9 distance 3000')
+cubit.cmd('split curve 10 distance 3000')
+cubit.cmd('split curve 11 distance 3000')
+cubit.cmd('split curve 12 distance 3000')
+
+# create surface for interface
+cubit.cmd('create surface vertex 9 10 12 11')
+
+cubit.cmd('section volume 1 with surface 7 keep normal')
+cubit.cmd('section volume 1 with surface 7 reverse')
+
+# create vertices for auxiliary interface to allow for refinement
+cubit.cmd('split curve 29 distance 9000')
+cubit.cmd('split curve 31 distance 9000')
+cubit.cmd('split curve 32 distance 9000')
+cubit.cmd('split curve 36 distance 9000')
+
+# create surface for buffer interface to refine BELOW the discontinuity
+cubit.cmd('create surface vertex 25 26 28 27')
+
+cubit.cmd('section volume 3 with surface 19 keep normal')
+cubit.cmd('section volume 3 with surface 19 reverse')
+
+cubit.cmd('delete volume 2 4')
+
+cubit.cmd('merge all')
+cubit.cmd('imprint all')
+
+# Meshing the volumes
+#elementsize = 1196.4 #hi-resolution
+#elementsize = 1500.0 # mid-resolution
+elementsize = 3000.0 # low-resolution
+
+cubit.cmd('volume 3 size '+str(elementsize))
+cubit.cmd('mesh volume 3')
+
+#cubit.cmd('refine surface 8 numsplit 1 bias 1.0 depth 1')
+cubit.cmd('volume 1 size '+str(elementsize))
+cubit.cmd('mesh volume 1')
+
+cubit.cmd('volume 5 size '+str(elementsize))
+cubit.cmd('mesh volume 5')
+
+#### 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 "acoustic" ') # acoustic material region
+cubit.cmd('block 1 attribute count 4')
+cubit.cmd('block 1 attribute index 1 1 ') # material 1
+cubit.cmd('block 1 attribute index 2 1480 ') # vp
+cubit.cmd('block 1 attribute index 3 0 ') # vs
+cubit.cmd('block 1 attribute index 4 1028 ') # rho (ocean salt water density:
+ # http://www.windows.ucar.edu/tour/link=/earth/Water/density.html
+
+cubit.cmd('block 2 name "acoustic" ') # acoustic material region
+cubit.cmd('block 2 attribute count 4')
+cubit.cmd('block 2 attribute index 1 2 ') # material 1
+cubit.cmd('block 2 attribute index 2 1480 ') # vp
+cubit.cmd('block 2 attribute index 3 0 ') # vs
+cubit.cmd('block 2 attribute index 4 1028 ') # rho (ocean salt water density:
+
+cubit.cmd('block 3 name "acoustic" ') # acoustic material region
+cubit.cmd('block 3 attribute count 4')
+cubit.cmd('block 3 attribute index 1 3 ') # material 1
+cubit.cmd('block 3 attribute index 2 1480 ') # vp
+cubit.cmd('block 3 attribute index 3 0 ') # vs
+cubit.cmd('block 3 attribute index 4 1028 ') # rho (ocean salt water density:
+
+cubit.cmd('export mesh "top.e" dimension 3 overwrite')
+cubit.cmd('save as "meshing.cub" overwrite')
+
+#### 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
+
+
+
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/waterlayered_halfspace/waterlayer_only-nodoubling.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list