Importing meshes from Gmsh - dbcavalcanti/porousLab GitHub Wiki

Using Gmsh to Generate Arbitrary 2D Meshes for PorousLab

To generate meshes for non-rectangular domains in PorousLab, we recommend using Gmsh, an open-source mesh generator with a built-in scripting language and Python API.

This guide demonstrates how to use the Gmsh Python API to:

Define custom 2D geometries
Configure transfinite and recombined meshing
Export meshes into .mat format
Load them into PorousLab

📦 Requirements

Install the required Python packages:

pip install gmsh numpy scipy

🛠 Auxiliary Export Function

To export a mesh in the correct format for PorousLab, we provide a Python utility function:

🔗 Gmsh2Porouslab.py

This function exports a .mat containing:

node: a [N x 2] array of node coordinates
elem: a {Nelem x 1} cell array of element connectivities

These can be loaded directly into the MATLAB workspace.

Example of mesh created using Gmsh Python API

Below is the script for generating the mesh for the Slope Stability Example.

import gmsh
from Gmsh2Porouslab import gmsh2porouslab

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)

# === Geometry Definition ===
dim = 2
Nel = 10
lc = 10.0 / Nel

# Define points
p1  = gmsh.model.geo.addPoint( 0.0,  0.0, 0.0)
p2  = gmsh.model.geo.addPoint(75.0,  0.0, 0.0)
p3  = gmsh.model.geo.addPoint(75.0, 30.0, 0.0)
p4  = gmsh.model.geo.addPoint(45.0, 30.0, 0.0)
p5  = gmsh.model.geo.addPoint(35.0, 40.0, 0.0)
p6  = gmsh.model.geo.addPoint( 0.0, 40.0, 0.0)
p7  = gmsh.model.geo.addPoint(20.0, 40.0, 0.0)
p8  = gmsh.model.geo.addPoint(20.0, 30.0, 0.0)
p9  = gmsh.model.geo.addPoint(20.0, 28.0, 0.0)
p10 = gmsh.model.geo.addPoint(45.0, 28.0, 0.0)
p11 = gmsh.model.geo.addPoint(0.0, 6.0, 0.0)

# Define lines
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p5)
l5 = gmsh.model.geo.addLine(p5, p7)
l6 = gmsh.model.geo.addLine(p7, p6)
l7 = gmsh.model.geo.addLine(p6, p11)
l15 = gmsh.model.geo.addLine(p11, p1)

l8 = gmsh.model.geo.addLine(p7, p8)
l9 = gmsh.model.geo.addLine(p8, p9)
l10 = gmsh.model.geo.addLine(p9, p10)
l11 = gmsh.model.geo.addLine(p10, p4)
l12 = gmsh.model.geo.addLine(p4, p8)
l13 = gmsh.model.geo.addLine(p10, p2)
l14 = gmsh.model.geo.addLine(p9, p1)
l16 = gmsh.model.geo.addLine(p11, p8)

# Define curve loops and surfaces
curv1 = gmsh.model.geo.addCurveLoop([l1, -l13, -l10, l14])
curv2 = gmsh.model.geo.addCurveLoop([l2, l3, -l11, l13])
curv3 = gmsh.model.geo.addCurveLoop([l4, l5, l8, -l12])
curv4 = gmsh.model.geo.addCurveLoop([l6, l7, l16, -l8])
curv5 = gmsh.model.geo.addCurveLoop([l15, -l14, -l9, -l16])
curv6 = gmsh.model.geo.addCurveLoop([l10, l11, l12, l9])

surf1 = gmsh.model.geo.addPlaneSurface([curv1])
surf2 = gmsh.model.geo.addPlaneSurface([curv2])
surf3 = gmsh.model.geo.addPlaneSurface([curv3])
surf4 = gmsh.model.geo.addPlaneSurface([curv4])
surf5 = gmsh.model.geo.addPlaneSurface([curv5])
surf6 = gmsh.model.geo.addPlaneSurface([curv6])

# Recombine triangles into quads
for s in [surf1, surf2, surf3, surf4, surf5, surf6]:
    gmsh.model.geo.mesh.setRecombine(2, s)

# Set transfinite curves and surfaces
# (see full script for detailed settings)

# Synchronize and define physical group
gmsh.model.geo.synchronize()
meshDomain = gmsh.model.addPhysicalGroup(dim, [surf1, surf2, surf3, surf4, surf5, surf6], name='ContinuumDomain')

# Mesh generation
gmsh.model.mesh.generate(2)

# Convert to second-order (serendipity) elements
gmsh.option.setNumber("Mesh.SecondOrderIncomplete", 1)
gmsh.model.mesh.setOrder(2)

# Renumber nodes for matrix efficiency
old, new = gmsh.model.mesh.computeRenumbering("RCMK")
gmsh.model.mesh.renumberNodes(old, new)

# Export to .mat file
gmsh2porouslab(gmsh.model, meshDomain, "MeshSlopeStabilityTransfiniteQuadratic.mat")

# (Optional) Launch GUI
gmsh.fltk.run()
gmsh.finalize()

📥 Loading the Mesh in MATLAB

The Slope Stability Example demonstrates how to use it.

load('MeshSlopeStabilityTransfiniteQuadratic.mat');  % Loads node and elem

These arrays can be passed directly into your PorousLab model setup.

% Set mesh to model
mdl.setMesh(node, elem);

🔗 Related Resources

Feel free to adapt the geometry and mesh configuration to your problem setup.