scalar_integral.cpp

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

#include <iostream>
#include <span>

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

int main() {
    // Start with a coarse 1D mesh over [0, 1].
    auto geometry = core::root_geometry(1, 2);
    auto cache = core::VertexCache<double>{};

    auto integrand = adaptive::simplex_integrand(
        cache,
        [](std::span<const double> point) {
            const double x = point[0];
            return x * x;
        },
        [](const core::Geometry& geometry, core::SimplexId simplex_id,
           const core::VertexCache<double>& vertex_cache) {
            const auto& simplex = geometry.simplices().simplex(simplex_id);
            double sum = 0.0;
            for (const auto vertex_id : simplex.vertex_ids) {
                sum += vertex_cache.get(vertex_id);
            }
            return simplex.volume * sum / static_cast<double>(simplex.vertex_ids.size());
        }
    );

    const auto result = adaptive::run(geometry, integrand, adaptive::Options{
        .target_error = 1e-4,
        .max_refinements = 128,
        .preview_depth = 1,
    });

    std::cout << "integral ~= " << result.integral << "\n";
    std::cout << "estimated error = " << result.stopping_error << "\n";
    return result.converged ? 0 : 1;
}