The top boundary is configured with CPML, yet the generated file is empty

Hang_Duan · August 5, 2025, 8:50am

I want to create a 3D cubic model with an internal inclined surface using a Python script in CUBIT software, apply CPML (Perfectly Matched Layer) boundary conditions to all boundaries, and generate all-hexahedral meshes. However, after running the code, the file free_or_absorbing_surface_file_zmax in the exported mesh files is empty. How can I modify the code to ensure all boundaries are correctly configured with CPML conditions?

The original 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')

# -------------------- Reset Cubit environment -------------------- #


# -------------------- Basic settings start -------------------- #
# cubit.cmd('view journal')   
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 right')    
# cubit.cmd('up 0 0 1')                           

# cubit.cmd('view top')    

# cubit.cmd('view bottom') 

# cubit.cmd('view front') 

cubit.cmd('view iso')
cubit.cmd('up 0 0 1') 
# -------------------- Basic settings end -------------------- # 


# -------------------- Mesh size & folder setup -------------------- #
elementsize = 4.0
path_folder = 'model_files_3D_20250803'


# -------------------- Define 3D geometry start -------------------- #
cubit.cmd('brick 	  x 340 y 340 z  300')
cubit.cmd('move vol 1 x 170 y 170 z -150')

cubit.cmd('create planar surface zplane')           
cubit.cmd('rotate volume 2 about Y angle -15.000000')
cubit.cmd('move volume 2 x 170 y 170 z -180')
cubit.cmd('webcut volume 1 with sheet extended from surface 7')    
cubit.cmd('delete volume 2')
# -------------------- Define 3D geometry end -------------------- #


cubit.cmd('view iso')
cubit.cmd('up 0 0 1') 


# -------------------- Set colors -------------------- #
cubit.cmd('color vol 1  id 4')		
cubit.cmd('color vol 3 id 60')		


# -------------------- Merge & mesh start -------------------- #
cubit.cmd('merge vol all')
cubit.cmd('imprint vol all')
cubit.cmd('volume all size ' + str(elementsize))
cubit.cmd('mesh volume all')
cubit.cmd('draw volume all')
# -------------------- Merge & mesh end -------------------- #


# -------------------- Define boundaries start -------------------- #
# cubit.cmd('block 1001 face in surface 1   68    69  70  72')    
# # cubit.cmd('block 1001 face in surface 1')                     
# cubit.cmd('block 1001 name "face_topo"')

cubit.cmd('block 1002 face in surface 2')
cubit.cmd('block 1002 name "face_bottom"')

cubit.cmd('block 1003 face in surface 10	17')
cubit.cmd('block 1003 name "face_abs_xmin"')

cubit.cmd('block 1004 face in surface 11	16')
cubit.cmd('block 1004 name "face_abs_ymin"')

cubit.cmd('block 1005 face in surface 12	15')
cubit.cmd('block 1005 name "face_abs_xmax"')

cubit.cmd('block 1006 face in surface 9		14')
cubit.cmd('block 1006 name "face_abs_ymax"')
# -------------------- Define boundaries end -------------------- #


# -------------------- Block properties start -------------------- #
cubit.cmd('block 1 hex in volume 1 ')
cubit.cmd('block 1 name "elastic 1 " ')        
cubit.cmd('block 1 attribute count 7')
cubit.cmd('block 1 attribute index 1 1  ')      
cubit.cmd('block 1 attribute index 2 1800 ') 		
cubit.cmd('block 1 attribute index 3 800 ')
cubit.cmd('block 1 attribute index 4 2000 ') 
cubit.cmd('block 1 attribute index 5 9999.0 ')
cubit.cmd('block 1 attribute index 6 9999.0 ')
cubit.cmd('block 1 attribute index 7 0 ')     

cubit.cmd('block 2 hex in volume 3')
cubit.cmd('block 2 name "elastic 2" ')        
cubit.cmd('block 2 attribute count 7')
cubit.cmd('block 2 attribute index 1 2  ')      
cubit.cmd('block 2 attribute index 2 2500 ') 
cubit.cmd('block 2 attribute index 3 1500 ')
cubit.cmd('block 2 attribute index 4 2200 ') 
cubit.cmd('block 2 attribute index 5 9999.0 ')
cubit.cmd('block 2 attribute index 6 9999.0 ')
cubit.cmd('block 2 attribute index 7 0 ')     
# -------------------- Block properties end -------------------- #


# -------------------- Import modules -------------------- #
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d


# -------------------- Export files start -------------------- #
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 = True)
# -------------------- Export files end -------------------- #

Norbert_Hofbauer · August 6, 2025, 8:59pm

Hi @Hang_Duan,
the file free_or_absorbing_surface_file_zmax seems to be created by the line

The repo you are using is not maintained by us. You would need to create an issue or start a discussion on github to get some help for specfem.

Hang_Duan · August 11, 2025, 12:54am

Thank you for your explanations and reminders. I initially presumed that the question I posed would be a common and fundamental one within the CUBIT software community, merely relating to the application realm of SPECFEM. In such a circumstance, it will be rather arduous for me to secure assistance. Should I submit inquiries on GitHub, it is likely that some will contend that this issue pertains to CUBIT software rather than SPECFEM, and thus advise seeking support within the CUBIT software community.

Hang_Duan · August 12, 2025, 5:24pm

To this day, I still haven’t been able to come up with a solution to this problem, and I’m quite at a loss as to what to do next. I originally thought it was a simple issue, but I never expected it to turn out like this.

Norbert_Hofbauer · August 13, 2025, 8:31am

Have you already started to debug the specfem code?

The function you call should be this one.

github.com

SPECFEM/specfem3d/blob/master/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py#L411

      
        
                    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(mesh_tools):
            
            
    def __init__(self, hex27=False, cpml=False, cpml_size=False,
                             top_absorbing=False):
                    super(mesh, self).__init__()
                    self.netcdf_db = False
                    self._netcdf_num_nod_hex = 8
                    self._netcdf_num_nod_quad = 4
                    self.abs_block_ids = [1001, 1002, 1003, 1004, 1005, 1006]
                    self._netcdf_len_string = 33
                    self._netcdf_four = 4

Hang_Duan · August 13, 2025, 4:56pm

Hi @Norbert_Hofbauer ,
I haven’t started debugging yet. The SPECFEM3D code I’m using is for the Linux system. For the time being, I don’t know how to debug it, but it’s definitely not as easy to debug as other common codes for the Windows version.