HeteQuadrature/algoim/organizer/primitive.hpp

#pragma once
#include <cassert>
#include <cstddef>
#include <vector>
#include "iostream"
#include "multiloop.hpp"
#include "polyset.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp"
#include "real.hpp"
#include "xarray.hpp"
#include "binomial.hpp"
#include "bernstein.hpp"

namespace algoim::Organizer
{

namespace detail
{
void compositePower(const std::vector<xarray<real, 3>>& powers,
                    int                                 powerIdx,
                    uvector<int, 3>                     powerSum,
                    real                                factor,
                    xarray<real, 3>&                    res)

{
    if (powerIdx == 0) {
        {
            uvector3 ext(1, 1, 1);
            for (auto& t : powers) { ext += t.ext() - 1; }
            assert(all(ext == res.ext()));
        }
        xarrayInit(res);
    }
    if (powerIdx == powers.size()) {
        res.m(powerSum) += factor;
        return;
    }
    auto& power = powers[powerIdx];
    for (auto i = power.loop(); ~i; ++i) {
        if (power.l(i) == 0) {
            factor = 0;
            continue;
        }
        compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
    }
    // int a = 1;
}

// void compositePower(const std::vector<xarray<real, 3>>& powers) {}

template <int N>
void power2BernsteinTensor(const xarray<real, N>& phiPower, xarray<real, N>& phiBernsetin)
{
    xarrayInit(phiBernsetin);
    for (auto i = phiPower.loop(); ~i; ++i) {
        // phi.l(i) = powerFactors.l(i);
        real factorBase = phiPower.l(i);
        if (factorBase == 0) continue;
        auto                           traverseRange = phiPower.ext() - i();
        std::vector<std::vector<real>> decompFactors(N, std::vector<real>(max(traverseRange), 0.));
        for (int dim = 0; dim < N; ++dim) {
            // Sigma
            size_t      nDim      = phiPower.ext()(dim) - 1;
            const real* binomNDim = Binomial::row(nDim);
            for (int j = i(dim); j <= nDim; ++j) {
                const real* binomJ             = Binomial::row(j);
                decompFactors[dim][j - i(dim)] = binomJ[i(dim)] / binomNDim[i(dim)];
            }
        }
        xarray<real, N> subgrid(nullptr, traverseRange);
        // algoim_spark_alloc(real, subgrid);

        for (auto ii = subgrid.loop(); ~ii; ++ii) {
            real factor = factorBase;
            for (int dim = 0; dim < N; ++dim) { factor *= decompFactors[dim][ii(dim)]; }
            phiBernsetin.m(i() + ii()) += factor;
        }
    }
}

void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
{
    std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
    const auto&                                 ext = phiPower.ext();
    for (int dim = 0; dim < 3; ++dim) {
        dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
        for (int degree = 0; degree < ext(dim); ++degree) {
            const real* binomDegree = Binomial::row(degree);
            dimOrderExpansion[dim][degree].reserve(degree + 1);
            // 根据二项定理展开
            for (int i = 0; i <= degree; ++i) {
                dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
            }
        }
    }
    for (auto i = phiPower.loop(); ~i; ++i) {
        // 迭代器必须按照坐标的升序进行访问，即，访问ijk时，(i-1)jk，i(j-1)k，ij(k-1)必须已经访问过
        real base     = phiPower.l(i);
        phiPower.l(i) = 0;
        auto exps     = i();
        for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
            real item = base;
            for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
            phiPower.m(j()) += item;
        }
    }
}
} // namespace detail

class Primitive
{
public:
    virtual void print() = 0;

    virtual real eval(uvector3) { return 0; }
};

template <int N>
real evalPower(const xarray<real, N>& phi, const uvector<real, N>& x)
{
    real res = 0;
    for (auto i = phi.loop(); ~i; ++i) {
        real item = phi.l(i);
        auto exps = i();
        for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), exps(dim)); }
        res += item;
    }
    return res;
}

template <int N>
real evalBernstein(const xarray<real, N>& phi, const uvector<real, N>& x)
{
    return bernstein::evalBernsteinPoly(phi, x);
}

class PowerTensor : public Primitive
{
public:
    xarray<real, 3> tensor;

    // SparkStack<real>* sparkStackPtr;

    void print() override { std::cout << "Power" << std::endl; }

    real eval(uvector3 p) override { return evalPower(tensor, p); }

    // PowerTensor() {}

    PowerTensor(uvector<int, 3> ext_)
    {
        tensor.ext_  = ext_;
        tensor.data_ = nullptr;
        algoim_spark_alloc(real, tensor);
        // sparkStackPtr = algoim_spark_alloc_heap(real, tensor);
        xarrayInit(tensor);
    }

    // PowerTensor(xarray<real, 3> t_) : tensor(t_) {}

    ~PowerTensor()
    {
        // if (sparkStackPtr) {
        //     algoim_spark_release_heap(sparkStackPtr);
        //     sparkStackPtr = nullptr;
        // }
    }
};

class PowerTensorComplex : public Primitive
{
public:
    xarray<real, 3>          compositeTensor; // 复合后的张量
    SparkStack<real>*        sparkStackPtr;
    // std::vector<xarray<real, 3>> tensors;         // 原始张量
    std::vector<PowerTensor> tensors; // 原始张量

    void print() override { std::cout << "PowerTensorComplex" << std::endl; }

    // PowerTensorComplex() {}

    // PowerTensorComplex(const std::vector<xarray<real, 3>>& ts_) : tensors(ts_)
    // {
    //     uvector3 ext(1, 1, 1);
    //     for (auto& t : ts_) { ext += t.ext() - 1; }
    //     compositeTensor.ext_ = ext;
    //     sparkStackPtr        = algoim_spark_alloc_heap(real, compositeTensor);
    //     // detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
    // }

    PowerTensorComplex(const std::vector<PowerTensor>& pts_)
    {
        uvector3 ext(1);
        for (auto& pt : pts_) { ext += pt.tensor.ext() - 1; }
        compositeTensor.ext_ = ext;
        sparkStackPtr        = algoim_spark_alloc_heap(real, compositeTensor);
        // detail::compositePower(tensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
    }

    PowerTensorComplex(const std::vector<PowerTensorComplex>& pcs_)
    {
        for (auto& pc : pcs_) {
            for (auto& t : pc.tensors) { tensors.push_back(t); }
        }
        std::vector<xarray<real, 3>> originCompositeTensors;
        uvector3                     ext(1, 1, 1);
        for (auto& pc : pcs_) {
            originCompositeTensors.push_back(pc.compositeTensor);
            ext += pc.compositeTensor.ext() - 1;
        }
        compositeTensor.ext_ = ext;
        sparkStackPtr        = algoim_spark_alloc_heap(real, compositeTensor);
        detail::compositePower(originCompositeTensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
    }

    PowerTensorComplex add(const PowerTensorComplex& pc)
    {
        std::vector<PowerTensorComplex> pcs;
        pcs.emplace_back(*this);
        pcs.emplace_back(pc);
        return PowerTensorComplex(pcs);
    }

    PowerTensorComplex add(const PowerTensor& pt)
    {
        std::vector<PowerTensorComplex> pcs;
        pcs.emplace_back(*this);
        // pcs.emplace_back(PowerTensorComplex({pt.tensor}));
        return PowerTensorComplex(pcs);
    }

    real eval(uvector3 p) override { return evalPower(compositeTensor, p); }

    bool isInside(uvector3 p)
    {
        for (auto& t : tensors) {
            // if (evalPower(t, p) >= 0) { return false; }
        }
        return true;
    }
};

PowerTensorComplex makeMesh(const std::vector<uvector3>& vertices, const std::vector<uvector<int, 3>>& indices)
{
    uvector3                 ext(1 + indices.size());
    // PowerTensorComplex pc(ext);
    std::vector<PowerTensor> pts;
    for (const auto& index : indices) {
        // xarray<real, 3> tensor(nullptr, ); // 最高1次
        // algoim_spark_alloc(real, tensor);
        PowerTensor pt(uvector3(2));

        uvector3 V01                    = vertices[index(1)] - vertices[index(0)];
        uvector3 V02                    = vertices[index(2)] - vertices[index(0)];
        uvector3 N                      = cross(V01, V02);
        N                              /= norm(N);
        real d                          = -dot(N, vertices[index(0)]);
        // 法线所指方向为>0区域
        pt.tensor.m(uvector3(0, 0, 0))  = d;
        pt.tensor.m(uvector3(1, 0, 0))  = N(0);
        pt.tensor.m(uvector3(0, 1, 0))  = N(1);
        pt.tensor.m(uvector3(0, 0, 1))  = N(2);
        pts.push_back(pt);
    }
    PowerTensorComplex pc(pts);
    pc.compositeTensor.ext_ = ext;
    algoim_spark_alloc(real, pc.compositeTensor);
    // detail::compositePower(pc.tensors, 0, 0, 1, pc.compositeTensor);
    return pc;
}

PowerTensor makeSphere(const real r, const uvector3& c = 0, const uvector3& a = 1)
{
    uvector<int, 3> ext = 3;
    PowerTensor     pt(ext);

    for (int dim = 0; dim < 3; ++dim) {
        uvector<int, 3> idx  = 0;
        pt.tensor.m(idx)    += c(dim) * c(dim);
        idx(dim)             = 1;
        pt.tensor.m(idx)     = 2 * a(dim) * (-c(dim));
        idx(dim)             = 2;
        pt.tensor.m(idx)     = a(dim) * a(dim);
    }
    pt.tensor.m(0) -= r * r;
    return pt;
}

PowerTensor makeCylinder(uvector3 startPt, uvector3 endPt, real r)
{
    // PowerTensor pt;
    // TODO:
    return PowerTensor({});
}

class BernsteinPrimitive : public Primitive
{
};

class BernsteinTensor : public BernsteinPrimitive
{
public:
    xarray<real, 3> tensor;

    void print() override { std::cout << "Bernstein" << std::endl; }

    real eval(uvector3 p) override { return evalBernstein(tensor, p); }

    BernsteinTensor(const PowerTensor& pt_)
    {
        auto v0 = xarray2StdVector(pt_.tensor);

        tensor.ext_ = pt_.tensor.ext();
        algoim_spark_alloc(real, tensor);
        auto v1 = xarray2StdVector(pt_.tensor);
        auto v2 = xarray2StdVector(tensor);
        detail::power2BernsteinTensor(pt_.tensor, tensor);

        uvector3 x(0.2, 0.5, 0.6);
        real     BernsteinValue = bernstein::evalBernsteinPoly(tensor, x);
        real     PowerValue     = evalPower(pt_.tensor, x);
        int      a              = 1;
    }

    bool isInside(uvector3 p) { return eval(p) < 0; }
};

class BernsteinTensorComplex : public BernsteinPrimitive
{
public:
    bool                         fromPower;
    xarray<real, 3>              compositeTensor; // 复合后的张量
    std::vector<xarray<real, 3>> tensors;         // 原始张量

    void print() override { std::cout << "Bernstein Complex" << std::endl; }

    real eval(uvector3 p) override { return evalBernstein(compositeTensor, p); }

    // BernsteinTensorComplex(const PowerTensorComplex& pc_) : fromPower(true), tensors(pc_.tensors)
    // {
    //     compositeTensor.ext_ = pc_.compositeTensor.ext();
    //     algoim_spark_alloc(real, compositeTensor);
    //     detail::power2BernsteinTensor(pc_.compositeTensor, compositeTensor);
    // };

    bool isInside(uvector3 p)
    {
        if (fromPower) {
            for (auto& t : tensors) {
                if (evalPower(t, p) >= 0) { return false; }
            }
            return true;
        } else {
            for (auto& t : tensors) {
                if (eval(p) >= 0) { return false; }
            }
            return true;
        }
        return true;
    };
};

class ParametricSurface
{
public:
    virtual uvector3 eval(uvector2 p) { return uvector3(0, 0, 0); }
};

class BezierSurface : public ParametricSurface
{
private:
    std::vector<std::vector<uvector3>> controlPoints;
};

class NURBSSurface : public ParametricSurface
{
    std::vector<std::vector<uvector3>> controlPoints;
    std::vector<std::vector<real>>     weights;
    std::vector<real>                  knotsU, knotsV;
    int                                degreeU, degreeV;
};

class ParametricCurve
{
public:
    virtual uvector3 eval(real p) { return uvector3(0, 0, 0); }
};

class BezierCurve : public ParametricCurve
{
private:
    std::vector<uvector3> controlPoints;
};

class NURBSCurve : public ParametricCurve
{
private:
    std::vector<uvector3> controlPoints;
    std::vector<real>     weights;
    std::vector<real>     knots;
    int                   degree;
};

class BRep : public Primitive
{
    std::vector<uvector3>          vertices;
    std::vector<ParametricCurve>   curves;
    std::vector<ParametricSurface> surfaces;

public:
    void print() override { std::cout << "FRep" << std::endl; }

    real eval(uvector3 p) override
    {
        // DOTO: the implicit conversion of Parametric BRep
        return 0;
    }

    BRep(const std::vector<uvector3>& vs_, const std::vector<ParametricCurve>& cs_, const std::vector<ParametricSurface>& ss_)
        : vertices(vs_), curves(cs_), surfaces(ss_)
    {
    }

    BRep(const BRep& brep)
    {
        vertices = brep.vertices;
        curves   = brep.curves;
        surfaces = brep.surfaces;
    }
};

class FRep : public Primitive
{
private:
    std::function<real(uvector3)> f;

public:
    void print() override { std::cout << "FRep" << std::endl; }

    real eval(uvector3 p) override { return f(p); }

    FRep(std::function<real(uvector3)> f_) : f(f_) {}
};
} // namespace algoim::Organizer