I’d like to create a mesh for a simple model representing a 3-layer geological section using Gmsh, to be imported into the Specfem software. For the interface between the first and second layers, I plan to use points imported from a file. This will allow for variations in elevation, with a total of 31 points defining the interface. The interface between the second and third layers will be identical to the previous one, with a constant value subtracted. I’m using the Python programming language, specifically the ‘pygmsh’ module.

The issue I’m encountering is the following error message: ‘Exception: Surface 2 is transfinite but has 8 corners’.

I appreciate any assistance. Below is the code I have written:

import pygmsh
import numpy as np
import pandas as pd

dat_pth = "Points_Uni_Elev.csv"
df = pd.read_csv(dat_pth)

xmin = min(df.distance)
xmax = max(df.distance)
topo_max = max(df.Elev_m)

lc = 10

base = df.Elev_m - 400

pontos = []

# Initialize empty geometry using the build in kernel in GMSH
with pygmsh.geo.Geometry() as model:
    # model top limit points
    pontos.append(model.add_point((xmin, 500, 0), lc))
    pontos.append(model.add_point((xmax, 500, 0), lc))

    interface1_pts = [model.add_point((x, z, 0), lc) for x, z in zip(
        df.distance, df.Elev_m)]  # interface 1 points
    interface1_pts_r = interface1_pts[::-1]  # reversing

    interface2_pts = [model.add_point((x, z, 0), lc)
                      for x, z in zip(df.distance, base)]  # interface 2 points
    interface2_pts_r = interface2_pts[::-1]  # reversing

    # model base limit points
    pontos.append(model.add_point((xmax, min(base)-200, 0), lc))
    pontos.append(model.add_point((xmin, min(base)-200, 0), lc))

    # lines
    l_top = model.add_line(pontos[0], pontos[1])  # first layer lines
    l_r = model.add_line(pontos[1], interface1_pts_r[0])
    interface1_line_r = model.add_spline(interface1_pts_r)
    l_l = model.add_line(interface1_pts_r[-1], pontos[0])

    interface1_line = model.add_spline(interface1_pts)  # second layer lines
    l_r2 = model.add_line(interface1_pts_r[0], interface2_pts_r[0])
    interface2_line_r = model.add_spline(interface2_pts_r)
    l_l2 = model.add_line(interface2_pts_r[-1], interface1_pts_r[-1])

    interface2_line = model.add_spline(interface2_pts)  # third layer lines
    l_r3 = model.add_line(interface2_pts_r[0], pontos[2])
    botton_l_r = model.add_line(pontos[2], pontos[3])
    l_l3 = model.add_line(pontos[3], interface2_pts_r[-1])

    # loops
    loop_layer1 = model.add_curve_loop([l_top, l_r, interface1_line_r, l_l])
    loop_layer2 = model.add_curve_loop(
        [-interface1_line_r, l_r2, interface2_line_r, l_l2])
    loop_layer3 = model.add_curve_loop(
        [-interface2_line_r, l_r3, botton_l_r, l_l3])

    # surfaces
    ps_layer1 = model.add_plane_surface(loop_layer1)
    ps_layer2 = model.add_plane_surface(loop_layer2)
    ps_layer3 = model.add_plane_surface(loop_layer3)

    # transfinite curves
    n_elm_pml = 6
    model.set_transfinite_curve(l_top, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_r, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(
        interface1_line_r, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_l, n_elm_pml, "Progression", 1.0)

    model.set_transfinite_curve(interface1_line, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_r2, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(
        interface2_line_r, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_l2, n_elm_pml, "Progression", 1.0)

    model.set_transfinite_curve(interface2_line, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_r3, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(botton_l_r, n_elm_pml, "Progression", 1.0)
    model.set_transfinite_curve(l_l3, n_elm_pml, "Progression", 1.0)

    # transfinite surfaces
    model.set_transfinite_surface(ps_layer1, "Left", corner_pts=[
                                  pontos[0], pontos[1], interface1_pts_r[0], interface1_pts_r[-1]])
    model.set_transfinite_surface(ps_layer2, "Left", corner_pts=[
                                  interface1_pts[0], interface1_pts[-1], interface2_pts_r[0], interface2_pts_r[-1]])
    model.set_transfinite_surface(ps_layer2, "Left", corner_pts=[
                                  interface2_pts[0], interface2_pts[-1], pontos[2], pontos[3]])

    model.set_recombined_surfaces([ps_layer1, ps_layer2, ps_layer3])

    # synchronize
    model.synchronize()

    # physical groups
    # boundary
    pg_top = model.add_physical(l_top, label="Top")
    pg_bot = model.add_physical(botton_l_r, label="Bottom")
    pg_left = model.add_physical([l_l, l_l2, l_l3], label="Left")
    pg_right = model.add_physical([l_r, l_r2, l_r3], label="Right")

    # material
    pg_m1 = model.add_physical(ps_layer1, label="M1")  # arenito
    pg_m2 = model.add_physical(ps_layer2, label="M2")  # basalto
    pg_m3 = model.add_physical(ps_layer3, label="M3")  # arenito

    # Generate mesh (meshio object)
    mesh = model.generate_mesh(dim=2, order=2, verbose=True)

    mesh.write("fe_mesh.msh")  # ,file_format='xdmf'