7 changed files with 638 additions and 25 deletions
@ -0,0 +1,89 @@ |
|||
#include <array> |
|||
#include <bitset> |
|||
#include <iostream> |
|||
#include <booluarray.hpp> |
|||
|
|||
|
|||
#include <cstddef> |
|||
#include <iostream> |
|||
#include <iomanip> |
|||
#include <fstream> |
|||
#include <vector> |
|||
#include "bernstein.hpp" |
|||
#include "multiloop.hpp" |
|||
#include "quadrature_multipoly.hpp" |
|||
#include "binomial.hpp" |
|||
|
|||
#include "real.hpp" |
|||
#include "uvector.hpp" |
|||
#include "vector" |
|||
#include "xarray.hpp" |
|||
#include <chrono> |
|||
#include <cmath> |
|||
#include <memory> |
|||
|
|||
#include "organizer/primitive.hpp" |
|||
|
|||
namespace algoim::Organizer |
|||
{ |
|||
class BasicTask |
|||
{ |
|||
public: |
|||
std::vector<BernsteinPrimitive> primitives; |
|||
|
|||
BasicTask(std::vector<std::shared_ptr<Primitive>> ps) {}; |
|||
|
|||
BasicTask(std::shared_ptr<Primitive> p) |
|||
{ |
|||
int q = 20; |
|||
real volume = 0; |
|||
real xmin = -1; |
|||
real xmax = 1; |
|||
auto integrand = [](const uvector<real, 3>& x) { return 1.0; }; |
|||
uvector3 range = xmax - xmin; |
|||
|
|||
|
|||
if (auto pt = std::dynamic_pointer_cast<PowerTensor>(p)) { |
|||
detail::powerTransformation(range, xmin, pt->tensor); |
|||
auto primitive = BernsteinTensor(*pt); |
|||
|
|||
uvector<real, 3> testX(0., 0., 0.25); |
|||
real testEvalBernstein = bernstein::evalBernsteinPoly(primitive.tensor, testX); |
|||
// auto vec1 = xarray2StdVector(phi);
|
|||
std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl; |
|||
|
|||
ImplicitPolyQuadrature<3> ipquad(primitive.tensor); |
|||
ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) { |
|||
if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin)); |
|||
}); |
|||
} else if (auto pc = std::dynamic_pointer_cast<PowerTensorComplex>(p)) { |
|||
detail::powerTransformation(range, xmin, pc->compositeTensor); |
|||
auto primitive = BernsteinTensorComplex(*pc); |
|||
ImplicitPolyQuadrature<3> ipquad(primitive.compositeTensor); |
|||
ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) { |
|||
if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin)); |
|||
}); |
|||
} |
|||
volume *= pow(xmax - xmin, 3); |
|||
std::cout << "Volume xxx: " << volume << std::endl; |
|||
}; |
|||
|
|||
// BasicTask(std::shared_ptr<PowerTensorComplex> pc)
|
|||
// {
|
|||
// int q = 10;
|
|||
// real volume;
|
|||
// uvector3 xmin = 0;
|
|||
// uvector3 xmax = 1;
|
|||
// auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
|
|||
// uvector3 range = xmax - xmin;
|
|||
|
|||
// detail::powerTransformation(range, xmin, pc->compositeTensor);
|
|||
// auto primitive = BernsteinTensorComplex(*pc);
|
|||
// ImplicitPolyQuadrature<3> ipquad(primitive.compositeTensor);
|
|||
// ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
|
|||
// if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin));
|
|||
// });
|
|||
// std::cout << "Volume xxx: " << volume << std::endl;
|
|||
// };
|
|||
}; |
|||
}; // namespace algoim::Organizer
|
@ -0,0 +1,415 @@ |
|||
#pragma once |
|||
#include <cassert> |
|||
#include <vector> |
|||
#include "iostream" |
|||
#include "multiloop.hpp" |
|||
#include "polyset.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; |
|||
|
|||
void print() override { std::cout << "Power" << std::endl; } |
|||
|
|||
real eval(uvector3 p) override { return evalPower(tensor, p); } |
|||
|
|||
PowerTensor() {} |
|||
|
|||
PowerTensor(xarray<real, 3> t_) : tensor(t_) {} |
|||
}; |
|||
|
|||
class PowerTensorComplex : public Primitive |
|||
{ |
|||
public: |
|||
xarray<real, 3> compositeTensor; // 复合后的张量
|
|||
std::vector<xarray<real, 3>> 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; |
|||
algoim_spark_alloc(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; |
|||
algoim_spark_alloc(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) |
|||
{ |
|||
PowerTensorComplex pc; |
|||
uvector3 ext(1, 1, 1); |
|||
for (const auto& index : indices) { |
|||
xarray<real, 3> tensor(nullptr, uvector3(2)); // 最高1次
|
|||
algoim_spark_alloc(real, tensor); |
|||
xarrayInit(tensor); |
|||
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区域
|
|||
tensor.m(uvector3(0, 0, 0)) = d; |
|||
tensor.m(uvector3(1, 0, 0)) = N(0); |
|||
tensor.m(uvector3(0, 1, 0)) = N(1); |
|||
tensor.m(uvector3(0, 0, 1)) = N(2); |
|||
pc.tensors.push_back(tensor); |
|||
ext += tensor.ext() - 1; |
|||
} |
|||
|
|||
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) |
|||
{ |
|||
PowerTensor pt; |
|||
uvector<int, 3> ext = 3; |
|||
pt.tensor.ext_ = ext; |
|||
algoim_spark_alloc(real, pt.tensor); |
|||
xarrayInit(pt.tensor); |
|||
|
|||
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 pt; |
|||
} |
|||
|
|||
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
|
@ -0,0 +1,74 @@ |
|||
#include <array> |
|||
#include <bitset> |
|||
#include <iostream> |
|||
#include <booluarray.hpp> |
|||
|
|||
|
|||
#include <cstddef> |
|||
#include <iostream> |
|||
#include <iomanip> |
|||
#include <fstream> |
|||
#include <vector> |
|||
#include "bernstein.hpp" |
|||
#include "multiloop.hpp" |
|||
#include "quadrature_multipoly.hpp" |
|||
#include "binomial.hpp" |
|||
|
|||
#include "real.hpp" |
|||
#include "uvector.hpp" |
|||
#include "vector" |
|||
#include "xarray.hpp" |
|||
#include <chrono> |
|||
#include <memory> |
|||
|
|||
#include "organizer/primitive.hpp" |
|||
#include "organizer/organizer.hpp" |
|||
|
|||
|
|||
using namespace algoim::Organizer; |
|||
using namespace algoim; |
|||
|
|||
void case0() |
|||
{ |
|||
std::vector<std::shared_ptr<Primitive>> ps; |
|||
|
|||
// mesh
|
|||
std::vector<uvector3> vertices = {uvector3(-0.8, -0.8, -0.8), |
|||
uvector3(-0.8, -0.8, 0.8), |
|||
uvector3(-0.8, 0.8, -0.8), |
|||
uvector3(-0.8, 0.8, 0.8), |
|||
uvector3(0.8, -0.8, -0.8), |
|||
uvector3(0.8, -0.8, 0.8), |
|||
uvector3(0.8, 0.8, -0.8), |
|||
uvector3(0.8, 0.8, 0.8) |
|||
|
|||
}; |
|||
std::vector<uvector<int, 3>> indixes = { |
|||
uvector<int, 3>(2, 1, 0), |
|||
uvector<int, 3>(1, 2, 3), |
|||
uvector<int, 3>(0, 1, 4), |
|||
uvector<int, 3>(4, 1, 5), |
|||
uvector<int, 3>(6, 2, 0), |
|||
uvector<int, 3>(6, 0, 4), |
|||
uvector<int, 3>(0, 3, 2), |
|||
uvector<int, 3>(0, 2, 6), |
|||
uvector<int, 3>(1, 3, 7), |
|||
uvector<int, 3>(1, 7, 5), |
|||
uvector<int, 3>(5, 7, 6), |
|||
uvector<int, 3>(4, 5, 6), |
|||
}; |
|||
|
|||
// ps.emplace_back(std::make_shared<PowerTensorComplex>(makeMesh(vertices, indixes)));
|
|||
// ps.emplace_back(std::make_shared<PowerTensor>(makeSphere(0.2)));
|
|||
auto phi0 = std::make_shared<PowerTensorComplex>(makeMesh(vertices, indixes)); |
|||
phi0->add(*std::make_shared<PowerTensor>(makeSphere(0.2))); |
|||
auto basicTask = BasicTask(phi0); |
|||
} |
|||
|
|||
void case1() |
|||
{ |
|||
auto phi0 = std::make_shared<PowerTensor>(makeSphere(0.2)); |
|||
auto basicTask = BasicTask(phi0); |
|||
} |
|||
|
|||
void testPrimitive() { case1(); } |
Loading…
Reference in new issue