fix bug of octree construction and blobtree propagation

algoim/organizer/blobtree.hpp

@@ -65,6 +65,10 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in)
 {
     const std::size_t rootIdx = tree.structure.size() - 1;
     auto&             node    = tree.structure[nodeIdx];
+    if (nodeIdx == 4) {
+        int aaa = 1;
+        int bbb = 1;
+    }
     if (node.inOut != NODE_IN_OUT_UNKNOWN) {
         int aaa = 1;
         int bb  = 1;
@@ -73,7 +77,7 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in)
     node.inOut = in ? NODE_IN : NODE_OUT;
     for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) {
         const auto&  nowNode  = tree.structure[nowIdx];
-        unsigned int nextIdx  = isLeft(nowNode) ? nodeIdx + nowNode.ancestor : nodeIdx + 1;
+        unsigned int nextIdx  = isLeft(nowNode) ? nowIdx + nowNode.ancestor : nowIdx + 1;
         auto&        nextNode = tree.structure[nextIdx]; // father
         if (nextNode.inOut != 0) { return; }
         if (nowNode.inOut == NODE_OUT) {

algoim/organizer/organizer.hpp

@@ -108,11 +108,7 @@ AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector<int, 3> side, int
 }
 } // namespace detail
 
-// 这里blobTree是拷贝传参
+bool keepQuadraturePoint(const Scene& scene, const OcTreeNode& ocTreeNode, const uvector3& originPt)
-bool keepQuadraturePoint(const Scene&        scene,
-                         organizer::BlobTree blobTree,
-                         const OcTreeNode&   ocTreeNode,
-                         const uvector3&     originPt)
 {
     // 只需要考虑intersect polys，不用考虑fully contained polys
     const auto&       polyIntersectIndices = ocTreeNode.polyIntersectIndices;
@@ -121,7 +117,9 @@ bool keepQuadraturePoint(const Scene&        scene,
         // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt);
         primitiveInOuts[i] = isInsideBernstein(scene.tensors[polyIntersectIndices[i]], originPt);
     }
-    int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts);
+    // 这里blobTree是拷贝传参
+    auto blobTree = ocTreeNode.blobTree;
+    int  res      = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts);
     assert(res == organizer::NODE_IN || res == organizer::NODE_OUT);
     return res == organizer::NODE_IN;
 }
@@ -285,13 +283,13 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
     }
 }
 
-int deepest = 0;
+static int deepest = 0;
 
-void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves, int depth)
+void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves, int depth, int& cnt)
 {
     if (deepest < depth) {
         deepest = depth;
-        std::cout << "depth: " << depth << std::endl;
+        std::cout << "depth: " << deepest << std::endl;
     }
     if (depth == 14300) {
         int aaa = 1;
@@ -302,6 +300,11 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
         leaves.emplace_back(node);
         return;
     }
+    int d = deepest;
+    if (d == 3 && depth == 1) {
+        int aaa = 1;
+        int bbb = 1;
+    }
     // std::array<OcTreeNode, CHILD_NUM> subNodes;
     std::vector<OcTreeNode> subNodes(CHILD_NUM);
     // intermediateNodes.resize(lastIdx + 8);
@@ -311,6 +314,15 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
         subNodes[subIdx].blobTree             = node.blobTree;
         subNodes[subIdx].polyIntersectIndices = std::vector<int>();
     }
+    if (node.blobTree.structure[4].inOut != NODE_IN_OUT_UNKNOWN
+        && std::find(polyIntersectIndices.begin(), polyIntersectIndices.end(), 3) != polyIntersectIndices.end()) {
+        int aaa = 1;
+        int bbb = 1;
+    }
+    if (cnt == 1) {
+        int aaa = 1;
+        int bbb = 1;
+    }
     uvector3 nodeMid = node.aabb.center();
     for (int i = 0; i < polyIntersectIndices.size(); ++i) {
         const int   polyIntersectIndex = polyIntersectIndices[i];
@@ -322,6 +334,10 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
             //     organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
             //     continue;
             // }
+            if (subIdx == 1) {
+                int aaa = 1;
+                int bbb = 1;
+            }
             tensor3 subcellPoly(nullptr, poly.ext());
             algoim_spark_alloc(real, subcellPoly);
             bernstein::deCasteljau(poly, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly);
@@ -340,7 +356,13 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
             }
         }
     }
-    for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { buildOcTreeV0(scene, subNodes[subIdx], leaves, depth + 1); }
+    if (subNodes[1].blobTree.structure[4].inOut != NODE_IN_OUT_UNKNOWN
+        & std::find(subNodes[1].polyIntersectIndices.begin(), subNodes[1].polyIntersectIndices.end(), 3)
+              != subNodes[1].polyIntersectIndices.end()) {
+        int aaa = 1;
+        int bbb = 1;
+    }
+    for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { buildOcTreeV0(scene, subNodes[subIdx], leaves, depth + 1, ++cnt); }
 }
 
 /** single object quadrature */
@@ -434,7 +456,7 @@ void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, in
     // for (auto& p : originTensorStacks) delete p;
 };
 
-BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 10)
+BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
 {
     auto                      integrand = [](const uvector<real, 3>& x) { return 1.0; };
     real                      volume = 0., surf = 0.;
@@ -442,10 +464,10 @@ BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree,
     ImplicitPolyQuadrature<3> ipquad(scene.tensors, node.polyIntersectIndices);
     ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
         auto realX = x * range + node.aabb.min;
-        if (keepQuadraturePoint(scene, blobTree, node, realX)) { volume += w * integrand(realX); }
+        if (keepQuadraturePoint(scene, node, realX)) { volume += w * integrand(realX); }
     });
     ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, 3>& x, real w, const uvector<real, 3>& wn) {
-        surf += w * integrand(x * range + node.aabb.min);
+        volume += w * integrand(x * range + node.aabb.min);
     });
     return {volume, surf};
 }
@@ -479,7 +501,7 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
         } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
             const int faceCount = pt->indexInclusiveScan.size();
             assert(faceCount > 1);
-            visiblePrimitiveReps[i].tensors = std::vector<tensor3>(faceCount, tensor3(nullptr, 2));
+            visiblePrimitiveReps[i].tensors = std::vector<tensor3>(faceCount, tensor3(nullptr, 3));
             auto& planeTensors              = visiblePrimitiveReps[i].tensors;
             algoimSparkAllocHeapVector(tensorStacks, planeTensors);
 
@@ -566,10 +588,11 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
     Scene      scene{realPrimitives, AABB(xmin, xmax)};
     OcTreeNode rootNode(0, 1, blobTree);
     for (int i = 0; i < realPrimitives.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
-    buildOcTreeV0(scene, rootNode, leaves, 1);
+    int cnt = 1;
+    buildOcTreeV0(scene, rootNode, leaves, 1, cnt);
 
     for (const auto& leaf : leaves) {
-        auto basicRes  = basicTask(scene, blobTree, leaf, q);
+        auto basicRes  = basicTask(scene, leaf, q);
         volume        += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min);
     }
 

algoim/organizer/primitive.hpp

@@ -699,12 +699,17 @@ struct VisiblePrimitiveRep {
     organizer::BlobTree  subBlobTree;
 };
 
+struct MinimalPrimitiveRep {
+    tensor3 tensor;
+    AABB    aabb;
+};
+
 struct Scene {
     // std::vector<CompleteTensorRep> polys;
     // std::vector<VisiblePrimitiveRep> visiblePrimitives;
     // organizer::BlobTree            blobTree;
-    std::vector<tensor3> tensors;
+    std::vector<MinimalPrimitiveRep> tensors;
-    AABB                 boundary;
+    AABB                             boundary;
 };
 
 void upwardGeneratingNodes(std::vector<int>&    levelLeftBook,

algoim/quadrature_multipoly.hpp

@@ -281,6 +281,9 @@ bool resultant_core(const xarray<real, N>& p, const xarray<real, N>* q, int k, x
     xarray<real, 2> mat(nullptr, uvector<int, 2>{M, M});
 
     /**gjj***/
+    // 当Q == 0，应该是对p求P类点，且p完全垂直降维方向的情形，这样的R类点可以直接丢弃
+    // 因为产生的R类点和Q类点是重合的（平行降维方向的多项式一定和上下边界交于同一位置，其降维后也是R类点降维后的结果）
+    // 当P和Q都<=1，要么是Q == 0，要么是两个垂直降维方向的多项式求R类点，此时交线退化为点，也丢弃
     if (Q == 0 || (P <= 1 && Q <= 1)) { return false; }
     if (!(P >= 2 && Q >= 1 && M >= 1)) {
         int aaa = 1;

gjj/primitiveDebug.hpp

@@ -140,8 +140,8 @@ void caseScene()
     const int                                   PRIMITIVE_CNT = 6;
     std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions(PRIMITIVE_CNT);
     primitiveDescriptions[0]                    = std::make_shared<CuboidDesc>(CuboidDesc(0., 1.6));
-    // primitiveDescriptions[1]                    = std::make_shared<CuboidDesc>(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.));
+    primitiveDescriptions[1]                    = std::make_shared<CuboidDesc>(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.));
-    primitiveDescriptions[1]                    = std::make_shared<SphereDesc>(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.));
+    // primitiveDescriptions[1]                    = std::make_shared<SphereDesc>(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.));
     primitiveDescriptions[2]                    = std::make_shared<SphereDesc>(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.));
     primitiveDescriptions[3]                    = std::make_shared<SphereDesc>(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.));
     std::vector<uvector3> pyramidBottomVertices = {