Add background mesh support (#62)

README.rst

@@ -123,6 +123,113 @@ Here is the cell quality as computed according to the minimum scaled jacobian.
 .. image:: https://github.com/pyvista/tetgen/raw/main/doc/images/sphere_qual.png
 
 
+Using a Background Mesh
+-----------------------
+A background mesh in TetGen is used to define a mesh sizing function for
+adaptive mesh refinement. This function informs TetGen of the desired element
+size throughout the domain, allowing for detailed refinement in specific areas
+without unnecessary densification of the entire mesh. Here's how to utilize a
+background mesh in your TetGen workflow:
+
+1. **Generate the Background Mesh**: Create a tetrahedral mesh that spans the
+   entirety of your input piecewise linear complex (PLC) domain. This mesh will
+   serve as the basis for your sizing function.
+
+2. **Define the Sizing Function**: At the nodes of your background mesh, define
+   the desired mesh sizes. This can be based on geometric features, proximity
+   to areas of interest, or any criterion relevant to your simulation needs.
+
+3. **Optional: Export the Background Mesh and Sizing Function**: Save your
+   background mesh in the TetGen-readable `.node` and `.ele` formats, and the
+   sizing function values in a `.mtr` file. These files will be used by TetGen
+   to guide the mesh generation process.
+
+4. **Run TetGen with the Background Mesh**: Invoke TetGen, specifying the
+   background mesh. TetGen will adjust the mesh according to the provided
+   sizing function, refining the mesh where smaller elements are desired.
+
+**Full Example**
+
+To illustrate, consider a scenario where you want to refine a mesh around a
+specific region with increased detail. The following steps and code snippets
+demonstrate how to accomplish this with TetGen and PyVista:
+
+1. **Prepare Your PLC and Background Mesh**:
+
+   .. code-block:: python
+
+      import pyvista as pv
+      import tetgen
+      import numpy as np
+
+      # Load or create your PLC
+      sphere = pv.Sphere(theta_resolution=10, phi_resolution=10)
+
+      # Generate a background mesh with desired resolution
+      def generate_background_mesh(bounds, resolution=20, eps=1e-6):
+          x_min, x_max, y_min, y_max, z_min, z_max = bounds
+          grid_x, grid_y, grid_z = np.meshgrid(
+              np.linspace(xmin - eps, xmax + eps, resolution),
+              np.linspace(ymin - eps, ymax + eps, resolution),
+              np.linspace(zmin - eps, zmax + eps, resolution),
+              indexing="ij",
+          )
+          return pv.StructuredGrid(grid_x, grid_y, grid_z).triangulate()
+
+      bg_mesh = generate_background_mesh(sphere.bounds)
+
+2. **Define the Sizing Function and Write to Disk**:
+
+   .. code-block:: python
+
+      # Define sizing function based on proximity to a point of interest
+      def sizing_function(points, focus_point=np.array([0, 0, 0]), max_size=1.0, min_size=0.1):
+          distances = np.linalg.norm(points - focus_point, axis=1)
+          return np.clip(max_size - distances, min_size, max_size)
+
+      bg_mesh.point_data['target_size'] = sizing_function(bg_mesh.points)
+
+      # Optionally write out the background mesh
+      def write_background_mesh(background_mesh, out_stem):
+          """Write a background mesh to a file.
+
+          This writes the mesh in tetgen format (X.b.node, X.b.ele) and a X.b.mtr file
+          containing the target size for each node in the background mesh.
+          """
+          mtr_content = [f"{background_mesh.n_points} 1"]
+          target_size = background_mesh.point_data["target_size"]
+          for i in range(background_mesh.n_points):
+              mtr_content.append(f"{target_size[i]:.8f}")
+
+      write_background_mesh(bg_mesh, 'bgmesh.b')
+
+3. **Use TetGen with the Background Mesh**:
+
+
+   Directly pass the background mesh from PyVista to ``tetgen``:
+
+   .. code-block:: python
+
+      tet_kwargs = dict(order=1, mindihedral=20, minratio=1.5)
+      tet = tetgen.TetGen(mesh)
+      tet.tetrahedralize(bgmesh=bgmesh, **tet_kwargs)
+      refined_mesh = tet.grid
+
+   Alternatively, use the background mesh files.
+
+   .. code-block:: python
+
+      tet = tetgen.TetGen(sphere)
+      tet.tetrahedralize(bgmeshfilename='bgmesh.b', **tet_kwargs)
+      refined_mesh = tet.grid
+
+
+This example demonstrates generating a background mesh, defining a spatially
+varying sizing function, and using this background mesh to guide TetGen in
+refining a PLC. By following these steps, you can achieve adaptive mesh
+refinement tailored to your specific simulation requirements.
+
+
 Acknowledgments
 ---------------
 Software was originally created by Hang Si based on work published in

pyproject.toml

@@ -19,12 +19,12 @@ filterwarnings = [
 [tool.cibuildwheel]
 archs = ["auto64"]  # 64-bit only
 skip = "pp* *musllinux* cp37-*"  # disable PyPy, musl-based wheels, and Python < 3.8
-test-requires = "pytest pymeshfix"
+before-test = "pip install -r requirements_test.txt"
 test-command = "pytest {project}/tests"
 
 [tool.cibuildwheel.macos]
 # https://cibuildwheel.readthedocs.io/en/stable/faq/#apple-silicon
-archs = ["x86_64", "universal2"]
+archs = ["universal2"]
 test-skip = ["*_arm64", "*_universal2:arm64"]
 
 [tool.codespell]

requirements_test.txt

@@ -1,3 +1,4 @@
+meshio
+pymeshfix>=0.13.4
 pytest
 pytest-cov
-pymeshfix>=0.13.4

tests/test_background_mesh.py

@@ -0,0 +1,121 @@
+"""Test the mesh resizing feature of tetgen with sizing function."""
+from pathlib import Path
+import tempfile
+
+import numpy as np
+import pyvista as pv
+import tetgen
+
+
+def sizing_function(points):
+    """Return the target size at a given point.
+
+    This is just an example. You can use any function you want.
+    """
+    x, y, z = points.T
+    return np.where(x < 0, 0.5, 0.1)
+
+
+def generate_background_mesh(boundary_mesh, resolution=20, eps=1e-6):
+    """Generate a new mesh with the same bounds as the boundary meshy.
+
+    We will use this as a background mesh for the sizing function.
+    """
+    xmin, xmax, ymin, ymax, zmin, zmax = boundary_mesh.bounds
+    new_vertices = np.meshgrid(
+        np.linspace(xmin - eps, xmax + eps, resolution),
+        np.linspace(ymin - eps, ymax + eps, resolution),
+        np.linspace(zmin - eps, zmax + eps, resolution),
+        indexing="ij",
+    )
+
+    # tetgen supports only tetrahedral meshes
+    return pv.StructuredGrid(*new_vertices).triangulate()
+
+
+def write_background_mesh(background_mesh, out_stem):
+    """Write a background mesh to a file.
+
+    This writes the mesh in tetgen format (X.b.node, X.b.ele) and a X.b.mtr file
+    containing the target size for each node in the background mesh.
+    """
+    mtr_content = [f"{background_mesh.n_points} 1"]
+    target_size = background_mesh.point_data["target_size"]
+    for i in range(background_mesh.n_points):
+        mtr_content.append(f"{target_size[i]:.8f}")
+
+    pv.save_meshio(f"{out_stem}.node", background_mesh)
+    mtr_file = f"{out_stem}.mtr"
+
+    with open(mtr_file, "w") as f:
+        f.write("\n".join(mtr_content))
+
+
+def mesh_resizing_with_bgmeshfilename(mesh, bgmesh, **kwargs):
+    """Performs mesh resizing with a background mesh file."""
+    with tempfile.TemporaryDirectory() as tmpdir:
+        tmpfile = Path(tmpdir).joinpath("bgmesh.b")
+        write_background_mesh(bgmesh, tmpfile)
+
+        # Pass the background mesh file to tetgen
+        tet = tetgen.TetGen(mesh)
+        tet.tetrahedralize(bgmeshfilename=str(tmpfile), metric=1, **kwargs)
+        grid = tet.grid
+
+    # extract cells below the 0 xy plane
+    # cell_center = grid.cell_centers().points
+    # subgrid = grid.extract_cells(cell_center[:, 2] < 0)
+
+    # debug: plot this
+    # plotter = pv.Plotter()
+    # plotter.set_background("w")
+    # plotter.add_mesh(subgrid, "lightgrey", lighting=True)
+    # plotter.add_mesh(grid, "r", "wireframe")
+    # plotter.add_legend([[" Input Mesh ", "r"], [" Tessellated Mesh ", "black"]])
+    # plotter.show() # Uncomment for visualisation of resized mesh
+
+    return grid
+
+
+def mesh_resizing_with_pyvista_bgmesh(mesh, bgmesh, **kwargs):
+    """Performs mesh resizing with a pyvista bgmesh."""
+    # Pass the background mesh to tetgen
+    tet = tetgen.TetGen(mesh)
+    tet.tetrahedralize(bgmesh=bgmesh, metric=1, **kwargs)
+    grid = tet.grid
+
+    # Extract cells below the 0 xy plane
+    # cell_center = grid.cell_centers().points
+    # subgrid = grid.extract_cells(cell_center[:, 2] < 0)
+
+    # Debug: uncomment for visualisation of resized mesh
+    # plotter = pv.Plotter()
+    # plotter.set_background("w")
+    # plotter.add_mesh(subgrid, "lightgrey", lighting=True)
+    # plotter.add_mesh(grid, "r", "wireframe")
+    # plotter.add_legend([[" Input Mesh ", "r"], [" Tessellated Mesh ", "black"]])
+    # plotter.show()
+    return grid
+
+
+def test_mesh_resizing():
+    """Test the mesh resizing feature of tetgen with sizing function."""
+    sphere = pv.Sphere(theta_resolution=10, phi_resolution=10)
+    tet_kwargs = dict(order=1, mindihedral=20, minratio=1.5)
+
+    # Vanilla tetgen for reference
+    tet = tetgen.TetGen(sphere)
+    tet.tetrahedralize(**tet_kwargs)
+    grid = tet.grid
+
+    # Generate background mesh
+    bgmesh = generate_background_mesh(sphere)
+    bgmesh.point_data["target_size"] = sizing_function(bgmesh.points)
+
+    resized_grid_file = mesh_resizing_with_bgmeshfilename(sphere, bgmesh, **tet_kwargs)
+    assert resized_grid_file.n_points >= grid.n_points
+
+    resized_grid_direct = mesh_resizing_with_pyvista_bgmesh(sphere, bgmesh, **tet_kwargs)
+    assert resized_grid_direct.n_points > grid.n_points
+
+    assert resized_grid_file == resized_grid_direct