Implement simple Laplace solver and some related utilities

.pre-commit-config.yaml

@@ -29,7 +29,7 @@ repos:
     hooks:
       - id: bandit
         name: Run bandit (security linter for python)
-        args: ["-c", "pyproject.toml", "-r", "-ll", "biomesh/", "test_biomesh/"]
+        args: ["-c", "pyproject.toml", "-r", "-ll", "src/biomesh/", "tests/"]
         additional_dependencies: ["bandit[toml]"]
   - repo: https://github.com/PyCQA/docformatter
     rev: v1.7.7
@@ -41,8 +41,8 @@ repos:
             "--in-place",
             "--config=./pyproject.toml",
             "-r",
-            "biomesh/",
-            "test_biomesh/",
+            "src/biomesh/",
+            "tests/",
           ]
   - repo: https://github.com/econchick/interrogate
     rev: 1.7.0
@@ -64,7 +64,7 @@ repos:
             "--check-untyped-defs",
             "--disallow-untyped-decorators",
           ]
-        exclude: "test_biomesh/"
+        exclude: "tests/"
   - repo: https://github.com/Lucas-C/pre-commit-hooks
     rev: v1.5.5
     hooks:

README.md

@@ -2,7 +2,7 @@
 
 [![pipeline](https://github.com/TUM-LNM/biomesh/actions/workflows/build_and_test.yml/badge.svg)](https://github.com/TUM-LNM/biomesh/actions/workflows/build_and_test.yml)
 
-## Installation
+## :rocket: Installation
 
 biomesh can be installed via pip as
 
@@ -16,27 +16,100 @@ If you want to generate a mesh with Gmsh, you also need to install the gmsh pyth
 pip install gmsh
 ```
 
-## Usage
+## :book: Usage
 
-To generate a finite element mesh from multiple colored stl-files, you simply need to do
+`biomesh` is composed of multiple utilities for working with complex biomechanical geometries. Below
+are some common workflows
+
+### Generate a mesh from colored STL files
+
+Colored STL-files (e.g., exported from Materialise 3-matic) can be used to generate a volume mesh.
+Although STL color encoding is not standardized, some software packages embed surface IDs in unused
+byte fields. `biomesh` leverages this to extract surface information.
+
+:point_right: Generating a mesh from stl-files requires the `gmsh` Python package.
 
 ```python
 import biomesh
 
-# note, this requires gmsh to be installed (pip install gmsh)
 mesh = biomesh.mesh_colored_stl_files(
     "path/to/part1.stl",
     "path/to/part2.stl",
     "path/to/part3.stl",
     mesh_size=2.0
 )
-# alternatively, you can also just use meshio to load a mesh: mesh = meshio.read("path/to/mesh.vtu")
+```
+
+Alternatively, you can load the meshes from any format supported by [meshio](https://github.com/nschloe/meshio):
+
+```python
+import meshio
+
+meshio.read("path/to/mesh.vtu")
+```
+
+### Convert linear to quadratic elements
 
-# make all elements quadratic
+Convert linear elements in your mesh to quadratic ones:
+
+```python
 mesh = biomesh.lin_to_quad(mesh)
+```
+
+### Reorder mesh nodes
+
+Finite element solvers often benefit from reducing bandwidth in the system matrix. `biomesh` provides
+a node reordering algorithm based on Cuthill-McKee's algorithm to improve efficiency:
 
-# reorder nodes to reduce bandwidth of matrix
+```python
 mesh = biomesh.reorder(mesh)
 ```
 
-The nodes of all stl-files are matched. Each stl-file will be considered as an own volume.
+### Merge multiple meshes
+
+Combine several meshes into a single mesh object:
+
+```python
+mesh_all = biomesh.merge(mesh1, mesh2, mesh3)
+```
+
+:warning: Overlapping points are **not** automatically merged.
+
+### Filter a mesh
+
+Extract a subset of a mesh using flexible filters. For example, filtering by cell type:
+
+```python
+# keep only hexahedral cells
+filter_hex = lambda block : block.type == 'hexahedron'
+mesh_filtered = biomesh.filter.by_cellblock(mesh, filter_hex)
+
+# Get point mapping from old -> new IDs
+point_mapping = biomesh.filter.points_map_by_cellblock(mesh, filter_hex)
+```
+
+### Solve simple Laplace problem
+
+`biomesh` can also solve basic Laplace problems, commonly used to estimate fiber directions with
+rule based methods.
+
+```python
+dbc_nodes = np.array([0, 1, 2, 4])
+dbc_values = np.array([0.0, 0.0, 1.0, 1.0])
+
+phi = biomesh.solve(mesh, dbc_nodes, dbc_values)
+
+# or equivalently
+phi = biomesh.solve_onezero(mesh, np.array([2, 4]), np.array([0, 1]))
+```
+
+The result `phi` is the solution vector of the Laplace problem.
+
+### Finite Element Utilities
+
+`biomesh.fe` provides helper functions for finite element analysis. For example, compute the nodal
+averaged gradient of a scalar field:
+
+```python
+grad_phi = biomesh.fe.grad(mesh, phi)
+```

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "biomesh"
-version = "0.4.1"
+version = "0.5.0"
 authors = [
   { name="The biomesh Authors" },
 ]
@@ -31,8 +31,12 @@ dependencies = [
     "PyYAML>=6.0.0, <7",
     "scipy",
     "lnmmeshio>=5.6.4",
+    "symfem==2025.6.0",
 ]
 
 [project.urls]
 "Homepage" = "https://github.com/TUM-LNM/biomesh"
 "Bug Tracker" = "https://github.com/TUM-LNM/biomesh/issues"
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/biomesh"]

requirements-dev.txt

@@ -6,7 +6,7 @@ pytest-cov==4.0.0
 pytest-isort==3.1.0
 pytest-sugar==0.9.6
 black==24.3.0
-pre-commit==3.0.2
+pre-commit==4.3.0
 parameterized==0.8.1
 coverage==5.2.1
 pdoc==14.5.1

requirements.txt

@@ -6,3 +6,4 @@ tqdm==4.66.5
 lnmmeshio==5.6.4
 scipy
 gmsh==4.14.0
+symfem==2025.6.0

biomesh/__init__.py → src/biomesh/__init__.py

@@ -3,13 +3,9 @@
 # See the LICENSE file in the top-level for license information.
 #
 # SPDX-License-Identifier: MIT
-"""
-biomesh
-=======
-
-biomesh is a Python package for working with 3D meshes, providing tools for mesh generation, manipulation, and analysis
-tailored for finite element simulations of biomechanical applications.
-"""
+"""Biomesh is a Python package for working with 3D meshes, providing tools for
+mesh generation, manipulation, and analysis tailored for finite element
+simulations of biomechanical applications."""
 from . import run_gmsh
 import pathlib
 from . import mesh
@@ -18,6 +14,15 @@
 import meshio
 from .reorder import reorder
 from .adapt import lin_to_quad
+from .merge import merge
+from .filter import (
+    filter_by_cellblock,
+    filter_by_block_ids,
+    filter_by_cellblock_point_mapping,
+)
+
+from . import utils
+from . import laplace
 
 
 def combine_colored_stl_files(*stl_files: pathlib.Path) -> meshio.Mesh:
@@ -29,17 +34,14 @@ def combine_colored_stl_files(*stl_files: pathlib.Path) -> meshio.Mesh:
 def mesh_colored_stl_files(*stl_files: pathlib.Path, mesh_size: float) -> meshio.Mesh:
     """Generate a mesh from multiple colored STL files.
 
-    Parameters
-    ----------
-    *stl_files : pathlib.Path
-        Paths to the STL files to be merged.
+    Args:
+        *stl_files:
+            Paths to the STL files to be merged.
 
-    mesh_size : float
-        The target size for the mesh elements.
+        mesh_size:
+            The target size for the mesh elements.
 
-    Returns
-    -------
-    meshio.Mesh
+    Returns:
         The generated mesh.
     """
     assert len(stl_files) > 0, "At least one STL file must be provided."

biomesh/adapt.py → src/biomesh/adapt.py

@@ -81,12 +81,11 @@ def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
     are converted into their corresponding quadratic elements.
 
     Args:
-    -------
-    mesh: meshio.Mesh
-        The input mesh containing linear elements.
+        mesh:
+            The input mesh containing linear elements.
 
     Returns:
-        meshio.Mesh: The modified mesh with quadratic elements.
+        The modified mesh with quadratic elements.
     """
 
     new_points = [coord for coord in mesh.points]
@@ -103,11 +102,7 @@ def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
         new_cells = []
         cell_type = cellblock.type
 
-        if cell_type in ["vertex"]:
-            # vertices do not need to be converted
-            continue
-
-        if cell_type in ["hexahedron27", "tetra10", "triangle6", "quad9"]:
+        if cell_type in ["vertex", "hexahedron27", "tetra10", "triangle6", "quad9"]:
             # this cell block is already quadratic, no need to convert
             new_cell_blocks.append(cellblock)
             for key, data in mesh.cell_data.items():