moves routines to save mesh as cubit to file save_mesh_as_cubit.f90 (for in-house mesher); code cleaning

Author: danielpeter

src/decompose_mesh/part_decompose_mesh.F90

@@ -613,8 +613,8 @@ end subroutine  write_material_props_database
   ! discontinuity interface
   !-------------------------------------------------
   subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
-                                       nb_wd, boundary_to_ispec_wd, side_wd, &
-                                                nspec, glob2loc_elmnts, part)
+                                                    nb_wd, boundary_to_ispec_wd, side_wd, &
+                                                    nspec, glob2loc_elmnts, part)
   use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_MESH, &
                        IFILE_WAVEFIELD_DISCONTINUITY, MAX_STRING_LEN
   implicit none
@@ -630,11 +630,13 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
                                         local_side_wd
   character(len=MAX_STRING_LEN), intent(in) :: 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)//&
             trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), &
        form='unformatted', action='write')
+
   local_nb_wd = 0
   do ib = 1, nb_wd
     ispec = boundary_to_ispec_wd(ib)
@@ -643,7 +645,9 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
       local_nb_wd = local_nb_wd + 1
     endif
   enddo
+
   print *, 'partition', iproc, ' has ', local_nb_wd, ' elements on wavefield discontinuity interface'
+
   allocate(local_boundary_to_ispec_wd(local_nb_wd))
   allocate(local_side_wd(local_nb_wd))
   local_nb_wd = 0
@@ -656,11 +660,15 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, &
       local_side_wd(local_nb_wd) = iside
     endif
   enddo
+
   write(IFILE_WAVEFIELD_DISCONTINUITY) local_nb_wd
   write(IFILE_WAVEFIELD_DISCONTINUITY) local_boundary_to_ispec_wd
   write(IFILE_WAVEFIELD_DISCONTINUITY) local_side_wd
+
   close(IFILE_WAVEFIELD_DISCONTINUITY)
+
   deallocate(local_boundary_to_ispec_wd, local_side_wd)
+
   end subroutine write_wavefield_discontinuity_database
 
 

src/decompose_mesh/write_mesh_databases.F90

@@ -202,12 +202,12 @@ subroutine write_mesh_databases()
   endif
 
 
-  !! setting up wavefield discontinuity boundary
+  ! setting up wavefield discontinuity boundary
   if (IS_WAVEFIELD_DISCONTINUITY) then
     do ipart = 0, nparts-1
       call write_wavefield_discontinuity_database(ipart, outputpath_name, &
-                                      nb_wd, boundary_to_ispec_wd, side_wd, &
-                                                nspec, glob2loc_elmnts, part)
+                                                  nb_wd, boundary_to_ispec_wd, side_wd, &
+                                                  nspec, glob2loc_elmnts, part)
     enddo
   endif
 

src/meshfem3D/create_CPML_regions.f90

@@ -36,7 +36,7 @@ subroutine create_CPML_regions(nglob)
   ! create the different regions of the mesh
   use constants, only: IMAIN,CUSTOM_REAL,SMALL_PERCENTAGE_TOLERANCE, &
     CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, &
-    PI,TINYVAL,HUGEVAL,MAX_STRING_LEN,NDIM,myrank
+    PI,HUGEVAL,MAX_STRING_LEN,myrank
 
   use constants_meshfem, only: NGLLX_M,NGLLY_M,NGLLZ_M
 
@@ -434,11 +434,11 @@ subroutine create_CPML_regions(nglob)
     write(IMAIN,*)
     ! converts to degrees
     if (.not. SUPPRESS_UTM_PROJECTION) then
-      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along X sides (in degree) = ',&
+      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along X sides (in degree) = ', &
                    sngl(elem_size_x_min*METERS_TO_DEGREES),'/',sngl(elem_size_x_max*METERS_TO_DEGREES)
-      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along Y sides (in degree) = ',&
+      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along Y sides (in degree) = ', &
                    sngl(elem_size_y_min*METERS_TO_DEGREES),'/',sngl(elem_size_y_max*METERS_TO_DEGREES)
-      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along Z sides (in degree) = ',&
+      write(IMAIN,'(a,f8.5,a,f8.5)') '                       min/max element size along Z sides (in degree) = ', &
                    sngl(elem_size_z_min*METERS_TO_DEGREES),'/',sngl(elem_size_z_max*METERS_TO_DEGREES)
       write(IMAIN,*)
     endif
@@ -525,9 +525,9 @@ subroutine add_CPML_region_as_extra_layers(nglob)
     UTM_X_MIN,UTM_X_MAX
 
   ! create the different regions of the mesh
-  use constants, only: IMAIN,CUSTOM_REAL,SMALL_PERCENTAGE_TOLERANCE, &
+  use constants, only: IMAIN, &
     CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, &
-    PI,TINYVAL,HUGEVAL,MAX_STRING_LEN,NDIM,myrank, &
+    MAX_STRING_LEN,NDIM,myrank, &
     GAUSSALPHA,GAUSSBETA
 
   use constants_meshfem, only: NGLLX_M,NGLLY_M,NGLLZ_M
@@ -617,7 +617,7 @@ subroutine add_CPML_region_as_extra_layers(nglob)
   integer :: elem_mapping(NGNOD)
   logical :: is_mapped(NGNOD)
   ! MPI Cartesian topology uses W for West (= XI_MIN), E for East (= XI_MAX), S for South (= ETA_MIN), N for North (= ETA_MAX)
-  integer, parameter :: W = 1,E = 2,S = 3,N = 4,B = 5,T = 6        ! B==Bottom, T==Top
+  integer, parameter :: W = 1,E = 2,S = 3,N = 4,B = 5,T = 6        ! B == Bottom, T == Top
 
   integer :: p1,p2,p3,p4,p5,p6,p7,p8,p9,ia,iglob
   integer :: factor_x,factor_y,factor_z

src/meshfem3D/create_meshfem_mesh.f90

@@ -83,9 +83,13 @@ subroutine create_meshfem_mesh()
   ! HDF5 file i/o
   use shared_parameters, only: HDF5_ENABLED
 
-  !! setting up wavefield discontinuity interface
+  ! setting up wavefield discontinuity interface
   use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY
 
+  ! CUBIT output
+  use shared_parameters, only: COUPLE_WITH_INJECTION_TECHNIQUE
+  use meshfem_par, only: SAVE_MESH_AS_CUBIT
+
   implicit none
 
   ! local parameters
@@ -194,6 +198,28 @@ subroutine create_meshfem_mesh()
                         ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top)
   endif
 
+  ! CUBIT output
+  if (SAVE_MESH_AS_CUBIT) then
+    ! add outputs as CUBIT
+    call save_mesh_files_as_cubit(nspec,nglob, &
+                                  nodes_coords, ispec_material_id, &
+                                  nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+                                  ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top)
+
+    ! output for AxiSEM coupling
+    if (COUPLE_WITH_INJECTION_TECHNIQUE) then
+      call save_mesh_files_for_coupled_model(nspec, &
+                                             nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+                                             ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+                                             xstore,ystore,zstore)
+    endif
+  endif
+
+  ! setting up wavefield discontinuity interface
+  if (IS_WAVEFIELD_DISCONTINUITY) then
+    call write_wavefield_discontinuity_database()
+  endif
+
   ! re-assign back actual number of elements and nodes
   NSPEC_AB = nspec
   NGLOB_AB = nglob

src/meshfem3D/get_wavefield_discontinuity.f90

@@ -27,51 +27,64 @@
 
 
   subroutine write_wavefield_discontinuity_database()
+
 ! write data related to wavefield discontinuity to proc*_Database
+
   use meshfem_par, only: nb_wd, boundary_to_ispec_wd, side_wd, prname
   use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_MESH, &
                        IFILE_WAVEFIELD_DISCONTINUITY
   implicit none
-  ! local variables
+
   open(unit=IFILE_WAVEFIELD_DISCONTINUITY, &
        file = prname(1:len_trim(prname))//&
             trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), &
        form='unformatted', action='write')
+
   write(IFILE_WAVEFIELD_DISCONTINUITY) nb_wd
   write(IFILE_WAVEFIELD_DISCONTINUITY) boundary_to_ispec_wd
   write(IFILE_WAVEFIELD_DISCONTINUITY) side_wd
+
   close(IFILE_WAVEFIELD_DISCONTINUITY)
+
   end subroutine write_wavefield_discontinuity_database
 
 !
 !-----------------------------------------------------------
 !
 
   subroutine write_wavefield_discontinuity_file()
+
 ! write wavefield discontinuity interfaces to an ascii file
 ! works only when NPROC = 1
+
   use meshfem_par, only: nb_wd, boundary_to_ispec_wd, side_wd
   use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, &
                        FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE
   implicit none
   ! local variables
   integer :: i
+
   open(unit=IFILE_WAVEFIELD_DISCONTINUITY, &
        file='MESH/'//trim(FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE), &
        form='formatted', action='write')
+
   do i = 1, nb_wd
     write(IFILE_WAVEFIELD_DISCONTINUITY, *) boundary_to_ispec_wd(i), side_wd(i)
   enddo
+
   close(IFILE_WAVEFIELD_DISCONTINUITY)
+
   end subroutine write_wavefield_discontinuity_file
 
 !
 !-----------------------------------------------------------
 !
 
   subroutine find_wavefield_discontinuity_elements()
+
 ! read the wavefield_discontinuity_box file
 ! find the wavefield discontinuity interface
+
   use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, &
                        FNAME_WAVEFIELD_DISCONTINUITY_BOX
   use meshfem_par, only: xstore,ystore,zstore,nspec
@@ -90,6 +103,7 @@ subroutine find_wavefield_discontinuity_elements()
   open(unit=IFILE_WAVEFIELD_DISCONTINUITY, &
        file=trim(FNAME_WAVEFIELD_DISCONTINUITY_BOX), &
        form='formatted', action='read')
+
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) IS_TOP_WAVEFIELD_DISCONTINUITY
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) IS_EXTRAPOLATION_MODE
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) x_min
@@ -98,7 +112,9 @@ subroutine find_wavefield_discontinuity_elements()
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) y_max
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) z_min
   read(IFILE_WAVEFIELD_DISCONTINUITY, *) z_max
+
   close(IFILE_WAVEFIELD_DISCONTINUITY)
+
   nb_wd = 0
   do ispec = 1, nspec
     covered(:) = .false.
@@ -473,9 +489,11 @@ subroutine find_wavefield_discontinuity_elements()
       endif
     endif
   enddo
+
   allocate(boundary_to_ispec_wd(nb_wd), side_wd(nb_wd))
   boundary_to_ispec_wd(1:nb_wd) = boundary_to_ispec_wd_temp(1:nb_wd)
   side_wd(1:nb_wd) = side_wd_temp(1:nb_wd)
+
   end subroutine find_wavefield_discontinuity_elements
 
 !
@@ -490,7 +508,9 @@ logical function is_boundary_wd(x, y, z, x_min, x_max, y_min, y_max, &
   double precision :: x_min, x_max, y_min, y_max, z_min, z_max
   double precision :: x, y, z, x_mid, y_mid, z_mid, dx, dy, dz
   logical :: IS_TOP_WAVEFIELD_DISCONTINUITY, IS_EXTRAPOLATION_MODE
+
   is_boundary_wd = .false.
+
   if (IS_EXTRAPOLATION_MODE) then
     if (((x > x_min - dx) .and. (x < x_max + dx) .and. &
          (y > y_min - dy) .and. (y < y_max + dy) .and. &
@@ -510,4 +530,5 @@ logical function is_boundary_wd(x, y, z, x_min, x_max, y_min, y_max, &
          (z_mid < z_max))) &
       is_boundary_wd = .true.
   endif
+
   end function is_boundary_wd

src/meshfem3D/read_mesh_parameter_file.f90

@@ -323,7 +323,7 @@ subroutine read_mesh_parameter_file()
     endif
   endif
   call synchronize_all()
-  
+
   ! set time step and radial distribution of elements
   ! right distribution is determined based upon maximum value of NEX
   NEX_MAX = max(NEX_XI,NEX_ETA)

src/meshfem3D/rules.mk

@@ -64,6 +64,7 @@ meshfem3D_OBJECTS = \
 	$O/read_value_mesh_parameters.mesh.o \
 	$O/save_databases.mesh.o \
 	$O/save_databases_hdf5.mesh_hdf5.o \
+	$O/save_mesh_as_cubit.mesh.o \
 	$O/store_boundaries.mesh.o \
 	$O/store_coords.mesh.o \
 	$(EMPTY_MACRO)