AdaptiveSimplex

AdaptiveSimplex is a header-only C++20 library for adaptive integration of discontinuous or localized non-smooth functions on simplex meshes.

It is built for high-performance computing workflows where function or vertex evaluation is expensive, and the important structure lives near an interface. A typical example is a level-set or occupation integral:

$$ I = \int_\Omega g(x),\mathbf{1}_{\phi(x)\le \ell},dx, $$

Why AdaptiveSimplex?

• ⚡ Low overhead: mesh, cache, and cut-simplex bookkeeping stay near the microsecond scale in the supplied benchmarks, so expensive application evaluations can dominate.

• 🎯 Adaptive by target: set an error target and refinement budget; the algorithm samples where the current estimator says more information is needed, instead of committing to a fixed tensor grid up front.

• 🔁 Strong vertex reuse: every evaluated vertex is stored in VertexCache<T>, so shared simplex vertices are reused instead of re-evaluated.

• 📈 Incremental updates: refinement updates the integral through local corrections, evaluating only newly introduced vertices and avoiding full recomputation.

• 🌐 Dimension-general: the adaptive simplex machinery is not tied to 2D or level sets; it works for general simplex integrals in any dimension.

• 🧵 Parallel vertex evaluation: missing vertices can be evaluated in parallel with OpenMP when available.

• ✂️ Interface-aware simplex integration: for level-set problems, each cut simplex contributes the clipped volume or barycentric moments implied by the vertex data.

When To Use This

Good fits:

• discontinuous or non-smooth integrands with localized structure
• level-set-defined regions, cut cells, and occupation functions
• Fermi-sea or zero-temperature occupation calculations
• expensive vertex evaluations where reuse matters
• problems with a meaningful local error or interface indicator

Poor fits:

• smooth integrals

For smooth integrands, use standard quadrature methods, sparse grids, or even uniform grids. They are usually simpler, faster, and more accurate for that problem class. The smooth two-peak comparison below is included as a reminder: AdaptiveSimplex is not a smooth-quadrature silver bullet, and tensor grids can be competitive or better when the function is cheap and smooth.

Quick Start: Level-Set Integral

This example integrates the occupied part of a rotated ellipse:

$$ \phi(x,y) \le 1. $$

The executable constructs a mesh, a vertex cache, and an integrand built from three domain-specific pieces:

#include "adaptivesimplex/adaptive/adaptive_loop.h"
#include "adaptivesimplex/adaptive/simplex_integrand.h"
#include "adaptivesimplex/cut.h"
#include "adaptivesimplex/core/root_mesh.h"
#include "adaptivesimplex/core/vertex_cache.h"
#include "adaptivesimplex_examples/level_set_function.h"

namespace adaptive = adaptivesimplex::adaptive;
namespace core = adaptivesimplex::core;
namespace examples = adaptivesimplex::examples;

auto function = examples::LevelSetIndicator2D{};

auto geometry = core::root_geometry(2, 2);
auto cache = core::VertexCache<examples::LevelSetVertexValue>{};

auto integrand = adaptive::simplex_integrand(
    cache,

    // 1. Evaluate one mesh vertex.
    [function](std::span<const double> point) {
        return examples::LevelSetVertexValue{
            .level_value = function.level_value(point[0], point[1]),
        };
    },

    // 2. Compute one simplex contribution from cached vertex data.
    [function](
        const core::Geometry& geometry,
        core::SimplexId simplex_id,
        const core::VertexCache<examples::LevelSetVertexValue>& vertex_cache
    ) {
        const auto moments = adaptivesimplex::cut::simplex_moments(
            geometry,
            simplex_id,
            [&](core::VertexId vertex_id) {
                return vertex_cache.get(vertex_id).level_value;
            },
            adaptivesimplex::cut::LevelOptions{.level = function.level()}
        );
        return moments.volume;
    }
);

auto options = adaptive::Options{
    .target_error = 5e-4,
    .max_refinements = 2048,
    .preview_depth = 2,
    .min_refinement_batch_size = 1,
    .max_refinement_batch_size = 4,
};

const auto result = adaptive::run(geometry, integrand, options);

For the same example with the helper from examples/common, the executable is shorter:

auto function = examples::LevelSetIndicator2D{};
auto geometry = examples::level_set_geometry();
auto cache = core::VertexCache<examples::LevelSetVertexValue>{};
auto integrand = examples::make_level_set_integrand(cache, function);

const auto result =
    adaptive::run(geometry, integrand, examples::level_set_options());

std::cout << result.integral << "\n";
std::cout << result.stopping_error << "\n";

The important data flow is still explicit:

1. adaptive::run(...) asks the mesh which vertices are missing.
2. The vertex evaluator computes data for each missing vertex.
3. VertexCache<T> stores that data by VertexId.
4. The simplex rule computes coarse and preview contributions.
5. Estimation policies turn local correction into global stopping error and local refinement score.
6. The adaptive loop refines high-score simplices and updates totals incrementally.

For scalar integrals, the default policies are:

adaptive::estimation_policies{
    adaptive::signed_correction_error{},
    adaptive::absolute_local_correction{},
};

signed_correction_error maintains sum(local.correction) incrementally and uses std::abs(sum) as the stopping error. Use adaptive::unsigned_correction_error{} when you want sum(std::abs(local.correction)) with no cancellation.

Set preview_depth to zero to evaluate only the active mesh, without speculative child simplices or preview-only vertices. In this mode, simplex_integrand(...) uses the coarse contribution as the preview and produces a shape-compatible zero correction. Refinement still bisects selected simplices by one level, and the returned integral represents the active mesh. Depth zero must be paired with stopping-error and refinement-score policies that use an intrinsic simplex estimate; the built-in correction-only policies otherwise report zero error immediately.

See examples/scalar_integral.cpp for the minimal scalar factory example and examples/level_set_integral.cpp for the runnable level-set example.

Advanced: Hand-Written Integrands

simplex_integrand(...) builds the ordinary integrand protocol for you. If you need full control, a hand-written integrand supplies:

• using vertex_value_type = ...
• using integral_value_type = ...
• cache()
• stopping_global_error()
• evaluate_vertex(geometry, vertex_id)
• estimate_simplices(geometry, simplex_ids, preview_depth)

estimate_simplices(...) returns one adaptive::SimplexEstimate<T> per simplex:

adaptive::SimplexEstimate<T>{
    .coarse = coarse,
    .preview = preview,
    .correction = preview - coarse,
    .refinement_score = score,
};

adaptive::IntegrationResult<T> contains the final accumulated values:

• integral: current integral estimate
• correction: accumulated signed correction state
• stopping_error: scalar stopping criterion
• evaluations: number of vertex evaluations
• refinements: number of selected parent refinements
• converged: whether the target was reached

Why Cut-Simplex Moments?

Adaptive refinement decides where to refine. Cut-simplex moments decide what a simplex contributes when a discontinuity crosses it.

For a level-set integral, simply averaging vertex indicator values is a poor local quadrature rule: the function jumps across an interface. Instead, AdaptiveSimplex linearly interpolates the level value on each simplex (T):

$$ \phi_h(x)=\sum_i\phi_i\lambda_i(x), $$

clips the simplex to the occupied region

$$ T^-={x\in T:\phi_h(x)\le \ell}, $$

and computes barycentric moments

$$ M_i=\int_{T^-}\lambda_i(x),dx. $$

For an indicator-only volume calculation, the moment routine already returns the occupied simplex volume. That is why the simplex rule above returns moments.volume.

This is the main reason adaptivesimplex::cut exists: it turns a level-set-defined discontinuity into a deterministic simplex contribution that can be reused by the adaptive algorithm.

Essential API Concepts

• core::Geometry

  Owns the dimension, dyadic vertices, active simplices, refinement, and simplex volumes. Create a basic domain with core::root_geometry(dimension, depth).

• core::VertexCache<T>

  Stores user-defined vertex data keyed by VertexId. Missing vertices are evaluated once; reused vertices are not counted again as evaluations.

• adaptive::simplex_integrand(...)

  Factory that connects a cache, vertex evaluator, simplex rule, and estimation policies to the adaptive loop.

• Vertex evaluator

  Computes the cache value for one vertex. In a physics code, this is where an expensive Hamiltonian diagonalization or band calculation would usually live.

• Simplex rule

  Computes the domain-specific integral contribution for one simplex. For level-set integrals this usually calls cut::simplex_moments(...).

• adaptive::estimation_policies

  Selects the global stopping-error accumulator and local refinement-score policy.

• adaptive::SimplexEstimate<T>

  Contains coarse, preview, correction, and refinement_score for one active simplex.

• adaptive::IntegrationResult<T>

  Contains integral, correction, stopping_error, evaluations, refinements, and converged for the full adaptive run.

• adaptive::Options

  Controls stopping and refinement behavior: target error, maximum refinements, preview depth, and refinement batch size.

• adaptive::run(...)

  Evaluates missing vertices, updates estimates, refines high-priority simplices, and returns the integration result.

Library Targets

The package exports:

• adaptivesimplex::core
• adaptivesimplex::cut
• adaptivesimplex::adaptive

The adaptive layer depends on core, but it does not depend on cut. Link adaptivesimplex::cut when your code directly includes or uses the cut-simplex utility.

Public headers are under include/adaptivesimplex:

#include <adaptivesimplex/core.h>
#include <adaptivesimplex/cut.h>
#include <adaptivesimplex/adaptive.h>
#include <adaptivesimplex/adaptivesimplex.h>

Build And Use

After installing AdaptiveSimplex:

find_package(adaptivesimplex CONFIG REQUIRED)

add_executable(my_integrator main.cpp)
target_compile_features(my_integrator PRIVATE cxx_std_20)
target_link_libraries(my_integrator PRIVATE adaptivesimplex::adaptive)

If your code uses cut-simplex moments directly:

target_link_libraries(
  my_integrator PRIVATE
  adaptivesimplex::cut
  adaptivesimplex::adaptive
)

With CMake FetchContent:

include(FetchContent)

FetchContent_Declare(
  adaptivesimplex
  GIT_REPOSITORY <repository-url>
  GIT_TAG v0.1.0
)
FetchContent_MakeAvailable(adaptivesimplex)

target_link_libraries(my_integrator PRIVATE adaptivesimplex::adaptive)

When used as a subproject, tests, examples, benchmarks, and documentation tools default to off. Main CMake options:

• ADAPTIVESIMPLEX_BUILD_TESTING
• ADAPTIVESIMPLEX_BUILD_EXAMPLES
• ADAPTIVESIMPLEX_BUILD_BENCHMARKS
• ADAPTIVESIMPLEX_BUILD_DOC_TOOLS
• ADAPTIVESIMPLEX_ENABLE_OPENMP

Examples

User-facing examples:

• examples/level_set_integral.cpp: adaptive integration of a level-set indicator using adaptivesimplex::cut
• examples/two_peak_integral.cpp: smooth two-peak integral, useful as a counterexample where tensor grids can be competitive
• examples/scalar_integral.cpp: minimal simplex_integrand setup
• examples/vector_integral.cpp: vector-valued result with explicit policies

Build and run examples:

cmake -S . -B build/examples -DADAPTIVESIMPLEX_BUILD_EXAMPLES=ON
cmake --build build/examples

build/examples/adaptivesimplex_level_set_integral
build/examples/adaptivesimplex_two_peak_integral

Development

cmake -S . -B build/local \
  -DADAPTIVESIMPLEX_BUILD_TESTING=ON \
  -DADAPTIVESIMPLEX_BUILD_EXAMPLES=ON \
  -DADAPTIVESIMPLEX_BUILD_BENCHMARKS=ON
cmake --build build/local
ctest --test-dir build/local --output-on-failure

CI runs the package check on Ubuntu GCC, Ubuntu Clang, and macOS Clang. The package check builds tests/examples/benchmarks, runs the tests, smoke-runs the benchmarks, installs the package, and builds a downstream package consumer.

OpenMP is optional. If CMake finds it and ADAPTIVESIMPLEX_ENABLE_OPENMP is on, adaptivesimplex::adaptive links OpenMP::OpenMP_CXX and the missing-vertex loop in adaptive/evaluation.h can run in parallel.