ExactBandIntegral.hpp

#pragma once
 
#include <Eigen/Core>
#include <Eigen/Geometry>
 
#include <wmtk/envelope/Envelope.hpp>
 
#include <array>
#include <cstdint>
#include <memory>
#include <mutex>
#include <utility>
#include <vector>
 
namespace wmtk::components::topological_offset {
 
class ExactBandIntegral
{
public:
    ExactBandIntegral(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
        : ExactBandIntegral(V, F, {}, 0)
    {}
    ExactBandIntegral(
        const Eigen::MatrixXd& V,
        const Eigen::MatrixXi& F,
        const std::vector<int8_t>& face_side,
        int sigma,
        size_t value_slots = 0,
        std::shared_ptr<const wmtk::Envelope> envelope = nullptr);
 
    struct Stats
    {
        int64_t cells = 0; 
        int64_t candidates = 0; 
        int64_t pieces = 0; 
        int64_t rule_a_pieces = 0;
        double vol = 0.; 
        double exact_vol = 0.; 
        double rule_a_vol = 0.; 
        double empty_vol = 0.; 
        int64_t clips = 0; 
        int64_t clips_closed_up = 0; 
        int64_t caps_split = 0; 
        bool check_closed = false;
        int64_t closed_checked = 0, open_results = 0;
        int64_t point_faces = 0; 
        int64_t line_faces_radial = 0; 
        int64_t line_faces_parallel = 0; 
        int64_t line_faces_series = 0; 
        int64_t closed_edges = 0; 
        int64_t series_short = 0; 
        int64_t series_long = 0; 
        int64_t side_clipped = 0, side_rule_a_split = 0, side_D = 0, side_D_pos = 0;
        bool timing = false;
        double t_total = 0., t_candidates = 0., t_clip = 0., t_closed = 0., t_series = 0.;
    };
 
    double value(const std::array<Eigen::Vector3d, 4>& t, double delta, Stats* stats = nullptr)
        const;
 
    struct Work;
    class Partition
    {
    public:
        Partition();
        ~Partition();
        Partition(Partition&&) noexcept;
        Partition& operator=(Partition&&) noexcept;
 
    private:
        friend class ExactBandIntegral;
        std::unique_ptr<Work> m_w;
        uint64_t m_owner = 0; 
        double m_key[12] = {};
        bool m_ok = false; 
    };
    double value(const std::array<Eigen::Vector3d, 4>& t, double delta, Partition& keep) const;
    void reserve_values(size_t slots) const;
 
    double value_grad_hess(
        const std::array<Eigen::Vector3d, 4>& t,
        int kv,
        double delta,
        Eigen::Vector3d& grad,
        Eigen::Matrix3d& hess,
        Stats* stats = nullptr) const;
    double value_grad_hess(
        const std::array<Eigen::Vector3d, 4>& t,
        int kv,
        double delta,
        Eigen::Vector3d& grad,
        Eigen::Matrix3d& hess,
        Partition& keep) const;
    void grad_hess(
        const std::array<Eigen::Vector3d, 4>& t,
        int kv,
        double delta,
        Eigen::Vector3d& grad,
        Eigen::Matrix3d& hess,
        Partition& keep) const;
 
    struct MaxInfo
    {
        Eigen::Vector3d at; 
        bool rule_a = false; 
        double U = 0.; 
        int candidates = 0, pieces = 0, rule_a_pieces = 0;
    };
 
    double max_distance(const std::array<Eigen::Vector3d, 3>& f, MaxInfo* info = nullptr) const;
 
    // --------------------------------------------------------------------------------------------
    // The front form (EXPERIMENTAL_band_volume_rule "front_flux"; Python prototype front_form.py,
    // d_vertex.py, e_hessian.py, scratchpad front_form/ of 2026-10-13). Per feature k, in its
    // local region, a field G_k with div G_k = d - delta, tangent to the region's planes and 0 on
    // the input (d = 0):
    //     face    G = (d^2/2 - delta d) W        d = W.x - c, W = +-n on each side of its plane
    //     edge    G = (d/3 - delta/2) (x - x_e)  x_e the projection of x on the edge's line
    //     vertex  G = (d/4 - delta/3) (x - v).
    // On a band B made of cells, the divergence theorem cell by cell gives
    //     E(B) = int_B (d - delta) dV = sum over the boundary faces f of B of int_f G.n dA
    //            + sum over the cells t of D_t,
    // D_t = int_t (d - delta) dV - oint_dt G.n dA: the jumps of G.n across the BISECTORS inside t
    // (the sets where two features are nearest), 0 for a cell no bisector crosses. Every piece of
    // the partition (as value()'s) is labelled with its nearest feature exactly: a member whose
    // squared distance is not below another's anywhere on the piece is dropped (the exact range
    // of the difference, a quadratic, over the piece); members that still cross must have a
    // planar bisector (face-face, vertex-vertex, parallel or intersecting edge-edge), and the
    // piece is split by it. A crossing pair with a curved bisector (face-edge, face-vertex,
    // edge-vertex, skew edges) THROWS std::runtime_error naming the pair and the position:
    // there is no approximation in this form (rule (a)'s bound is value()'s only).
    //
    // The offset side (class doc), as value(): a labelled part counts only on the offset side of
    // its label's plane. On a labelled part the label is the nearest feature, whose plane meets it
    // only where d = 0 (the pseudonormal theorem), so the side drops whole parts: a face's far
    // side, and the edge and vertex parts beyond the input (MEASURED, the side cases' cells of
    // the tests: 47 edge parts clipped by their side plane, every one dropped whole, none cut).
    // Where G = 0 (on the input) the
    // dropped parts carry no flux, so the identity holds with the side as without it. The
    // Python prototype used the unsigned d (both sides of every face counted); the two agree on
    // every part on the offset side, and constructed without a side this class is the prototype.
 
    struct FaceWork;
    class FaceParts
    {
    public:
        FaceParts();
        ~FaceParts();
        FaceParts(FaceParts&&) noexcept;
        FaceParts& operator=(FaceParts&&) noexcept;
 
    private:
        friend class ExactBandIntegral;
        std::unique_ptr<FaceWork> m_w;
        uint64_t m_owner = 0; 
        double m_key[9] = {};
    };
    double front_face(
        const std::array<Eigen::Vector3d, 3>& f,
        double delta,
        FaceParts* keep = nullptr) const;
    void front_face_derivs(
        const std::array<Eigen::Vector3d, 3>& f,
        double delta,
        int kv,
        Eigen::Vector3d& grad,
        Eigen::Matrix3d& hess,
        FaceParts* keep = nullptr) const;
    class QuietCurvedBisector
    {
    public:
        QuietCurvedBisector();
        ~QuietCurvedBisector();
        QuietCurvedBisector(const QuietCurvedBisector&) = delete;
        QuietCurvedBisector& operator=(const QuietCurvedBisector&) = delete;
 
    private:
        bool m_prev = false;
    };
    double bisector_jump(
        const std::array<Eigen::Vector3d, 4>& t,
        double delta,
        Partition* keep = nullptr) const;
    void prepare_front_flux(double delta);
    double front_flux_reach() const { return m_ff_reach; }
    size_t front_flux_polytopes() const { return m_ff_boxes.size(); }
    bool front_flux_flag(const std::array<Eigen::Vector3d, 4>& t) const
    {
        return front_flux_flag(t.data(), 4);
    }
    bool front_flux_flag(const Eigen::Vector3d* t, int n) const;
    double distance(const Eigen::Vector3d& x) const { return upper_bound(&x, 1); }
    double distance_bound(const std::array<Eigen::Vector3d, 4>& t) const
    {
        return upper_bound(t.data(), 4);
    }
    double distance_bound(const Eigen::Vector3d* t, const int n) const { return upper_bound(t, n); }
 
    size_t num_faces() const { return m_A.size(); }
 
private:
    struct Feature
    {
        int kind = 0;
        Eigen::Vector3d p, dir;
        int h_begin = 0, h_end = 0;
        Eigen::Vector3d lo, hi; 
    };
    std::vector<Feature> m_feat; 
    std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> m_hs;
    int m_nF = 0, m_nE = 0;
    std::vector<Eigen::Vector3d> m_A, m_B, m_C; 
    std::vector<std::array<int, 6>> m_face_feats; 
    int m_sigma = 0;
    std::vector<char> m_side_has;
    std::vector<Eigen::Vector3d> m_side_n, m_side_q;
 
    std::shared_ptr<const wmtk::Envelope> m_env;
    double m_pad = 0.;
    struct Node
    {
        Eigen::Vector3d lo, hi;
        int left = -1, right = -1, begin = 0, end = 0;
    };
 
    struct Local;
    struct TriWork;
    void localize(const Eigen::Vector3d* t, int n, Local& w, Stats* stats) const;
    bool partition(const std::array<Eigen::Vector3d, 4>& t, Work& w, Stats* stats) const;
    const Work* cell_partition(
        const std::array<Eigen::Vector3d, 4>& t,
        Stats* stats,
        Work* own,
        Partition* keep) const;
    uint64_t m_id = 0;
    double evaluate(
        const std::array<Eigen::Vector3d, 4>& t,
        double delta,
        Stats* stats,
        Partition* keep) const;
    double value_tabled(const std::array<Eigen::Vector3d, 4>& t, double delta, Partition* keep)
        const;
    double value_grad_hess_impl(
        const std::array<Eigen::Vector3d, 4>& t,
        int kv,
        double delta,
        Eigen::Vector3d& grad,
        Eigen::Matrix3d& hess,
        Stats* stats,
        Partition* keep,
        bool want_value) const;
    struct ValueSlot
    {
        double in[13];
        double out = 0.;
        bool full = false;
    };
    static constexpr size_t kValueStripes = 1024;
    mutable std::vector<ValueSlot> m_values;
    mutable std::unique_ptr<std::mutex[]> m_value_locks;
    void tri_partition(const std::array<Eigen::Vector3d, 3>& f, TriWork& w) const;
    double upper_bound(const Eigen::Vector3d* t, int n) const;
    void check_cell_pieces(
        const std::array<Eigen::Vector3d, 4>& t,
        double delta,
        const Work& w,
        double vol_sum,
        double empty_vol) const;
    void check_face_pieces(
        const std::array<Eigen::Vector3d, 3>& f,
        const TriWork& w,
        const char* who) const;
    void label_face(const std::array<Eigen::Vector3d, 3>& f, FaceWork& fw) const;
    const FaceWork& face_parts(const std::array<Eigen::Vector3d, 3>& f, FaceParts* keep) const;
    struct FFBox
    {
        Eigen::Vector3d lo, hi;
        std::vector<std::pair<Eigen::Vector3d, double>> planes;
    };
    std::vector<FFBox> m_ff_boxes;
    std::vector<Node> m_ff_nodes;
    std::vector<int> m_ff_order;
    double m_ff_reach = 0.;
    int build_ff_node(int begin, int end);
    void candidates(
        const Eigen::Vector3d& lo,
        const Eigen::Vector3d& hi,
        double r,
        std::vector<int>& out) const;
};
 
} // namespace wmtk::components::topological_offset