solve the inf of volume

algoim/organizer/organizer.hpp

@@ -108,14 +108,19 @@ AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector<int, 3> side, int
 }
 } // namespace detail
 
-bool keepQuadraturePoint(const Scene& scene, const OcTreeNode& ocTreeNode, const uvector3& originPt)
+bool keepQuadraturePoint(const std::vector<tensor3>& tensors, const OcTreeNode& ocTreeNode, const uvector3& point)
 {
     // 只需要考虑intersect polys，不用考虑fully contained polys
-    const auto&       polyIntersectIndices = ocTreeNode.polyIntersectIndices;
+    const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices;
+    assert(tensors.size() == polyIntersectIndices.size());
+    if (polyIntersectIndices.size() == 0) {
+        assert(ocTreeNode.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
+        return ocTreeNode.blobTree.structure.back().inOut;
+    }
     std::vector<bool> primitiveInOuts(polyIntersectIndices.size());
     for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) {
         // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt);
-        primitiveInOuts[i] = isInsideBernstein(scene.tensors[polyIntersectIndices[i]], originPt);
+        primitiveInOuts[i] = isInsideBernstein(tensors[i], point);
     }
     // 这里blobTree是拷贝传参
     auto blobTree = ocTreeNode.blobTree;
@@ -188,7 +193,7 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
     uvector3 nodeMid = node.aabb.center();
     for (int i = 0; i < polyIntersectIndices.size(); ++i) {
         const int       polyIntersectIndex = polyIntersectIndices[i];
-        const auto&     poly               = scene.tensors[polyIntersectIndex];
+        const auto&     poly               = scene.minimalReps[polyIntersectIndex].tensor;
         uvector<int, 3> mark(0, 0, 0);
         for (int faceAxis = 0; faceAxis < 3; ++faceAxis) {
             real    centerPlane = nodeMid(faceAxis);
@@ -296,7 +301,7 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
         int bbb = 1;
     }
     const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
-    if (polyIntersectIndices.size() <= 5) {
+    if (polyIntersectIndices.size() <= 3) {
         leaves.emplace_back(node);
         return;
     }
@@ -326,21 +331,21 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
     uvector3 nodeMid = node.aabb.center();
     for (int i = 0; i < polyIntersectIndices.size(); ++i) {
         const int   polyIntersectIndex = polyIntersectIndices[i];
-        const auto& poly               = scene.tensors[polyIntersectIndex];
+        const auto& minimalRep         = scene.minimalReps[polyIntersectIndex];
         subIdx                         = 0;
         for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
-            // if (!poly.aabb.intersect(subNodes[subIdx].aabb)) {
-            //     // out of the subcell
-            //     organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
-            //     continue;
-            // }
+            if (!minimalRep.aabb.intersect(subNodes[subIdx].aabb)) {
+                // out of the subcell
+                organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
+                continue;
+            }
             if (subIdx == 1) {
                 int aaa = 1;
                 int bbb = 1;
             }
-            tensor3 subcellPoly(nullptr, poly.ext());
+            tensor3 subcellPoly(nullptr, minimalRep.tensor.ext());
             algoim_spark_alloc(real, subcellPoly);
-            bernstein::deCasteljau(poly, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly);
+            bernstein::deCasteljau(minimalRep.tensor, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly);
             int sign = bernstein::uniformSign(subcellPoly);
             if (sign == 1) {
                 organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
@@ -458,18 +463,36 @@ void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, in
 
 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.;
-    auto                      range = node.aabb.size();
-    ImplicitPolyQuadrature<3> ipquad(scene.tensors, node.polyIntersectIndices);
+    if (node.polyIntersectIndices.size() == 0) {
+        assert(node.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN);
+        if (node.blobTree.structure.back().inOut == NODE_IN) {
+            return {node.aabb.volume(), 0};
+        } else {
+            return {0, 0};
+        }
+    }
+
+    auto                           integrand = [](const uvector<real, 3>& x) { return 1.0; };
+    real                           volume = 0., surf = 0.;
+    auto                           range = node.aabb.size();
+    std::vector<tensor3>           phis; // phis using subcell as [0,1]^3
+    std::vector<SparkStack<real>*> phiStacks;
+    for (int i = 0; i < node.polyIntersectIndices.size(); i++) {
+        const auto& polyWithinScene = scene.minimalReps[node.polyIntersectIndices[i]].tensor;
+        phis.emplace_back(tensor3(nullptr, polyWithinScene.ext()));
+        phiStacks.emplace_back(algoim_spark_alloc_heap(real, phis.back()));
+        bernstein::deCasteljau(polyWithinScene, node.aabb.min, node.aabb.max, phis.back());
+    }
+    // ImplicitPolyQuadrature<3> ipquad(scene.minimalReps, node.polyIntersectIndices);
+    ImplicitPolyQuadrature<3> ipquad(phis);
+
     ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
-        auto realX = x * range + node.aabb.min;
-        if (keepQuadraturePoint(scene, node, realX)) { volume += w * integrand(realX); }
+        auto realX = x * range + node.aabb.min; // 这里realX应该是最原始空间下的点？不过因为算体积，所以不影响
+        if (keepQuadraturePoint(phis, node, x)) { volume += w * integrand(realX); }
     });
-    ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, 3>& x, real w, const uvector<real, 3>& wn) {
-        volume += w * integrand(x * range + node.aabb.min);
-    });
-    return {volume, surf};
+
+    for (auto& p : phiStacks) delete p;
+    return {volume * node.aabb.volume(), surf};
 }
 
 void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives,
@@ -513,7 +536,7 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
             }
             visiblePrimitiveReps[i].aabb.normalize(range, xmin);
         } else if (auto pt = std::dynamic_pointer_cast<CylinderDesc>(primitives[i])) {
-            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 2), tensor3(nullptr, 2)};
+            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)};
             auto& tensors                   = visiblePrimitiveReps[i].tensors;
             tensors[0].ext_(pt->alignAxis)  = 1;
             algoimSparkAllocHeapVector(tensorStacks, tensors);
@@ -526,7 +549,7 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
             }
             visiblePrimitiveReps[i].aabb.normalize(range, xmin);
         } else if (auto pt = std::dynamic_pointer_cast<ConeDesc>(primitives[i])) {
-            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 2), tensor3(nullptr, 2)};
+            visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)};
             auto& tensors                   = visiblePrimitiveReps[i].tensors;
             algoimSparkAllocHeapVector(tensorStacks, tensors);
 
@@ -542,9 +565,9 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
         }
     }
 
-    std::vector<tensor3> realPrimitives;
+    std::vector<MinimalPrimitiveRep> minimalReps;
     /*** merge subtrees to main tree ***/
-    std::vector<int>     realLeafIndices;
+    std::vector<int>                 realLeafIndices;
     for (int i = 0; i < visiblePrimitiveReps.size() - 1; ++i) {
         // 必须以自左向右的顺序访问所有叶节点
         assert(blobTree.primitiveNodeIdx[i] < blobTree.primitiveNodeIdx[i + 1]);
@@ -574,26 +597,33 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
             for (auto primitiveIdx : visiblePrimitiveReps[i].subBlobTree.primitiveNodeIdx) {
                 realLeafIndices.push_back(primitiveIdx + originLeafIdx);
             }
-
-            realPrimitives.insert(realPrimitives.end(), visiblePrimitiveReps[i].tensors.begin(),
-                                  visiblePrimitiveReps[i].tensors.end());
+            minimalReps.reserve(minimalReps.size() + visiblePrimitiveReps[i].tensors.size());
+            const auto& aabb = visiblePrimitiveReps[i].aabb;
+            for (const auto& tensor : visiblePrimitiveReps[i].tensors) {
+                minimalReps.emplace_back(MinimalPrimitiveRep{tensor, aabb});
+            }
         } else {
             blobTree.structure[originLeafIdx].isPrimitive = true;
             realLeafIndices.push_back(originLeafIdx);
-            realPrimitives.emplace_back(visiblePrimitiveReps[i].tensors[0]);
+            minimalReps.emplace_back(MinimalPrimitiveRep{visiblePrimitiveReps[i].tensors[0], visiblePrimitiveReps[i].aabb});
         }
     }
     blobTree.primitiveNodeIdx = realLeafIndices;
     /*** merge subtrees to main tree ***/
-    Scene      scene{realPrimitives, AABB(xmin, xmax)};
+    Scene      scene{minimalReps, AABB(xmin, xmax)};
     OcTreeNode rootNode(0, 1, blobTree);
-    for (int i = 0; i < realPrimitives.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
+    for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
     int cnt = 1;
     buildOcTreeV0(scene, rootNode, leaves, 1, cnt);
 
+    basicTask(scene, leaves[14], q);
+
+    int i = 0;
     for (const auto& leaf : leaves) {
-        auto basicRes  = basicTask(scene, leaf, q);
-        volume        += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min);
+        auto basicRes = basicTask(scene, leaf, q);
+        if (std::isinf(basicRes.volume)) { std::cout << "inf volume when solving leaf: " << i << std::endl; }
+        volume += basicRes.volume;
+        std::cout << "Solved leaves: " << ++i << "/" << leaves.size() << std::endl;
     }
 
     volume *= prod(xmax - xmin);

algoim/organizer/primitive.hpp

@@ -635,9 +635,9 @@ public:
         uvector3 halfSize = size / 2;
         for (int i = 0; i < 8; ++i) {
             vertices[i]     = center;
-            vertices[i](0) += (i & 1) ? halfSize(0) : -halfSize(0);
+            vertices[i](0) += (i & 4) ? halfSize(0) : -halfSize(0);
             vertices[i](1) += (i & 2) ? halfSize(1) : -halfSize(1);
-            vertices[i](2) += (i & 4) ? halfSize(2) : -halfSize(2);
+            vertices[i](2) += (i & 1) ? halfSize(2) : -halfSize(2);
         }
         indices            = {3, 2, 0, 1,  // left
                               4, 6, 7, 5,  // right
@@ -678,6 +678,8 @@ public:
 
     uvector3 size() const { return max - min; }
 
+    real volume() const { return prod(size()); }
+
     bool intersect(const AABB& aabb) const
     {
         for (int i = 0; i < 3; ++i) {
@@ -708,7 +710,7 @@ struct Scene {
     // std::vector<CompleteTensorRep> polys;
     // std::vector<VisiblePrimitiveRep> visiblePrimitives;
     // organizer::BlobTree            blobTree;
-    std::vector<MinimalPrimitiveRep> tensors;
+    std::vector<MinimalPrimitiveRep> minimalReps;
     AABB                             boundary;
 };
 

algoim/quadrature_multipoly.hpp

@@ -18,6 +18,7 @@
 // a good choice, and this is specified by the following macro def.
 #include <cassert>
 #include <cmath>
+#include <iterator>
 #include <vector>
 #include "organizer/primitive.hpp"
 #define ALGOIM_M 8
@@ -614,6 +615,15 @@ struct ImplicitPolyQuadrature {
         build(true, false);
     }
 
+    ImplicitPolyQuadrature(const std::vector<organizer::MinimalPrimitiveRep>& minimalReps, const std::vector<int>& polyIndices)
+    {
+        for (auto i : polyIndices) {
+            auto mask = detail::nonzeroMask(minimalReps[i].tensor, booluarray<N, ALGOIM_M>(true));
+            if (!detail::maskEmpty(mask)) phi.push_back(minimalReps[i].tensor, mask);
+        }
+        build(true, false);
+    }
+
     // Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks
     ImplicitPolyQuadrature(const xarray<real, N>&         p,
                            const booluarray<N, ALGOIM_M>& pmask,
@@ -654,12 +664,13 @@ struct ImplicitPolyQuadrature {
             // Compute score; penalise any directions which likely contain vertical tangents
             uvector<bool, N> has_disc;
             uvector<real, N> score = detail::score_estimate(phi, has_disc);
-            // if (max(abs(score)) == 0)
-            //     score(0) = 1.0;
+            /**gjj */
+            if (max(abs(score)) == 0) score(0) = 1.0;
+            /**gjj */
             assert(max(abs(score)) > 0);
             score /= 2 * max(abs(score));
             for (int i = 0; i < N; ++i)
-                if (!has_disc(i)) score(i) += 1.0; // 1个与k相切的多项式都没有，才+1
+                if (!has_disc(i)) score(i) += 1.0; // 一个与k相切的多项式都没有时，才+1
 
             // Choose height direction and form base polynomials; if tanh-sinh is being used at this
             // level, suggest the same all the way down; moreover, suggest tanh-sinh if a non-empty
@@ -748,9 +759,20 @@ struct ImplicitPolyQuadrature {
                     if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j];
                 }
             };
-
+            /**gjj */
+            auto newEnd = std::remove_if(nodes, nodes + count, [](real x) { return std::isnan(x); });
+            count       = std::distance(nodes, newEnd);
+            // for (int i = 0; i < count; ++i) {
+            //     if (std::isnan(nodes[i])) { nanCount++; }
+            // }
+            // count -= nanCount;
+            /**gjj */
             // Sort the nodes
             std::sort(nodes, nodes + count);
+            if (!(nodes[0] == real(0) && nodes[count - 1] == real(1))) {
+                int aaa = 1;
+                int bbb = 1;
+            }
             assert(nodes[0] == real(0) && nodes[count - 1] == real(1));
 
             // Force nearly-degenerate sub-intervals to be exactly degenerate

gjj/primitiveDebug.hpp

@@ -151,8 +151,8 @@ void caseScene()
         uvector3{-1, 1,  0}
     };
     // primitiveDescriptions[4] = std::make_shared<PyramidDesc>(pyramidBottomVertices, uvector3{0, 0, 1});
-    primitiveDescriptions[4] = std::make_shared<SphereDesc>(SphereDesc(0.2, uvector3(0., 0., -1), 1.));
-    primitiveDescriptions[5] = std::make_shared<ConeDesc>(ConeDesc(uvector3(0., -0.2, 0.), -0.7, 0.4, 1));
+    primitiveDescriptions[4] = std::make_shared<SphereDesc>(SphereDesc(0.2, uvector3(0.2, -0.7, 0.), 1.));
+    primitiveDescriptions[5] = std::make_shared<ConeDesc>(ConeDesc(uvector3(0., -0.2, 0.), 0.4, -0.7, 1));
     // primitiveDescriptions[6] = std::make_shared<CylinderDesc>(CylinderDesc(uvector3(-0.3, 0.3, 2.3), 0.4, 3.6, 1));
     organizer::BlobTree blobTree;
     blobTree.structure.resize(PRIMITIVE_CNT * 2 - 1);
@@ -164,16 +164,16 @@ void caseScene()
     blobTree.structure[4]     = {1, 0, 0, 0, 0, 0};       // sphere2
     blobTree.structure[5]     = {0, 0, 0, 0, 0, 0};       // Union of spheres = opNode2
     blobTree.structure[6]     = {0, 2, 0, 0, 1, 10 - 6};  // Difference of opNode1 and opNode2 = opNode3
-    blobTree.structure[7]     = {1, 0, 0, 0, 1, 9 - 7};   // Pyramid
+    blobTree.structure[7]     = {1, 0, 0, 0, 1, 9 - 7};   // Pyramid (sphere3)
     blobTree.structure[8]     = {1, 0, 0, 0, 0, 0};       // Cone
-    blobTree.structure[9]     = {0, 1, 0, 0, 0, 0};       // Intersection of Pyramid and Cone = opNode4
-    blobTree.structure[10]    = {0, 0, 0, 0, 1, 12 - 10}; // Union of opNode3 and opNode4 = opNode5
+    blobTree.structure[9]     = {0, 0, 0, 0, 0, 0};       // UNION of Pyramid (sphere3) and Cone = opNode4
+    blobTree.structure[10]    = {0, 2, 0, 0, 1, 12 - 10}; // Difference of opNode3 and opNode4 = opNode5
                                                           // blobTree.structure[11]    = {1, 0, 0, 0, 0, 0};  // Cylinder
     // blobTree.structure[12]    = {0, 2, 0, 0, 1, 0};  // Difference of opNode5 and Cylinder
     // blobTree.primitiveNodeIdx = {0, 1, 3, 4, 7, 8, 11};
     // quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 2.3), blobTree);
     blobTree.primitiveNodeIdx = {0, 1, 3, 4, 7, 8};
-    quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree);
+    quadratureScene(primitiveDescriptions, uvector3(-0.8, -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree);
 }
 
 void testDeCasteljau()
@@ -279,6 +279,30 @@ void testPlaneUniformSign()
     std::cout << "sign = " << sign << std::endl;
 }
 
+void testPlaneWithinSubcell()
+{
+    AABB sceneAABB{uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 1.5)};
+    AABB subCell{
+        uvector3{0.076923076923076983, 0.17241379310344829, 0.77419354838709675},
+        uvector3{0.076923076923077094, 0.1724137931034484,  0.77419354838709686}
+    };
+    tensor3 planeTensor(nullptr, 3), subPlaneTensor(nullptr, 3);
+    algoim_spark_alloc(real, planeTensor, subPlaneTensor);
+    xarrayInit(planeTensor);
+    planeTensor.m(0)                 = -0.8;
+    planeTensor.m(uvector3{0, 0, 1}) = 1;
+    organizer::detail::powerTransformation(sceneAABB.size(), sceneAABB.min, planeTensor);
+    organizer::detail::power2BernsteinTensor(planeTensor);
+    bernstein::deCasteljau(planeTensor, subCell.min, subCell.max, subPlaneTensor);
+    int sign = bernstein::uniformSign(subPlaneTensor);
+    std::cout << sign << std::endl;
+    ImplicitPolyQuadrature<3> ipquad(subPlaneTensor);
+    real                      volume = 0;
+    ipquad.integrate(AutoMixed, 10, [&](const uvector<real, 3> &x, real w) {
+        if (evalBernstein(subPlaneTensor, x)) volume += w * 1;
+    });
+}
+
 void testBlob()
 {
     organizer::Blob blob = {5, 3, 4, 5, 6};
@@ -320,6 +344,7 @@ void testPrimitive()
     // testSubDivideWithDeCasteljau();
     // testBlob();
     caseScene();
+    // testPlaneWithinSubcell();
     // caseCubeSphere();
     // caseTwoCube();
     // testTensorInverse();