subdivision algorithm design

algoim/bernstein.hpp

@@ -322,7 +322,7 @@ template <int N>
 void deCasteljau(const xarray<real, N>& alpha, const uvector<real, N>& a, const uvector<real, N>& b, xarray<real, N>& out)
 {
     assert(all(out.ext() == alpha.ext()));
-    out = alpha;
+    out = alpha; // 这里是一个深拷贝，alpha不会被修改
     deCasteljau(out, a.data(), b.data());
 }
 

algoim/organizer/organizer.hpp

@@ -1,5 +1,6 @@
 #include <array>
 #include <bitset>
+#include <cassert>
 #include <iostream>
 #include <booluarray.hpp>
 
@@ -28,6 +29,11 @@
 namespace algoim::Organizer
 {
 
+namespace detail
+{
+void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {}
+} // namespace detail
+
 bool keepQuadraturePoint(std::vector<tensor3>& originTensors, const uvector3& originPt)
 {
     // TODO: using blobtree to filter tensors
@@ -40,7 +46,22 @@ bool keepQuadraturePoint(std::vector<tensor3>& originTensors, const uvector3& or
 // std::vector<std::shared_ptr<PrimitiveDesc>> primitives;
 
 // BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> 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<PrimitiveDesc>& 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<Node>& nodes, const uvector3& min, const uvector3& max) {}
 
+template <int N>
+int numOfZero(const uvector<int, N>& v)
+{
+    int res = 0;
+    for (int i = 0; i < N; ++i) {
+        if (v(i) == 0) { res++; }
+    }
+    return res;
+}
+
+template <int N>
+int binaryToDecimal(const uvector<int, N>& 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 N>
+int replaceFirst(uvector<int, N>& 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<Node>& nodes, int nowNodeIdx, const uvector3& nowNodeMin, const uvector3& nowNodeMax)
 {
     const std::vector<int>& polyIntersectIndices = nodes[nowNodeIdx].polyIntersectIndices;
     if (polyIntersectIndices.size() == 1) {
         // TODO 相交很少时
     }
-    buildOctree(scene, nodes, nowNodeMin, nowNodeMax);
+    int lastIdx = nodes.size();
-    nodes.resize(nodes.size() + 8);
+    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<int, 3> 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<int, 3> 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

algoim/organizer/primitive.hpp

@@ -591,8 +591,8 @@ void makeCylinder(xarray<real, 3>& tensor, uvector3 startPt, uvector3 endPt, rea
 
 struct Scene {
     std::vector<tensor3> polys;
-    real*                boolDescription; // blobtree
+    int*                 boolDescription; // blobtree
-    uvector3             min, max;
+    // uvector3             min, max;
 };
 
 struct Node {

algoim/xarray.hpp

@@ -379,6 +379,7 @@ public:
 };
 
 using tensor3 = xarray<real, 3>;
+using tensor2 = xarray<real, 2>;
 
 template <typename T, int N>
 void xarrayInit(xarray<T, N>& arr)