#pragma once #include #include #include #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>& powers, int powerIdx, uvector powerSum, real factor, xarray& 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>& powers) {} template void power2BernsteinTensor(const xarray& phiPower, xarray& 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> decompFactors(N, std::vector(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 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& scale, const uvector& bias, xarray& phiPower) { std::vector>> dimOrderExpansion; const auto& ext = phiPower.ext(); for (int dim = 0; dim < 3; ++dim) { dimOrderExpansion.push_back(std::vector>(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 real evalPower(const xarray& phi, const uvector& 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 real evalBernstein(const xarray& phi, const uvector& x) { return bernstein::evalBernsteinPoly(phi, x); } class PowerTensor : public Primitive { public: xarray tensor; SparkStack* sparkStackPtr; void print() override { std::cout << "Power" << std::endl; } real eval(uvector3 p) override { return evalPower(tensor, p); } // PowerTensor() {} PowerTensor(uvector ext_) { tensor.ext_ = ext_; tensor.data_ = nullptr; sparkStackPtr = algoim_spark_alloc_heap(real, tensor); xarrayInit(tensor); } // PowerTensor(xarray t_) : tensor(t_) {} ~PowerTensor() { if (sparkStackPtr) { sparkStackPtr = nullptr; algoim_spark_release_heap(sparkStackPtr); } } }; class PowerTensorComplex : public Primitive { public: xarray compositeTensor; // 复合后的张量 SparkStack* sparkStackPtr; // std::vector> tensors; // 原始张量 std::vector tensors; // 原始张量 void print() override { std::cout << "PowerTensorComplex" << std::endl; } // PowerTensorComplex() {} PowerTensorComplex(const std::vector>& 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& 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& pcs_) { for (auto& pc : pcs_) { for (auto& t : pc.tensors) { tensors.push_back(t); } } std::vector> 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 pcs; pcs.emplace_back(*this); pcs.emplace_back(pc); return PowerTensorComplex(pcs); } PowerTensorComplex add(const PowerTensor& pt) { std::vector 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& vertices, const std::vector>& indices) { uvector3 ext(1 + indices.size()); // PowerTensorComplex pc(ext); std::vector pts; for (const auto& index : indices) { // xarray 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 ext = 3; PowerTensor pt(ext); for (int dim = 0; dim < 3; ++dim) { uvector 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 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 compositeTensor; // 复合后的张量 std::vector> 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> controlPoints; }; class NURBSSurface : public ParametricSurface { std::vector> controlPoints; std::vector> weights; std::vector 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 controlPoints; }; class NURBSCurve : public ParametricCurve { private: std::vector controlPoints; std::vector weights; std::vector knots; int degree; }; class BRep : public Primitive { std::vector vertices; std::vector curves; std::vector 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& vs_, const std::vector& cs_, const std::vector& 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 f; public: void print() override { std::cout << "FRep" << std::endl; } real eval(uvector3 p) override { return f(p); } FRep(std::function f_) : f(f_) {} }; } // namespace algoim::Organizer