Hi @kgmerkley,
Thank you for your assistance. Inspired by the code you provided, I successfully implemented the mesh generation for a model featuring a cube containing an upper hemisphere. However, as you noted, the SPECFEM3D code requires 8-node hexahedral (HEX8) meshes as input. Using the following Python script, I generated 11 model files and used them as inputs for SPECFEM3D. During execution, an error occurred: “Error while attempting to read 8 element data values from the mesh file” , indicating that the mesh file does not contain HEX8 elements. In this situation, how can I ensure that the output mesh consists of 8-node hexahedral elements?
My Python code is as follows:
#!/usr/bin/env python
from future import print_function
import cubit
cubit.init([""])
try:
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d
except:
import boundary_definition
import cubit2specfem3d
import os
import sys
cubit.cmd(‘reset’)
cubit.cmd(‘Timer Start’)
cubit.cmd(‘clear’)
cubit.cmd(‘Color Background white’)
cubit.cmd(‘Color Axis Labels red’)
cubit.cmd(‘graphics perspective off’)
cubit.cmd(‘graphics parallel scale 303.79899’)
cubit.cmd(‘from -2209.0335 -2399.8038 840.56353’)
cubit.cmd(‘at 0 0 0’)
cubit.cmd(‘up 0.15018128 0.20082845 0.96804624’)
cubit.cmd(‘set duplicate block elements on’)
cubit.cmd(‘view iso’)
cubit.cmd(‘up 0 0 1’)
elementsize = 4.0
path_folder = ‘model_files_3D_20250812’
cubit.cmd(‘brick x 200 y 200 z 200’)
cubit.cmd(‘move vol 1 x 100 y 100 z -100’)
cubit.cmd(‘create sphere radius 30 Zpositive’)
cubit.cmd(‘move volume 2 x 100 y 100 z -150’)
cubit.cmd(‘subtract volume 2 from 1 keep_tool’)
cubit.cmd(‘volume 2 name “hemi”’)
cubit.cmd(‘volume 3 name “outer”’)
cubit.cmd(‘webcut volume all with cylinder radius 20 axis z center 100 100 -150’)
cubit.cmd(‘volume with name “outer@*” name “cyl_cut”’)
cubit.cmd(‘volume with name “hemi@*” name “cyl_cut”’)
cubit.cmd(‘webcut volume all except with name “cyl_cut*” with plane xplane offset 100’)
cubit.cmd(‘webcut volume all except with name “cyl_cut*” with plane yplane offset 100’)
cubit.cmd(‘webcut volume all with plane zplane offset -150’)
cubit.cmd(‘webcut volume all with cylinder radius 40 axis z center 100 100 -150’)
cubit.cmd(‘webcut volume 21 to 24 with plane zplane offset -132’)
cubit.cmd(‘color volume with name “hemi*” blue’)
cubit.cmd(‘color volume with name “outer*” id 4’)
cubit.cmd(‘draw volume all’)
cubit.cmd(‘imprint all’)
cubit.cmd(‘merge all’)
cubit.cmd('volume all size ’ + str(elementsize))
cubit.cmd(‘curve with radius < 4 interval 6’)
cubit.cmd(‘surface with z_coord < -199 scheme pave’)
cubit.cmd(‘surface with area < 100 with area > 20 scheme map’)
cubit.cmd(‘mesh volume all’)
cubit.cmd(‘Geometry visibility ON’)
cubit.cmd(‘Mesh visibility ON’)
cubit.cmd(‘Surface 25 144 151 154 162 164 172 174 182 Name “Zmax”’)
cubit.cmd(‘Surface 26 185 192 195 202 207 209 215 222 Name “Zmin”’)
cubit.cmd(‘Surface 108 115 129 130 Name “Xmax”’)
cubit.cmd(‘Surface 117 123 136 142 Name “Xmin”’)
cubit.cmd(‘Surface 126 133 135 141 Name “Ymax”’)
cubit.cmd(‘Surface 111 112 118 124 Name “Ymin”’)
cubit.cmd(‘block 1001 face in surface with name “Zmax@*”’)
cubit.cmd(‘block 1001 name “face_topo”’)
cubit.cmd(‘block 1002 face in surface with name “Zmin@*”’)
cubit.cmd(‘block 1002 name “face_bottom”’)
cubit.cmd(‘block 1003 face in surface with name “Xmin@*”’)
cubit.cmd(‘block 1003 name “face_abs_xmin”’)
cubit.cmd(‘block 1004 face in surface with name “Ymin@*”’)
cubit.cmd(‘block 1004 name “face_abs_ymin”’)
cubit.cmd(‘block 1005 face in surface with name “Xmax@*”’)
cubit.cmd(‘block 1005 name “face_abs_xmax”’)
cubit.cmd(‘block 1006 face in surface with name “Ymax@*”’)
cubit.cmd(‘block 1006 name “face_abs_ymax”’)
cubit.cmd(‘block 1 hex in volume with name “hemi@*”’)
cubit.cmd('block 1 name "elastic 1 " ') # material region
cubit.cmd(‘block 1 attribute count 7’)
cubit.cmd('block 1 attribute index 1 1 ') # volume 2
cubit.cmd('block 1 attribute index 2 3800 ') # vp
cubit.cmd('block 1 attribute index 3 2220 ')# vs
cubit.cmd('block 1 attribute index 4 2250 ') # rho
cubit.cmd('block 1 attribute index 5 9999.0 ')# Qkappa
cubit.cmd('block 1 attribute index 6 9999.0 ')# Qmu
cubit.cmd('block 1 attribute index 7 0 ') # anisotropy_flag
cubit.cmd(‘block 2 hex in volume with name “outer@*”’)
cubit.cmd('block 2 name “elastic 2” ') # material region
cubit.cmd(‘block 2 attribute count 7’)
cubit.cmd('block 2 attribute index 1 2 ') # volume 2
cubit.cmd('block 2 attribute index 2 4000 ') # vp
cubit.cmd('block 2 attribute index 3 2350 ')# vs
cubit.cmd('block 2 attribute index 4 2300 ') # rho
cubit.cmd('block 2 attribute index 5 9999.0 ')# Qkappa
cubit.cmd('block 2 attribute index 6 9999.0 ')# Qmu
cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d
os.system('mkdir ’ + path_folder)
cubit.cmd(‘export mesh "’ + path_folder + ‘/top.e" dimension 3 overwrite’)
cubit.cmd(‘save as "’ + path_folder + ‘/meshing.cub" overwrite’)
cubit.cmd(‘save as "’ + path_folder + ‘/top.cub" overwrite’)
cubit2specfem3d.export2SPECFEM3D(path_folder, cpml = True, cpml_size = 20, top_absorbing = False)