level_set_integral.cpp

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

#include <cmath>
#include <iomanip>
#include <iostream>

int main() {
    namespace adaptive = adaptivesimplex::adaptive;
    namespace core = adaptivesimplex::core;
    namespace examples = adaptivesimplex::examples;

    // This example integrates the occupied part of a rotated ellipse:
    //     level_value(x, y) <= 1
    // The integrand uses cut-simplex moments, so each simplex contributes the
    // clipped area implied by linearly interpolated vertex level values.
    auto function = examples::LevelSetIndicator2D{};

    // Geometry owns the adaptive simplex mesh over [0, 1]^2. The first
    // argument is dimension; the second is the initial root subdivision depth.
    auto geometry = core::root_geometry(2, 2);

    // VertexCache stores only newly evaluated vertices. Reused vertices are not
    // counted again as function evaluations.
    auto cache = core::VertexCache<examples::LevelSetVertexValue>{};

    // The integrand connects vertex evaluation, cut-simplex contribution, and
    // estimation policies to the generic adaptive loop.
    auto integrand = examples::make_level_set_integrand(cache, function);

    // Options control stopping and refinement. preview_depth controls how far
    // each active simplex is temporarily refined when estimating local error.
    auto options = adaptive::Options{
        .target_error = 5e-4,
        .max_refinements = 2048,
        .preview_depth = 2,
        .min_refinement_batch_size = 1,
        .max_refinement_batch_size = 4,
    };

    // adaptive::run evaluates missing vertices, estimates active simplices, and
    // refines the simplices with the largest local scores.
    const auto result = adaptive::run(geometry, integrand, options);
    const double exact = function.exact_integral();

    std::cout << std::setprecision(17);
    std::cout << "integral ~= " << result.integral << "\n";
    std::cout << "exact = " << exact << "\n";
    std::cout << "absolute error = " << std::abs(result.integral - exact) << "\n";
    std::cout << "estimated error = " << result.stopping_error << "\n";
    std::cout << "evaluations = " << result.evaluations << "\n";
    std::cout << "refinements = " << result.refinements << "\n";
    std::cout << "active simplices = " << geometry.simplices().n_active() << "\n";
    std::cout << "cached vertices = " << cache.size() << "\n";

    return result.converged ? 0 : 1;
}