module;

#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <memory>
#include <span>
#include <stdexcept>
#include <utility>
#include <vector>

#include <mfem.hpp>
#include <mpi.h>

// Outward enclosure depends on ordinary floating-point operations. Optimizers
// that discard NaNs or reassociate expressions invalidate the certificate.
#if defined(__FAST_MATH__) || (defined(__FINITE_MATH_ONLY__) && __FINITE_MATH_ONLY__ > 0)
#error "Polynomial enthalpy certification does not support fast-math or finite-math-only compilation."
#endif

export module mean_field:solver.enthalpy_admissibility;
export import :field.mfem;

export namespace mean_field::solver {
    enum class PolynomialSignStatus { certified_nonnegative, unresolved, negative, non_finite };

    struct EnthalpyAdmissibilityOptions final {
        int maximumSubdivisionDepth{24};
        int maximumSubcells{8192};
    };

    struct PolynomialSignReport final {
        PolynomialSignStatus status{PolynomialSignStatus::certified_nonnegative};
        double lowerBound{std::numeric_limits<double>::infinity()};
        double upperBound{-std::numeric_limits<double>::infinity()};
        int visitedSubcells{0};
        int maximumDepth{0};
        int witnessElement{-1};
        std::array<double, 3> witnessReferencePoint{};
    };

    enum class EnthalpyDirectionRepairStatus {
        feasible, invalid_base, invalid_anchor, invalid_direction, unresolved, no_descent
    };

    struct EnthalpyDirectionRepairReport final {
        EnthalpyDirectionRepairStatus status{EnthalpyDirectionRepairStatus::unresolved};
        double minimumSafeFraction{1.0};
        double chosenFraction{0.0};
        double meritDirectionalDerivative{std::numeric_limits<double>::quiet_NaN()};
        double linearModelResidualNorm{std::numeric_limits<double>::quiet_NaN()};
        int jacobianApplications{0};
        // Early rejection can leave later checks unvisited. Such a report is
        // never evidence of feasibility, even though the empty-domain check is.
        PolynomialSignReport base{.status = PolynomialSignStatus::unresolved};
        PolynomialSignReport anchor{.status = PolynomialSignStatus::unresolved};
        PolynomialSignReport candidate{.status = PolynomialSignStatus::unresolved};
    };
}

namespace mean_field::solver::enthalpy_detail {
    static_assert(std::numeric_limits<double>::is_iec559 && std::numeric_limits<double>::radix == 2 &&
                  std::numeric_limits<double>::digits == 53 && std::numeric_limits<double>::max_exponent == 1024 &&
                  std::numeric_limits<double>::denorm_min() > 0.0 &&
                  std::numeric_limits<double>::denorm_min() < std::numeric_limits<double>::min(),
                  "Polynomial enthalpy certification requires IEEE 754 binary64 with subnormal support.");

    // Bounds concern the represented polynomial, not a tolerance-relaxed EOS.
    // Exact zero is retained algebraically so a prescribed zero trace does not
    // become an unresolved negative interval merely through outward rounding.
    struct Interval final {
        double lo{0.0}, hi{0.0};
        explicit Interval(double value = 0.0) : lo(value), hi(value) {}
        Interval(double lower, double upper) : lo(lower), hi(upper) {}
        bool zero() const { return lo == 0.0 && hi == 0.0; }
        bool finite() const { return std::isfinite(lo) && std::isfinite(hi); }
    };
    inline double Down(double value) { return std::nextafter(value, -std::numeric_limits<double>::infinity()); }
    inline double Up(double value) { return std::nextafter(value, std::numeric_limits<double>::infinity()); }
    inline Interval Add(Interval a, Interval b) {
        if (a.zero()) return b;
        if (b.zero()) return a;
        return {Down(a.lo + b.lo), Up(a.hi + b.hi)};
    }
    inline Interval Negate(Interval a) { return {-a.hi, -a.lo}; }
    inline Interval Multiply(Interval a, Interval b) {
        if (a.zero() || b.zero()) return Interval{};
        if (a.lo == 1.0 && a.hi == 1.0) return b;
        if (b.lo == 1.0 && b.hi == 1.0) return a;
        const std::array products{a.lo*b.lo, a.lo*b.hi, a.hi*b.lo, a.hi*b.hi};
        return {Down(*std::min_element(products.begin(), products.end())),
                Up(*std::max_element(products.begin(), products.end()))};
    }
    inline Interval Divide(Interval a, Interval b) {
        if (b.lo <= 0.0 && b.hi >= 0.0) throw std::invalid_argument("A tensor nodal basis has repeated nodes.");
        if (a.zero()) return Interval{};
        return Multiply(a, {Down(1.0/b.hi), Up(1.0/b.lo)});
    }
    inline void ValidateOptions(const EnthalpyAdmissibilityOptions &options) {
        if (options.maximumSubdivisionDepth < 0 || options.maximumSubdivisionDepth > 60 ||
            options.maximumSubcells < 1)
            throw std::invalid_argument("Polynomial certification requires a finite positive work budget.");
    }
    inline int TensorSize(int points, int dimension) {
        std::size_t count = 1;
        for (int d = 0; d < dimension; ++d) {
            count *= static_cast<std::size_t>(points);
            if (count > static_cast<std::size_t>(std::numeric_limits<int>::max()))
                throw std::invalid_argument("Tensor polynomial is too large.");
        }
        return static_cast<int>(count);
    }
    inline int CollectiveMaximum(int value, MPI_Comm comm) {
        int result = 0;
        if (MPI_Allreduce(&value, &result, 1, MPI_INT, MPI_MAX, comm) != MPI_SUCCESS)
            throw std::runtime_error("Enthalpy admissibility could not combine rank-local status.");
        return result;
    }
}

export namespace mean_field::solver {
    /*
     * Dimension belongs to the type; polynomial order belongs to the mesh.
     * Coefficients use lexicographic ordering with the first coordinate fastest.
     * Only tensor Lagrange polynomials on [0,1]^Dimension are represented.
     */
    template <int Dimension> requires (Dimension >= 1 && Dimension <= 3)
    class TensorBernsteinPolynomial final {
        using Interval = enthalpy_detail::Interval;
    public:
        explicit TensorBernsteinPolynomial(std::span<const double> nodes)
            : m_nodes(nodes.begin(), nodes.end()), m_points(static_cast<int>(nodes.size())),
              m_size(enthalpy_detail::TensorSize(m_points, Dimension)) {
            if (m_points < 2 || m_nodes.front() != 0.0 || m_nodes.back() != 1.0)
                throw std::invalid_argument("Tensor certification requires closed nodal points including 0 and 1.");
            for (int i = 0; i < m_points; ++i)
                if (!std::isfinite(m_nodes[i]) || (i && m_nodes[i] <= m_nodes[i-1]))
                    throw std::invalid_argument("Tensor nodes must be finite and strictly increasing.");
            BuildTransform();
        }

        [[nodiscard]] int Size() const noexcept { return m_size; }
        [[nodiscard]] int Order() const noexcept { return m_points-1; }

        [[nodiscard]] PolynomialSignReport Check(std::span<const double> nodalValues,
                                                 EnthalpyAdmissibilityOptions options = {}) const {
            enthalpy_detail::ValidateOptions(options);
            auto coefficients = Transform(nodalValues);
            double rootLower = std::numeric_limits<double>::infinity(), rootUpper = -rootLower;
            for (const auto &coefficient : coefficients) {
                rootLower = std::min(rootLower,coefficient.lo);
                rootUpper = std::max(rootUpper,coefficient.hi);
            }
            struct Cell { std::vector<Interval> coefficients; int depth; std::array<double,3> lo, hi; };
            std::vector<Cell> pending;
            pending.push_back({std::move(coefficients), 0, {0,0,0}, {1,1,1}});
            PolynomialSignReport report;
            bool unresolved = false;
            while (!pending.empty()) {
                Cell cell = std::move(pending.back());
                pending.pop_back();
                ++report.visitedSubcells;
                report.maximumDepth = std::max(report.maximumDepth, cell.depth);
                double lower = std::numeric_limits<double>::infinity(), upper = -lower;
                bool finite = true;
                for (const auto &coefficient : cell.coefficients) {
                    finite = finite && coefficient.finite();
                    lower = std::min(lower, coefficient.lo);
                    upper = std::max(upper, coefficient.hi);
                }
                if (!finite) { report.status = PolynomialSignStatus::non_finite; return report; }
                // Bernstein coefficients at tensor corners are actual values;
                // an arbitrary negative control coefficient is not a witness.
                for (int corner = 0; corner < (1<<Dimension); ++corner) {
                    int index = 0, stride = 1;
                    for (int axis = 0; axis < Dimension; ++axis) {
                        const bool high = (corner & (1<<axis)) != 0;
                        index += high ? (m_points-1)*stride : 0;
                        stride *= m_points;
                    }
                    if (cell.coefficients[index].hi < 0.0) {
                        report.status = PolynomialSignStatus::negative;
                        report.lowerBound = rootLower;
                        report.upperBound = rootUpper;
                        for (int axis = 0; axis < Dimension; ++axis)
                            report.witnessReferencePoint[axis] = (corner & (1<<axis)) ? cell.hi[axis] : cell.lo[axis];
                        return report;
                    }
                }
                if (lower >= 0.0) {
                    report.lowerBound = std::min(report.lowerBound, lower);
                    report.upperBound = std::max(report.upperBound, upper);
                    continue;
                }
                if (cell.depth >= options.maximumSubdivisionDepth ||
                    report.visitedSubcells + static_cast<int>(pending.size()) + 2 > options.maximumSubcells) {
                    unresolved = true;
                    report.lowerBound = std::min(report.lowerBound, lower);
                    report.upperBound = std::max(report.upperBound, upper);
                    continue;
                }
                const int axis = cell.depth % Dimension;
                auto children = Split(cell.coefficients, axis);
                const double midpoint = 0.5*(cell.lo[axis]+cell.hi[axis]);
                auto leftHi = cell.hi, rightLo = cell.lo;
                leftHi[axis] = midpoint;
                rightLo[axis] = midpoint;
                pending.push_back({std::move(children.second), cell.depth+1, rightLo, cell.hi});
                pending.push_back({std::move(children.first), cell.depth+1, cell.lo, leftHi});
            }
            report.status = unresolved ? PolynomialSignStatus::unresolved : PolynomialSignStatus::certified_nonnegative;
            return report;
        }

        // This computes a sufficient fraction in a convex inner approximation,
        // not the minimum over all nonnegative polynomials or a Newton optimum.
        [[nodiscard]] double MinimumRepairFraction(std::span<const double> anchor,
                                                   std::span<const double> candidate,
                                                   EnthalpyAdmissibilityOptions options = {}) const {
            enthalpy_detail::ValidateOptions(options);
            struct Cell { std::vector<Interval> anchor, candidate; int depth; };
            std::vector<Cell> pending;
            pending.push_back({Transform(anchor), Transform(candidate), 0});
            double fraction = 0.0;
            int visited = 0;
            while (!pending.empty()) {
                Cell cell = std::move(pending.back()); pending.pop_back(); ++visited;
                bool anchorCertified = true;
                for (const auto &a : cell.anchor) {
                    if (!a.finite()) return std::numeric_limits<double>::quiet_NaN();
                    anchorCertified = anchorCertified && a.lo >= 0.0;
                }
                if (!anchorCertified) {
                    if (cell.depth >= options.maximumSubdivisionDepth ||
                        visited + static_cast<int>(pending.size()) + 2 > options.maximumSubcells)
                        return std::numeric_limits<double>::quiet_NaN();
                    const int axis = cell.depth % Dimension;
                    auto a = Split(cell.anchor, axis), c = Split(cell.candidate, axis);
                    pending.push_back({std::move(a.second), std::move(c.second), cell.depth+1});
                    pending.push_back({std::move(a.first), std::move(c.first), cell.depth+1});
                    continue;
                }
                for (int i = 0; i < m_size; ++i) {
                    const double c = cell.candidate[i].lo, a = cell.anchor[i].lo;
                    if (!cell.candidate[i].finite()) return std::numeric_limits<double>::quiet_NaN();
                    if (c < 0.0) {
                        // Outward rounding covers both subtraction and division.
                        const double denominator = enthalpy_detail::Down(a-c);
                        const double needed = a == 0.0 ? 1.0 : enthalpy_detail::Up((-c)/denominator);
                        fraction = std::max(fraction, std::min(1.0, needed));
                    }
                }
            }
            return fraction;
        }

    private:
        void BuildTransform() {
            using namespace enthalpy_detail;
            m_transform.resize(static_cast<std::size_t>(m_points*m_points));
            const int degree = m_points-1;
            for (int node = 0; node < m_points; ++node) {
                std::vector<Interval> power(1, Interval(1.0));
                for (int other = 0; other < m_points; ++other) {
                    if (other == node) continue;
                    const Interval denominator = Add(Interval(m_nodes[node]), Negate(Interval(m_nodes[other])));
                    std::vector<Interval> next(power.size()+1);
                    for (std::size_t k = 0; k < power.size(); ++k) {
                        next[k] = Add(next[k], Divide(Multiply(power[k], Interval(-m_nodes[other])), denominator));
                        next[k+1] = Add(next[k+1], Divide(power[k], denominator));
                    }
                    power = std::move(next);
                }
                for (int coefficient = 0; coefficient < m_points; ++coefficient) {
                    Interval value;
                    Interval ratio(1.0);
                    for (int k = 0; k <= coefficient; ++k) {
                        value = Add(value, Multiply(power[k], ratio));
                        if (k < coefficient)
                            ratio = Multiply(ratio, Divide(Interval(coefficient-k), Interval(degree-k)));
                    }
                    m_transform[coefficient*m_points+node] = value;
                }
            }
            // Endpoint interpolation is exact and is also the algebraic reason
            // a face whose nodal values vanish has an exactly zero trace.
            for (int node = 0; node < m_points; ++node) {
                m_transform[node] = Interval(node == 0 ? 1.0 : 0.0);
                m_transform[(m_points-1)*m_points+node] = Interval(node == m_points-1 ? 1.0 : 0.0);
            }
        }

        [[nodiscard]] std::vector<Interval> Transform(std::span<const double> values) const {
            using namespace enthalpy_detail;
            if (values.size() != static_cast<std::size_t>(m_size))
                throw std::invalid_argument("Tensor polynomial has incompatible nodal coefficients.");
            std::vector<Interval> current; current.reserve(values.size());
            for (double value : values) current.emplace_back(value);
            int stride = 1;
            for (int axis = 0; axis < Dimension; ++axis) {
                std::vector<Interval> next(m_size);
                for (int index = 0; index < m_size; ++index) {
                    const int row = (index/stride)%m_points, line = index-row*stride;
                    Interval sum;
                    for (int column = 0; column < m_points; ++column)
                        sum = Add(sum, Multiply(m_transform[row*m_points+column], current[line+column*stride]));
                    next[index] = sum;
                }
                current = std::move(next); stride *= m_points;
            }
            return current;
        }

        [[nodiscard]] std::pair<std::vector<Interval>,std::vector<Interval>>
        Split(const std::vector<Interval> &values, int axis) const {
            using namespace enthalpy_detail;
            std::vector<Interval> left(m_size), right(m_size), row(m_points);
            int stride = 1;
            for (int d = 0; d < axis; ++d) stride *= m_points;
            for (int base = 0; base < m_size; ++base) {
                if ((base/stride)%m_points != 0) continue;
                for (int j = 0; j < m_points; ++j) row[j] = values[base+j*stride];
                left[base] = row[0]; right[base+(m_points-1)*stride] = row[m_points-1];
                for (int level = 1; level < m_points; ++level) {
                    for (int j = 0; j < m_points-level; ++j)
                        row[j] = Multiply(Add(row[j], row[j+1]), Interval(0.5));
                    left[base+level*stride] = row[0];
                    right[base+(m_points-1-level)*stride] = row[m_points-1-level];
                }
            }
            return {std::move(left),std::move(right)};
        }

        std::vector<double> m_nodes;
        int m_points, m_size;
        std::vector<Interval> m_transform;
    };

    /* Prepared against an immutable conforming scalar H1 space. All mutations
     * remain in reduced true-DOF coordinates; elementwise coefficient clipping
     * would break continuity and is deliberately absent. Collective methods
     * must be called in the same order on all ranks of the supplied space. */
    class PreparedEnthalpyAdmissibility final {
        struct ElementData final {
            int element;
            const TensorBernsteinPolynomial<3> *basis;
            std::vector<int> lexicographicToNative;
        };
    public:
        PreparedEnthalpyAdmissibility(mfem::ParFiniteElementSpace &space, field::FieldDofMap dofMap,
                                     std::span<const int> supportedElementIds,
                                     std::span<const int> zeroTraceReducedDofs = {},
                                     EnthalpyAdmissibilityOptions options = {})
            : m_space(&space), m_adapter(std::move(dofMap), space),
              m_zeroTrace(zeroTraceReducedDofs.begin(),zeroTraceReducedDofs.end()), m_options(options) {
            const bool validOptions = options.maximumSubdivisionDepth >= 0 &&
                options.maximumSubdivisionDepth <= 60 && options.maximumSubcells >= 1;
            if (enthalpy_detail::CollectiveMaximum(validOptions ? 0 : 1,space.GetComm()))
                throw std::invalid_argument("Polynomial certification requires a finite positive work budget.");
            bool eligible = space.GetVDim() == 1 && !space.Nonconforming() &&
                dynamic_cast<const mfem::H1_FECollection *>(space.FEColl()) != nullptr;
            for (int e : supportedElementIds) {
                if (e < 0 || e >= space.GetNE()) { eligible = false; continue; }
                const auto *fe = dynamic_cast<const mfem::NodalTensorFiniteElement *>(space.GetFE(e));
                eligible = eligible && fe && fe->GetDim() == 3 && fe->GetGeomType() == mfem::Geometry::CUBE &&
                    fe->GetMapType() == mfem::FiniteElement::VALUE &&
                    (fe->GetBasisType() == mfem::BasisType::GaussLobatto ||
                     fe->GetBasisType() == mfem::BasisType::ClosedUniform ||
                     fe->GetBasisType() == mfem::BasisType::ClosedGL) && fe->GetDof() ==
                    enthalpy_detail::TensorSize(fe->GetOrder()+1,3);
            }
            for (int row : m_zeroTrace) eligible = eligible && row >= 0 && row < m_adapter.dof_map().reduced_size();
            if (enthalpy_detail::CollectiveMaximum(eligible ? 0 : 1, space.GetComm()))
                throw std::invalid_argument("Enthalpy certification requires conforming scalar tensor H1 hexahedra.");
            const std::array localOptions{options.maximumSubdivisionDepth,options.maximumSubcells};
            std::array<int,2> minimum{}, maximum{};
            if (MPI_Allreduce(localOptions.data(),minimum.data(),2,MPI_INT,MPI_MIN,space.GetComm()) != MPI_SUCCESS ||
                MPI_Allreduce(localOptions.data(),maximum.data(),2,MPI_INT,MPI_MAX,space.GetComm()) != MPI_SUCCESS)
                throw std::runtime_error("Could not validate collective polynomial-certification options.");
            if (minimum != maximum) throw std::invalid_argument("Polynomial-certification budgets must agree on all ranks.");
            for (int e : supportedElementIds) {
                const auto &fe = dynamic_cast<const mfem::NodalTensorFiniteElement &>(*space.GetFE(e));
                const auto key = std::pair{fe.GetOrder(),fe.GetBasisType()};
                const auto &map = fe.GetDofMap();
                std::vector<int> ordering(fe.GetDof());
                for (int i = 0; i < fe.GetDof(); ++i) ordering[i] = map.Size() ? map[i] : i;
                if (!m_bases.contains(key)) {
                    std::vector<double> nodes(fe.GetOrder()+1);
                    for (int i = 0; i <= fe.GetOrder(); ++i) nodes[i] = fe.GetNodes().IntPoint(ordering[i]).x;
                    m_bases.emplace(key,std::make_unique<TensorBernsteinPolynomial<3>>(nodes));
                }
                m_elements.push_back({e,m_bases.at(key).get(),std::move(ordering)});
            }
        }

        [[nodiscard]] int Size() const noexcept { return m_adapter.dof_map().reduced_size(); }

        [[nodiscard]] PolynomialSignReport Check(const mfem::Vector &reducedEnthalpy) const {
            RequireSize(reducedEnthalpy);
            PolynomialSignReport local;
            bool finite = true;
            for (int i = 0; i < reducedEnthalpy.Size(); ++i) finite = finite && std::isfinite(reducedEnthalpy(i));
            if (enthalpy_detail::CollectiveMaximum(finite ? 0 : 1,m_space->GetComm())) {
                local.status = PolynomialSignStatus::non_finite; return local;
            }
            mfem::ParGridFunction grid(m_space); m_adapter.scatter(reducedEnthalpy,grid);
            for (const auto &element : m_elements) {
                const auto values = ElementValues(grid,element);
                auto report = element.basis->Check(values,m_options);
                local.visitedSubcells += report.visitedSubcells;
                local.maximumDepth = std::max(local.maximumDepth,report.maximumDepth);
                local.lowerBound = std::min(local.lowerBound,report.lowerBound);
                local.upperBound = std::max(local.upperBound,report.upperBound);
                if (static_cast<int>(report.status) > static_cast<int>(local.status)) {
                    local.status = report.status;
                    local.witnessElement = element.element;
                    local.witnessReferencePoint = report.witnessReferencePoint;
                }
            }
            local.status = static_cast<PolynomialSignStatus>(enthalpy_detail::CollectiveMaximum(
                static_cast<int>(local.status),m_space->GetComm()));
            double lower = local.lowerBound, upper = local.upperBound;
            if (MPI_Allreduce(&lower,&local.lowerBound,1,MPI_DOUBLE,MPI_MIN,m_space->GetComm()) != MPI_SUCCESS ||
                MPI_Allreduce(&upper,&local.upperBound,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
                throw std::runtime_error("Could not reduce distributed polynomial bounds.");
            // Witness indices and work counts intentionally remain rank-local.
            return local;
        }

        [[nodiscard]] EnthalpyDirectionRepairReport MinimumRepairFraction(
            const mfem::Vector &base, const mfem::Vector &anchor, const mfem::Vector &direction) const {
            RequireSize(base); RequireSize(anchor); RequireSize(direction);
            EnthalpyDirectionRepairReport report;
            report.base = Check(base);
            if (report.base.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(base)) {
                report.status = EnthalpyDirectionRepairStatus::invalid_base; return report;
            }
            // The anchor can be reused only for identical coefficients. All
            // ranks agree on the cache path before any later scatter/reduction.
            bool cached = m_cachedAnchor.Size() == anchor.Size() && m_hasCachedAnchor;
            if (cached)
                for (int i = 0; i < anchor.Size(); ++i) cached = cached && m_cachedAnchor(i) == anchor(i);
            if (enthalpy_detail::CollectiveMaximum(cached ? 0 : 1,m_space->GetComm())) {
                m_cachedAnchorReport = Check(anchor);
                m_cachedAnchor = anchor;
                m_hasCachedAnchor = true;
            }
            report.anchor = m_cachedAnchorReport;
            if (report.anchor.status != PolynomialSignStatus::certified_nonnegative || !ZeroTrace(anchor)) {
                report.status = EnthalpyDirectionRepairStatus::invalid_anchor; return report;
            }
            bool finiteDirection = true;
            for (int i = 0; i < direction.Size(); ++i)
                finiteDirection = finiteDirection && std::isfinite(direction(i));
            if (enthalpy_detail::CollectiveMaximum(finiteDirection ? 0 : 1,m_space->GetComm())) {
                report.status = EnthalpyDirectionRepairStatus::invalid_direction;
                report.candidate.status = PolynomialSignStatus::non_finite;
                return report;
            }
            mfem::Vector candidate(base); candidate += direction;
            for (int row : m_zeroTrace) candidate(row) = 0.0;
            report.candidate = Check(candidate);
            if (report.candidate.status == PolynomialSignStatus::non_finite) {
                report.status = EnthalpyDirectionRepairStatus::invalid_direction; return report;
            }
            if (report.candidate.status == PolynomialSignStatus::certified_nonnegative) {
                report.minimumSafeFraction = 0.0;
                report.status = EnthalpyDirectionRepairStatus::feasible; return report;
            }
            mfem::ParGridFunction anchorGrid(m_space), candidateGrid(m_space);
            m_adapter.scatter(anchor,anchorGrid); m_adapter.scatter(candidate,candidateGrid);
            double localFraction = 0.0;
            bool resolved = true;
            for (const auto &element : m_elements) {
                const double fraction = element.basis->MinimumRepairFraction(
                    ElementValues(anchorGrid,element),ElementValues(candidateGrid,element),m_options);
                resolved = resolved && std::isfinite(fraction);
                if (std::isfinite(fraction)) localFraction = std::max(localFraction,fraction);
            }
            if (enthalpy_detail::CollectiveMaximum(resolved ? 0 : 1,m_space->GetComm())) {
                report.status = EnthalpyDirectionRepairStatus::unresolved; return report;
            }
            if (MPI_Allreduce(&localFraction,&report.minimumSafeFraction,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
                throw std::runtime_error("Could not combine feasible enthalpy repair fractions.");
            // A small inward move covers rounding in the later state addition;
            // ApplyRepairFraction still checks the actual represented candidate.
            if (report.minimumSafeFraction > 0.0 && report.minimumSafeFraction < 1.0)
                report.minimumSafeFraction = std::min(1.0, report.minimumSafeFraction +
                    64.0*std::numeric_limits<double>::epsilon()*(1.0-report.minimumSafeFraction));
            report.status = EnthalpyDirectionRepairStatus::feasible;
            return report;
        }

        [[nodiscard]] PolynomialSignReport ApplyRepairFraction(
            const mfem::Vector &base, const mfem::Vector &anchor, const mfem::Vector &originalDirection,
            double fraction, mfem::Vector &repairedDirection) const {
            RequireSize(base); RequireSize(anchor); RequireSize(originalDirection);
            const bool validFraction = std::isfinite(fraction) && fraction >= 0.0 && fraction <= 1.0;
            if (enthalpy_detail::CollectiveMaximum(validFraction ? 0 : 1,m_space->GetComm()))
                throw std::invalid_argument("An enthalpy repair fraction must lie in [0,1].");
            double minimum = 0.0, maximum = 0.0;
            if (MPI_Allreduce(&fraction,&minimum,1,MPI_DOUBLE,MPI_MIN,m_space->GetComm()) != MPI_SUCCESS ||
                MPI_Allreduce(&fraction,&maximum,1,MPI_DOUBLE,MPI_MAX,m_space->GetComm()) != MPI_SUCCESS)
                throw std::runtime_error("Could not validate the collective enthalpy repair fraction.");
            if (minimum != maximum) throw std::invalid_argument("The enthalpy repair fraction must agree across ranks.");
            bool finite = true;
            for (int i = 0; i < Size(); ++i)
                finite = finite && std::isfinite(base(i)) && std::isfinite(anchor(i)) &&
                    std::isfinite(originalDirection(i));
            if (enthalpy_detail::CollectiveMaximum(finite ? 0 : 1,m_space->GetComm())) {
                PolynomialSignReport report;
                report.status = PolynomialSignStatus::non_finite;
                return report;
            }
            mfem::Vector repaired(Size()), candidate(Size());
            for (int i = 0; i < Size(); ++i) {
                repaired(i) = fraction == 0.0 ? originalDirection(i) :
                    fraction == 1.0 ? anchor(i)-base(i) :
                    (1.0-fraction)*originalDirection(i)+fraction*(anchor(i)-base(i));
            }
            for (int row : m_zeroTrace) repaired(row) = 0.0;
            candidate = base; candidate += repaired;
            auto report = Check(candidate);
            if (report.status == PolynomialSignStatus::certified_nonnegative) repairedDirection = repaired;
            return report;
        }

    private:
        void RequireSize(const mfem::Vector &vector) const {
            if (enthalpy_detail::CollectiveMaximum(vector.Size() == Size() ? 0 : 1,m_space->GetComm()))
                throw std::invalid_argument("Enthalpy feasibility received incompatible reduced true-DOF vectors.");
        }
        [[nodiscard]] bool ZeroTrace(const mfem::Vector &values) const {
            bool zero = true;
            for (int row : m_zeroTrace) zero = zero && values(row) == 0.0;
            return enthalpy_detail::CollectiveMaximum(zero ? 0 : 1,m_space->GetComm()) == 0;
        }
        [[nodiscard]] std::vector<double> ElementValues(const mfem::ParGridFunction &grid,
                                                       const ElementData &element) const {
            mfem::Array<int> dofs; mfem::Vector native;
            auto *transformation = m_space->GetElementDofs(element.element,dofs);
            grid.GetSubVector(dofs,native);
            if (transformation) transformation->InvTransformPrimal(native);
            std::vector<double> values(element.lexicographicToNative.size());
            for (std::size_t i = 0; i < values.size(); ++i) values[i] = native(element.lexicographicToNative[i]);
            return values;
        }
        mfem::ParFiniteElementSpace *m_space;
        field::FieldDofGridFunctionAdapter m_adapter;
        std::vector<int> m_zeroTrace;
        EnthalpyAdmissibilityOptions m_options;
        std::map<std::pair<int,int>,std::unique_ptr<TensorBernsteinPolynomial<3>>> m_bases;
        std::vector<ElementData> m_elements;
        mutable mfem::Vector m_cachedAnchor;
        mutable PolynomialSignReport m_cachedAnchorReport;
        mutable bool m_hasCachedAnchor{false};
    };
}