search

SPECFEM / specfem3d

SPECFEM3D_Cartesian simulates acoustic (fluid), elastic (solid), coupled acoustic/elastic, poroelastic or seismic wave propagation in any type of conforming mesh of hexahedra (structured & unstructured).
https://specfem.org
GNU General Public License v3.0
411 stars 229 forks source link

STOP Inconsistent fault mesh: corresponding node in the other fault face was not found #1614

Open TangSEU opened 1 year ago

TangSEU commented 1 year ago

Hi all： When I tried to run a model that includes faults, topograpy data, and different material layers in specmem3d, I received this message:STOP Inconsistent fault mesh: corresponding node in the other fault face was not found.

jpampuero commented 1 year ago

This sounds like a user error in the mesh generation. Please pay special attention to the instructions about meshing faults with split nodes on chapter 6 of the manual (https://specfem3d.readthedocs.io/en/latest/06_fault_sources/). If you need further help from developers, please add sufficient information about your problem so that someone can reproduce it.

TangSEU commented 1 year ago

This sounds like a user error in the mesh generation. Please pay special attention to the instructions about meshing faults with split nodes on chapter 6 of the manual (https://specfem3d.readthedocs.io/en/latest/06_fault_sources/). If you need further help from developers, please add sufficient information about your problem so that someone can reproduce it.

Thanks for your reply,the picture blow is my model created by cubit 2022.11.I have separated the two faults surface by a small distance. 07b1bc12da3a4f4108f825114ea0927

TangSEU commented 1 year ago 
#############################################################
#
# script uses ACIS surface formats
#
# note: this script seems to work with CUBIT version > 12.2
#
# created by Tang.S
#
# date: 2023/05
#############################################################
from __future__ import print_function

import os
import math
from math import sin
from math import pi
from math import tan
from math import cos
import sys
import os.path
import time

sys.path.append('/home/tang/specfem3d/specfem3d/CUBIT_GEOCUBIT/geocubitlib')
sys.path.append('/opt/Coreform-Cubit-2022.11/bin')
sys.path.append('/home/tang/specfem3d/specfem3d/CUBIT_GEOCUBIT')

import cubit
cubit.init([""])
import geocubitlib
from geocubitlib import boundary_definition
from geocubitlib import cubit2specfem3d
from geocubitlib import save_fault_nodes_elements
from geocubitlib import absorbing_boundary

# time stamp
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))

# working directory
cwd = os.getcwd()
print("#current working directory: " + str(cwd))

cubit.cmd('version')
cubit.cmd('reset')

# uses developer commands
cubit.cmd('set developer commands on')
cubit.cmd('set import mesh tolerance 1')

#############################################################
#
# 0. step: loading topography surface
#
#############################################################
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#loading topo surface...")

# topography surface
if os.path.exists("topo.cub"):
  print("opening existing topography surface")
  # topography surface
  # previously run, just reopen the cubit file
  cubit.cmd('open "topo.cub"')

# healing the surface...
cubit.cmd('Auto_clean volume 1 small_surfaces small_curve_size 10')
cubit.cmd('regularize volume 1')

#############################################################
#
# 1. step: creates temporary brick volume and fault
#
#############################################################
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#creating brick...")

# initializing coordinates x,y,z
km = 1000
z_surf = 40*km

x = []   # fault
y = []
z = []

xbulk = []   #bulk
ybulk = []
zbulk = []

xbulk.append(-30*km)   #x1
xbulk.append(30*km)    #x2
xbulk.append(30*km)    #x3
xbulk.append(-30*km)   #x4

ybulk.append(-30*km)  #y1
ybulk.append(-30*km)  #y2
ybulk.append(30*km)   #y3
ybulk.append(30*km)   #y4

zbulk = [z_surf]*4

a1 = -27*km
a2 = 0*km
a3 = 27*km
a4 = 0*km

b1 = 0
b2 = 0.1
b3 = 0
b4 = -0.1

x.append(a1)   #x5
x.append(a2)   #x6
x.append(a3)  #x7
x.append(a4)   #x8

y.append(b1)       #y5
y.append(b2)    #y6
y.append(b3)       #y7
y.append(b4)   #y8

z_surf2 = 44*km
z = [z_surf2]*4

# bulk
for i in range(len(xbulk)):
  vert = "create vertex x "+str(xbulk[i])+" y "+str(ybulk[i])+" z "+str(zbulk[i])
  cubit.cmd(vert)

# loading fault points profile
for i in range(len(x)):
  vert="create vertex x "+str(x[i])+" y "+str(y[i])+" z "+str(z[i])
  cubit.cmd(vert)

cubit.cmd("create vertex x -27000 y 0 z -24000")
cubit.cmd("create vertex x 27000 y 0 z -24000")

# create fault domains
bulk1="create curve vertex 40406 40407"
bulk2="create curve vertex 40407 40408"
bulk3="create curve vertex 40408 40409"
bulk4="create curve vertex 40409 40406"

fault_up="create curve spline vertex 40410 40411 40412"
fault_down="create curve spline vertex 40410 40413 40412"

cubit.cmd(bulk1)
cubit.cmd(bulk2)
cubit.cmd(bulk3)
cubit.cmd(bulk4)

cubit.cmd(fault_up)
cubit.cmd(fault_down)

cubit.cmd("create curve vertex 40410 40414")
cubit.cmd("create curve vertex 40412 40415")
cubit.cmd("create curve vertex 40414 40415")

cubit.cmd("create surface 210 212 213 214")
cubit.cmd("create surface 211 212 213 214")

# rotate fault plane
cubit.cmd("rotate surface 2 about x 0 angle 6")
cubit.cmd("rotate surface 3 about x 0 angle 6")

cubit.cmd("rotate surface 2 about z 0 angle 341")
cubit.cmd("rotate surface 3 about z 0 angle 341")

# cut fault plane
cubit.cmd('merge all')
cubit.cmd('webcut surface 2 with plane zplane offset 40000')
cubit.cmd('webcut surface 3 with plane zplane offset 40000')

cubit.cmd('delete surface 4')
cubit.cmd('delete surface 6')

cubit.cmd("sweep curve 206 vector 0 0 -1 distance "+str(80*km))
cubit.cmd("sweep curve 207 vector 0 0 -1 distance "+str(80*km))
cubit.cmd("sweep curve 208 vector 0 0 -1 distance "+str(80*km))
cubit.cmd("sweep curve 209 vector 0 0 -1 distance "+str(80*km))

cubit.cmd("create surface curve 221 227 206 207 208 209")
cubit.cmd("create surface curve 230 233 236 239")
cubit.cmd("create volume surface 5 7 8 9 10 11 12 13")

# moves volume to UTM coordinates of topography surface
cubit.cmd('volume 11 move x 701538 y 4189202.3 z 0 ')

# temporary backup
cubit.cmd('save as "topo_1.cub" overwrite')

#############################################################
#
# 2. step: creates volume with topography surface
#
#############################################################
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#creating volume with topography...")

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

# exports only surfaces which will create single volume
cubit.cmd('export acis "topo_2.acis" surface 13 20 21 25 28 17 23 14 ascii overwrite')

# backup
cubit.cmd('save as "topo_2.cub" overwrite')

#############################################################
#
# 3. step: manipulate ACIS file to create a single volume
#
#############################################################
# checks if new file available
if not os.path.exists("topo_2.acis"):
  print("")
  print("error creating new volume, please check manually...")
  print("")
  cubit.cmd('pause')
# clears workspace
cubit.cmd('reset')

# imports surfaces and merges to single volume
cubit.cmd('import acis "topo_2.acis" ascii merge_globally')

print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#creating new volume, this will take another 2 min...")
cubit.cmd('create volume surface all heal')

# backup
cubit.cmd('save as "topo_3.cub" overwrite')

#############################################################
#
# 4. step: create mesh
#
#############################################################
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#initial meshing...")
print("#(will take around 100 min)")

# cubit.cmd('webcut volume 1 with plane zplane offset -2000')
# cubit.cmd('webcut volume 2 with plane zplane offset -7000')
# cubit.cmd('webcut volume 3 with plane zplane offset -25000')

# create top volume
cubit.cmd('split curve 45  distance 38000 from vertex 37')
cubit.cmd('split curve 43  distance 38000 from vertex 36')
cubit.cmd('split curve 40  distance 38000 from vertex 35')
cubit.cmd('split curve 41  distance 38000 from vertex 38')

cubit.cmd('create surface vertex 47 48 49 50')
cubit.cmd('section volume 1 with surface 17 keep normal')
cubit.cmd('section volume 1 with surface 17 reverse')

# create mid volume 1
cubit.cmd('split curve 70  distance 5000 from vertex 59')
cubit.cmd('split curve 72  distance 5000 from vertex 60')
cubit.cmd('split curve 80  distance 5000 from vertex 58')
cubit.cmd('split curve 73  distance 5000 from vertex 57')

cubit.cmd('create surface vertex 69 70 71 72')

cubit.cmd('section volume 3 with surface 33 keep normal')
cubit.cmd('section volume 3 with surface 33 reverse')

# create mid volume 2
cubit.cmd('split curve 107  distance 18000 from vertex 81')
cubit.cmd('split curve 109  distance 18000 from vertex 82')
cubit.cmd('split curve 115  distance 18000 from vertex 80')
cubit.cmd('split curve 110  distance 18000 from vertex 79')

cubit.cmd('create surface vertex 91 92 93 94')

cubit.cmd('section volume 5 with surface 49 keep normal')
cubit.cmd('section volume 5 with surface 49 reverse')

cubit.cmd('delete volume 2 4 6')

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

cubit.cmd("surface 18 size 600")
cubit.cmd("surface 18 scheme pave")
cubit.cmd("mesh surface 18")
cubit.cmd("volume 3 size 600")
cubit.cmd("mesh volume 3")
cubit.cmd('refine surface 18 numsplit 1 bias 1.0 depth 1')

cubit.cmd("surface 34 size 600")
cubit.cmd("volume 5 size 600")
cubit.cmd("surface 34 scheme pave")
cubit.cmd("mesh surface 34")
cubit.cmd("mesh volume 5")

cubit.cmd("surface 50 size 600")
cubit.cmd("volume 7 size 600")
cubit.cmd("surface 50 scheme pave")
cubit.cmd("mesh surface 50")
cubit.cmd("mesh volume 7")

cubit.cmd("volume 1 size 200")
cubit.cmd("mesh volume 1")

# fault output
cubit.cmd("unmerge surface 29 32")
cubit.cmd("unmerge surface 48 45")
cubit.cmd("unmerge surface 37 40")

os.system('mkdir -p MESH')

# time stamp
print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime()))
print("#done meshing...")

# backup
cubit.cmd('save as "topo_4.cub" overwrite')

########### Fault elements and nodes ###############
Au = [29, 48, 37]
Ad = [32, 45, 40]
faultA = save_fault_nodes_elements.fault_input(1,Au,Ad)

########### boundary definition #####################
entities=['face']
absorbing_boundary.define_parallel_bc(entities) # in absorbing_boundary.py

#### Define material properties for the 2 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')

cubit.cmd('block 1 name "elastic 1" ')        # material region
cubit.cmd('block 1 attribute count 6')
cubit.cmd('block 1 attribute index 1 1')      # flag for fault side 1
cubit.cmd('block 1 attribute index 2 4800')   # vp
cubit.cmd('block 1 attribute index 3 2800')    # vs
cubit.cmd('block 1 attribute index 4 3186')   # rho
cubit.cmd('block 1 attribute index 5 200')     # q flag (see constants.h: iattenuation_ ... )
cubit.cmd('block 1 attribute index 6 0')     # q flag (see constants.h: iattenuation_ ... )

cubit.cmd('block 2 name "elastic 2" ')        # material region
cubit.cmd('block 2 attribute count 6')
cubit.cmd('block 2 attribute index 1 2')      # flag for fault side 1
cubit.cmd('block 2 attribute index 2 5350')   # vp
cubit.cmd('block 2 attribute index 3 3200')    # vs
cubit.cmd('block 2 attribute index 4 3350')   # rho
cubit.cmd('block 2 attribute index 5 500')     # q flag (see constants.h: iattenuation_ ... )
cubit.cmd('block 2 attribute index 6 0')     # q flag (see constants.h: iattenuation_ ... )

cubit.cmd('block 3 name "elastic 3" ')        # material region
cubit.cmd('block 3 attribute count 6')
cubit.cmd('block 3 attribute index 1 3')      # flag for fault side 1
cubit.cmd('block 3 attribute index 2 6000')   # vp
cubit.cmd('block 3 attribute index 3 3500')    # vs
cubit.cmd('block 3 attribute index 4 3544')   # rho
cubit.cmd('block 3 attribute index 5 9000')     # q flag (see constants.h: iattenuation_ ... )
cubit.cmd('block 3 attribute index 6 0')     # q flag (see constants.h: iattenuation_ ... )

cubit.cmd('block 4 name "elastic 4" ')        # material region
cubit.cmd('block 4 attribute count 6')
cubit.cmd('block 4 attribute index 1 4')      # flag for fault side 1
cubit.cmd('block 4 attribute index 2 6150')   # vp
cubit.cmd('block 4 attribute index 3 3600')    # vs
cubit.cmd('block 4 attribute index 4 3588')   # rho
cubit.cmd('block 4 attribute index 5 9000')     # q flag (see constants.h: iattenuation_ ... )
cubit.cmd('block 4 attribute index 6 0')     # q flag (see constants.h: iattenuation_ ... )
#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT

cubit2specfem3d.export2SPECFEM3D('MESH')
TangSEU commented 1 year ago

This sounds like a user error in the mesh generation. Please pay special attention to the instructions about meshing faults with split nodes on chapter 6 of the manual (https://specfem3d.readthedocs.io/en/latest/06_fault_sources/). If you need further help from developers, please add sufficient information about your problem so that someone can reproduce it.

It would be appreciated that if you can also take a look at my Python file above and point out my mistake.Thank you very much.

jxhxb commented 2 months ago

As of now, have you done it successfully? I have also faced the same problem, if you have fixed the problem ,can you told me how to do it ?Thank you very much. 01

jxhxb commented 2 months ago

fault.txt

jxhxb commented 2 months ago

And have you succeed setting the fault's strike angle and dip angle ? I want to know how it work in your example,thank you a lot!