The Honeycomb User Guide

Honeycomb

Honeycomb aims to provide a safe, efficient and scalable implementation of combinatorial maps for meshing applications. More specifically, the goal is to converge towards a (or multiple) structure(s) adapted to algorithms exploiting GPUs and many-core architectures.

The current objective is to

• write a first implementation in Rust
• improve the structure without having to deal with data races and similar issues, thanks to the Rust's guarantees
• implement basic meshing algorithms to evaluate the viability of the implementation & improve our structure using Rust's framework to streamline the refactoring and parallelization process
• Benchmark and/or profile and/or parallelize our first algorithm, grisubal
• Ship a first stable version of the library (see this issue)
• Work on efficient parallelism

Core Requirements

• Rust stable release - The MSRV may not be the latest stable release, but we do not give any guarantees for older versions compatibility

Quickstart

Rust

You can add honeycomb as a dependency of your project by adding the following lines to its Cargo.toml:

# [dependencies]
honeycomb = { git = "https://github.com/LIHPC-Computational-Geometry/honeycomb", tag = "0.6.0"} # remove tag for master branch build

Alternatively, you can add the sub-crates that are currently published on crates.io:

# [dependencies]
honeycomb-core = "0.6.0"
honeycomb-kernels = "0.6.0"
honeycomb-render = "0.6.0"

Note that if you want to access the latest changes and documentation, you may have to specify a commit instead of a version, and use the GitHub Pages documentation instead of the one hosted on docs.rs.

Documentation

You can generate this book and the Rust documentation locally using respectively mdbook and cargo doc:

mdbook serve --open user-guide/

cargo +nightly doc --all --all-features --no-deps

Note that generating the doc using a stable toolchain is possible, the features just won't be documented as clearly.

Contributing

Contributions are welcome and accepted as pull requests on GitHub. Feel free to use issues to report bugs, missing documentation or suggest improvements of the project.

Note that a most of the code possess documentation, including private modules / items / sections. You can generate the complete documentation by using the instructions above and passing the option --document-private-items to cargo doc.

License

Licensed under either of

• Apache License, Version 2.0 (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or http://opensource.org/licenses/MIT)

at your preference.

The SPDX license identifier for this project is MIT OR Apache-2.0.

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

References

Combinatorial Maps

• Damiand and Lienhardt. Combinatorial Maps: Efficient Data Structures for Computer Graphics and Image Processing. Chapman&Hall/CRC, 2015.
  • Provides an in-depth presentation of the structure and its variants
  • Link
• The CGAL Project. CGAL User and Reference Manual. CGAL Editorial Board, 5.6.1 edition, 2024.
  • Provides concrete examples as well as code snippets of the CGAL implementation of the structure. The CGAL implementation uses a different approach than ours, & support N-dimensionnal map.
  • Link

Algorithms

• Staten, Noble, and Wilson. Constructing Tetrahedral Meshes No Matter How Ugly. SIAM, 2024
  • Describes the logic behind an overlay grid algorithm.
  • Link
• Rangarajan and Lew. Provably Robust Directional Vertex Relaxation for geometric mesh optimization. SIAM, 2017
  • usage TBD
  • Link

Integration

• The repository structure and workspace system is heavily inspired by the wgpu and bevy repositories.

Workspace

The project root is organized using Cargo workspaces at the moment. This may change when other languages are introduced to the project.

The repository hosts both published crates (usable content) as well as complementary content such as benchmarks, examples or this guide.

Published crates

• honeycomb Main crate
• honeycomb-core Core definitions and tools for combinatorial map implementation
• honeycomb-kernels Meshing kernel implementations using combinatorial maps
• honeycomb-render Visualization tool for combinatorial maps

Other content

• honeycomb-benches Benchmarks of the main map structures and methods
• honeycomb-examples Examples of usage of the project's features
• user guide Source files of the user guide

Published crates

Several crates of this project are published on the registry crates.io: the main crate, honeycomb (not yet published), as well as specialized crates honeycomb-core, honeycomb-kernels, and honeycomb-render.

honeycomb

honeycomb is the main crate provided to users and serve as the entrypoint for combinatorial map usage. It is exclusively made up of re-exports from the core, kernels and render crate to provide a clean, all-in-one dependency.

Usage

At the moment, the honeycomb name is not available on crates.io; this means that using this crate requires adding the dependency using the git repository:

# [dependencies]
honeycomb = { git = "https://github.com/LIHPC-Computational-Geometry/honeycomb" }

honeycomb-core

honeycomb-core is a Rust crate that provides basic structures and operations for combinatorial map manipulation. This includes map structures, methods implementation, type aliases and geometric modeling for mesh representation.

Implemented features

• a builder structure to handle map creation: CMapBuilder.
• 2D combinatorial maps, usable in concurrent contexts: CMap2. this includes:
  • all regular operations (sew, unsew, beta images, ...),
  • a custom embedding logic to associate vertices and attributes to darts.
• 2D orbit computation: Orbit2.
• abstractions over attributes, to allow arbitrary items binding to the map using the same embedding logic as vertices:
  • AttributeBind & AttributeUpdate traits,
  • AttrSparseVec as a predefined collection for attributes,
  • additional traits to describe custom collection structures.
• geometry primitives:
  • Vector2 / Vertex2,
  • Vector3 / Vertex3.

honeycomb-kernels

honeycomb-kernels is a Rust crate that provides implementations of meshing kernels using the core crate's combinatorial maps. These implementations have multiple purposes:

1. Writing code using n-maps from a user's perspective
2. Covering a wide range of operations, with routines that are more topology-heavy / geometry-heavy / balanced
3. Stressing the data structure, to identify its advantages and its pitfalls in a meshing context
4. Testing for more unwanted behaviors / bugs

Explanations provided in this guide focus on the overall workflow of algorithms; Implementation-specific details and hypothesis are listed in the documentation of the crate.

Implemented algorithms

• Directional Vertex Relaxation
• Grisubal
• Polygon triangulations
• Cell splits

honeycomb-render

honeycomb-render is a Rust crate that provides a simple visualization framework to allow the user to render their combinatorial map. It is designed to be used directly in the code by reading data through a reference to the map (as opposed to a binary that would read serialized data). This render tool can be used to debug algorithm results in a significantly easier way than reading internal data would.

Usage

Use the App structure to render a given combinatorial map. You may need to run the program in release mode to render large maps. All items used to build that tool are kept public to allow users to customize the render logic (e.g. to render a specific attribute).

Other content

Benchmarks

honeycomb-benches is the crate used to group benchmarking routines of the Rust code. It contains:

• binaries used to profile code and kernels
• benchmarks implemented using the criterion crate
• scritps used to aggregate and (partially) process results of both

Examples

honeycomb-examples is the crate used to group examples & snippets illustrating possible usages. It also contains sample VTK files to ensure examples and benchmarks can run out-of-the-box.

User guide

The user guide is the documentation you are currently reading right now. It is generated using mdbook. Its content mainly focuses on definition and feature-listing rather than technical details on implementation. The latter can be found in the code documentation.

Building

You can generate this documentation locally using mdbook:

mdbook serve --open user-guide/

Additional Information

A few observations on writing documentation using mdbook:

• If you edit the user guide's content, you will have to generate the rust doc again as mdbook remove all files of its target folder at each update.
• Linking to html files (and not markdown) has a varying level of success when working locally. Your browser may or may not like links toward folders instead of explicit index.html.

Benchmarks

Binaries

cargo run --profile=profiling --bin=<BINARY> -- <ARGS>

builder

Generate a grid using GridDescriptor w.r.t. the specified arguments. Arguments are the following: <N_CELL_X> <N_CELL_Y> <SPLIT>, with:

• <N_CELL_X>: number of cell of the grid along the X axis
• <N_CELL_Y>: number of cell of the grid along the Y axis
• <SPLIT>: if non-empty, the program will build a grid of tris instead of quads

grisubal

Run the grisubal algorithm w.r.t. the specified arguments. Arguments are the following: <INPUT_FILE> <LEN_CELL_X> <LEN_CELL_Y> <CLIP>, with:

• <INPUT_FILE>: input geometry, passed as a VTK file input
• <LEN_CELL_X>: length of the cells of the overlay grid along the X axis
• <LEN_CELL_X>: length of the cells of the overlay grid along the Y axis
• <CLIP>: left to clip the left side of the boundary, right the right side, anything else is a no-op

Benchmarks

cargo bench --bench <BENCHMARK>

The following benchmarks are available:

A detailed explanation about the purpose of each benchmark is provided at the beginning of their respective source files.

Scripts

Benchmarking

Both scripts provide an interactive menu with four options:

(0) all
(1) fixed-size profiling
(2) grid size scaling
(3) thread number scaling (not yet implmented)

builder.py

The script aggregates metrics about the grid building routines used by CMapBuilder. Data is collected from Criterion (runtime benchmarks), perf (CPU profiling), flamegraph (visualization), heaptrack (memory analysis), and across different grid sizes (from 128x128 to 8192x8192).

grisubal.py

The script aggregates metrics about the grisubal kernel. Data is collected from Criterion (runtime benchmarks), perf (CPU profiling), flamegraph (visualization), and heaptrack (memory analysis). Additionally, internal timers are implemented to measure time-per-sections of the algorithm. Measurements are done across different grid granularities (from 0.1 to 1.0 in 0.1 increments).

Plotting

grisubal_plot.py

This script generates a pie chart from per-section timings sampled by the grisubal.py script.

Examples

# Run a specific example
cargo run --example <EXAMPLE>

The following examples are available:

Note that the memory_usage example is deprecated and should be updated along with its plotting script.

Directional Vertices Relaxation

WIP: https://github.com/LIHPC-Computational-Geometry/honeycomb/pull/212

honeycomb-kernels - GRISUBAL

Grisubal is a mesh generation algorithm inspired by Morph. The mesh is built by capturing the input geometry in an overlapping grid, by first computing intersection vertices and then rebuilding boundaries from the captured vertices.

The algorithm can be called using this function.

Input

The algorithm expects a (2D) geometry specified via a path to a VTK legacy file. The geometry's boundaries should be described via segments & vertices; Segments should be consistently oriented.

Some vertices can be explicitly listed as cells in order for the algorithm to interpret those as Points of Interests. This can be used to ensure irregular geometries are captured correctly by the algorithm, which uses a uniform grid size.

Geometry Pre-processing

Before running the main algorithm, a few steps are followed to ensure correctness of later computations.

First, we compute the characteristics of a grid overlapping the entire geometry. The origin is automatically computed according to the input geometry and the desired cell sizes; It is then, if necessary, adjusted to avoid a few edge cases that would create issues in the main algorithm.

After grid characteristics are obtained, they're used along the geometry to detect and remove redundant Points of Interest. These are defined as PoI that land exactly on the grid; Their removal is necessary to avoid creating duplicate vertices in the main algorithm (since there would be both an intersection and a PoI at this location).

As a last step before calling the main kernel, we check for trivial orientation issues by ensuring no vertex is the start (resp. the end) of two different segments. This detects inconsistencies per-boundary, not overall consistency.

Grid Submersion Algorithm

Step 1 - Intersect Grid & Geometry

The goal of this step is to edit and complete the segment list to obtain a list of non-dividable segments that can be used for reconstruction at a later step. Consider the previous geometry, submerged in an overlapping grid:

We check each segment of the geometry for intersection(s) with the grid, and replace the original segment with new ones given the result of the check. A dedicated section goes over the method we use, these are the rough cases:

• Both vertices belong to the same cell: the new segment is the same as the original.
• Vertices belong to neighboring cells: there are two new segments, first from start vertex to intersection, second from intersection to end vertex.
• Vertices belong to different non-neighboring cells: there are d new segments, first from start vertex to intersection, then between intersections, last from intersection to end vertex.

Intersection information is collected and returned along with the list of non-dividable segments. The information is made up of a dart identifier (the intersected dart) as well as the relative position of the intersection on this dart (a floating-point t between 0 and 1).

At the same time, vertices are labeled as one of three types: Regular, PoI, or Intersec. This is used by the processing logic of the next steps.

Step 2 - Transform data for step 3 and 4

This is an intermediate step which enables better implementation of the next two step. It roughly results in:

• the parallelization of step 3
• the decorrelation of step 3 from step 4, making their concurrent execution posssible

This is achieved by:

1. grouping intersection data (obtained from step 1) per edge
2. pre-allocating darts for step 2
3. computing each intersection's corresponding new dart

Step 3 - Insert Intersection Vertices

Given intersection information per edge (step 2 1.), and pre-allocated darts (step 2 2.), we iterate through edges, building intersections into the map. Thanks to the splitn_edge implementation and dart pre-allocation, processing edges should be an embarassingly parallel section.

Step 4 - Filter & Rebuild Segment Data

Given information computed at step 2 3., we can rebuild new segments where both ends are either points of interest or intersections. This corresponds to building segments using the following vertices:

This can be done in two substeps:

1. Filter segments by starting vertex to only keep intersections;
2. For each of these segments, follow through until landing on another intersection; While searching through, keep track of any PoI encountered.

Using a set of data made up of starting intersection, ending intersection, and (optional) intermediates, we can build edges into the final 2-map.

Step 5 - Insert Segments

Given the data built up at the last step, we can proceed with insertion into the map. At this point, only darts linking the first intersection to the following vertex need to be added to the map.

The main work consist of fetching correct dart identifiers and update the topology by using link and sew methods. By following this process recursively for intermediates, we can build the final map, capturing the geometry's boundaries:

Clipping

Optionally, some cells of the resulting combinatorial map can be removed. These cell would correspond to the inside or the outside of the geometry. We choose to simply consider the left side and the right side of the boundary, to minimize orientation assumptions and avoid confusion.

During the last step of the main algorithm, darts of the boundary are marked according to their respective side. From this, we can retrieve faces of a given side, and use them as a starting point for a coloring-like algorithm.

Faces are searched and marked using a BFS; only adjacent faces with an unmarked dart are considered. If, at any point, a face with a dart of the other side of the boundary is reached it means that:

• the geometry was open, or
• nested boundaries are inconsistently oriented.

After this, all darts of the marked faces are deleted. The attribute used to mark boundary darts is then removed from the map before it is returned.

Appendices

Intersection Computation

Consider a given segment of the geometry. Each of the segment follow one of three cases:

• A - Vertices belong to neighboring cells.
• B - Vertices belong to different non-neighboring cells.
• C - Both vertices belong to the same cell.

We can identify which case we're in by computing the Manhattan distance between grid cells of the respective vertices:

• A - Distance is equal to 1.
• B - Distance is strictly superior to 1.
• C - Distance is 0.

The A case is pretty straight forward: By considering that the segment is oriented from the start to its end, we can deduce which side of the cell it's intersecting, and compute the relative position of the intersection t:

The B case is trickier to handle because we need to compute multiple intersections, and know their order to be able to build the segment back into the map.

By considering the minimal subgrid containing both vertices, as well as the direction of the segment, we can list all edges that are potentially intersected.

We can compute coefficients s and t for each of the potentially intersected segments. If the segment is actually intersected, both coefficients will have a value between 0 and 1. s being the relative position of the intersection along the original segment, we can use its value to reorder all valid intersections for segment building.

Insertion Logic

For a given instance of data containing:

• a starting dart: d_start
• (optional) one or more intermediate darts: { d_i } _i
• an ending dart: d_end

These steps are followed:

• 1-unsew d_start
• 0-unsew d_end, i.e. 1-unsew β₀ (d_end)
• create a pair of dart & link it via β₂
• 1-sew d_start to this pair, as well as the dart that was deconnected from d_start
• if there's no intermediate, 1-sew the pair to the ending dart, as well as the dart that was deconnected from d_end
• if there are intermediates:
  • 1-sew the pair to the first intermediate, as well as the dart that was deconnected from d_start
  • successively 1-sew intermediates
  • 1-sew the last intermediate to d_end, as well as the dart that was deconnected from d_end

Polygon triangulation

We implement two functions for polygon triangulation. These are not meshing functions; our goal with these is to cut existing cells of an irregular mesh into triangular cells.

With consideration to the above, we implement two polygon triangulation methods: fanning, and ear-clipping. Both implementations are designed to operate in parallel; they take pre-allocated darts as argument and do not create any contention over data as long as two calls are not made on the same face (2-cell).

Fanning

Two versions of this algorithm are implemented:

The first implementation is a defensive one where the function actively search for a valid vertex to fan from.

The second implementation assume the cell is convex; it fans the polygon from its first vertex. Convexity is not checked, so use this only if you know all your cells fit the requirements!

Ear-clipping

This method isn't algorithmically efficient, but we operate on small cells, and it covers our needs: it is a potential fallback for non-fannable polygons without holes.

Splits

We refer to "splits" or "splitting functions" when talking about higher-level routines involving the split of a geometrical cell (e.g. an edge/1-cell). They are implemented in this module.

No-allocation variants

To avoid small, repetitive call to system allocator, we implement no-allocation variants of these routines. Instead of adding darts to the map in the function, the variants take an additional argument: a slice of free darts. These darts are then used to apply modifications to the map.

This, coupled with grouped pre-allocation, allows to significantly reduce the number of allocator calls in programs.

Functions

Edges

Two splitting functions are implemented over edges, each having a no-allocation variant:

• split_edge / split_edge_noalloc: split an edge in two, inserting one new vertex in the map
• splitn_edge / splitn_edge_noalloc: split an edge into n+1 segments, n being the number of new intermediate vertices

Introduction

Scope

Because our implementation has very practical goals, the definitions presented here use a more intuitive approach rather than strict mathematical concepts.

Our main interest at the moment being 2D and 3D combinatorial maps, some operations are defined with consideration to these restrictions. This is especially useful for orbits and sewing operations.

Combinatorial Maps

N-dimensional combinatorial maps, noted N-map, are objects made up of two main elements:

• A set of darts, darts being the smallest elements making up the map
• N beta functions, linking the darts of the map

Additionally, we can define embedded data as spatial anchoring of the darts making up the map. While the first two elements hold topological information, embedded data gives information about the "shape" of the map (e.g. vertices position in a spatial domain).

With these elements, we can represent and operate on meshes.

Example

Operations on a combinatorial map can affect its topology, shape or both:

The specifics on how data is encoded is detailed in attribute-specific sections.

Darts

Darts are the finest grain composing combinatorial maps. The structure of the map is given by the relationship between darts, defined through beta functions. Additionally, a null dart is defined, we denote it ∅.

In our implementation, darts exist implicitly through indexing and their associated data. There are no dart objects in a strict sense, there is only a given number of dart, their associated data ordered by an array-like logic, and a record of "unused" slots that can be used for dart insertion. Because of this, we assimilate dart and dart index.

Beta Functions

Each combinatorial map of dimension N defines N beta functions linking the set of darts together (e.g. a 2-map contains β₁ and β₂). These functions model the topology of the map, giving information about connections of the different cells of the map / mesh. In our case, we mostly use:

• β₁, a (partial) permutation,
• β₂, β₃, two (partial) involutions

Additionally, we define β₀ as the inverse of β₁, i.e. β₀(β₁(d)) = d. This comes from a practical consideration for performances and efficiency of the implementation.

The β_i functions can be interpreted as navigation functions along the i-th dimensions: β₁ makes you navigate along the edges, β₂ along the faces, etc. This can be generalized to N dimensions, but we are only interested in 2D and 3D at the moment.

Properties

For a given dart d, we define two properties:

• d is i-free if β_i(d) = ∅, ∅ being the null dart
• d is free if it is i-free for all i

Construction Example

Orbits

We define orbits as a set of darts that are accessible from a given dart, using a certain set of beta functions. For example:

• ⟨β₁⟩(d) refers to all darts accessible from d using β₁ recursively any number of times.
• ⟨β₁, β₃⟩(d) refers to all darts accessible from d using any combination of β₁ and β₃.

i-cells

A specific subset of orbits, referred to as i-cells are defined and often used in algorithms. The general definition is the following:

• if i = 0: 0-cell(d) = ⟨{ β_j o β_k with 1 ≤ j < k ≤ N }⟩(d)
• else: i-cell(d) = ⟨β₁, β₂, ..., β_i-1, β_i+1, ..., β_N⟩(d)

In our case, we can use specialized definitions for our dimensions:

Examples

Embedding

Embedding, or embedded data refers to the association of topological entities (darts, i-cells) to geometrical data (spatial positions, vertices, faces, volumes).

The embedding of geometrical data has implication for operations on the map. This is detailed along operation specificities in their dedicated sections.

Implementation

Because embedding is defined per-application or per-needs, our combinatorial map implementation uses a generic system to handle this; The storage follows our own specific logic, which is detailed in this documentation entry.

Examples

Note that darts and vertex IDs are not actually stored. Formers exist implictly and latters are computed on the fly.

In the above example, data would be organized like this:

In the above example, data would be organized like this:

Basic Operations

A number of basic operations are defined on combinatorial maps. In this section, we go over their definition, a simple example and eventual implementation specificities. Qualifies as a basic operation:

• changing the dimension of the map
• adding and removing free darts from the map
• sewing and unsewing darts

Other operations might be defined, but these are the useful core for writing higher-level abstractions over the structure.

Add / Remove a dimension

Adding or removing a dimension on a given combinatorial maps effectively corresponds, respectively, to adding or removing a beta function. In the case of decreasing the dimension, this operation can result in two disjoint dart set in the same map.

Because the current implementation only covers 2D combinatiorial maps, this operation is not implemented. When 3D maps are implemented, it would be possible to implement this using the From trait provided by the Rust language.

Add / Remove a free dart

Adding darts

As explained in the Darts section of this guide, these entities exist implictly through indexing of the internal storage structures of the map. Because of this, adding darts translates to extending internal vectors and storages in our implementation.

An internal counter is incremented at each dart addition. This, coupled with an unused dart tracking mechanism, constitutes a way to keep track of attributed darts.

Removing darts

Removing a dart would technically require us to remove an entry inside storage structures, which are often ordered, contiguous vectors. There are two way to approach this problem:

• Actually remove the entry
  • requires adjustments on all the structure to keep consistent indices
  • keeps the storage compact, i.e. all allocated slots are used
• "Forget" the entry
  • does not require any re-arrangements besides making sure no beta functions lands on the dart
  • creates "holes" in the storage

Our implementation uses the second solution, along with a structure used to store unused slots. In turns, we can use these "holes" in the storage to reinsert darts or collapse the structure at a later point during execution.

Sewing and unsewing operation

Sewing

The sew operation can be divided into two parts:

• a topological update, which corresponds to linking the darts
• a geometrical update, which corresponds to an update of the affected embedded data, called attributes in our code

Note that the implementation is not as simple as doing one and then the other for consistency reasons: changing the topology affects our ability to retrieve the embedded data, therefore the result is highly sensitive to operation order.

Topology

The i-link operation corresponds to the aforementionned topological update. Given two darts d_a and d_b, and a given beta function β_i, a link operation corresponds to the update of the β_i function in order to have β_i(d_a) = d_b and/or β_i(d_b) = d_a depending on darts order and i. For example:

• 1-link(d_a,d_b) results in:
  • β₁(d_a) = d_b
  • if β₀ is defined, β₀(d_b) = d_a
• 1-link(d_b,d_a) results in:
  • β₁(d_b) = d_a
  • if β₀ is defined, β₀(d_a) = d_b
• 2-link(d_a,d_b) results in:
  • β₂(d_a) = d_b
  • β₂(d_b) = d_a
• 2-link(d_b,d_a) results in the same changes as 2-link(d_a,d_b)

Exact properties of the link operation directly depends on the property of the modified beta function.

Geometry

The i-sew operation corresponds to an i-link operation, coupled with an update of the affected attributes. How the attributes are updated is defined through trait implementation in the Rust crate (see AttributeUpdate, AttributeBind). Which attributes are updated can be deduced from the dimension i of the sewing operation. This is summarized in the following table:

Unsewing

The unsew operation is the complementary to the sew operation. It behaves according to similar properties, but is used to remove links between darts. It does so by replacing values of the beta functions by the null dart. Geometrical updates are handled and defined in the same way as for the sew operation.