PyMeshIt Headless Batch Workflow¶
This notebook shows how to run the new PyMeshIt headless workflow without opening the GUI.
Use it for simple batch cases where the GUI workflow would normally be:
Load surface and optional well files.
Compute hulls, segments, initial triangulations, and intersections.
In the Refine & Mesh tab, choose select intersections.
Generate conforming surface meshes.
Run TetGen and export the final mesh.
The equivalent API option is constraint_mode="intersections". With preserve_boundary_hulls=True, PyMeshIt keeps each surface hull as the PLC boundary and adds selected intersection lines inside it, matching the GUI benchmark workflow.
1. Install or Import PyMeshIt¶
If you installed PyMeshIt from PyPI after the headless API release, the normal import should work. If you are running this notebook from the GitHub repository, use the local repository copy or install the repository in editable mode first.
Important: inside Jupyter, use %pip install .... Do not type pip3 install ... as plain Python; that command only works in a terminal. After installing into the notebook kernel, restart the kernel and run the notebook again.
If you edit Pymeshit/headless.py while this notebook is already open, rerun the import cell below or restart the kernel. Jupyter otherwise keeps the old module in memory.
# Run ONE of these only if the imports below fail in a fresh environment.
# Keep the leading % when running inside Jupyter.
# If this notebook is opened from the cloned GitHub repository:
# %pip install -e ..
# If you want the PyPI release instead:
# %pip install pymeshit
# Terminal-only equivalent, not valid as plain Python in a notebook cell:
# pip3 install pymeshit
from pathlib import Path
import sys
import numpy as np
def find_repo_root(start: Path):
"""Find a local MeshIt checkout that contains Pymeshit/headless.py."""
for candidate in [start, *start.parents]:
if (candidate / "Pymeshit" / "headless.py").exists():
return candidate
return None
def path_is_under(path: Path, root: Path) -> bool:
try:
path.resolve().relative_to(root.resolve())
return True
except ValueError:
return False
# Prefer this source checkout over an older installed Pymeshit package.
# Also force-reload Pymeshit.headless from disk, because Jupyter caches modules.
REPO_ROOT = find_repo_root(Path.cwd())
if REPO_ROOT is not None:
sys.path.insert(0, str(REPO_ROOT))
existing = sys.modules.get("Pymeshit")
if existing is not None and getattr(existing, "__file__", None):
if not path_is_under(Path(existing.__file__), REPO_ROOT):
for module_name in list(sys.modules):
if module_name == "Pymeshit" or module_name.startswith("Pymeshit."):
del sys.modules[module_name]
else:
sys.modules.pop("Pymeshit.headless", None)
else:
sys.modules.pop("Pymeshit.headless", None)
try:
import Pymeshit.headless as headless_module
from Pymeshit.headless import (
MeshCase,
MeshOptions,
SurfaceSpec,
WellSpec,
MaterialSpec,
read_points,
run_mesh_case,
)
except ImportError as exc:
import Pymeshit
raise ImportError(
"This Python kernel is importing an older PyMeshIt package from:\n"
f" {getattr(Pymeshit, '__file__', 'unknown')}\n\n"
"Install the current repository into this notebook kernel with:\n"
" %pip install -e ..\n\n"
"Then restart the kernel and run the notebook again."
) from exc
print("PyMeshIt headless API imported successfully")
print("Using headless module:", Path(headless_module.__file__).resolve())
if REPO_ROOT is not None:
print("Using source checkout:", REPO_ROOT)
2. Put Input Files in One Folder¶
Accepted point inputs are text/CSV-style files with columns x y z, or .vtu/.vtk/.vtp files readable by PyVista.
Recommended naming convention:
outer/boundary surfaces:
bottom.txt,top.txt,xmin.txt, etc.geological layer/internal surfaces:
layer_1.txt,layer_2.txt, etc.faults:
fault_1.txt,fault_2.txt, etc.wells:
well_1.txt,well_2.txt, etc.
For a volume mesh, the border surfaces must define a closed outer domain. Internal layers should usually be unit; fault surfaces should be fault.
# Option A: local Jupyter / VS Code / JupyterLab
# This repo includes a small example dataset in examples/data.
# If you run the notebook from the repository root, this path works directly.
DEFAULT_DATA_DIR = Path("examples/data")
DATA_DIR = DEFAULT_DATA_DIR if DEFAULT_DATA_DIR.exists() else Path("data")
DATA_DIR.mkdir(exist_ok=True)
# Option B: Google Colab upload
# Uncomment this block in Colab, then upload your files.
# from google.colab import files
# uploaded = files.upload()
# for filename, content in uploaded.items():
# (DATA_DIR / filename).write_bytes(content)
print("Data folder:", DATA_DIR.resolve())
print("Files:")
for path in sorted(DATA_DIR.glob("*")):
print(" -", path.name)
3. Define Surfaces and Roles¶
The example below uses the files in examples/data. You can add any number of files to each group.
Roles:
border: outer boundary of the 3D domain. These are kept even in intersection-only mode.unit: internal stratigraphic/layer surfaces. You can have multiple unit surfaces.fault: internal 2D fault/fracture surfaces. You can have multiple fault surfaces. These become surface constraints, not volumetric material regions.
Each group has one target_size. Smaller values make denser meshes. You can also split files into separate groups if different surfaces need different sizes.
# Example input using examples/data.
# Add/remove filenames in each group for your own model.
surface_groups = [
{
"role": "border",
"target_size": 15.0,
"files": [
"bottom.dat",
"top.dat",
"border1.txt",
"border2.txt",
"border3.txt",
"border4.txt",
],
},
{
"role": "unit",
"target_size": 15.0,
"files": [
"middle.dat",
# Add more unit/layer files here, for example: "layer_2.dat"
],
},
{
"role": "fault",
"target_size": 15.0,
"files": [
"fault1.txt",
"fault2.txt",
],
},
]
def build_surfaces_from_groups(surface_groups, data_dir: Path):
surfaces = []
for group in surface_groups:
role = group["role"]
target_size = float(group.get("target_size", 15.0))
for filename in group["files"]:
path = data_dir / filename
if not path.exists():
print(f"WARNING: missing {path}. Edit surface_groups or upload the file.")
continue
pts = read_points(path)
name = path.stem
print(f"{name:20s} role={role:7s} target={target_size:6.2f} points={len(pts):5d} path={path}")
surfaces.append(SurfaceSpec(name=name, points=pts, role=role, target_size=target_size))
return surfaces
surfaces = build_surfaces_from_groups(surface_groups, DATA_DIR)
print("Loaded surfaces:", len(surfaces))
4. Optional Wells¶
Wells are 1D polylines. If you do not have wells, leave well_inputs empty.
well_inputs = [
# name, filename, target_size
# ("well_1", "well_1.txt", 20.0),
]
wells = []
for name, filename, target_size in well_inputs:
path = DATA_DIR / filename
if not path.exists():
print(f"WARNING: missing {path}. Edit well_inputs or upload the file.")
continue
pts = read_points(path)
print(f"{name:20s} points={len(pts)} path={path}")
wells.append(WellSpec(name=name, points=pts, target_size=target_size))
print("Loaded wells:", len(wells))
5. Define Material Seed Points¶
Each material seed point must be inside the volume region it represents. For a two-layer model, add one seed point in each layer.
If you define no materials, PyMeshIt uses one default region seed near the PLC center. For geological models, explicit seed points are safer.
# Example material seed points for examples/data.
# Put each seed inside the volume region it should label.
# With middle.dat near z=0, these two seeds label lower and upper domains.
materials = [
MaterialSpec("lower_unit", [0.0, 0.0, -25.0], attribute=1),
MaterialSpec("upper_unit", [0.0, 0.0, 25.0], attribute=2),
]
for material in materials:
print(material)
6. Run the Workflow and Write Exodus¶
constraint_mode="intersections" is the notebook equivalent of the GUI Select intersections button: PyMeshIt keeps each surface hull as the PLC boundary and adds selected intersection lines inside it.
Keep preserve_boundary_hulls=True for volume meshes. If you only want conforming surface meshes, set generate_volume=False.
For Golem/MOOSE, pass an .exo or .e path to output_path. The headless API then uses the same Exodus writer as the GUI instead of PyVista’s generic .save() path. The Exodus file contains tetrahedral material element blocks, surface/fault sidesets, and node sets derived from those sidesets. Well polylines are exported as node sets by default.
For the benchmark data, the fault surfaces should report constraint_source="selected". hull_fallback is only a safety fallback for cases where selected constraints still cannot be triangulated.
if not surfaces:
raise RuntimeError("No surfaces loaded. Upload files and edit surface_groups first.")
OUTPUT_DIR = Path("outputs")
OUTPUT_DIR.mkdir(exist_ok=True)
options = MeshOptions(
target_size=15.0,
min_angle=20.0,
gradient=2.0,
hull_method="delaunay",
interpolation="Thin Plate Spline (TPS)",
smoothing=0.0,
constraint_mode="intersections",
preserve_boundary_hulls=True,
include_hull_fallback_faults_in_volume=False,
tetgen_switches="pq1.414aAY",
generate_volume=True,
)
case = MeshCase(
surfaces=surfaces,
wells=wells,
materials=materials,
options=options,
)
result = run_mesh_case(case, output_path=OUTPUT_DIR / "example_data_case.exo")
if result.tetra_mesh is not None:
result.export_mesh(OUTPUT_DIR / "example_data_case.vtu")
print("Failures:", result.failures)
if result.tetra_mesh is not None:
print("Tetra mesh points:", result.tetra_mesh.n_points)
print("Tetra mesh cells:", result.tetra_mesh.n_cells)
print("Saved Exodus:", (OUTPUT_DIR / "example_data_case.exo").resolve())
print("Saved VTU:", (OUTPUT_DIR / "example_data_case.vtu").resolve())
else:
print("No tetra mesh was generated. Check failures and logs above.")
7. Exodus Sidesets, Node Sets, and Names¶
PyMeshIt assigns deterministic Exodus markers from the order of surfaces:
borderandunitsurfaces use markersurface_index + 1.faultsurfaces use marker1000 + surface_index.well node sets use
WellSpec.markerwhen supplied; otherwise they uselen(surfaces) + local_well_index + 2.boundary/fault node sets are created automatically from every sideset and named
nodes_<sideset_name>.
Use custom_sideset_names to control the names used by MOOSE boundary conditions and fault generators. Boundary node set names are derived from these sideset names. Use custom_block_names to name material element blocks and custom_well_names to name well node sets. Exodus names are sanitized and truncated to the Exodus name length limit.
For fault workflows such as LowerDBlockFromSidesetGenerator, use the fault sideset names, for example ss_fault1 and ss_fault2 below.
# Build the marker/name dictionaries used by the Exodus writer.
surface_marker_by_name = {}
custom_sideset_names = {}
for surface_idx, surface in enumerate(surfaces):
marker = 1000 + surface_idx if surface.role == "fault" else surface_idx + 1
surface_marker_by_name[surface.name] = marker
prefix = "ss_fault" if surface.role == "fault" else "ss_boundary"
custom_sideset_names[marker] = f"{prefix}_{surface.name}"
custom_block_names = {material.attribute: material.name for material in materials}
custom_well_names = {
(well.marker if well.marker is not None else len(surfaces) + well_idx + 2): well.name
for well_idx, well in enumerate(wells)
}
print("Surface markers:")
for name, marker in surface_marker_by_name.items():
print(f" {marker:5d} {name} -> {custom_sideset_names[marker]}")
print("Material block names:", custom_block_names)
print("Well node set names:", custom_well_names)
# Re-export an existing result with these names, or pass the same dictionaries
# directly to run_mesh_case(..., output_path=..., custom_sideset_names=...).
named_exodus_path = OUTPUT_DIR / "example_data_case_named.exo"
if result.tetra_mesh is not None:
ok = result.export_mesh(
named_exodus_path,
custom_block_names=custom_block_names,
custom_sideset_names=custom_sideset_names,
custom_well_names=custom_well_names,
well_export_type="Node Sets",
)
print("Named Exodus export:", ok, named_exodus_path.resolve())
# Optional: inspect Exodus block, sideset, and nodeset names.
# Requires netCDF4, which is installed with PyMeshIt from pyproject.toml.
try:
import netCDF4 as nc
except ImportError:
nc = None
print("Install netCDF4 to inspect Exodus metadata in Python: %pip install netCDF4")
def exodus_names(exodus_path: Path, variable_name: str):
if nc is None or not exodus_path.exists():
return []
with nc.Dataset(exodus_path) as dataset:
if variable_name not in dataset.variables:
return []
names = []
for row in dataset.variables[variable_name][:]:
raw = b"".join(row).split(b"\x00", 1)[0]
names.append(raw.decode("ascii", "ignore"))
return names
exodus_to_check = named_exodus_path if "named_exodus_path" in globals() else OUTPUT_DIR / "example_data_case.exo"
print("Element blocks:", exodus_names(exodus_to_check, "eb_names"))
print("Sidesets:", exodus_names(exodus_to_check, "ss_names"))
print("Node sets:", exodus_names(exodus_to_check, "ns_names"))
if nc is not None and exodus_to_check.exists():
with nc.Dataset(exodus_to_check) as dataset:
print("num_el_blk:", len(dataset.dimensions.get("num_el_blk", [])))
print("num_side_sets:", len(dataset.dimensions.get("num_side_sets", [])))
print("num_node_sets:", len(dataset.dimensions.get("num_node_sets", [])))
To write the named Exodus file in one pass, pass the same dictionaries to run_mesh_case:
result = run_mesh_case(
case,
output_path=OUTPUT_DIR / "case_named.exo",
custom_block_names=custom_block_names,
custom_sideset_names=custom_sideset_names,
custom_well_names=custom_well_names,
well_export_type="Node Sets",
)
Use well_export_type="Element Blocks" only when you want wells as BAR2 element blocks instead of well node sets. Boundary and fault node sets are still derived from the sidesets.
8. Inspect or Save Intermediate Surface Meshes¶
The result also contains conforming surface meshes. These are useful for checking whether the intersection-only constraint selection worked before TetGen.
for surface_idx, data in result.conforming_surface_data.items():
print(
surface_idx,
data["name"],
"source=", data.get("constraint_source", "selected"),
"vertices=", len(data["vertices"]),
"triangles=", len(data["triangles"]),
)
Optional PyVista visualization. This works best in local Jupyter. In some remote notebooks you may need off-screen rendering or skip plotting.
# Optional visualization
# import pyvista as pv
# if result.tetra_mesh is not None:
# scalars = "MaterialID" if "MaterialID" in result.tetra_mesh.cell_data else None
# result.tetra_mesh.plot(scalars=scalars, show_edges=True)
9. Batch Run Multiple Cases to Exodus¶
For parameter studies, put one complete geometry case in each folder and export one Exodus file per case. Each case folder should contain the filenames listed in surface_groups, and optionally the filenames listed in well_inputs.
Example layout:
data/batch_cases/
dip_30/
bottom.dat
top.dat
border1.txt
middle.dat
fault1.txt
fault2.txt
dip_45/
bottom.dat
top.dat
border1.txt
middle.dat
fault1.txt
fault2.txt
If your cases arrive as a zip file, extract it once into data/batch_cases and then run the batch cell. The output files are written to outputs/batch_exodus/<case_name>.exo.
import zipfile
def expected_surface_filenames(surface_groups):
names = []
for group in surface_groups:
names.extend(group["files"])
return names
def build_wells_from_inputs(well_inputs, data_dir: Path):
local_wells = []
for name, filename, target_size in well_inputs:
path = data_dir / filename
if not path.exists():
print(f" WARNING: missing well file {path}")
continue
local_wells.append(WellSpec(name=name, points=read_points(path), target_size=target_size))
return local_wells
def build_exodus_name_maps(local_surfaces, local_materials, local_wells):
custom_sideset_names = {}
for surface_idx, surface in enumerate(local_surfaces):
marker = 1000 + surface_idx if surface.role == "fault" else surface_idx + 1
prefix = "ss_fault" if surface.role == "fault" else "ss_boundary"
custom_sideset_names[marker] = f"{prefix}_{surface.name}"
custom_block_names = {material.attribute: material.name for material in local_materials}
custom_well_names = {
(well.marker if well.marker is not None else len(local_surfaces) + well_idx + 2): well.name
for well_idx, well in enumerate(local_wells)
}
return custom_block_names, custom_sideset_names, custom_well_names
def discover_case_folders(batch_root: Path, required_files):
if not batch_root.exists():
return []
cases = []
for case_dir in sorted(path for path in batch_root.iterdir() if path.is_dir()):
present = [filename for filename in required_files if (case_dir / filename).exists()]
if present:
cases.append((case_dir.name, case_dir))
return cases
# Optional: extract a zip containing one folder per case.
# The zip should contain folders such as dip_30/, dip_45/, ... with the same filenames.
BATCH_ROOT = DATA_DIR / "batch_cases"
BATCH_ROOT.mkdir(parents=True, exist_ok=True)
# BATCH_ZIP = Path("/path/to/batch_cases.zip")
# if BATCH_ZIP.exists():
# with zipfile.ZipFile(BATCH_ZIP) as archive:
# archive.extractall(BATCH_ROOT)
# Option A: explicit cases.
batch_cases = [
# ("dip_30", BATCH_ROOT / "dip_30"),
# ("dip_45", BATCH_ROOT / "dip_45"),
# ("dip_60", BATCH_ROOT / "dip_60"),
]
# Option B: auto-discover every subfolder that contains at least one expected surface file.
if not batch_cases:
batch_cases = discover_case_folders(BATCH_ROOT, expected_surface_filenames(surface_groups))
BATCH_OUTPUT_DIR = OUTPUT_DIR / "batch_exodus"
BATCH_OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
batch_summary = []
for case_name, case_dir in batch_cases:
print("Running", case_name, "from", case_dir)
local_surfaces = build_surfaces_from_groups(surface_groups, case_dir)
if not local_surfaces:
batch_summary.append({"case": case_name, "ok": False, "reason": "no surfaces loaded"})
print(" skipped: no surfaces loaded")
continue
local_wells = build_wells_from_inputs(well_inputs, case_dir)
local_case = MeshCase(
surfaces=local_surfaces,
wells=local_wells,
materials=materials,
options=options,
)
custom_block_names, custom_sideset_names, custom_well_names = build_exodus_name_maps(
local_surfaces,
materials,
local_wells,
)
output_exo = BATCH_OUTPUT_DIR / f"{case_name}.exo"
batch_result = run_mesh_case(
local_case,
output_path=output_exo,
custom_block_names=custom_block_names,
custom_sideset_names=custom_sideset_names,
custom_well_names=custom_well_names,
well_export_type="Node Sets",
)
ok = batch_result.tetra_mesh is not None and not batch_result.failures and output_exo.exists()
batch_summary.append({
"case": case_name,
"ok": ok,
"cells": None if batch_result.tetra_mesh is None else batch_result.tetra_mesh.n_cells,
"surfaces": len(local_surfaces),
"failures": batch_result.failures,
"output": output_exo,
})
print(" ok:", ok)
print(" failures:", batch_result.failures)
if batch_result.tetra_mesh is not None:
print(" cells:", batch_result.tetra_mesh.n_cells)
print(" exodus:", output_exo.resolve())
print("\nBatch summary")
for item in batch_summary:
print(item)
Troubleshooting¶
No surfaces loaded: upload files or fixsurface_inputsfilenames.No conforming surface meshes: the selected constraints do not form enough PLC geometry. Tryconstraint_mode="all"first, or keeppreserve_boundary_hulls=True.TetGen fails: check that
bordersurfaces form a closed domain and material seed points are inside valid regions.Exodus export fails: install
netCDF4in the active environment, then rerun the export cell.Wrong material IDs: move material seed points away from faults, holes, or boundaries.
Too coarse/fine mesh: change
target_sizeglobally or per surface.Only need surface meshes: set
generate_volume=FalseinMeshOptions.