The Honeycomb Book

Statement of need

Honeycomb aims to provide a safe, efficient and scalable implementation of combinatorial maps for meshing applications.

The goal is to converge towards a (or multiple) structure(s) which could be used to easily experiment with parallel meshing algorithm, specifically targeting many-core architectures.

More extensive explanations regarding our needs and design choices of this solution is included in one of our paper.

Requirements

Core

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

Optional

• CUDA
  • Used in one of the application, as a PoC for GPU usage
  • nvcc and libcudart may be sufficient, but we recommend installing the full toolkit
• hwloc
  • The library is used by the application binaries to bind threads to physical cores;
  • Consider this a core requirement if you are running benchmarks
  • Disable usage by compiling application binaries with the --no-default-features option

Quickstart

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.10.2" # it is highly encouraged to pin version using a tag or a revision
}

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

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

Note that the documentation hosted on GitHub corresponds to the master branch. Versioned documentation is available on docs.rs.

Combinatorial maps

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

Definition

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.

Components

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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

Cell representation

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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 implicitly 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:

Element insertion and deletion

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

Dart insertion

As explained in the Darts section of this guide, these entities exist implicitly 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.

Dart deletion

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.

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.

Sewing operation

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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

Topology

The i-link operation corresponds to the aforementioned 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.

Transactional memory

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

Needs

Rust’s ownership semantics require us to add synchronization mechanism to our structure if we want to use it in concurrent contexts. Using primitives such as atomics and mutexes would be enough to get programs to compile, but it would respectively yield an incorrect or impractical implementation:

• Atomics give guarantees on instructions interleaving for a single given variable, they do not give any guarantees for instructions affecting different atomic variables.
• Mutexes (and similar locks, e.g. RWLocks) can be used to create guarantees when accessing multiple variable: for example, we can write an operation that does not progress until all of the used data is locked. However, locks are error-prone, have very poor composability.

The nature of meshing operations makes both mechanisms very unpractical. They are complex, access many variables, and are often comprised of multiple steps. For example, the following operation is executed on all affected attributes of a sew:

Because the map can go through invalid intermediate states during a single operation, we need to ensure another thread will not use one of these as the starting point for another operation. This rules out fine-grained atomics.

The sew operation is the main method used to create new connectivities in the map. This means that most high-level meshing operations will call this method multiple times. If these meshing operations require locking all of the variables to ensure correct execution, the locks must be returned or exposed to the user so that he can unlock them at the right time. Manual lock management is error-prone, and becomes impossible in practice for complex meshing operations.

Software Transactional Memory

We choose to use Software Transactional Memory (STM) to handle high-level synchronization of the structure. Unlike locks, STM has great composability and allows users of the crate to easily define pseudo-atomic segments in their own algorithms.

Exposing an API that allows users to handle synchronization also means that the implementation isn’t bound to a given parallelization framework. Instead of relying on predefinite parallel routines (e.g. a provided parallel_for on given cells), the structure can be used to implement existing algorithms regardless of their approach (data-oriented, task-based, …).

Contributing

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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.

Environment

The repository contains a Nix flake to easily setup a development environment:

nix develop

Most notably, it handles hwloc install on both MacOs and Linux, as well as Linux-specific dependencies for our visualizer (used by bevy).

Checks

Nix

The flake also defines checks. They are identical to those of the CI, so use this rather than the pre-commit if possible.

nix flake check

Pre-commit hook

The repository contains a pre-commit hook config file. To use it:

pip install pre-commit # or whichever package manager
pre-commit install
pre-commit run # test it!

While it is not identical to the CI (most notably, it excludes honeycomb-render due to compile time), it is fine for core and kernel crates development.

The hook can be bypassed by using the --no-verify option to git commit.

Documentation

Note that a most of the code possess documentation, including private modules / items / sections. You can generate the complete documentation by using the following instructions:

mdbook serve --open user-guide/

cargo +nightly doc --all --all-features --no-deps --document-private-items

Project structure

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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.

The following libraries are available:

• honeycomb Main crate, which re-exports items from the three subcrates below
• 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

The repository also hosts:

• The applications crate, which contains a collection of algorithms which are used as benchmarks and/or examples
• This book’s source files, available in the user-guide directory

Libraries

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

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.

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. The following features are available:

• a builder structure to handle map creation: CMapBuilder.
• 2D and 3D combinatorial maps, usable in concurrent contexts: CMap2/CMap3. this includes:
  • all regular operations (sew, unsew, beta images, …),
  • a custom embedding logic to associate vertices and attributes to darts.
• 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.

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 quickly debug algorithm results by looking at the resulting structure instead of reading hard-to-interpret numerical data.

Use the exported functions 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).

(De)Serialization

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

Our crate implements two different serialization logics. We use a custom format for combinatorial map representation, while we use the VTK format for mesh representation. The latter can also be hijacked to initialize maps by restraining the range of supported data.

Serialization is available directly as map structures methods, while deserialization is available through the usage of the builder structure.

Usage

#![allow(unused)]
fn main() {
use std::{fs::File, io::Write};

use honeycomb_core::cmap::{CMap2, CMapBuilder};

// Init from a VTK file; only implemented for 2D map
let map: CMap2<f32> = match CMapBuilder::<2>::from_vtk_file("path/to/file.vtk").build() {
    Ok(cmap) => cmap,
    Err(e) => panic!("Error while building map: {e:?}"),
};
// Init from serialized data; implemented for 2D and 3
let map: CMap2<f32> = match CMapBuilder::<2>::from_cmap_file("path/to/file.cmap").build() {
    Ok(cmap) => cmap,
    Err(e) => panic!("Error while building map: {e:?}"),
};

// Save to VTK file; only implemented for 2D map
let file = File::create_new("out.vtk").unwrap();
map.to_vtk_binary(file);
// Serialize map data
let mut file = std::fs::File::create("out.cmap").unwrap();
let mut out = String::new();
map.serialize(&mut out);
file.write_all(out.as_bytes()).unwrap();

}

Custom serialization

Format specification

The file should contain 4 sections:

• [META]: header with miscellaneous data to help parsing / checking
• [BETAS]: values of beta functions
• [UNUSED]: unused darts
• [VERTICES]: vertex identifiers and values

Single line comments are supported using the # character.

META

Single line section specifying:

• format version,
• map dimension,
• number of darts.

The format version is the same as the crate’s version; it serves as a hint to use the correct crate version with a file you did not generate.

BETAS

Values of beta functions organized as one line per beta functions. The number of line should be equal to dimension plus one.

All lines should:

• have the same length
• have a length equal the number of darts specified in the header plus one (for the null dart)

Values on a single line are separated by a space, and the first value of each line should be 0 as they corresponds to beta images of the null dart.

UNUSED (optional)

Single line, optional section listing all unused darts. Unused darts should be free.

VERTICES

List of identifiers and corresponding vertex values. Vertices should have correct dimension (e.g. x, y, and z coordinates for a 3-map).

Examples

2D

# unit square
[META] 
0.8.0 2 4 # <VERSION> <DIM> <N_DARTS_EXCL_0>

[BETAS] # 3 beta functions for 2D, 5 columns = 4 darts + the null dart
0 4 1 2 3 # b0 
0 1 3 4 1 # b1
0 0 0 0 0 # b2

[VERTICES]
1 0.0 0.0 # <ID> <X> <Y>
2 1.0 0.0
3 1.0 1.0
4 0.0 1.0

3D

# simple tetrahedron
[META]
0.8.0 3 14

[BETAS] 
0 3 1 2  6 4  5 9 7 8  12 10 11 0 0 # columns don't have to be aligned,
0 2 3 1  5 6  4 8 9 7  11 12 10 0 0 # though our routines will align items
0 4 7 10 1 12 8 2 6 11 3  9  5  0 0
0 0 0 0  0 0  0 0 0 0  0  0  0  0 0

[UNUSED]
13 14

[VERTICES]
1 0.0 0.0 0.0 # <ID> <X> <Y> <Z>
2 1.0 0.0 0.0
3 1.0 1.0 0.0
6 0.5 0.5 1.0

VTK serialization

We use the vtkio crate to handle file IO. Only the legacy format is supported, in both its binary or ASCII form.

Expected input for deserialization

Using a VTK file to initialize a map can fail for two main reason:

• The file contains general inconsistencies:
  • the number of coordinates cannot be divided by 3, meaning a tuple is incomplete
  • the number of Cells and CellTypes isn’t equal,
  • a given cell has an inconsistent number of vertices with its specified cell type.
• The file contains unsupported data:
  • file format isn’t Legacy,
  • data set is something other than UnstructuredGrid,
  • coordinate representation type isn’t float or double,
  • mesh contains unsupported cell types (PolyVertex, PolyLine, …, or anything 3D for a 2D map for example).

Generic attribute system

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

The attributes module of the core crate provides the necessary tools for to add custom attributes to given orbits of the map. Each attribute should be uniquely typed (i.e. to type aliases) as the maps’ internal storages use std::any::TypeId for identification.

An attribute struct should implement both AttributeBind and AttributeUpdate. It can then be added to the map using the dedicated CMapBuilder method. This is showcased in n example below, where we add

Implementation example

#![allow(unused)]
fn main() {
use honeycomb_core::{
    attributes::{AttrSparseVec, AttributeBind, AttributeError, AttributeUpdate},
    cmap::{OrbitPolicy, VertexIdType},
};

#[derive(Debug, Clone, Copy, Default, PartialEq)]
struct Weight(pub u32);

impl AttributeUpdate for Weight {
    // when merging two weights, we add them
    fn merge(attr1: Self, attr2: Self) -> Result<Self, AttributeError> {
        Ok(Self(attr1.0 + attr2.0))
    }

    // when splitting, we do an approximate 50/50
    fn split(attr: Self) -> Result<(Self, Self), AttributeError> {
        // adding the % to keep things conservative
        Ok((Self(attr.0 / 2 + attr.0 % 2), Self(attr.0 / 2)))
    }

    // if we have to merge from a single value, we assume the "other" is 0
    fn merge_incomplete(attr: Self) -> Result<Self, AttributeError> {
        Ok(attr)
    }
}

impl AttributeBind for Weight {
    // Weight values will be stored in an `AttrSparseVec`
    type StorageType = AttrSparseVec<Self>;
    // Weights bind to vertices
    type IdentifierType = VertexIdType;
    const BIND_POLICY: OrbitPolicy = OrbitPolicy::Vertex;
}
}

Usage example

use honeycomb_core::cmap::{CMap2, CMapBuilder};

fn main() {
    let map: CMap2<_> = CMapBuilder::<2, f64>::from_n_darts(4)
        .add_attribute::<Weight>()
        .build()
        .unwrap();

    let _ = map.link::<2>(1, 2);
    let _ = map.link::<2>(3, 4);
    map.write_vertex(2, (0.0, 1.0));
    map.write_vertex(3, (1.0, 1.0));
    map.write_attribute::<Weight>(2, Weight(5));
    map.write_attribute::<Weight>(3, Weight(6));

    let _ = map.sew::<1>(1, 3);

    assert_eq!(map.read_attribute::<Weight>(2), Some(Weight(11)));

    let _ = map.unsew::<1>(1);

    assert_eq!(map.read_attribute::<Weight>(2), Some(Weight(6)));
    assert_eq!(map.read_attribute::<Weight>(3), Some(Weight(5)));
}

Parallel Laplace smoothing

This content has been copy-pasted from the previous guide. It is up-to-date but should be improved at some point.

In the following routine, we shift each vertex that’s not on a boundary to the average of its neighbors positions. In this case, transactions allow us to ensure we won’t compute a new position from a value that has been replaced since the start of the computation.

Code

use honeycomb_core::{
    cmap::{CMap2, DartIdType, OrbitPolicy, VertexIdType, NULL_DART_ID},
    geometry::Vertex2,
    stm::atomically,
};
use honeycomb_kernels::grid_generation::GridBuilder;
use rayon::prelude::*;

const N_SQUARES: usize = 256;
const N_ROUNDS: usize = 100;

fn main() {
    let map: CMap2<f64> = GridBuilder::<2, _>::unit_triangles(N_SQUARES);

    // fetch all vertices that are not on the boundary of the map
    let nodes: Vec<(VertexIdType, Vec<VertexIdType>)> = map
        .iter_vertices()
        .filter_map(|v| {
            // the condition detects if we're on the boundary
            if map
                .orbit(OrbitPolicy::Vertex, v as DartIdType)
                .any(|d| map.beta::<2>(d) == NULL_DART_ID)
            {
                None
            } else {
                // the orbit transformation yields neighbor IDs
                Some((
                    v,
                    map.orbit(OrbitPolicy::Vertex, v as DartIdType)
                        .map(|d| map.vertex_id(map.beta::<2>(d)))
                        .collect(),
                ))
            }
        })
        .collect();

    // main loop
    let mut round = 0;
    loop {
        // process nodes in parallel
        nodes.par_iter().for_each(|(vid, neigh)| {
            // we need a transaction here to avoid UBs, since there's
            // no guarantee we won't process neighbor nodes concurrently
            //
            // the transaction will ensure that we do not validate an operation
            // where inputs have changed due to instruction interleaving between threads
            // here, it will retry the transaction until it can be validated
            atomically(|trans| {
                let mut new_val = Vertex2::default();
                for v in neigh {
                    let vertex = map.read_vertex(trans, *v)?.unwrap();
                    new_val.0 += vertex.0;
                    new_val.1 += vertex.1;
                }
                new_val.0 /= neigh.len() as f64;
                new_val.1 /= neigh.len() as f64;
                map.write_vertex(trans, *vid, new_val)
            });
        });

        round += 1;
        if round >= N_ROUNDS {
            break;
        }
    }

    std::hint::black_box(map);
}

Breakdown

The main map structure, CMap2, can be edited in parallel using transactions to ensure algorithm correctness.

In the main computation loop, we use a transaction to ensure each new vertex value is computed from the current neighbor’s values. The errors generated by read_vertex and write_vertex are used to (early) detect any changes to the data used in the transaction, here, the list of neigh vertices.

At the end of the transaction block, the commit routine will check again if any used data has been altered. If not, results of the transaction will be validated and written to memory.

License

The following notice only applies to this book. For the library’s license, refer to the repository.

This work, including all text and figures, is licensed under Creative Commons Attribution 4.0 International (CC BY 4.0). You are free to share and adapt the material for any purpose, provided you give appropriate credit, link to the license, and indicate if changes were made.