HeteQuadrature/algoim/organizer/primitive.hpp

#pragma once
#include <cassert>
#include <cmath>
#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) {}

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(const 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(const uvector3& p) override { return evalPower(tensor, p); }

//     // PowerTensor() {}

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

//     // const xarray<real, 3>& getTensor() { return tensor; }

//     ~PowerTensor() = default;
// };

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

//     // std::vector<PowerTensor> powerTensors; // 原始张量

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

//     static void compositePower(const std::vector<PowerTensor>& powers,
//                                int                             powerIdx,
//                                uvector<int, 3>                 powerSum,
//                                real                            factor,
//                                xarray<real, 3>&                res)

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

//     static 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);
//         }
//     }

//     // PowerTensorComplex() {}

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

//     PowerTensorComplex(const std::vector<PowerTensor>& pts_, xarray<real, 3>& ct_) : compositeTensor(ct_)
//     {
//         uvector3 ext(1);
//         tensors.resize(pts_.size());
//         for (int i = 0; i < pts_.size(); ++i) {
//             tensors[i]  = pts_[i].tensor;
//             ext        += pts_[i].tensor.ext() - 1;
//         }
//         assert(all(ext == compositeTensor.ext()));
//         compositePower(tensors, 0, uvector3(0), 1, compositeTensor);
//     }

//     PowerTensorComplex(const std::vector<PowerTensorComplex>& ptcs_, xarray<real, 3>& ct_) : compositeTensor(ct_)
//     {
//         std::vector<xarray<real, 3>> originCompositeTensors;
//         uvector3                     ext(1);

//         for (auto& ptc : ptcs_) {
//             for (auto& t : ptc.tensors) { tensors.emplace_back(t); }
//             originCompositeTensors.push_back(ptc.compositeTensor);
//             ext += ptc.compositeTensor.ext() - 1;
//         }
//         compositePower(originCompositeTensors, 0, uvector3(0, 0, 0), 1, compositeTensor);
//     }

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

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

// class BernsteinPrimitive : public Primitive
// {
// };

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

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

//     real eval(const 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(const 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;
//     };
// };

bool isInsidePowers(const std::vector<tensor3>& tensors, const uvector3& p)
{
    for (auto& t : tensors) {
        if (evalPower(t, p) >= 0) { return false; }
    }
    return true;
};

bool isInsideBernstein(const tensor3& t, const uvector3& p) { return evalBernstein(t, p) < 0; }

bool isInsidePower(const tensor3& t, const uvector3& p) { return evalPower(t, p) < 0; }

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

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

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

    FRep(std::function<real(uvector3)> f_) : f(f_) {}
};

enum PrimitiveType { Sphere, Cylinder, Cone, Mesh, BRep };

class PrimitiveDesc
{
public:
    // const static PrimitiveType type;
    PrimitiveDesc() = default;

    virtual void print() {} // 空定义也可以，但是一定要有定义
};

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 BRepDesc : virtual public PrimitiveDesc
{
    const static PrimitiveType     type = BRep;
    std::vector<uvector3>          vertices;
    std::vector<ParametricCurve>   curves;
    std::vector<ParametricSurface> surfaces;

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

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

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

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

class SphereDesc : virtual public PrimitiveDesc
{
public:
    const static PrimitiveType type = Sphere;
    real                       radius;
    uvector3                   center;
    uvector3                   amplitude;

    SphereDesc(real r_, const uvector3& c_, const uvector3& a_) : PrimitiveDesc(), radius(r_), center(c_), amplitude(a_) {}

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

class CylinderDesc : virtual public PrimitiveDesc
{
    const static PrimitiveType type = Cylinder;
    uvector3                   node1;
    uvector3                   node2;
    real                       radius;

    CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {}

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

class ConeDesc : virtual public PrimitiveDesc
{
    const static PrimitiveType type = Cone;
    uvector3                   node1;
    uvector3                   node2;
    real                       radius;

    ConeDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {}

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

class MeshDesc : virtual public PrimitiveDesc
{
public:
    const static PrimitiveType type = Mesh;
    std::vector<uvector3>      vertices;
    std::vector<int>           indices;
    std::vector<int>           indexInclusiveScan;

    MeshDesc(const std::vector<uvector3>& vertices_,
             const std::vector<int>&      indices_,
             const std::vector<int>&      indexInclusiveScan_)
        : PrimitiveDesc(), vertices(vertices_), indices(indices_), indexInclusiveScan(indexInclusiveScan_)
    {
    }

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

void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<real, 3>>& planeTensors)
{
    uvector3 ext(1 + mesh.indexInclusiveScan.size());
    assert(all(ext == tensor.ext()));
    assert(planeTensors.size() == mesh.indexInclusiveScan.size());
    // for (const auto& index : indices) {
    for (int i = 0; i < mesh.indexInclusiveScan.size(); ++i) {
        const int indexBeg  = i == 0 ? 0 : mesh.indexInclusiveScan[i - 1];
        const int indexSize = mesh.indexInclusiveScan[i] - indexBeg;
        assert(indexSize >= 3);
        auto& planeTensor = planeTensors[i];
        xarrayInit(planeTensor);
        auto&    vertices                 = mesh.vertices;
        auto&    indices                  = mesh.indices;
        uvector3 V01                      = vertices[indices[indexBeg + 1]] - mesh.vertices[indices[indexBeg]];
        uvector3 V02                      = vertices[indices[indexBeg + 2]] - mesh.vertices[indices[indexBeg]];
        uvector3 N                        = cross(V01, V02);
        N                                /= norm(N);
        real d                            = -dot(N, vertices[indices[indexBeg]]);
        // 法线所指方向为>0区域
        planeTensor.m(uvector3(0, 0, 0))  = d;
        planeTensor.m(uvector3(1, 0, 0))  = N(0);
        planeTensor.m(uvector3(0, 1, 0))  = N(1);
        planeTensor.m(uvector3(0, 0, 1))  = N(2);
        // test other vertices
        for (int j = indexBeg + 3; j < mesh.indexInclusiveScan[i]; ++j) {
            assert(dot(N, vertices[indices[j]]) + d < std::numeric_limits<real>::epsilon());
        }
    }
    // compositePower(planeTensors, 0, 0, 1, tensor);
    // return PowerTensorComplex(planeTensors, tensor);
};

void makeSphere(const SphereDesc& sphereDesc, xarray<real, 3>& tensor)
{
    uvector<int, 3> ext = 3;
    assert(all(ext == tensor.ext()));

    for (int dim = 0; dim < 3; ++dim) {
        uvector<int, 3> idx  = 0;
        tensor.m(idx)       += sphereDesc.center(dim) * sphereDesc.center(dim);
        idx(dim)             = 1;
        tensor.m(idx)        = 2 * sphereDesc.amplitude(dim) * (-sphereDesc.center(dim));
        idx(dim)             = 2;
        tensor.m(idx)        = sphereDesc.amplitude(dim) * sphereDesc.amplitude(dim);
    }
    tensor.m(0) -= sphereDesc.radius * sphereDesc.radius;
    // return PowerTensor(tensor);
};

void makeCylinder(xarray<real, 3>& tensor, uvector3 startPt, uvector3 endPt, real r)
{
    // PowerTensor pt;
    // TODO:
}


} // namespace algoim::Organizer