From 273cf707ef21d89ef826fdfcf00719c1656dddd6 Mon Sep 17 00:00:00 2001 From: gjj Date: Fri, 30 Aug 2024 20:41:55 +0800 Subject: [PATCH] subdivision algorithm design --- algoim/bernstein.hpp | 2 +- algoim/organizer/organizer.hpp | 161 ++++++++++++++++++++++++++++++++- algoim/organizer/primitive.hpp | 4 +- algoim/xarray.hpp | 1 + 4 files changed, 161 insertions(+), 7 deletions(-) diff --git a/algoim/bernstein.hpp b/algoim/bernstein.hpp index bd9e7b1..f760c09 100644 --- a/algoim/bernstein.hpp +++ b/algoim/bernstein.hpp @@ -322,7 +322,7 @@ template void deCasteljau(const xarray& alpha, const uvector& a, const uvector& b, xarray& out) { assert(all(out.ext() == alpha.ext())); - out = alpha; + out = alpha; // 这里是一个深拷贝，alpha不会被修改 deCasteljau(out, a.data(), b.data()); } diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index ae363aa..d6ad4f6 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -1,5 +1,6 @@ #include #include +#include #include #include @@ -28,6 +29,11 @@ namespace algoim::Organizer { +namespace detail +{ +void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {} +} // namespace detail + bool keepQuadraturePoint(std::vector& originTensors, const uvector3& originPt) { // TODO: using blobtree to filter tensors @@ -40,7 +46,22 @@ bool keepQuadraturePoint(std::vector& originTensors, const uvector3& or // std::vector> primitives; // BasicTask(std::vector> ps) {}; - +void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) +{ + assert(0 <= k && k < 3); + assert(all(res.ext() == remove_component(phi.ext(), k))); + assert(0 < facePos && facePos < 1); + real* basisAlongK; + int P = phi.ext(k); + algoim_spark_alloc(real, &basisAlongK, P); + bernstein::evalBernsteinBasis(facePos, P, basisAlongK); + + for (auto i = res.loop(); ~i; ++i) { + real evalAlongKAtFacePos = 0; + for (int j = 0; j < P; ++j) { evalAlongKAtFacePos += basisAlongK[j] * phi.m(add_component(i(), k, j)); } + res.l(i) = evalAlongKAtFacePos; + } +} void basicTask(const std::shared_ptr& p, int q = 20, real xmin = -1, real xmax = 1) { @@ -179,16 +200,148 @@ void quadratureTask(const Scene& scene) {} void buildOctree(const Scene& scene, std::vector& nodes, const uvector3& min, const uvector3& max) {} +template +int numOfZero(const uvector& v) +{ + int res = 0; + for (int i = 0; i < N; ++i) { + if (v(i) == 0) { res++; } + } + return res; +} + +template +int binaryToDecimal(const uvector& v, int zeroRep = 0) +{ + int res = 0; + int base = 1; + for (int i = N - 1; i >= 0; --i) { + if (v(i) != zeroRep) { res += base; } + base *= 2; + } + return res; +} + +template +int replaceFirst(uvector& v, int oldVal, int newVal) +{ + for (int i = 0; i < N; ++i) { + if (v(i) == oldVal) { + v(i) = newVal; + return i; + } + } + return -1; +} + void build(const Scene& scene, std::vector& nodes, int nowNodeIdx, const uvector3& nowNodeMin, const uvector3& nowNodeMax) { const std::vector& polyIntersectIndices = nodes[nowNodeIdx].polyIntersectIndices; if (polyIntersectIndices.size() == 1) { // TODO 相交很少时 } - buildOctree(scene, nodes, nowNodeMin, nowNodeMax); - nodes.resize(nodes.size() + 8); + int lastIdx = nodes.size(); + nodes.resize(lastIdx + 8); + uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; for (int i = 0; i < polyIntersectIndices.size(); ++i) { - for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { restrictToFace } + auto poly = scene.polys[polyIntersectIndices[i]]; + uvector mark(0, 0, 0); + for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { + real centerPlane = nodeMid(faceAxis); + tensor2 restrictToCenterPlane(nullptr, remove_component(poly.ext(), faceAxis)); + algoim_spark_alloc(real, restrictToCenterPlane); + restrictToInnerFace(poly, faceAxis, centerPlane, restrictToCenterPlane); + if (bernstein::uniformSign(restrictToCenterPlane) == -1) { + // primitive contain the conterface + mark(faceAxis) = 2; + } else { + // primitive is on one side of the centerface + // deCasteljau to transformation + uvector3 halfCellMin = 0, halfCellMax = 1; + halfCellMax(faceAxis) = 0.5; + tensor3 negativeHalfCell(nullptr, poly.ext()); + algoim_spark_alloc(real, negativeHalfCell); + bernstein::deCasteljau(poly, halfCellMin, halfCellMax, negativeHalfCell); + if (bernstein::uniformSign(negativeHalfCell) == -1) { + // 不在负半空间，此时一定在正半空间，可以assert一下正半空间的sign，不能为-1 + mark(faceAxis) = 1; + } else { + mark(faceAxis) = -1; + } + } + } + if (any(mark == uvector3(2, 2, 2))) { + int subIdx = 0; + for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { + nodes[lastIdx + subIdx].polyIntersectIndices.emplace_back(poly); + nodes[lastIdx + subIdx].min = nowNodeMin; + nodes[lastIdx + subIdx].max = nowNodeMax; + for (int k = 0; k < 3; ++k) { + if (j(k) == 0) + nodes[lastIdx + subIdx].max(k) = nodeMid(k); + else + nodes[lastIdx + subIdx].min(k) = nodeMid(k); + } + } + continue; + } + int zeroCount = numOfZero<3>(mark); + if (zeroCount == 0) { + // mark has -1 or 1 only + nodes[lastIdx + binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + } else if (zeroCount == 1) { + // poly related to 2 subcells + uvector sideMark = mark; + int zeroIdx = replaceFirst(sideMark, 0, 1); + nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + sideMark(zeroIdx) = -1; + nodes[lastIdx + binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndices[i]); + } else { + // zeroCount == 2 or 3 + } } + // TODO: 考虑fully contain 问题 } }; // namespace algoim::Organizer + +// clang-format off +/** + 最一开始的cell应当包含所有的primitive + func(given cell) { + 得到8个subcell + for each primitive { + mark = [0,0,0] + for k in 0,1,2 { + 测试primitive和k轴向的centerface有无交 + if 无交 { + if primitive包裹centerface { + mark[k] = 2 + } else { + 判断primitive位于哪一侧 + if k正方向 mark[k] = 1 + else mark[k] = -1 + } + } + } + if mark 含2 { + // TODO: 可能需要优化：在极端情况下， + // 例如primitive的边刚好与centerface重合， + // 此时不是所有subcell都关联primitive + 将primitive放入所有subcell + } elif mark 没有0 { + 将primitive放入对应的subcell + } elif mark 有1个0 { + 将primitive放入两个对应的subcell + } else { // >= 2个0 + 测试每个可能有交subcell是否与primitive相交 + 相交的放入对于subcell + } + } + for each subcell { + func(subcell) + } + } + + */ + +// clang-format on \ No newline at end of file diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 642ef74..66aca3c 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -591,8 +591,8 @@ void makeCylinder(xarray& tensor, uvector3 startPt, uvector3 endPt, rea struct Scene { std::vector polys; - real* boolDescription; // blobtree - uvector3 min, max; + int* boolDescription; // blobtree + // uvector3 min, max; }; struct Node { diff --git a/algoim/xarray.hpp b/algoim/xarray.hpp index 78861ac..a7b876c 100644 --- a/algoim/xarray.hpp +++ b/algoim/xarray.hpp @@ -379,6 +379,7 @@ public: }; using tensor3 = xarray; +using tensor2 = xarray; template void xarrayInit(xarray& arr)