Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Question about creating external models with cubit, can't generate absorption boundary files #1710

Closed
zuoxuan-de opened this issue Jul 4, 2024 · 2 comments
Closed

Question about creating external models with cubit, can't generate absorption boundary files #1710

zuoxuan-de opened this issue Jul 4, 2024 · 2 comments

Comments

@zuoxuan-de
Copy link

zuoxuan-de commented Jul 4, 2024 •
edited by danielpeter
Loading

Hello guys, I'm creating a three-layer model containing an alleyway, where the alleyway is located in the middle layer, in addition, I manually set his free surface and other absorbing surfaces. But when I run cubit2SPECFEM3D, I have a problem with my absorbing_surface_file_bottom,free_or_absorbing_surface_file_zmax,absorbing_surface_file_ymin and other similar absorbing surface files are all empty. I saw your replies related to this in my previous question, I hope you can be generous with your suggestions, thanks a lot! @casarotti @homnath
image
This is the specific case of the lane in the second layer, which is empty in the center, without any medium, and has all free surfaces.
image
This is the free surface of the lane I set manually.
image
In order to set the absorbing layer on the upper surface, I set both the upper and lower interfaces (zmax,zmin) to face_bottom again.

Next is my specific python script：

#!/usr/bin/env python

import cubit
import cubit2specfem3d

import os
import sys


SEMoutput='MESH03'
CUBIToutput='MESH_GEOCUBIT03'

os.system('mkdir -p '+ SEMoutput)
os.system('mkdir -p '+ CUBIToutput)

# Simulating a coal seam with a roadway
# This geological structure is mainly divided into three layers, the upper and lower peripheral rocks and the coal seam,
# which contains two roadways, of which the roadway portion is not defined by any medium and is all free surface
cubit.cmd('reset')
# underlying rock
cubit.cmd('brick x 50 y 40 z 10')
cubit.cmd('vol 1 move x 25 y 20 z 5')

#################################################  intermediate seam   ##############################################

#Leftmost part of the coal seam in the y-direction, also the left wall of the left roadway
cubit.cmd('brick x 50 y 15 z 5')
cubit.cmd('vol 2 move x 25 y 7.5 z -2.5')

# #Rightmost part of the coal seam in the y-direction, also the right wall of the right roadway
cubit.cmd('brick x 50 y 15 z 5')
cubit.cmd('volume 3 move x 25 y 32.5 z -2.5')

# The top of the lane in the y-direction
cubit.cmd('brick x 5 y 10 z 5')
cubit.cmd('volume 4 move x 2.5 y 20 z -2.5')

# The end of the lane in the y-direction
cubit.cmd('brick x 5 y 10 z 5')
cubit.cmd('volume 5 move x 47.5 y 20 z -2.5')


####################################################  end ####################################################
# top layer of surrounding rock
cubit.cmd('brick x 50 y 40 z 10')
cubit.cmd('volume 6 move x 25 y 20 z -10')

cubit.cmd('merge all')
cubit.cmd('imprint all')

# Meshing the volumes
elementsize = 1.0

cubit.cmd('volume all size 1')
cubit.cmd('mesh volume all')

cubit.cmd('block 1001 add surface 24 28 38 42 46 49 ')
cubit.cmd('block 1001 name "face_topo"')

cubit.cmd('block 1002 add surface 1 32 ')
cubit.cmd('block 1002 name "face_bottom"')

cubit.cmd('block 1003 add surface 4 10 16 22 34 ')
cubit.cmd('block 1003 name "face_abs_xmin"')

cubit.cmd('block 1004 add surface 3 9 33')
cubit.cmd('block 1004 name "face_abs_ymin"')

cubit.cmd('block 1005 add surface 6 12 18 36  30 ')
cubit.cmd('block 1005 name "face_abs_xmax"')

cubit.cmd('block 1006 add surface 5 17 35')
cubit.cmd('block 1006 name "face_abs_ymax"')


from geocubitlib import cubit2specfem3d


#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
cubit.cmd('block 1 hex in volume 1')
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')      # flag for material: 1 for 1. material
cubit.cmd('block 1 attribute index 2 3500')   # vp
cubit.cmd('block 1 attribute index 3 2000')   # vs
cubit.cmd('block 1 attribute index 4 2400')   # rho
cubit.cmd('block 1 attribute index 5 9999.0') # Q_kappa
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 2')
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')      # flag for material: 1 for 1. material
cubit.cmd('block 2 attribute index 2 1900')   # vp
cubit.cmd('block 2 attribute index 3 1100')   # vs
cubit.cmd('block 2 attribute index 4 1300')   # rho
cubit.cmd('block 2 attribute index 5 9999.0')  # Q_mu
cubit.cmd('block 2 attribute index 6 9999.0')  # Qmu
cubit.cmd('block 2 attribute index 7 0 ')      # anisotropy_flag

cubit.cmd('block 3 hex in volume 3')
cubit.cmd('block 3 name "elastic 3"')
cubit.cmd('block 3 attribute count 7')
cubit.cmd('block 3 attribute index 1 3 ')
cubit.cmd('block 3 attribute index 2  1900 ')
cubit.cmd('block 3 attribute index 3  1100 ')
cubit.cmd('block 3 attribute index 4  1300 ')
cubit.cmd('block 3 attribute index 5 9999')
cubit.cmd('block 3 attribute index 6 9999')
cubit.cmd('block 3 attribute index 7 0')

cubit.cmd('block 4 hex in volume 4')
cubit.cmd('block 4 name "elastic 4"')
cubit.cmd('block 4 attribute count 7')
cubit.cmd('block 4 attribute index 1 4 ')
cubit.cmd('block 4 attribute index 2 1900 ')
cubit.cmd('block 4 attribute index 3 1100 ')
cubit.cmd('block 4 attribute index 4 1300 ')
cubit.cmd('block 4 attribute index 5 9999')
cubit.cmd('block 4 attribute index 6 9999')
cubit.cmd('block 4 attribute index 7 0')

cubit.cmd('block 5 hex in volume 5')
cubit.cmd('block 5 name "elastic 5"')
cubit.cmd('block 5 attribute count 7')
cubit.cmd('block 5 attribute index 1 5 ')
cubit.cmd('block 5 attribute index 2 1900 ')
cubit.cmd('block 5 attribute index 3 1100 ')
cubit.cmd('block 5 attribute index 4 1300 ')
cubit.cmd('block 5 attribute index 5 9999')
cubit.cmd('block 5 attribute index 6 9999')
cubit.cmd('block 5 attribute index 7 0')


cubit.cmd('block 6 hex in volume 6')
cubit.cmd('block 6 name "elastic 6"')
cubit.cmd('block 6 attribute count 7')
cubit.cmd('block 6 attribute index 1 6 ')
cubit.cmd('block 6 attribute index 2 3500 ')
cubit.cmd('block 6 attribute index 3 2000 ')
cubit.cmd('block 6 attribute index 4 2400 ')
cubit.cmd('block 6 attribute index 5 9999')
cubit.cmd('block 6 attribute index 6 9999')
cubit.cmd('block 6 attribute index 7 0')


cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')

#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT

os.system('mkdir -p MESH0703')
cubit2specfem3d.export2SPECFEM3D('MESH0703')

# all files needed by SCOTCH are now in directory MESH

One last question, can cubit determine the dimensions of the three dimensions of a hexahedral cell, similar to 110.5, instead of a square with side lengths of 1.
Thank you again and I hope to hear back from you all sooner rather than later!
@homnath
Copy link

homnath commented Jul 5, 2024 via email

@danielpeter
Copy link
Contributor

danielpeter commented Jul 8, 2024 •
edited
Loading

there are a few problems with your meshing script:

-- use an order of "imprint & merge":
that is, use first imprint, then merge command. this will connect your volumes correctly, otherwise your volumes and mesh won't be conforming, but still be separated. thus,
replace:

cubit.cmd('merge all')
cubit.cmd('imprint all')

with

cubit.cmd('imprint all')
cubit.cmd('merge all')

-- the blocks defined for the surfaces need to have faces, not surfaces added:
thus, replace any occurrences of add surface with face in surface, like:

cubit.cmd('block 1002 add surface 1 32 ')

with

cubit.cmd('block 1002 face in surface 1 32 ')

this will then find all faces for those blocks.

-- replace your block name face_topo to something like face_free and use a different block number than 1001:
the face_topo is used for a free surface that is also a topography surface at the top of a mesh. when outputting such topography surfaces, it will include a vertical normal check to make sure that the surface has a correct face point numbering.

since your "free" surface has vertical sides, you should avoid this check and define it as a more general free surface. thus, use a block name like face_free or free_surface (anything with a free added will do) and use a block number other than 1001 which is reserved for topography surfaces again. thus,
replace:

cubit.cmd('block 1001 add surface 24 28 38 42 46 49 ')
cubit.cmd('block 1001 name "face_topo"')

with

cubit.cmd('block 1000 face in surface 24 28 38 42 46 49 ')
cubit.cmd('block 1000 name "face_free"')

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
3 participants
@homnath @danielpeter @zuoxuan-de