diff --git a/README.md b/README.md
index 3df8c00..0f99802 100644
--- a/README.md
+++ b/README.md
@@ -2,19 +2,95 @@
[](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`.
-
-
-
-
-
+**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.
+
+
+
+
+ 
+ 2D mesh cut by a circular level set
+ |
+
+ 
+ Hybrid mesh (quads+tris) cut by flower level set
+ |
+
+ 
+ Quad cases with disambiguation panels
+ |
+
+
+
+ 
+ Cut prism (wedge) demo
+ |
+
+ 
+ Cut pyramid demo
+ |
+
+
+ 
+ Cut tetrahedron demo
+ |
+
+
+
+ 
+ Tetrahedral mesh cut by a popcorn level set
+ |
+
+ 
+ Hexahedral mesh cut by popcorn level set
+ |
+
+ 
+ Hexahedral mesh cut by gyroid level set
+ |
+
+
# 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.
diff --git a/cpp/src/CMakeLists.txt b/cpp/src/CMakeLists.txt
index f33c25f..8885064 100644
--- a/cpp/src/CMakeLists.txt
+++ b/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
+ $
+ $)
+
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
diff --git a/cpp/src/cell_topology.h b/cpp/src/cell_topology.h
new file mode 100644
index 0000000..266cbff
--- /dev/null
+++ b/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
+#include
+#include
+
+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, 1> interval_edges = {{
+ {0, 1}
+}};
+
+/// Triangle: 3 edges
+/// VTK edge order: (0,1), (1,2), (2,0)
+inline constexpr std::array, 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, 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, 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, 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, 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, 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> 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, 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, 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, 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 hexahedron_face_sizes = {4, 4, 4, 4, 4, 4};
+inline constexpr std::array prism_face_sizes = {3, 3, 4, 4, 4};
+inline constexpr std::array 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
diff --git a/cpp/src/cell_types.h b/cpp/src/cell_types.h
index e345a2f..71c9450 100644
--- a/cpp/src/cell_types.h
+++ b/cpp/src/cell_types.h
@@ -218,8 +218,18 @@ namespace cutcells
cell_types[1]= type::quadrilateral;
}
break;}
- case type::quadrilateral: {throw std::invalid_argument("quadrilateral yet not implemented");
- break;}
+ case type::quadrilateral: {if(ls_part=="phi=0")
+ {
+ cell_types.resize(1);
+ cell_types[0]= type::interval;
+ }
+ else
+ {
+ cell_types.resize(2);
+ cell_types[0]= type::triangle;
+ cell_types[1]= type::quadrilateral;
+ }
+ break;}
case type::tetrahedron: {if(ls_part=="phi=0")
{
cell_types.resize(2);
@@ -233,8 +243,19 @@ namespace cutcells
cell_types[1]= type::prism;
}
break;}
- case type::hexahedron: {throw std::invalid_argument("hexahedron yet not implemented");
- break;}
+ case type::hexahedron: {if(ls_part=="phi=0")
+ {
+ cell_types.resize(2);
+ cell_types[0]= type::triangle;
+ cell_types[1]= type::quadrilateral;
+ }
+ else
+ {
+ cell_types.resize(2);
+ cell_types[0]= type::tetrahedron;
+ cell_types[1]= type::prism;
+ }
+ break;}
case type::prism: {throw std::invalid_argument("prism yet not implemented");
break;}
case type::pyramid: {throw std::invalid_argument("pyramid yet not implemented");
diff --git a/cpp/src/cut_cell.cpp b/cpp/src/cut_cell.cpp
index 007f1e8..01ba9de 100644
--- a/cpp/src/cut_cell.cpp
+++ b/cpp/src/cut_cell.cpp
@@ -7,6 +7,10 @@
#include "cut_cell.h"
#include "cut_tetrahedron.h"
#include "cut_triangle.h"
+#include "cut_quadrilateral.h"
+#include "cut_hexahedron.h"
+#include "cut_prism.h"
+#include "cut_pyramid.h"
#include "cut_interval.h"
#include "cell_flags.h"
#include "cell_subdivision.h"
@@ -16,6 +20,7 @@
#include
#include
#include