Merge pull request #1711 from danielpeter/devel

code cleaning & cubit script updates
SPECFEM · Jul 8, 2024 · 1b0d98c · 1b0d98c
2 parents acd52e3 + b17f716
commit 1b0d98c
Show file tree

Hide file tree

Showing 20 changed files with 142 additions and 89 deletions.
diff --git a/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py b/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py
@@ -652,14 +652,16 @@ def block_definition(self):
             self.bc = bc
             print('HEX Blocks:')
             for m, f in zip(self.block_mat, self.block_flag):
-                print('block ', m, 'material flag ', f)
+                print('  block ', m, 'material flag ', f)
             print('Absorbing Boundary Conditions:')
             for m, f in zip(self.block_bc, self.block_bc_flag):
-                print('bc ', m, 'bc flag ', f)
-            print('Topography (free surface)')
-            print(self.topography)
-            print('Free surface')
-            print(self.free)
+                print('  bc ', m, 'bc flag ', f)
+            if self.topography:
+                print('Topography (free surface):')
+                print('  ',self.topography)
+            if self.free:
+                print('Free surface:')
+                print('  ',self.free)
         except:
             print('****************************************')
             print('sorry, no blocks or blocks not properly defined')
@@ -765,7 +767,7 @@ def mat_parameter(self, properties):
         else:
             raise RuntimeError('Error: material id must be strictly positive or negative, but not equal to 0')
         # info output
-        print("material: ",txt)
+        print("  material: ",txt)
         return txt
 
     def nummaterial_write(self, nummaterial_name, placeholder=True):
@@ -853,8 +855,7 @@ def create_hexnode_string(self, hexa, hexnode_string=True):
         else:
             map(int, txt.split())
 
-    def create_facenode_string(self, hexa, face, normal=None, cknormal=True,
-                               facenode_string=True):
+    def create_facenode_string(self, hexa, face, normal=None, cknormal=True, facenode_string=True):
         nodes = self.get_face_connectivity(face)
         if cknormal:
             nodes_ok = self.normal_check(nodes[0:4], normal)
@@ -886,7 +887,7 @@ def mesh_write(self, mesh_name):
         meshfile.write(str(num_elems) + '\n')
         for block, flag in zip(self.block_mat, self.block_flag):
             hexes = cubit.get_block_hexes(block)
-            print('block ', block, ' hexes ', len(hexes))
+            print('  block ', block, ' hexes ', len(hexes))
             for hexa in hexes:
                 txt = self.create_hexnode_string(hexa)
                 meshfile.write(txt)
@@ -897,7 +898,7 @@ 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):
-            print('block ', block, 'flag ', flag)
+            print('  block ', block, 'flag ', flag)
             hexes = cubit.get_block_hexes(block)
             for hexa in hexes:
                 mat.write(('%10i %10i\n') % (hexa, flag))
@@ -949,8 +950,7 @@ def free_write(self, freename=None):
         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('free surface (topography) block name:', \
-                      name, 'id:', block)
+                print('  free surface (topography): block name:', name, 'id:', block)
                 quads_all = cubit.get_block_faces(block)
                 print('  number of faces = ', len(quads_all))
                 dic_quads_all = dict(zip(quads_all, quads_all))
@@ -968,13 +968,13 @@ def free_write(self, freename=None):
                     faces = cubit.get_sub_elements('hex', h, 2)
                     for f in faces:
                         if f in s:
-                            txt = self.create_facenode_string(
-                                h, f, normal, cknormal=True)
+                            # topography surface: checks face orientation against vertical normal (0,0,1)
+                            txt = self.create_facenode_string(h, f, normal, cknormal=True)
                             freehex.write(txt)
                 print('    100 %')
             elif block == self.free:
                 name = cubit.get_exodus_entity_name('block', block)
-                print('free surface block name:', name, 'id:', block)
+                print('  free surface: block name:', name, 'id:', block)
                 quads_all = cubit.get_block_faces(block)
                 print('  number of faces = ', len(quads_all))
                 dic_quads_all = dict(zip(quads_all, quads_all))
@@ -992,8 +992,8 @@ def free_write(self, freename=None):
                     faces = cubit.get_sub_elements('hex', h, 2)
                     for f in faces:
                         if f in s:
-                            txt = self.create_facenode_string(
-                                h, f, normal, cknormal=False)
+                            # free surface: no face orientation check against vertical normal
+                            txt = self.create_facenode_string(h, f, normal, cknormal=False)
                             freehex.write(txt)
                 print('    100 %')
         freehex.close()
@@ -1122,60 +1122,60 @@ def abs_write(self, absname=None):
             # loops through all block definitions
             list_hex = cubit.parse_cubit_list('hex', 'all')
             for block, flag in zip(self.block_bc, self.block_bc_flag):
-                if block != self.topography:
+                if block != self.topography and block != self.free:
                     name = cubit.get_exodus_entity_name('block', block)
                     print('  block name:', name, 'id:', block)
                     cknormal = True
                     abshex_local = False
                     # opens file
                     if re.search('xmin', name):
-                        print("xmin")
+                        print("  xmin")
                         abshex_local = open(absname + '_xmin', 'w')
                         normal = (-1, 0, 0)
                     elif re.search('xmax', name):
-                        print("xmax")
+                        print("  xmax")
                         abshex_local = open(absname + '_xmax', 'w')
                         normal = (1, 0, 0)
                     elif re.search('ymin', name):
-                        print("ymin")
+                        print("  ymin")
                         abshex_local = open(absname + '_ymin', 'w')
                         normal = (0, -1, 0)
                     elif re.search('ymax', name):
-                        print("ymax")
+                        print("  ymax")
                         abshex_local = open(absname + '_ymax', 'w')
                         normal = (0, 1, 0)
                     elif re.search('bottom', name):
-                        print("bottom")
+                        print("  bottom")
                         abshex_local = open(absname + '_bottom', 'w')
                         normal = (0, 0, -1)
                     elif re.search('abs', name):
-                        print("abs all - experimental, check the output")
+                        print("  abs all - experimental, check the output")
                         cknormal = False
                         abshex_local = open(absname, 'w')
                         normal = None
                     else:
                         if block == 1003:
-                            print("xmin")
+                            print("  xmin")
                             abshex_local = open(absname + '_xmin', 'w')
                             normal = (-1, 0, 0)
                         elif block == 1004:
-                            print("ymin")
+                            print("  ymin")
                             abshex_local = open(absname + '_ymin', 'w')
                             normal = (0, -1, 0)
                         elif block == 1005:
-                            print("xmax")
+                            print("  xmax")
                             abshex_local = open(absname + '_xmax', 'w')
                             normal = (1, 0, 0)
                         elif block == 1006:
-                            print("ymax")
+                            print("  ymax")
                             abshex_local = open(absname + '_ymax', 'w')
                             normal = (0, 1, 0)
                         elif block == 1002:
-                            print("bottom")
+                            print("  bottom")
                             abshex_local = open(absname + '_bottom', 'w')
                             normal = (0, 0, -1)
                         elif block == 1000:
-                            print("custom")
+                            print("  custom")
                             abshex_local = open(absname, 'w')
                             cknormal = False
                             normal = None
@@ -1199,8 +1199,7 @@ def abs_write(self, absname=None):
                             faces = cubit.get_sub_elements('hex', h, 2)
                             for f in faces:
                                 if f in s:
-                                    txt = self.create_facenode_string(
-                                        h, f, normal=normal, cknormal=cknormal)
+                                    txt = self.create_facenode_string(h, f, normal=normal, cknormal=cknormal)
                                     abshex_local.write(txt)
                         abshex_local.close()
                         print('    100 %')
@@ -1215,19 +1214,23 @@ def surface_write(self, pathdir=None):
         #  > block 10 name 'moho_surface'
         import re
         for block in self.block_bc:
-            if block != self.topography:
+            if block != self.topography and block != self.free:
                 name = cubit.get_exodus_entity_name('block', block)
-                # skips block names like face_abs**, face_topo**
+                # skips block names like face_abs**, face_topo**, face_free**
                 if re.search('abs', name):
                     continue
                 elif re.search('topo', name):
                     continue
+                elif re.search('free', name):
+                    continue
                 elif re.search('surface', name):
+                    # e.g. moho_surface
                     filename = pathdir + name + '_file'
                 else:
                     continue
                 # gets face elements
-                print('  surface block name: ', name, 'id: ', block)
+                print('Optional surface:')
+                print('  surface: block name: ', name, 'id: ', block)
                 quads_all = cubit.get_block_faces(block)
                 print('  face = ', len(quads_all))
                 if len(quads_all) == 0:
@@ -1243,17 +1246,16 @@ def surface_write(self, pathdir=None):
                 Nhex = len(list_hex)
                 # Create a set to speed up "if f in dic_quads_all.keys()" big time
                 s = set(dic_quads_all.keys())
-                percentageDone = -1
-                for h in list_hex:
+                lastDisplayed = -1
+                for ih, h in enumerate(list_hex):
                     percentageDone = float(ih)/Nhex*100.0
                     if int(percentageDone)%10 == 0 and lastDisplayed != int(percentageDone):
                         print('    ', int(percentageDone), '%')
                         lastDisplayed = int(percentageDone)
                     faces = cubit.get_sub_elements('hex', h, 2)
                     for f in faces:
                         if f in s:
-                            txt = self.create_facenode_string(
-                                h, f, cknormal=False)
+                            txt = self.create_facenode_string(h, f, cknormal=False)
                             surfhex_local.write(txt)
                 # closes file
                 surfhex_local.close()
@@ -1388,7 +1390,7 @@ def _get_bc_flag(self, block):
             normal = (0, 0, -1)
         elif re.search('abs', name):
             label = "all"
-            print("abs all - experimental, check the output")
+            print("  abs all - experimental, check the output")
             cknormal = False
         else:
             if block == 1003:
@@ -1424,9 +1426,7 @@ def _get_bc_flag(self, block):
 #            faces = cubit.get_sub_elements('hex', h, 2)
 #            for f in faces:
 #                if dic_quads_all.has_key(f):
-#                    abs_array.append(self.create_facenode_string(
-#                        h, f, normal=normal, cknormal=cknormal,
-#                        facenode_string=False))
+#                    abs_array.append(self.create_facenode_string(h, f, normal=normal, cknormal=cknormal, facenode_string=False))
 #        print('Ok')
 #        cubit.cd('set info on')
 #        cubit.cmd('set echo on')
@@ -1474,9 +1474,7 @@ def _get_bc_flag(self, block):
 #                    for f in faces:
 #                        if dic_quads_all.has_key(f):
 #                            # print(f)
-#                            free_array.append(self.create_facenode_string(
-#                                h, f, normal, cknormal=True,
-#                                facenode_string=False))
+#                            free_array.append(self.create_facenode_string(h, f, normal, cknormal=True, facenode_string=False))
 #            elif block == self.free:
 #                name = cubit.get_exodus_entity_name('block', block)
 #                print('free surface block name:', name, 'id:', block)
@@ -1488,9 +1486,7 @@ def _get_bc_flag(self, block):
 #                    faces = cubit.get_sub_elements('hex', h, 2)
 #                    for f in faces:
 #                        if dic_quads_all.has_key(f):
-#                            free_array.append(self.create_facenode_string(
-#                                h, f, normal, cknormal=False,
-#                                facenode_string=False))
+#                            free_array.append(self.create_facenode_string(h, f, normal, cknormal=False, facenode_string=False))
 #        print('Ok')
 #        cubit.cmd('set info on')
 #        cubit.cmd('set echo on')

diff --git a/EXAMPLES/applications/Mount_StHelens/mesh_mount.py b/EXAMPLES/applications/Mount_StHelens/mesh_mount.py
@@ -158,8 +158,8 @@
 else:
     print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
     print("# imprinting volume, this will take around 1 min, please be patience...")
-    cubit.cmd('merge all')
     cubit.cmd('imprint all')
+    cubit.cmd('merge all')
 
     # exports only surfaces which will create single volume
     cubit.cmd('export acis "topo_2.acis" surface 3 10 12 14 15 9 ascii overwrite')

diff --git a/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py b/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py
@@ -132,7 +132,7 @@
 # imprints topography surfaces into brick volume, creates connected surfaces
 # note: this is a develop feature
 cubit.cmd('imprint all')
-
+cubit.cmd('merge all')
 
 cubit.cmd('version')
 

diff --git a/EXAMPLES/applications/fault_examples/tpv5/TPV5.py b/EXAMPLES/applications/fault_examples/tpv5/TPV5.py
@@ -131,6 +131,7 @@
 
 cubit.cmd("imprint all")
 cubit.cmd("merge all")
+
 cubit.cmd("surface 1 size "+str(elementsize))
 cubit.cmd("volume 1 size "+str(elementsize))
 cubit.cmd("surface 1 scheme pave")

diff --git a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py
@@ -94,8 +94,8 @@
 
 cubit.cmd('delete volume 2 4')
 
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 cubit.cmd('volume 3 size '+str(elementsize))
 cubit.cmd('mesh volume 3')

diff --git a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py
@@ -94,8 +94,8 @@
 
 cubit.cmd('delete volume 2 4')
 
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 # Meshing the volumes
 ## middle volume

diff --git a/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py b/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py
@@ -88,8 +88,8 @@
 
 cubit.cmd('delete volume 2 4')
 
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 # Meshing the volumes
 cubit.cmd('volume 3 size '+str(elementsize))

diff --git a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py
@@ -68,8 +68,8 @@
 
 cubit.cmd('delete volume 2 4')
 
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 # Meshing the volumes
 cubit.cmd('volume 3 size 3589.2')

diff --git a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py
@@ -64,8 +64,8 @@
 
 cubit.cmd('delete volume 2 4')
 
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 # Meshing the volumes
 #elementsize = 1196.4 #hi-resolution

diff --git a/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py b/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py
@@ -95,8 +95,8 @@
     cubit.cmd('delete volume 4')
 
 # merges surfaces
-cubit.cmd('merge all')
 cubit.cmd('imprint all')
+cubit.cmd('merge all')
 
 # resets numbering
 cubit.cmd('compress all')

diff --git a/src/decompose_mesh/part_decompose_mesh.F90 b/src/decompose_mesh/part_decompose_mesh.F90
@@ -620,7 +620,7 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
   implicit none
   integer, intent(in)  :: iproc
   integer, intent(in)  :: nspec
-  integer, intent(in)  :: nb_wd 
+  integer, intent(in)  :: nb_wd
   integer, dimension(:), pointer :: glob2loc_elmnts
   integer, dimension(1:nspec)  :: part
   integer, dimension(1:nb_wd), intent(in)  :: boundary_to_ispec_wd, side_wd
@@ -632,7 +632,7 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
   character(len=MAX_STRING_LEN) :: prname
   write(prname, "('proc',i6.6,'_')") iproc
   open(unit=IFILE_WAVEFIELD_DISCONTINUITY, &
-       file=trim(outputpath_name)//'/'//trim(prname)//&
+       file = trim(outputpath_name)//'/'//trim(prname)//&
             trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), &
        form='unformatted', action='write')
   local_nb_wd = 0

diff --git a/src/decompose_mesh/read_mesh_files.F90 b/src/decompose_mesh/read_mesh_files.F90
@@ -39,7 +39,7 @@ subroutine read_mesh_files()
   use fault_scotch, only: ANY_FAULT,read_fault_files,save_nodes_coords,close_faults
 
   use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE
-  
+
   implicit none
 
   ! local parameters
@@ -592,7 +592,7 @@ subroutine read_mesh_files()
       read(IIN_DB, *, iostat=ier) boundary_to_ispec_wd(ib), side_wd(ib)
     enddo
     close(IIN_DB)
-  endif  
+  endif
 
   ! reads in absorbing boundary files
   open(unit=IIN_DB, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', &