Implement lin_to_quad for wedges and to_serendipity

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "biomesh"
-version = "0.5.1"
+version = "0.6.0"
 authors = [
   { name="The biomesh Authors" },
 ]

src/biomesh/adapt.py

@@ -6,6 +6,8 @@
 """Small utilities for mesh adaption (e.g. converting linear elements to
 quadratic elements)"""
 
+from typing import Callable, Optional
+
 import meshio
 import numpy as np
 
@@ -54,6 +56,23 @@
         [1, 3],
         [2, 3],
     ],
+    "wedge": [
+        [0],
+        [1],
+        [2],
+        [3],
+        [4],
+        [5],
+        [0, 1],
+        [1, 2],
+        [2, 0],
+        [3, 4],
+        [4, 5],
+        [5, 3],
+        [0, 3],
+        [1, 4],
+        [2, 5],
+    ],
     "triangle": [
         [0],
         [1],
@@ -68,26 +87,69 @@
 _lin_to_quad_cell_type = {
     "hexahedron": "hexahedron27",
     "tetra": "tetra10",
+    "wedge": "wedge15",
     "triangle": "triangle6",
     "quad": "quad9",
     "line": "line3",
 }
 
+_hex8_to_hex20 = [
+    [0],
+    [1],
+    [2],
+    [3],
+    [4],
+    [5],
+    [6],
+    [7],
+    [0, 1],
+    [1, 2],
+    [2, 3],
+    [3, 0],
+    [4, 5],
+    [5, 6],
+    [6, 7],
+    [7, 4],
+    [0, 4],
+    [1, 5],
+    [2, 6],
+    [3, 7],
+]
 
-def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
-    """Convert linear elements to quadratic elements in a mesh.
 
-    This function returns a new mesh in which all linear elements (triangles, quadrilaterals, tetrahedra, and hexahedra)
-    are converted into their corresponding quadratic elements.
+def adapt_mesh(
+    mesh: meshio.Mesh,
+    adapting_scheme: Callable[[str], Optional[tuple[str, list[list[int]]]]],
+) -> meshio.Mesh:
+    """Adapt a mesh to a higher-order (quadratic) representation using a
+    provided adapting scheme.
 
-    Args:
-        mesh:
-            The input mesh containing linear elements.
+    This function takes a mesh and, for each cell block, determines whether it should be adapted
+    to a new cell type. If so, it generates new mid-edge nodes and updates
+    the mesh accordingly. The adapting scheme determines the mapping from the original cell type
+    to the new cell type and specifies how to construct new nodes.
+
+    Parameters:
+        mesh: The input mesh to be adapted.
+        adapting_scheme:
+            A function that takes a cell type string and returns either:
+                - None, if the cell type should not be adapted,
+                - or a tuple of (new_cell_type, node_mapping), where:
+                    new_cell_type : str
+                        The name of the quadratic cell type (e.g., "triangle6").
+                    node_mapping : list
+                        A list specifying how to map old nodes and create new mid-edge nodes.
+                        Each entry is either a tuple with one index (existing node) or multiple
+                        indices (to create a new node as the mean of those nodes).
 
     Returns:
-        The modified mesh with quadratic elements.
-    """
+        The adapted mesh with quadratic cells and updated points, point data, and cell data.
 
+    Notes
+        - Cell blocks that are already quadratic or are not to be adapted are left unchanged.
+        - New mid-edge nodes are created only once per unique edge and reused across cells.
+        - Point data is interpolated for new nodes using the mean of the corresponding data.
+    """
     new_points = [coord for coord in mesh.points]
     new_point_data = {key: [d for d in data] for key, data in mesh.point_data.items()}
     new_cell_blocks = []
@@ -102,16 +164,17 @@ def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
         new_cells = []
         cell_type = cellblock.type
 
-        if cell_type in ["vertex", "hexahedron27", "tetra10", "triangle6", "quad9"]:
-            # this cell block is already quadratic, no need to convert
+        # determine the mapping of the nodes to the new quadratic celltype
+        scheme = adapting_scheme(cell_type)
+
+        if not scheme:
+            # don't adapt this block
             new_cell_blocks.append(cellblock)
             for key, data in mesh.cell_data.items():
                 new_cell_data[key].append(data[i])
             continue
 
-        # determine the mapping of the nodes to the new quadratic celltype
-        node_mapping = _lin_to_quad_nodes[cell_type]
-        new_cell_type = _lin_to_quad_cell_type[cell_type]
+        new_cell_type, node_mapping = scheme
 
         for cell in cellblock.data:
             cell_nodes = []
@@ -135,6 +198,7 @@ def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
                                 np.mean(
                                     [mesh.point_data[name][cell[n]] for n in node],
                                     axis=0,
+                                    dtype=mesh.point_data[name].dtype,
                                 )
                             )
 
@@ -159,3 +223,53 @@ def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
             key: [np.array(d) for d in data] for key, data in new_cell_data.items()
         },
     )
+
+
+def lin_to_quad(mesh: meshio.Mesh) -> meshio.Mesh:
+    """Convert linear elements to quadratic elements in a mesh.
+
+    This function returns a new mesh in which all linear elements (triangles, quadrilaterals, tetrahedra, and hexahedra)
+    are converted into their corresponding quadratic elements.
+
+    Args:
+        mesh:
+            The input mesh containing linear elements.
+
+    Returns:
+        The modified mesh with quadratic elements.
+    """
+
+    def lin_to_quad(cell_type: str) -> Optional[tuple[str, list[list[int]]]]:
+        """Returns the node mapping for quadratic elements."""
+        if cell_type in ["vertex", "hexahedron27", "tetra10", "triangle6", "quad9"]:
+            return None
+
+        return _lin_to_quad_cell_type[cell_type], _lin_to_quad_nodes[cell_type]
+
+    return adapt_mesh(mesh, lin_to_quad)
+
+
+def to_serendipity(mesh: meshio.Mesh) -> meshio.Mesh:
+    """Converts elements to serendipity elements in a mesh.
+
+    Args:
+        mesh:
+            The input mesh containing linear elements.
+
+    Returns:
+        The modified mesh with serendipity elements.
+    """
+
+    def _to_serendipity(cell_type: str) -> Optional[tuple[str, list[list[int]]]]:
+        """Returns the node mapping for serendipity elements."""
+        if cell_type in ["hexahedron20"]:
+            return None
+
+        if cell_type not in ["hexahedron"]:
+            raise ValueError(
+                f"Cell type {cell_type} not supported for serendipity conversion."
+            )
+
+        return "hexahedron20", _hex8_to_hex20
+
+    return adapt_mesh(mesh, _to_serendipity)