models may not match

[liaohaiyang1534]

Hi,

i use a python scipt to make a model (given at the end), and it looks like this in the software CoreForm Cubit 2022.11:

when I set

# save AVS or OpenDX mesh files to check the mesh
SAVE_MESH_FILES                 = .true.

and check part_array.vtk from OUTPUT_FILES/DATABASES_MPI, it looks like this:

and I think this is alright, it matches the size of my model.

but when I check proc000010_vs.vtk:

It looks a little strange.

Should I ignore it? Or is there some potential bug?

the python script that I use:

> `
> from __future__ import print_function
> 
> import os
> import sys
> 
> elementsize = 4
> 
> # default directories
> SEMoutput='MESH'
> CUBIToutput='MESH_GEOCUBIT'
> 
> os.system('mkdir -p '+ SEMoutput)
> os.system('mkdir -p '+ CUBIToutput)
> 
> import cubit
> try:
>     #cubit.init([""])
>     cubit.init(["-noecho","-nojournal"])
> except:
>     pass
> 
> version = cubit.get_version()
> version_major = int(version.split(".")[0])
> version_minor = int(version.split(".")[1])
> print("cubit version: ",version)
> 
> cubit.cmd('reset')
> cubit.cmd('brick x 300 y 150 z 100')
> 
> # This seems to conflict with boundary_definition.py
> 
> cubit.cmd('volume 1 move x 150 y 75 z -50')
> 
> # create vertices for discontinuity
> cubit.cmd('split curve 9  distance 20')
> cubit.cmd('split curve 10  distance 20')
> cubit.cmd('split curve 11  distance 20')
> cubit.cmd('split curve 12  distance 20')
> 
> # create surface for interface
> cubit.cmd('create surface vertex 9 10 12 11')
> 
> cubit.cmd('section volume 1 with surface 7 keep normal')
> cubit.cmd('section volume 1 with surface 7 reverse')
> 
> # create vertices for auxiliary interface to allow for refinement
> cubit.cmd('split curve 29  distance 30')
> cubit.cmd('split curve 31  distance 30')
> cubit.cmd('split curve 32  distance 30')
> cubit.cmd('split curve 36  distance 30')
> 
> # create surface for buffer interface to refine BELOW the discontinuity
> cubit.cmd('create surface vertex 25 26 28 27')
> 
> cubit.cmd('section volume 3 with surface 19 keep normal')
> cubit.cmd('section volume 3 with surface 19 reverse')
> 
> cubit.cmd('delete volume 2 4')
> 
> cubit.cmd('merge all')
> cubit.cmd('imprint all')
> 
> # Meshing the volumes
> cubit.cmd('volume 3 size '+str(elementsize))
> cubit.cmd('mesh volume 3')
> 
> cubit.cmd('volume 1 size '+str(elementsize))
> cubit.cmd('mesh volume 1')
> 
> cubit.cmd('volume 5 size '+str(elementsize))
> cubit.cmd('mesh volume 5')
> 
> #### End of meshing
> 
> ## obsolete:
> #import boundary_definition
> #import cubit2specfem3d
> ## new:
> from geocubitlib import boundary_definition
> from geocubitlib import cubit2specfem3d
> 
> ###### This is boundary_definition.py of GEOCUBIT
> #..... which extracts the bounding faces and defines them into blocks
> boundary_definition.entities=['face']
> boundary_definition.define_bc(boundary_definition.entities,parallel=True)
> 
> #### Define material properties for the 3 volumes ################
> cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
> cubit.cmd('block 1 name "elastic 1" ')        # elastic material region
> cubit.cmd('block 1 attribute count 7')
> cubit.cmd('block 1 attribute index 1 1  ')     # material 1
> cubit.cmd('block 1 attribute index 2 1200 ')  # vp
> cubit.cmd('block 1 attribute index 3 400 ')      # vs
> cubit.cmd('block 1 attribute index 4 2000 ')  # rho 
> cubit.cmd('block 2 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 2 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 2 attribute index 7 0 ')     # anisotropy_flag
> 
> cubit.cmd('block 2 name "elastic 2" ')        # elastic material region
> cubit.cmd('block 2 attribute count 7')
> cubit.cmd('block 2 attribute index 1 2  ')     # material 2
> cubit.cmd('block 2 attribute index 2 2000 ')  # vp
> cubit.cmd('block 2 attribute index 3 700 ')  # vs
> cubit.cmd('block 2 attribute index 4 3200 ')  # rho
> cubit.cmd('block 2 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 2 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 2 attribute index 7 0 ')     # anisotropy_flag
> 
> cubit.cmd('block 3 name "elastic 3" ')        # elastic material region
> cubit.cmd('block 3 attribute count 7')
> cubit.cmd('block 3 attribute index 1 3  ')     
> cubit.cmd('block 3 attribute index 2 2500 ')  # vp
> cubit.cmd('block 3 attribute index 3 800 ')  # vs
> cubit.cmd('block 3 attribute index 4 3600 ')  # rho
> cubit.cmd('block 3 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 3 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 3 attribute index 7 0 ')     # anisotropy_flag
> 
> cubit.cmd('export mesh "' + CUBIToutput + '/top.e" dimension 3 overwrite')
> cubit.cmd('save as "' + CUBIToutput + '/top.cub" overwrite')
> 
> #### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT
> 
> cubit2specfem3d.export2SPECFEM3D(SEMoutput)
> 
> # screen shot
> # (could crash version < 16.4)
> if version_major >= 16 and version_minor >= 4:
>     cubit.cmd('view top')
>     cubit.cmd('hardcopy "' + CUBIToutput + '/waterlayered.png" png')
> 
> # all files needed by SCOTCH are now in directory MESH`
> 

I'm very appreciated to have any advice and suggestions for this.
(If you need any more clarifications, please let me know)

Regards

Haiyang Liao
haiyangliao@smail.nju.edu.cn

[homnath]

Hi Haiyang Liao, How does it look like when you visualize all processor files, proc000000_vs.vtk, proc000001_vs.vtk, proc000002_vs.vtk, etc., altogether? The individual mesh may look irregular due to partitioning but altogether, they should match your original model. Best, Hom Nath

________________________________ From: liaohaiyang1534 ***@***.***> Sent: Saturday, May 27, 2023 11:15 PM To: SPECFEM/specfem3d ***@***.***> Cc: Subscribed ***@***.***> Subject: [SPECFEM/specfem3d] models may not match (Issue #1609) Hi, i use a python scipt to make a model (given at the end), and it looks like this in the software CoreForm Cubit 2022.11: [)EP3~DHA8G8%TPP{AN03@{6]<https://user-images.githubusercontent.com/104192749/241474009-81f180e1-06e3-4825-8374-0acb5c092df9.png> when I set # save AVS or OpenDX mesh files to check the mesh SAVE_MESH_FILES = .true. and check part_array.vtk from OUTPUT_FILES/DATABASES_MPI, it looks like this: [image]<https://user-images.githubusercontent.com/104192749/241476203-a3c7f16f-a220-43b3-8e37-65cd37947672.png> and I think this is alright, it matches the size of my model. but when I check proc000010_vs.vtk: [image]<https://user-images.githubusercontent.com/104192749/241476548-8c2f60ac-5cca-47b2-9b18-f12a0b4afcfc.png> It looks a little strange. Should I ignore it? Or is there some potential bug? the python script that I use:

` from __future__ import print_function import os import sys elementsize = 4 # default directories SEMoutput='MESH' CUBIToutput='MESH_GEOCUBIT' os.system('mkdir -p '+ SEMoutput) os.system('mkdir -p '+ CUBIToutput) import cubit try: #cubit.init([""]) cubit.init(["-noecho","-nojournal"]) except: pass version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) print("cubit version: ",version) cubit.cmd('reset') cubit.cmd('brick x 300 y 150 z 100') # This seems to conflict with boundary_definition.py cubit.cmd('volume 1 move x 150 y 75 z -50') # create vertices for discontinuity cubit.cmd('split curve 9 distance 20') cubit.cmd('split curve 10 distance 20') cubit.cmd('split curve 11 distance 20') cubit.cmd('split curve 12 distance 20') # create surface for interface cubit.cmd('create surface vertex 9 10 12 11') cubit.cmd('section volume 1 with surface 7 keep normal') cubit.cmd('section volume 1 with surface 7 reverse') # create vertices for auxiliary interface to allow for refinement cubit.cmd('split curve 29 distance 30') cubit.cmd('split curve 31 distance 30') cubit.cmd('split curve 32 distance 30') cubit.cmd('split curve 36 distance 30') # create surface for buffer interface to refine BELOW the discontinuity cubit.cmd('create surface vertex 25 26 28 27') cubit.cmd('section volume 3 with surface 19 keep normal') cubit.cmd('section volume 3 with surface 19 reverse') cubit.cmd('delete volume 2 4') cubit.cmd('merge all') cubit.cmd('imprint all') # Meshing the volumes cubit.cmd('volume 3 size '+str(elementsize)) cubit.cmd('mesh volume 3') cubit.cmd('volume 1 size '+str(elementsize)) cubit.cmd('mesh volume 1') cubit.cmd('volume 5 size '+str(elementsize)) cubit.cmd('mesh volume 5') #### End of meshing ## obsolete: #import boundary_definition #import cubit2specfem3d ## new: from geocubitlib import boundary_definition from geocubitlib import cubit2specfem3d ###### This is boundary_definition.py of GEOCUBIT #..... which extracts the bounding faces and defines them into blocks boundary_definition.entities=['face'] boundary_definition.define_bc(boundary_definition.entities,parallel=True) #### Define material properties for the 3 volumes ################ cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################') cubit.cmd('block 1 name "elastic 1" ') # elastic material region cubit.cmd('block 1 attribute count 7') cubit.cmd('block 1 attribute index 1 1 ') # material 1 cubit.cmd('block 1 attribute index 2 1200 ') # vp cubit.cmd('block 1 attribute index 3 400 ') # vs cubit.cmd('block 1 attribute index 4 2000 ') # rho cubit.cmd('block 2 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 2 attribute index 6 9999.0') # Q_mu cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag cubit.cmd('block 2 name "elastic 2" ') # elastic material region cubit.cmd('block 2 attribute count 7') cubit.cmd('block 2 attribute index 1 2 ') # material 2 cubit.cmd('block 2 attribute index 2 2000 ') # vp cubit.cmd('block 2 attribute index 3 700 ') # vs cubit.cmd('block 2 attribute index 4 3200 ') # rho cubit.cmd('block 2 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 2 attribute index 6 9999.0') # Q_mu cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag cubit.cmd('block 3 name "elastic 3" ') # elastic material region cubit.cmd('block 3 attribute count 7') cubit.cmd('block 3 attribute index 1 3 ') cubit.cmd('block 3 attribute index 2 2500 ') # vp cubit.cmd('block 3 attribute index 3 800 ') # vs cubit.cmd('block 3 attribute index 4 3600 ') # rho cubit.cmd('block 3 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 3 attribute index 6 9999.0') # Q_mu cubit.cmd('block 3 attribute index 7 0 ') # anisotropy_flag cubit.cmd('export mesh "' + CUBIToutput + '/top.e" dimension 3 overwrite') cubit.cmd('save as "' + CUBIToutput + '/top.cub" overwrite') #### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT cubit2specfem3d.export2SPECFEM3D(SEMoutput) # screen shot # (could crash version < 16.4) if version_major >= 16 and version_minor >= 4: cubit.cmd('view top') cubit.cmd('hardcopy "' + CUBIToutput + '/waterlayered.png" png') # all files needed by SCOTCH are now in directory MESH`

I'm very appreciated to have any advice and suggestions for this. (If you need any more clarifications, please let me know) Regards Haiyang Liao ***@***.******@***.***> — Reply to this email directly, view it on GitHub<#1609>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABMCQ4XI3Q3GTV3BTOBUE5DXIK7MZANCNFSM6AAAAAAYRR4LDE>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>