Because your fault surface touches the top surface, you need to embed the surface trace of the fault in the top surface and embed the fault surface in the volume.
For reference, here is the complete script:
#!/usr/bin/python3
# Code to generate PyLith compatible 3D mesh using Gmsh for
# with a fault surface embedded in it
import gmsh
from pylith.meshio.gmsh_utils import (VertexGroup, MaterialGroup, GenerateMesh)
import utm
import numpy as np
gmsh.initialize()
gmsh.clear()
filename = "3d_mesh_rev.msh"
# Define the coordinates of the fault from the slip model
faultXComp = np.array([69.510483, 69.420603, 69.363187, 69.453003])
faultYComp = np.array([33.061283, 32.917904, 32.943429, 33.086843])
# Define the depth of the fault
Depth = np.array([0, 0, -10, -10])
numberCode = np.empty(len(faultXComp), dtype=str)
letterCode = np.empty(len(faultXComp), dtype=str)
# Convert WGS84 coordinates to UTM
# Set as the center of the model
for i in range(len(faultXComp)):
faultXComp[i], faultYComp[i], numberCode[i], letterCode[i] \
= utm.from_latlon(faultYComp[i], faultXComp[i])
# Make the coordinates relative to the center of the model
# Select the hexahedral mesh
faultXCompCenter = np.mean(faultXComp)
faultYCompCenter = np.mean(faultYComp)
faultXComp = faultXComp - faultXCompCenter
faultYComp = faultYComp - faultYCompCenter
cell = "hex"
print("FaultXComp: ",faultXComp)
print("FaultYComp: ",faultYComp)
print("Depth: ",Depth)
# Define the corner of the box
xmin = np.min(faultXComp) * 2
ymin = np.min(faultYComp) * 2
zmax = Depth[0]
xmax = np.max(faultXComp) * 2
ymax = np.max(faultYComp) * 2
zmin = float(Depth[2] *1000) -5000
# Create a 3D model of the fault
# Create points for depth = 0 km for a rectangular prism
p1 = gmsh.model.geo.addPoint(xmin, ymin, zmax)
p2 = gmsh.model.geo.addPoint(xmin, ymax, zmax)
p3 = gmsh.model.geo.addPoint(xmax, ymax, zmax)
p4 = gmsh.model.geo.addPoint(xmax, ymin, zmax)
# Create points for depth = -12 km for a rectangular prism
p5 = gmsh.model.geo.addPoint(xmin, ymin, zmin)
p6 = gmsh.model.geo.addPoint(xmin, ymax, zmin)
p7 = gmsh.model.geo.addPoint(xmax, ymax, zmin)
p8 = gmsh.model.geo.addPoint(xmax, ymin, zmin)
# Create 6surfaces of the rectangular prism boundary surface
# from 12 lines
boundaryCondition1 = gmsh.model.geo.addLine(p1, p2)
boundaryCondition2 = gmsh.model.geo.addLine(p2, p3)
boundaryCondition3 = gmsh.model.geo.addLine(p3, p4)
boundaryCondition4 = gmsh.model.geo.addLine(p4, p1)
boundaryCondition5 = gmsh.model.geo.addLine(p5, p6)
boundaryCondition6 = gmsh.model.geo.addLine(p6, p7)
boundaryCondition7 = gmsh.model.geo.addLine(p7, p8)
boundaryCondition8 = gmsh.model.geo.addLine(p8, p5)
boundaryCondition9 = gmsh.model.geo.addLine(p1, p5)
boundaryCondition10 = gmsh.model.geo.addLine(p4, p8)
boundaryCondition11 = gmsh.model.geo.addLine(p3, p7)
boundaryCondition12 = gmsh.model.geo.addLine(p2, p6)
boundaryCurve1 = gmsh.model.geo.addCurveLoop([boundaryCondition1, boundaryCondition2,\
boundaryCondition3, boundaryCondition4])
boundarySurface1 = gmsh.model.geo.addPlaneSurface([boundaryCurve1])
boundaryCurve2 = gmsh.model.geo.addCurveLoop([boundaryCondition5, boundaryCondition6,\
boundaryCondition7, boundaryCondition8])
boundarySurface2 = gmsh.model.geo.addPlaneSurface([boundaryCurve2])
boundaryCurve3 = gmsh.model.geo.addCurveLoop([-boundaryCondition4, boundaryCondition10,\
boundaryCondition8, -boundaryCondition9])
boundarySurface3 = gmsh.model.geo.addPlaneSurface([boundaryCurve3])
boundaryCurve4 = gmsh.model.geo.addCurveLoop([boundaryCondition2, boundaryCondition11,\
-boundaryCondition6, -boundaryCondition12])
boundarySurface4 = gmsh.model.geo.addPlaneSurface([boundaryCurve4])
boundaryCurve5 = gmsh.model.geo.addCurveLoop([boundaryCondition1, boundaryCondition12,\
-boundaryCondition5, -boundaryCondition9])
boundarySurface5 = gmsh.model.geo.addPlaneSurface([boundaryCurve5])
boundaryCurve6 = gmsh.model.geo.addCurveLoop([-boundaryCondition3, boundaryCondition11,\
boundaryCondition7, -boundaryCondition10])
boundarySurface6 = gmsh.model.geo.addPlaneSurface([boundaryCurve6])
# Create the fault surface
p9 = gmsh.model.geo.addPoint(faultXComp[0], faultYComp[0], Depth[0]*1000)
p10 = gmsh.model.geo.addPoint(faultXComp[1], faultYComp[1], Depth[1]*1000)
p11 = gmsh.model.geo.addPoint(faultXComp[2], faultYComp[2], Depth[2]*1000)
p12 = gmsh.model.geo.addPoint(faultXComp[3], faultYComp[3], Depth[3]*1000)
faultline1 = gmsh.model.geo.addLine(p9, p10)
faultline2 = gmsh.model.geo.addLine(p10, p11)
faultline3 = gmsh.model.geo.addLine(p11, p12)
faultline4 = gmsh.model.geo.addLine(p12, p9)
faultCurve = gmsh.model.geo.addCurveLoop([faultline1, faultline2, faultline3, faultline4])
faultSurface = gmsh.model.geo.addPlaneSurface([faultCurve])
#Synchronize the geometry in Gmsh
gmsh.model.geo.synchronize()
# Embed the surface trace in the top boundary surface
gmsh.model.mesh.embed(1, [faultline1], 2, boundarySurface1)
#Create volume from the boundary surfaces
materialVolume = gmsh.model.geo.addSurfaceLoop([boundarySurface5, boundarySurface1,\
boundarySurface6, boundarySurface4, boundarySurface2, boundarySurface3])
v1 = gmsh.model.geo.addVolume([materialVolume])
#Synchronize the geometry in Gmsh
gmsh.model.geo.synchronize()
# Embed the fault surface in the volumne
gmsh.model.mesh.embed(2, [faultSurface], 3, v1)
gmsh.model.geo.synchronize()
#Create a physical group for the homogeneous material
materials = (
MaterialGroup(tag=1, entities=[v1]),)
for material in materials:
material.create_physical_group()
#Create vertex groups for the boundary conditions and fault
vertex_groups = (
VertexGroup(name="boundary_xneg", tag=2, dim=2, entities=[boundarySurface5]),
VertexGroup(name="boundary_xpos", tag=3, dim=2, entities=[boundarySurface6]),
VertexGroup(name="boundary_yneg", tag=4, dim=2, entities=[boundarySurface3]),
VertexGroup(name="boundary_ypos", tag=5, dim=2, entities=[boundarySurface4]),
VertexGroup(name="boundary_zneg", tag=6, dim=2, entities=[boundarySurface2]),
VertexGroup(name="boundary_zpos", tag=7, dim=2, entities=[boundarySurface1]),
VertexGroup(name="fault", tag=8, dim=2, entities=[faultSurface]),
VertexGroup(name="fault_edges", tag=9, dim=1, entities=[faultline1, faultline2,\
faultline3, faultline4])
)
for group in vertex_groups:
group.create_physical_group()
print("Vertex Groups: ",vertex_groups)
gmsh.model.geo.synchronize()
# We turn off the default sizing methods.
gmsh.option.set_number("Mesh.MeshSizeFromPoints", 0)
gmsh.option.set_number("Mesh.MeshSizeFromCurvature", 0)
gmsh.option.set_number("Mesh.MeshSizeExtendFromBoundary", 0)
field_distance = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(field_distance, "CurvesList", [faultSurface])
field_size = gmsh.model.mesh.field.add("MathEval")
math_exp = GenerateMesh.get_math_progression(field_distance, min_dx=4.0e+3, bias=1.05)
gmsh.model.mesh.field.setString(field_size, "F", math_exp)
gmsh.model.mesh.field.setAsBackgroundMesh(field_size)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Laplace2D")
# Save the geometry to a file
gmsh.write(filename)
gmsh.finalize()