I want to create a model with skew faults, but I get an error after I modify the example in SPECFEM3D: STOP INCONSISTENT FAULT MESH:Corresponding node in the other fault face was not found.
Perhaps you could take a look at my Python files and give me some help. I’d be grateful！
#!/usr/bin/env python
from __future__ import print_function
import cubit
import boundary_definition
import cubit2specfem3d
import math
import os
import sys
import numarray
from save_fault_nodes_elements import *
cubit.cmd('reset')
km = 1000
y_vert=0*km
h = 0.1#(h=0.1)
### initializing coordinates x,y,z ###
x=[] # fault
y=[]
z=[]
xbulk=[] # bulk
ybulk=[]
zbulk=[]
### Bulk points ###
xbulk.append(-80*km) #x1
xbulk.append(80*km) #x2
xbulk.append(80*km) #x3
xbulk.append(-80*km) #x4
zbulk.append(-40*km) #y1
zbulk.append(-40*km) #y2
zbulk.append(40*km) #y3
zbulk.append(40*km) #y4
ybulk=[y_vert]*4
### CRACKS ##########
x.append(80*km) #x5
x.append(80*km) #x6
x.append(-80*km) #x7
x.append(-80*km) #x8
z.append(-20*km-h) #y5
z.append(-20*km+h) #y6
z.append(10*km+h) #y7
z.append(10*km-h) #y8
y=[y_vert]*9
#################### 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)
################ creating fault domains #################################
bulk1="create curve vertex 1 2" #c1
bulk2="create curve vertex 2 5" #c2
bulk3="create curve vertex 6 3" #c3
bulk4="create curve vertex 3 4" #c4
bulk5="create curve vertex 4 7" #c5
bulk6="create curve vertex 8 1" #c6
cubit.cmd(bulk1)
cubit.cmd(bulk2)
cubit.cmd(bulk3)
cubit.cmd(bulk4)
cubit.cmd(bulk5)
cubit.cmd(bulk6)
fault_up="create curve spline vertex 6 7" #c7
fault_down="create curve spline vertex 5 8" #c8
cubit.cmd(fault_up)
cubit.cmd(fault_down)
surface1="create surface curve 7 8"
cubit.cmd(surface1)
surface2="create surface curve 1 2 3 4 5 6"
cubit.cmd(surface2)
cubit.cmd("sweep surface 1 vector 0 1 0 distance "+str(160*km))
cubit.cmd("sweep surface 2 vector 0 1 0 distance "+str(160*km))
#####################################################
#####################################################
elementsize = 2000 #(2500)
# IMPRINTING
cubit.cmd("imprint all")
# MERGING
cubit.cmd("merge all")
# Meshing coincide fault_A upper boundaries .
cubit.cmd('curve 7 8 size 2000')
cubit.cmd('curve 7 8 scheme equal')
cubit.cmd('mesh curve 7 8')
cubit.cmd("surface 12 size "+str(elementsize))
cubit.cmd("volume 1 2 size "+str(elementsize))
cubit.cmd("surface 12 scheme pave")
cubit.cmd("mesh surface 12 ")
cubit.cmd("mesh volume 1 2 ")
#cubit.cmd("unmerge surface 1")
########## loading cracks ##########
# SAVING FAULT NODES AND ELEMENTS. #
os.system('mkdir -p MESH')
########## FAULT A ##########
Au = [12] #face_up
Ad = [2] #face_down
faultA = fault_input(1,Au,Ad)
###### 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 2 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
# Material properties in concordance with tpv5 benchmark.
cubit.cmd('block 1 name "elastic 1" ') # material region
cubit.cmd('block 1 attribute count 5')
cubit.cmd('block 1 attribute index 1 1') # volume 1
cubit.cmd('block 1 attribute index 2 5477.2') # vp
cubit.cmd('block 1 attribute index 3 3162.3') # vs
cubit.cmd('block 1 attribute index 4 3000') # rho
cubit.cmd('block 1 attribute index 5 600') # Q flag (see constants.h: IATTENUATION_ ... )
# Material properties in concordance with tpv5 benchmark.
cubit.cmd('block 2 name "elastic 2" ') # material region
cubit.cmd('block 2 attribute count 5')
cubit.cmd('block 2 attribute index 1 1') # volume 2
cubit.cmd('block 2 attribute index 2 5477.2') # vp
cubit.cmd('block 2 attribute index 3 3162.3') # vs
cubit.cmd('block 2 attribute index 4 3000') # rho
cubit.cmd('block 2 attribute index 5 600') # Q flag (see constants.h: IATTENUATION_ ... )
#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT ####
cubit2specfem3d.export2SPECFEM3D('MESH')
# all files needed by SCOTCH are now in directory MESH