Error when ingesting mesh in PyLith

I’m encountering an error when ingesting the mesh I created in Gmsh to PyLith.

Here is the error.
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: No support for this operation for this object type
[0]PETSC ERROR: Could not determine Plex facet for Gmsh element 558 (Plex cell 1253)
[0]PETSC ERROR: See FAQ — PETSc 3.21.1 documentation for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.20.2-242-gc574b33dbb6 GIT Date: 2023-12-05 20:46:15 +0000

I’m also attaching the code to generate the meshes using Gmsh. (2.4 KB)

You added the fault geometry but it is not part of the volume, so the volume mesh is completely independent of your fault mesh. You need to use the embed function (api/ · gmsh_4_12_2 · gmsh / gmsh · GitLab) to embed the fault geometry in the volume geometry. Refer to the Gmsh t15 tutorial for how to embed geometry (tutorials/python/ · gmsh_4_12_2 · gmsh / gmsh · GitLab).

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:


# 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

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

# 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

# Embed the fault surface in the volumne
gmsh.model.mesh.embed(2, [faultSurface], 3, v1)

#Create a physical group for the homogeneous material
materials = (
     MaterialGroup(tag=1, entities=[v1]),)

for material in materials:

#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:
print("Vertex Groups: ",vertex_groups)


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



# Save the geometry to a file

It worked. Thank you @baagaard!