Tablegen

README.md

@@ -2,19 +2,95 @@
 
 [![CutCells CI](https://github.com/sclaus2/cutcells/actions/workflows/python-app.yml/badge.svg)](https://github.com/sclaus2/cutcells/actions/workflows/python-app.yml)
 
-a library to compute intersections between mesh cells and implicit (level set) functions using marching cubes/tetrahedra algorithms.
 
-The current version supports intervals, triangles and tetrahedra. Future versions will include quadrilaterals, hexahedra, pyramids and prisms.
-
-The library contains utility functions to cut through vtk meshes using pyvista. It is also used as a basis to generate run time quadrature rules in immersed boundary methods/cut finite element methods (see CutFEMx).
-
-The images below show a cut through a vtk mesh using pyvista and CutCells.
-The corresponding demonstration code is located in `python/demo/cut_vtk_mesh`.
-
-<p float="left">
-<img src="./img/mesh2D.png" alt="Regular mesh cut by zero contour line of a circular level set function" width="500"/>
-<img src="./img/popcorn.png" alt="Regular mesh cut by zero contour line of a popcorn level set function" width="500"/>
-</p>
+**CutCells** is a lightweight **C++20 geometry kernel with Python bindings** to compute **cut-cell decompositions** of standard mesh elements by an **implicit surface / level set** `φ(x)`.
+
+Given one parent cell (or a whole mesh) and the level set values on its vertices, CutCells builds a **local sub-mesh** that describes:
+
+* the **inside** region `φ < 0`,
+* the **outside** region `φ > 0`,
+* and/or the **interface** `φ = 0` (as facets/segments embedded in the cell),
+
+with **explicit connectivity** and **element types**.
+
+It is designed as a building block for **unfitted / immersed FEM** workflows (e.g. **CutFEM**), where robust cut geometry and stable parent mappings are required for runtime quadrature.
+
+## Supported cell types
+
+Current parent cell support:
+
+* 1D: **interval**
+* 2D: **triangle**, **quadrilateral**
+* 3D: **tetrahedron**, **hexahedron**, **prism**, **pyramid**
+
+Higher-order (currently pragmatic) support:
+
+* **P2 triangle (6 nodes)** and **P2 tetrahedron (10 nodes)** via subdivision into linear sub-cells, cut, then merged back.
+
+## Performance and dependencies
+
+* **Few dependencies**: the computational core is plain C++.
+* **Python bindings via nanobind**: CutCells exposes arrays in NumPy-friendly form and is designed to avoid unnecessary overhead.
+* **Zero-copy where possible**: large buffers (coords, connectivity/offsets, parent IDs, provenance tags) are exposed as **views** when layout permits, minimizing memory traffic between C++ and Python.
+
+## Applications
+
+CutCells is intended as a reusable geometry backend for:
+
+* **CutFEM / unfitted FEM**: robust sub-cell and interface extraction for integration on `Ω∩K` and `Γ∩K`
+* **runtime quadrature** generation 
+* general embedded geometry workflows where remeshing is undesirable
+
+## Example gallery 
+
+
+Below is a gallery of example outputs generated by CutCells and the provided Python demos. See `python/demo/` for scripts to reproduce these images.
+
+<table>
+	<tr>
+		<td align="center">
+			<img src="./img/mesh2D.png" alt="2D mesh cut by circle" width="320"/><br/>
+			<sub>2D mesh cut by a circular level set</sub>
+		</td>
+        <td align="center">
+			<img src="./img/hybrid_flower.png" alt="Hybrid mesh flower cut" width="320"/><br/>
+			<sub>Hybrid mesh (quads+tris) cut by flower level set</sub>
+		</td>
+		<td align="center">
+			<img src="./img/quad_cases_with_disambiguation_panels_quads.png" alt="Quad cases with disambiguation" width="320"/><br/>
+			<sub>Quad cases with disambiguation panels</sub>
+		</td>
+	</tr>
+	<tr>
+		<td align="center">
+			<img src="./img/prism_n0_demo.png" alt="Cut prism demo" width="320"/><br/>
+			<sub>Cut prism (wedge) demo</sub>
+		</td>
+		<td align="center">
+			<img src="./img/pyramid_n0_demo.png" alt="Cut pyramid demo" width="320"/><br/>
+			<sub>Cut pyramid demo</sub>
+		</td>
+        <td align="center">
+			<img src="./img/cut_tetra_quad.png" alt="Cut pyramid demo" width="180"/>
+            <img src="./img/cut_tetra_tri.png" alt="Cut pyramid demo" width="180"/><br/>
+			<sub>Cut tetrahedron demo</sub>
+		</td>
+	</tr>
+	<tr>
+		<td align="center">
+			<img src="./img/popcorn.png" alt="2D mesh cut by popcorn function" width="320"/><br/>
+			<sub>Tetrahedral mesh cut by a popcorn level set</sub>
+		</td>
+		<td align="center">
+			<img src="./img/popcorn_hex_mesh3D.png" alt="3D hex mesh cut by popcorn" width="320"/><br/>
+			<sub>Hexahedral mesh cut by popcorn level set</sub>
+		</td>
+        		<td align="center">
+			<img src="./img/cut_hex_gyroid_grid.png" alt="Hexahedron gyroid grid" width="320"/><br/>
+			<sub>Hexahedral mesh cut by gyroid level set</sub>
+		</td>
+	</tr>
+</table>
 
 # Installation
 
@@ -64,6 +140,11 @@ in the corresponding demo folder, e.g. `cpp/demo/cut_triangle`
 
 The python demos are located in `python/demo`.
 
+## Third-party notices
+
+Some generated clip/cut case tables (hexahedron/prism/pyramid) are derived from VTK's `vtkTableBasedClipCases.h`.
+VTK's BSD-3-Clause license text is included in `third_party/VTK-Copyright.txt`.
+
 ## Dependencies
 
 CutCells requires a C++20 compiler and depends on the C++ standard template library.

cpp/src/CMakeLists.txt

@@ -8,12 +8,22 @@ add_library(cutcells)
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/version.h.in version.h)
 include_directories(${CMAKE_CURRENT_BINARY_DIR})
 
+# Make headers in this directory discoverable for generated table headers
+# (e.g. generated/*_tables.h includes "cell_types.h").
+target_include_directories(cutcells PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>)
+
 set(HEADERS
   ${CMAKE_CURRENT_SOURCE_DIR}/cell_types.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cell_flags.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cell_subdivision.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_triangle.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_interval.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_quadrilateral.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_hexahedron.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_prism.h
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_pyramid.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_tetrahedron.h
   ${CMAKE_CURRENT_SOURCE_DIR}/triangulation.h
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_cell.h
@@ -27,6 +37,10 @@ set(HEADERS
 target_sources(cutcells PRIVATE
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_interval.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_triangle.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_quadrilateral.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_hexahedron.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_prism.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/cut_pyramid.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_tetrahedron.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_cell.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/cut_mesh.cpp

cpp/src/cell_topology.h

@@ -0,0 +1,230 @@
+// Copyright (c) 2022 ONERA
+// Authors: Susanne Claus
+// This file is part of CutCells
+//
+// SPDX-License-Identifier:    MIT
+#pragma once
+
+#include "cell_types.h"
+#include <array>
+#include <span>
+#include <stdexcept>
+
+namespace cutcells::cell
+{
+/// @brief Cell topology data for edge-based clipping algorithms.
+/// Edge and face numbering follows VTK conventions.
+
+//-----------------------------------------------------------------------------
+// Edge definitions per cell type (VTK ordering)
+//-----------------------------------------------------------------------------
+
+/// Interval: 1 edge (0-1)
+inline constexpr std::array<std::array<int, 2>, 1> interval_edges = {{
+    {0, 1}
+}};
+
+/// Triangle: 3 edges
+/// VTK edge order: (0,1), (1,2), (2,0)
+inline constexpr std::array<std::array<int, 2>, 3> triangle_edges = {{
+    {0, 1}, {1, 2}, {2, 0}
+}};
+
+/// Quadrilateral: 4 edges
+/// VTK edge order: (0,1), (1,2), (2,3), (3,0)
+///
+///   3 --- 2
+///   |     |
+///   0 --- 1
+inline constexpr std::array<std::array<int, 2>, 4> quadrilateral_edges = {{
+    {0, 1}, {1, 2}, {2, 3}, {3, 0}
+}};
+
+/// Tetrahedron: 6 edges
+/// VTK edge order: (0,1), (1,2), (2,0), (0,3), (1,3), (2,3)
+inline constexpr std::array<std::array<int, 2>, 6> tetrahedron_edges = {{
+    {0, 1}, {1, 2}, {2, 0}, {0, 3}, {1, 3}, {2, 3}
+}};
+
+/// Hexahedron: 12 edges
+/// VTK vertex numbering:
+///        7 -------- 6
+///       /|         /|
+///      / |        / |
+///     4 -------- 5  |
+///     |  3 ------|-- 2
+///     | /        | /
+///     |/         |/
+///     0 -------- 1
+///
+/// VTK edge order:
+///   Base: (0,1), (1,2), (2,3), (3,0)
+///   Top:  (4,5), (5,6), (6,7), (7,4)
+///   Verticals: (0,4), (1,5), (2,6), (3,7)
+inline constexpr std::array<std::array<int, 2>, 12> hexahedron_edges = {{
+    // Base edges (0-3)
+    {0, 1}, {1, 2}, {2, 3}, {3, 0},
+    // Top edges (4-7)
+    {4, 5}, {5, 6}, {6, 7}, {7, 4},
+    // Vertical edges (8-11)
+    {0, 4}, {1, 5}, {2, 6}, {3, 7}
+}};
+
+/// Prism (Wedge): 9 edges
+/// VTK vertex numbering:
+///       2
+///      /|\
+///     / | \
+///    0-----1
+///    |  |  |
+///    |  5  |
+///    | /|\ |
+///    |/ | \|
+///    3-----4
+///
+/// VTK edge order:
+///   Bottom triangle: (0,1), (1,2), (2,0)
+///   Top triangle: (3,4), (4,5), (5,3)
+///   Verticals: (0,3), (1,4), (2,5)
+inline constexpr std::array<std::array<int, 2>, 9> prism_edges = {{
+    // Bottom triangle (0-2)
+    {0, 1}, {1, 2}, {2, 0},
+    // Top triangle (3-5)
+    {3, 4}, {4, 5}, {5, 3},
+    // Vertical edges (6-8)
+    {0, 3}, {1, 4}, {2, 5}
+}};
+
+/// Pyramid: 8 edges
+/// VTK vertex numbering:
+///         4
+///        /|\
+///       / | \
+///      /  |  \
+///     /   |   \
+///    3----|----2
+///    |    |    |
+///    0---------1
+///
+/// VTK edge order:
+///   Base: (0,1), (1,2), (2,3), (3,0)
+///   Apex connections: (0,4), (1,4), (2,4), (3,4)
+inline constexpr std::array<std::array<int, 2>, 8> pyramid_edges = {{
+    // Base edges (0-3)
+    {0, 1}, {1, 2}, {2, 3}, {3, 0},
+    // Apex edges (4-7)
+    {0, 4}, {1, 4}, {2, 4}, {3, 4}
+}};
+
+//-----------------------------------------------------------------------------
+// Edge count per cell type
+//-----------------------------------------------------------------------------
+
+inline constexpr int num_edges(type cell_type)
+{
+    switch (cell_type)
+    {
+    case type::point:         return 0;
+    case type::interval:      return 1;
+    case type::triangle:      return 3;
+    case type::quadrilateral: return 4;
+    case type::tetrahedron:   return 6;
+    case type::hexahedron:    return 12;
+    case type::prism:         return 9;
+    case type::pyramid:       return 8;
+    default:
+        throw std::invalid_argument("Unknown cell type in num_edges");
+    }
+}
+
+//-----------------------------------------------------------------------------
+// Edge accessor (returns span to edge array for cell type)
+//-----------------------------------------------------------------------------
+
+/// Get edges for a given cell type.
+/// @param cell_type The cell type
+/// @return Span of edge vertex pairs [v0, v1] for each edge
+inline std::span<const std::array<int, 2>> edges(type cell_type)
+{
+    switch (cell_type)
+    {
+    case type::interval:
+        return std::span(interval_edges);
+    case type::triangle:
+        return std::span(triangle_edges);
+    case type::quadrilateral:
+        return std::span(quadrilateral_edges);
+    case type::tetrahedron:
+        return std::span(tetrahedron_edges);
+    case type::hexahedron:
+        return std::span(hexahedron_edges);
+    case type::prism:
+        return std::span(prism_edges);
+    case type::pyramid:
+        return std::span(pyramid_edges);
+    default:
+        throw std::invalid_argument("Unknown cell type in edges()");
+    }
+}
+
+//-----------------------------------------------------------------------------
+// Face definitions (for interface stitching, optional future use)
+//-----------------------------------------------------------------------------
+
+/// Hexahedron: 6 faces (quads)
+/// Face order: -X, +X, -Y, +Y, -Z, +Z (left, right, front, back, bottom, top)
+inline constexpr std::array<std::array<int, 4>, 6> hexahedron_faces = {{
+    {0, 3, 7, 4},  // -X (left)
+    {1, 2, 6, 5},  // +X (right)
+    {0, 1, 5, 4},  // -Y (front)
+    {3, 2, 6, 7},  // +Y (back)
+    {0, 1, 2, 3},  // -Z (bottom)
+    {4, 5, 6, 7}   // +Z (top)
+}};
+
+/// Prism: 5 faces (2 triangles + 3 quads)
+/// Represented as max 4 vertices per face, with -1 padding for triangles
+inline constexpr std::array<std::array<int, 4>, 5> prism_faces = {{
+    {0, 1, 2, -1},  // bottom triangle
+    {3, 4, 5, -1},  // top triangle
+    {0, 1, 4, 3},   // quad face
+    {1, 2, 5, 4},   // quad face
+    {2, 0, 3, 5}    // quad face
+}};
+
+/// Pyramid: 5 faces (1 quad base + 4 triangles)
+inline constexpr std::array<std::array<int, 4>, 5> pyramid_faces = {{
+    {0, 1, 2, 3},   // base quad
+    {0, 1, 4, -1},  // triangle
+    {1, 2, 4, -1},  // triangle
+    {2, 3, 4, -1},  // triangle
+    {3, 0, 4, -1}   // triangle
+}};
+
+//-----------------------------------------------------------------------------
+// Face vertex count per face (for mixed cells)
+//-----------------------------------------------------------------------------
+
+inline constexpr std::array<int, 6> hexahedron_face_sizes = {4, 4, 4, 4, 4, 4};
+inline constexpr std::array<int, 5> prism_face_sizes = {3, 3, 4, 4, 4};
+inline constexpr std::array<int, 5> pyramid_face_sizes = {4, 3, 3, 3, 3};
+
+/// Get number of faces for a cell type
+inline constexpr int num_faces(type cell_type)
+{
+    switch (cell_type)
+    {
+    case type::point:         return 0;
+    case type::interval:      return 0;
+    case type::triangle:      return 1;
+    case type::quadrilateral: return 1;
+    case type::tetrahedron:   return 4;
+    case type::hexahedron:    return 6;
+    case type::prism:         return 5;
+    case type::pyramid:       return 5;
+    default:
+        throw std::invalid_argument("Unknown cell type in num_faces");
+    }
+}
+
+} // namespace cutcells::cell