Skip to content

Commit

Permalink
fixes mesh exporting script for poroelastic domains
Browse files Browse the repository at this point in the history
  • Loading branch information
danielpeter committed Apr 14, 2020
1 parent adc47a7 commit 0caec10
Show file tree
Hide file tree
Showing 18 changed files with 179 additions and 25 deletions.
10 changes: 5 additions & 5 deletions CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py
Original file line number Diff line number Diff line change
Expand Up @@ -502,11 +502,11 @@ def block_definition(self):
qmu = 9999.0
ani = 0
# material domain id
if "acoustic" in name:
if name.startswith('acoustic'):
imaterial = 1
elif "elastic" in name:
elif name.startswith('elastic'):
imaterial = 2
elif "poroelastic" in name:
elif name.startswith('poroelastic'):
imaterial = 3
else:
imaterial = 0
Expand Down Expand Up @@ -761,7 +761,7 @@ def nummaterial_write(self, nummaterial_name, placeholder=True):
! #(1)domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_k #(7)Q_mu #(8)ani
!
! where
! domain_id : 1=acoustic / 2=elastic
! domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : POSITIVE integer identifier of material block
! rho : density
! vp : P-velocity
Expand All @@ -778,7 +778,7 @@ def nummaterial_write(self, nummaterial_name, placeholder=True):
! #(1)domain_id #(2)material_id tomography elastic #(3)filename #(4)positive
!
! where
! domain_id : 1=acoustic / 2=elastic
! domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : NEGATIVE integer identifier of material block
! filename : filename of the tomography file
! positive : a positive unique identifier
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 3
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 2000 2000 1500 9999. 9999. 0 2
2 2000 2000 0 9999. 9999. 0 1
3 2000 2000 0.001 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 3
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 2000 5000 3000 9999. 9999. 0 2
2 2000 2000 0 9999. 9999. 0 1
3 2000 2000 0.001 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
! #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_kappa #(7)Q_mu #(8)anisotropy_flag
!
! where
! material_domain_id : 1=acoustic / 2=elastic
! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : POSITIVE integer identifier corresponding to the identifier of material block
! rho : density
! vp : P-velocity
Expand All @@ -24,7 +24,7 @@
! #(1)material_domain_id #(2)material_id tomography elastic #(3)tomography_filename #(4)positive_unique_number
!
! where
! material_domain_id : 1=acoustic / 2=elastic
! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : NEGATIVE integer identifier corresponding to the identifier of material block
! tomography_filename: filename of the tomography file
! positive_unique_number: a positive unique identifier
Expand Down
154 changes: 154 additions & 0 deletions EXAMPLES/homogeneous_poroelastic/create_mesh.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,154 @@
#!/usr/bin/env python
#
# "create_mesh.py" is a script that generates mesh specific to homogenous halfspace example
# i.e., a uniform mesh of 134 km x 134 km x 60 km with an element size 3.75 km.
# It is not applicable to other examples.
from __future__ import print_function

import os
import sys

# to run this script from command line, python must have its PATH environment set such that it
# includes the path to CUBIT/Trelis (cubit.py).
#
# you can also explicitly set it here e.g. like:
#sys.path.append('/opt/Trelis-15.0/bin/')

try:
import cubit
except ImportError:
print("Error: Importing cubit as python module failed")
print("could not import cubit, please check your PYTHONPATH settings...")
print("")
print("current path: ")
print(sys.path)
print("")
print("try to include path to directory which includes file cubit.py, e.g. /opt/Trelis-15.0/bin/")
print("")
sys.exit("Import cubit failed")

print(sys.path)

cubit.init([""])

# gets version string
cubit_version = cubit.get_version()
print("version: ",cubit_version)

# extracts major number
v = cubit_version.split('.')
cubit_version_major = int(v[0])
print("major version number: ",cubit_version_major)

# current work directory
cubit.cmd('pwd')

# Creating the volumes
cubit.cmd('reset')

# single volume
cubit.cmd('brick x 134000 y 134000 z 60000')
cubit.cmd('volume 1 move x 67000 y 67000 z -30000')

# two merged volumes
#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 = 3750.0

cubit.cmd('volume all size '+str(elementsize))
cubit.cmd('mesh volume all')

# End of meshing

#
# GEOCUBIT
#
# adds path to geocubit (if not setup yet)
sys.path.append('../../CUBIT_GEOCUBIT/')

print("path: ")
print(sys.path)
print("")

# avoids assigning empty blocks
cubit.cmd('set duplicate block elements on')

# creates MESH/ directory for file output
os.system('mkdir -p MESH')

# conversion to specfem-format
# use_explicit will explicitly assign material properties as block attributes
use_explicit=1

if use_explicit == 1:
from geocubitlib import boundary_definition
# bounding faces
boundary_definition.entities=['face']
boundary_definition.define_bc(boundary_definition.entities,parallel=True)
from geocubitlib import cubit2specfem3d
print("")
print("material properties: assigned as block attributes")
print("")
# sets the id of the volume block
# (volume block starts at id 4)
id_block = 1
print("cubit block:")
print(" volume block id = " + str(id_block))
print("")
# Define material properties
print("#### DEFINE MATERIAL PROPERTIES #######################")
# elastic material
cubit.cmd('block '+str(id_block)+' name "poroelastic 1" ') # poroelastic material region
cubit.cmd('block '+str(id_block)+' attribute count 7')
cubit.cmd('block '+str(id_block)+' attribute index 1 1') # flag for material: 1 for 1. material
cubit.cmd('block '+str(id_block)+' attribute index 2 3371.0') # vp
cubit.cmd('block '+str(id_block)+' attribute index 3 2128.0') # vs
cubit.cmd('block '+str(id_block)+' attribute index 4 2473.0') # rho
cubit.cmd('block '+str(id_block)+' attribute index 5 9999.0') # Qkappa
cubit.cmd('block '+str(id_block)+' attribute index 6 9999.0') # Qmu
cubit.cmd('block '+str(id_block)+' attribute index 7 0') # anisotropy_flag
print("")
print("exporting to SPECFEM3D-format:")
print("")
# Export to SPECFEM3D format
cubit2specfem3d.export2SPECFEM3D('MESH/')
# backup cubit
cubit.cmd('export mesh "MESH/top.e" dimension 3 overwrite')
cubit.cmd('save as "MESH/meshing.cub" overwrite')
else:
from geocubitlib import exportlib
print("")
print("exporting to SPECFEM3D-format:")
print("")
# Export to SPECFEM3D format
# note: exportlib-commands will overwrite material properties
exportlib.define_blocks(outdir='MESH/',save_cubfile=True,outfilename='top')
exportlib.e2SEM(outdir='MESH/')
# Define material properties
print("#### DEFINE MATERIAL PROPERTIES #######################")
# elastic material
material_cfg=[{'material region':'3','id_block':'1','vp':'3371.0','vs':'2128.0','rho':'2473.0','Qkappa':'9999.0','Qmu':'9999.0','anisotropy_flag':'0'}]
# modifies material file
nummaterial_velocity_file='MESH/nummaterial_velocity_file'
f=open(nummaterial_velocity_file,'w')
for block in material_cfg:
print(block)
s=block['material region']+' '
s=s+block['id_block']+' '
s=s+block['rho']+' '
s=s+block['vp']+' '
s=s+block['vs']+' '
s=s+block['Qkappa']+' '
s=s+block['Qmu']+' '
s=s+block['anisotropy_flag']
f.write(s+'\n')
f.close()

# backup cubit

# all files needed by xdecompose_mesh are now in directory MESH/
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 1
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 1450 391 230 9999. 9999. 0 2

#-----------------------------------------------------------
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 1200 1500 750 9999. 40 0 2
2 1100 1600 800 9999. 60 0 2
3 1000 1500 0 9999. 9999. 0 1
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 2
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 1000 2500 0 9999. 50.0 0 1
2 2000 4500 2500 9999. 40.0 0 2

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 1200 1500 750 9999. 40.0 0 2
2 1100 1600 0 9999. 50.0 0 1
3 1000 1500 700 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 1000.0 0 2
2 2800 6700 3870 9999. 500.0 0 2
3 2670 6300 3640 9999. 300.0 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ NMATERIALS = 4
# Q_Kappa : Q_Kappa attenuation quality factor
# Q_mu : Q_mu attenuation quality factor
# anisotropy_flag : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
# domain_id : 1 = acoustic / 2 = elastic
# domain_id : 1 = acoustic / 2 = elastic / 3 = poroelastic
1 3000 7800 4500 9999. 9999. 0 2
2 2800 6700 3870 9999. 9999. 0 2
3 2670 6300 3640 9999. 9999. 0 2
Expand Down
6 changes: 3 additions & 3 deletions EXAMPLES/tomographic_model/tomoblock_mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,11 @@
print("")

try:
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d
except:
import boundary_definition
import cubit2specfem3d
import cubit2specfem3d

# bounding faces
boundary_definition.entities=['face']
Expand Down
4 changes: 2 additions & 2 deletions utils/CPML/add_CPML_layers_to_an_existing_mesh.f90
Original file line number Diff line number Diff line change
Expand Up @@ -260,7 +260,7 @@ program add_CPML_layers_to_a_given_mesh
! #(1)material_domain_id #(2)material_id #(3)rho #(4)vp #(5)vs #(6)Q_kappa #(7)Q_mu #(8)anisotropy_flag
!
! where
! material_domain_id : 1=acoustic / 2=elastic
! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : POSITIVE integer identifier corresponding to the identifier of material block
! rho : density
! vp : P-velocity
Expand All @@ -281,7 +281,7 @@ program add_CPML_layers_to_a_given_mesh
! #(1)material_domain_id #(2)material_id tomography elastic #(3)tomography_filename #(4)positive_unique_number
!
! where
! material_domain_id : 1=acoustic / 2=elastic
! material_domain_id : 1=acoustic / 2=elastic / 3=poroelastic
! material_id : NEGATIVE integer identifier corresponding to the identifier of material block
! tomography_filename: filename of the tomography file
! positive_unique_number: a positive unique identifier
Expand Down

0 comments on commit 0caec10

Please sign in to comment.