concave polyhedron

algoim/organizer/loader.hpp

@@ -1,8 +1,10 @@
 #pragma once
+#include <array>
 #include <sparkstack.hpp>
 #include <uvector.hpp>
 #include "organizer.hpp"
 #include "primitive.hpp"
+#include "timer.hpp"
 
 class Vector3D
 {
@@ -47,6 +49,21 @@ public:
         return node;
     }
 
+    void rotation(const std::array<std::array<double, 3>, 3>& rotationMatrix)
+    {
+        for (int i = 0; i < 3; i++) {
+            double temp = 0;
+            for (int j = 0; j < 3; j++) { temp += rotationMatrix[i][j] * this->getUVector3Data()(j); }
+            if (i == 0) {
+                this->m_x = temp;
+            } else if (i == 1) {
+                this->m_y = temp;
+            } else {
+                this->m_z = temp;
+            }
+        }
+    }
+
     double m_x, m_y, m_z;
 };
 
@@ -536,6 +553,7 @@ public:
      */
     std::pair<double, double> getAreaAndVolume(const BodyTag& tag)
     {
+        Timer timer("Building scene octree");
         auto& rep = this->m_allVisible[tag];
         assert(rep.tensors.size() == rep.aabbs.size());
 
@@ -565,23 +583,31 @@ public:
         std::vector<algoim::organizer::OcTreeNode> leaves;
         algoim::organizer::buildOcTreeV0(scene, rootNode, leaves, 1, cnt, 0, -1);
         std::cout << "octree built over" << std::endl;
+        timer.stop();
         // basicTask(scene, leaves[14], q);
 
-        int    i      = 0;
         int    q      = 10;
         double volume = 0;
         double area   = 0;
-        for (const auto& leaf : leaves) {
-            auto basicRes = algoim::organizer::basicTask(scene, leaf, q);
+        // #pragma omp parallel for reduction(+ : volume, area)
+        timer.rename("Quadrature");
+        std::vector<double> volumes(leaves.size(), 0);
+#pragma omp parallel for
+        for (int i = 0; i < leaves.size(); i++) {
+            const auto& leaf     = leaves[i];
+            auto        basicRes = algoim::organizer::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()
-                      << ", primitive cnt: " << leaf.polyIntersectIndices.size() << ", volume: " << volume << std::endl;
+            // volume += basicRes.volume;
+            volumes[i] = basicRes.volume;
+            std::cout << "Solved leaves: " << i + 1 << "/" << leaves.size()
+                      << ", primitive cnt: " << leaf.polyIntersectIndices.size() << ", leaf volume: " << basicRes.volume
+                      << std::endl;
         }
 
+        for (auto& v : volumes) { volume += v; }
         volume *= algoim::prod(rep.aabb.max - rep.aabb.min);
         std::cout << "Volume xxx: " << volume << std::endl;
-
+        timer.stop();
         return std::make_pair(area, volume);
     }
 
@@ -594,6 +620,8 @@ public:
         std::cout << std::endl;
     }
 
+    algoim::organizer::VisiblePrimitiveRep& getVisible(BodyTag tag) { return this->m_allVisible[tag]; }
+
 protected:
     /**
      * @brief Compute the barycentric coordinates of polygon

algoim/organizer/organizer.hpp

@@ -344,7 +344,7 @@ void buildOcTreeV0(const Scene&             scene,
         int bbb = 1;
     }
     const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
-    if (polyIntersectIndices.size() <= 7) {
+    if (polyIntersectIndices.size() <= 6) {
         leaves.emplace_back(node);
         return;
     }
@@ -454,7 +454,7 @@ void basicTask(const std::shared_ptr<PrimitiveDesc>& p, int q = 20, real xmin =
             if (isInsideBernstein(tensor, x)) volume += w * integrand(xmin + x * range);
         });
     } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
-        const int faceCount = pt->indexInclusiveScan.size();
+        const int faceCount = pt->faces.size();
         assert(faceCount > 1);
         std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 2));
         algoim_spark_alloc(real, planeTensors);
@@ -496,7 +496,7 @@ void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, in
             detail::power2BernsteinTensor(tensor);
             phis.emplace_back(tensor);
         } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
-            const int faceCount = pt->indexInclusiveScan.size();
+            const int faceCount = pt->faces.size();
             assert(faceCount > 1);
             std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 3));
             algoimSparkAllocHeapVector(phiStacks, planeTensors);
@@ -543,9 +543,9 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
     std::vector<tensor3>           phis; // phis using subcell as [0,1]^3
     std::vector<SparkStack<real>*> phiStacks;
 
-    uvector3i numSamples(50, 50, 50);
+    uvector3i numSamples(20, 20, 20);
     // if (false) {
-    if (node.polyIntersectIndices.size() > 8) {
+    if (node.polyIntersectIndices.size() >= 5) {
         // TODO: 蒙特卡洛
         std::cout << "Monte Carlo 1" << std::endl;
         ;
@@ -573,7 +573,6 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
     }
     // ImplicitPolyQuadrature<3> ipquad(scene.minimalReps, node.polyIntersectIndices);
     ImplicitPolyQuadrature<3> ipquad(phis);
-    int                       aaa = 1;
 #if STOP_WHEN_BLOCKED
     // if (stopCurrentQuadrature
     //     || std::chrono::duration_cast<std::chrono::seconds>(std::chrono::high_resolution_clock::now() -
@@ -585,9 +584,8 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
         stopCurrentQuadrature = false;
         // std::cout << "blocked" << std::endl;
         // 注意这里tensor已经经过deCasteljau变换，所以采样点属于[0,1]^3
-        uvector3i numSamples(20, 20, 20);
-        uvector3f interval = 1.f / (numSamples);
-        int       inCnt    = 0;
+        uvector3f interval    = 1.f / (numSamples);
+        int       inCnt       = 0;
         for (auto i = MultiLoop<3>(0, numSamples); ~i; ++i) {
             uvector3f x = i() * interval + interval * 0.5;
             if (keepQuadraturePoint(phis, node, x)) ++inCnt;
@@ -612,9 +610,8 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
         stopCurrentQuadrature = false;
         // std::cout << "blocked" << std::endl;
         // 注意这里tensor已经经过deCasteljau变换，所以采样点属于[0,1]^3
-        uvector3i numSamples(20, 20, 20);
-        uvector3f interval = 1.f / (numSamples);
-        int       inCnt    = 0;
+        uvector3f interval    = 1.f / (numSamples);
+        int       inCnt       = 0;
         for (auto i = MultiLoop<3>(0, numSamples); ~i; ++i) {
             uvector3f x = i() * interval + interval * 0.5;
             if (keepQuadraturePoint(phis, node, x)) ++inCnt;
@@ -626,6 +623,7 @@ BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10)
 #endif
 
     for (auto& p : phiStacks) delete p;
+    double aaa = volume * node.aabb.volume();
     return {volume * node.aabb.volume(), surf};
 }
 
@@ -656,7 +654,7 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
             detail::power2BernsteinTensor(phi);
             visiblePrimitiveReps[i].aabb.normalize(range, xmin);
         } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
-            const int faceCount = pt->indexInclusiveScan.size();
+            const int faceCount = pt->faces.size();
             assert(faceCount > 1);
             visiblePrimitiveReps[i].tensors = std::vector<tensor3>(faceCount, tensor3(nullptr, 3));
             auto& planeTensors              = visiblePrimitiveReps[i].tensors;
@@ -785,12 +783,14 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
     // basicTask(scene, leaves[14], q);
 
     int i = 0;
-    for (const auto& leaf : leaves) {
-        auto basicRes = basicTask(scene, leaf, q);
+
+    for (int i = 0; i < leaves.size(); ++i) {
+        const auto& leaf     = leaves[i];
+        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() << ", primitive cnt: " << leaf.polyIntersectIndices.size()
-                  << ", volume: " << volume << std::endl;
+                  << ", volume: " << basicRes.volume << std::endl;
     }
 
     volume *= prod(xmax - xmin);

algoim/organizer/primitive.hpp

@@ -684,13 +684,19 @@ public:
     const static PrimitiveType type = Mesh;
     std::vector<uvector3>      vertices;
     std::vector<int>           indices;
-    std::vector<int>           indexInclusiveScan;
+    // std::vector<int>           indexInclusiveScan;
+    std::vector<uvector2i>     faces;
 
     MeshDesc(const std::vector<uvector3>& vertices_,
              const std::vector<int>&      indices_,
              const std::vector<int>&      indexInclusiveScan_)
-        : PrimitiveDesc(), vertices(vertices_), indices(indices_), indexInclusiveScan(indexInclusiveScan_)
+        : PrimitiveDesc(), vertices(vertices_), indices(indices_)
     {
+        faces.resize(indexInclusiveScan_.size());
+        for (int i = 0; i < indexInclusiveScan_.size(); ++i) {
+            faces[i](0) = i == 0 ? 0 : indexInclusiveScan_[i];
+            faces[i](1) = indexInclusiveScan_[i] - faces[i](0);
+        }
     }
 
     /**
@@ -699,11 +705,11 @@ public:
              @param faces_ 面索引 {顶点索引其实偏移（exclusive scan），顶点数}
      */
     MeshDesc(const std::vector<uvector3>& vertices_, const std::vector<int>& indices_, const std::vector<uvector2i>& faces_)
-        : PrimitiveDesc(), vertices(vertices_), indices(indices_)
+        : PrimitiveDesc(), vertices(vertices_), indices(indices_), faces(faces_)
     {
-        indexInclusiveScan.resize(faces_.size());
-        for (int i = 0; i < faces_.size(); ++i) { indexInclusiveScan[i] = faces_[i](0) + faces_[i](1); }
-        int aaa = 1;
+        // indexInclusiveScan.resize(faces_.size());
+        // for (int i = 0; i < faces_.size(); ++i) { indexInclusiveScan[i] = faces_[i](0) + faces_[i](1); }
+        // int aaa = 1;
     }
 
     MeshDesc() = default;
@@ -720,18 +726,24 @@ public:
         int bottomSize = bottomVertices.size();
         vertices.emplace_back(topVertex);
         indices.reserve(bottomSize * 3 + bottomSize);
-        indexInclusiveScan.reserve(bottomSize + 1);
+        faces.reserve(bottomSize + 1);
+        // indexInclusiveScan.reserve(bottomSize + 1);
         for (int i = 0; i < bottomSize; ++i) {
+            // 所有三角侧面
             indices.emplace_back(i);
             indices.emplace_back(i == bottomSize - 1 ? 0 : indices.emplace_back(i + 1));
             indices.emplace_back(bottomSize);
-            indexInclusiveScan.emplace_back(indexInclusiveScan.empty() ? 3 : indexInclusiveScan.back() + 3);
+            if (indices.empty()) {
+                faces.emplace_back(0, 3);
+            } else {
+                faces.emplace_back(faces.back()(0) + faces.back()(1), 3);
+            }
         }
         for (int i = 0; i < bottomSize; ++i) {
             indices.emplace_back(i);
             indices.emplace_back(i);
         }
-        indexInclusiveScan.emplace_back(indexInclusiveScan.back() + bottomSize);
+        faces.emplace_back(faces.back()(0) + faces.back()(1), bottomSize);
     }
 };
 
@@ -748,13 +760,20 @@ public:
             vertices[i](1) += (i & 2) ? halfSize(1) : -halfSize(1);
             vertices[i](2) += (i & 1) ? halfSize(2) : -halfSize(2);
         }
-        indices            = {3, 2, 0, 1,  // left
-                              4, 6, 7, 5,  // right
-                              6, 2, 3, 7,  // top
-                              1, 0, 4, 5,  // bottom
-                              7, 3, 1, 5,  // front
-                              2, 6, 4, 0}; // back
-        indexInclusiveScan = {4, 8, 12, 16, 20, 24};
+        indices = {3, 2, 0, 1,  // left
+                   4, 6, 7, 5,  // right
+                   6, 2, 3, 7,  // top
+                   1, 0, 4, 5,  // bottom
+                   7, 3, 1, 5,  // front
+                   2, 6, 4, 0}; // back
+        faces   = {
+            uvector2i{0,  4},
+            uvector2i{4,  4},
+            uvector2i{8,  4},
+            uvector2i{12, 4},
+            uvector2i{16, 4},
+            uvector2i{20, 4}
+        };
     }
 };
 
@@ -825,10 +844,10 @@ uvector3 getFaceNorm(const std::vector<uvector3>& points, const std::vector<int>
 
 void makeMesh(const MeshDesc& mesh, VisiblePrimitiveRep& visiblePrimitive)
 {
-    assert(visiblePrimitive.tensors.size() == mesh.indexInclusiveScan.size());
-    for (int i = 0; i < mesh.indexInclusiveScan.size(); ++i) {
-        const int indexBeg  = i == 0 ? 0 : mesh.indexInclusiveScan[i - 1];
-        const int indexSize = mesh.indexInclusiveScan[i] - indexBeg;
+    assert(visiblePrimitive.tensors.size() == mesh.faces.size());
+    for (int i = 0; i < mesh.faces.size(); ++i) {
+        const int indexBeg  = mesh.faces[i](0);
+        const int indexSize = mesh.faces[i](1);
         assert(indexSize >= 3);
         auto& planeTensor = visiblePrimitive.tensors[i];
         xarrayInit(planeTensor);
@@ -849,7 +868,7 @@ void makeMesh(const MeshDesc& mesh, VisiblePrimitiveRep& visiblePrimitive)
         // planeTensor.m(uvector3(0, 1, 0))  = -N(1);
         // planeTensor.m(uvector3(0, 0, 1))  = -N(2);
         // test other vertices
-        for (int j = indexBeg + 3; j < mesh.indexInclusiveScan[i]; ++j) {
+        for (int j = indexBeg + 3; j < indexBeg + indexSize; ++j) {
             auto tmp = dot(N, vertices[indices[j]]) + d;
             if (std::abs(tmp) > std::numeric_limits<real>::epsilon() * 1e6) {
                 std::cerr << "Points are not on the same plane! " << tmp << std::endl;
@@ -861,7 +880,7 @@ void makeMesh(const MeshDesc& mesh, VisiblePrimitiveRep& visiblePrimitive)
     for (const auto& v : mesh.vertices) { visiblePrimitive.aabb.extend(v); }
     // subBlobTree
     visiblePrimitive.subBlobTree.clear();
-    buildNearBalancedBlobTree(visiblePrimitive.subBlobTree, mesh.indexInclusiveScan.size());
+    buildNearBalancedBlobTree(visiblePrimitive.subBlobTree, mesh.faces.size());
 };
 
 void makeSphere(const SphereDesc& sphereDesc, VisiblePrimitiveRep& visiblePrimitive)

algoim/organizer/timer.hpp

@@ -0,0 +1,30 @@
+#pragma once
+#include <chrono>
+#include <string>
+#include <iostream>
+
+class Timer
+{
+private:
+    std::chrono::steady_clock::time_point m_begin;
+    std::chrono::steady_clock::time_point m_end;
+    std::string                           m_name;
+
+public:
+    Timer(std::string name) : m_name(std::move(name)) { m_begin = std::chrono::steady_clock::now(); }
+
+    void restart() { m_begin = std::chrono::steady_clock::now(); }
+
+    void rename(std::string name)
+    {
+        m_name  = std::move(name);
+        m_begin = std::chrono::steady_clock::now();
+    }
+
+    void stop()
+    {
+        m_end           = std::chrono::steady_clock::now();
+        double duration = std::chrono::duration<double>(m_end - m_begin).count();
+        std::cout << m_name << " duration: " << duration << "s" << std::endl;
+    }
+};

algoim/quadrature_multipoly.hpp

@@ -36,13 +36,13 @@
 #include "tanhsinh.hpp"
 #include "bernstein.hpp"
 
-#define STOP_WHEN_BLOCKED false
+#define STOP_WHEN_BLOCKED true
 
 #if STOP_WHEN_BLOCKED
 #include <chrono>
-auto         timerStart            = std::chrono::high_resolution_clock::now();
-bool         stopCurrentQuadrature = false;
-const double MAX_DURATION          = 10; // seconds
+static thread_local auto timerStart            = std::chrono::high_resolution_clock::now();
+static thread_local bool stopCurrentQuadrature = false;
+const double             MAX_DURATION          = 3; // seconds
 #endif
 
 // #define

examples/example_loader.cpp

@@ -321,12 +321,12 @@ void loaderTest1()
 
     after = loader.getAreaAndVolume(tag1);
 
-    /* ��ȡ ��1 ����ǰ��� ����� �� ����� */
+    /* 锟斤拷取 锟斤拷1 锟斤拷锟斤拷前锟斤拷锟? 锟斤拷锟斤拷锟? 锟斤拷 锟斤拷锟斤拷锟? */
     areaDifference   = after.first - before.first;
     volumeDifference = after.second - before.second;
     std::cout << areaDifference << ", " << volumeDifference << std::endl;
 
-    /* ������� �ۼ���1 */
+    /* 锟斤拷锟斤拷锟斤拷锟? 锟桔硷拷锟斤拷1 */
     auto subTag1 = loader.addEmpty();
     auto cycle   = [&loader, subTag1](const Point3D& point) {
         std::vector<Point3D> points;
@@ -441,7 +441,7 @@ void loaderTest1()
     loader.differentNode(tag1, subTag1);
     after = loader.getAreaAndVolume(tag1);
 
-    /* ��ȡ ��1 ����ǰ��� ����� �� ����� */
+    /* 锟斤拷取 锟斤拷1 锟斤拷锟斤拷前锟斤拷锟? 锟斤拷锟斤拷锟? 锟斤拷 锟斤拷锟斤拷锟? */
     areaDifference   = after.first - before.first;
     volumeDifference = after.second - before.second;
     std::cout << areaDifference << ", " << volumeDifference << std::endl;
@@ -475,7 +475,7 @@ void loaderTest1()
     loader.differentNode(tag1, tag16);
     after = loader.getAreaAndVolume(tag1);
 
-    /* ��ȡ ��1 ����ǰ��� ����� �� ����� */
+    /* 锟斤拷取 锟斤拷1 锟斤拷锟斤拷前锟斤拷锟? 锟斤拷锟斤拷锟? 锟斤拷 锟斤拷锟斤拷锟? */
     areaDifference   = after.first - before.first;
     volumeDifference = after.second - before.second;
     std::cout << areaDifference << ", " << volumeDifference << std::endl;
@@ -621,7 +621,7 @@ void loaderTest2Simplified()
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
 
-    /* ������� �ۼ���1 */
+    /* 锟斤拷锟斤拷锟斤拷锟? 锟桔硷拷锟斤拷1 */
     // auto subTag1 = loader.addEmpty();
     // auto cycle   = [&loader, subTag1](const Point3D& point) {
     //     std::vector<Point3D> points;
@@ -687,7 +687,7 @@ void loaderTest2Simplified()
     // loader.differentNode(tag1, subTag1);
     // after = loader.getAreaAndVolume(tag1);
 
-    /* ��ȡ ��1 ����ǰ��� ����� �� ����� */
+    /* 锟斤拷取 锟斤拷1 锟斤拷锟斤拷前锟斤拷锟? 锟斤拷锟斤拷锟? 锟斤拷 锟斤拷锟斤拷锟? */
     // areaDifference = after.first - before.first;
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
@@ -721,6 +721,7 @@ void loaderTest2Simplified()
 
 void loaderTest2()
 {
+    Timer  timer("loading and blob construction");
     Loader loader;
 
     std::vector<Point3D> points;
@@ -1031,7 +1032,7 @@ void loaderTest2()
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
 
-    /* ������� �ۼ���1 */
+    /* 锟斤拷锟斤拷锟斤拷锟? 锟桔硷拷锟斤拷1 */
     auto subTag1 = loader.addEmpty();
     auto cycle   = [&loader, subTag1](const Point3D& point) {
         std::vector<Point3D> points;
@@ -1146,7 +1147,7 @@ void loaderTest2()
     loader.differentNode(tag1, subTag1);
     // after = loader.getAreaAndVolume(tag1);
 
-    /* ��ȡ ��1 ����ǰ��� ����� �� ����� */
+    /* 锟斤拷取 锟斤拷1 锟斤拷锟斤拷前锟斤拷锟? 锟斤拷锟斤拷锟? 锟斤拷 锟斤拷锟斤拷锟? */
     // areaDifference = after.first - before.first;
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
@@ -1346,7 +1347,7 @@ void loaderTest2()
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
 
-    /* ����Ϊһ��ѭ����͹����̶��������鲻һ������������һ��
+    /* 锟斤拷锟斤拷为一锟斤拷循锟斤拷锟斤拷凸锟斤拷锟斤拷潭锟斤拷锟斤拷锟斤拷锟斤拷椴灰伙拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟揭伙拷锟?
      */
     auto cycle1 = [&loader, tag1](const Point3D& point1, const Point3D point2) {
         /* ѭ���� */
@@ -1530,7 +1531,7 @@ void loaderTest2()
     points.push_back(Point3D{-22821.612136159150, -11781.067532595984, 0.0000000000000000});
     for (size_t i = 0; i < points.size(); i += 2) { cycle1(points[i], points[i + 1]); }
 
-    /* ����Ϊһ��ѭ����͹����̶��������鲻һ������������һ��
+    /* 锟斤拷锟斤拷为一锟斤拷循锟斤拷锟斤拷凸锟斤拷锟斤拷潭锟斤拷锟斤拷锟斤拷锟斤拷椴灰伙拷锟斤拷锟斤拷锟斤拷锟斤拷锟斤拷锟揭伙拷锟?
      */
     auto cycle2 = [&loader, tag1](const Point3D& point1, const Point3D point2, const Point3D point3, const Point3D point4) {
         /* ѭ���� */
@@ -1576,6 +1577,7 @@ void loaderTest2()
     // volumeDifference = after.second - before.second;
     // std::cout << areaDifference << ", " << volumeDifference << std::endl;
 
+    timer.stop();
     auto result = loader.getAreaAndVolume(tag1);
     // std::cout << "result:" << result.first << ", " << result.second << std::endl;
 
@@ -1612,7 +1614,7 @@ int main(int argc, char** argv)
     // loaderEntBool();
     // loaderTest2Simplified();
     // loaderTest2();
-    // loaderPolyTilted();
-    loaderPolyThickPlane();
+    loaderPolyTilted();
+    // loaderPolyThickPlane();
     return 0;
 }

examples/example_loader2.hpp

@@ -1,11 +1,14 @@
 #pragma once
+#include <chrono>
+#include "organizer/blobtree.hpp"
 #include "organizer/loader.hpp"
+#include "organizer/timer.hpp"
 
 void loaderPolyTilted()
 {
-    Loader loader;
-
-    std::vector<Point3D> points = {
+    Loader               loader;
+    Timer                timer("loading and blob construction");
+    std::vector<Point3D> pointsA = {
         {11748.40944836227, -4325.95442631630,  21549.59021205360},
         {11748.40944836227, -4325.95442631630,  21688.15427665911},
         {12484.47483358811, -3577.50000026381,  21256.03391230463},
@@ -295,7 +298,7 @@ void loaderPolyTilted()
         {17624.99999998905, -8632.93491411901,  18198.71313931236},
         {17624.99999998906, -9996.24203115674,  18211.84978571097}
     };
-    std::vector<int> indices = {
+    std::vector<int> indicesA = {
         0,   1,   2,   3,   4,   1,   0,   5,   6,   7,   8,   2,   1,   4,   9,   3,   2,   8,   10,  11,  10,  8,   9,
         12,  6,   5,   11,  12,  11,  5,   0,   3,   10,  7,   6,   12,  9,   4,   7,   12,  85,  86,  87,  88,  89,  90,
         91,  92,  86,  85,  91,  93,  94,  95,  96,  97,  87,  86,  92,  98,  99,  100, 101, 102, 103, 96,  95,  104, 105,
@@ -306,7 +309,7 @@ void loaderPolyTilted()
         197, 188, 187, 186, 148, 197, 187, 271, 272, 273, 274, 241, 266, 275, 276, 277, 272, 271, 278, 279, 280, 281, 282,
         283, 266, 241, 284, 280, 279, 285, 286, 273, 285, 279, 284, 274, 273, 272, 277, 287, 286, 285, 287, 277, 276, 282,
         281, 280, 286, 276, 275, 283, 282, 266, 283, 275, 284, 241, 278, 274, 278, 271, 274};
-    std::vector<std::pair<int, int>> faces = {
+    std::vector<std::pair<int, int>> facesA = {
         {0,   4 },
         {4,   6 },
         {10,  5 },
@@ -353,33 +356,58 @@ void loaderPolyTilted()
         {217, 4 },
         {221, 3 }
     };
+
+    std::vector<int>                              objExclusiveScan = {0, 9, 24, 34};
+    std::vector<std::vector<std::pair<int, int>>> objFaces(objExclusiveScan.size());
+    for (int i = 0; i < objExclusiveScan.size(); i++) {
+        objFaces[i].assign(facesA.begin() + objExclusiveScan[i],
+                           i == objExclusiveScan.size() - 1 ? facesA.end() : facesA.begin() + objExclusiveScan[i + 1]);
+    }
+
     // debug strange face
     std::vector<Point3D> strangeFace;
-    for (int index = faces[13].first; index < faces[13].first + faces[13].second; index++) {
-        strangeFace.push_back(points[indices[index]]);
+    for (int index = facesA[13].first; index < facesA[13].first + facesA[13].second; index++) {
+        strangeFace.push_back(pointsA[indicesA[index]]);
     }
 
-    auto tag1 = loader.addMesh(points, indices, faces);
     // std::cout << "Des volume:" << std::endl;
     // loader.getAreaAndVolume(tag1);
 
-    auto tag2 = loader.addMesh(
-        {
-            {17624.99999998905, -8510.00401819123, 17830.90112057791},
-            {17374.99999998905, -8510.00401819123, 17830.90112057791},
-            {17374.99999998905, -8510.00401819123, 18408.25138976753},
-            {17624.99999998905, -8510.00401819123, 18408.25138976753},
-            {17624.99999998906, -3577.50000024273, 20678.68364311939},
-            {17374.99999998906, -3577.50000024273, 20678.68364311939},
-            {17374.99999998906, -3577.50000024273, 21256.03391230901},
-            {17624.99999998906, -3577.50000024273, 21256.03391230901}
-    },
-        {0, 1, 2, 3, 0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7, 4, 5, 6, 7},
-        {{0, 4}, {4, 4}, {8, 4}, {12, 4}, {16, 4}, {20, 4}});
-    // std::cout << "Src volume:" << std::endl;
+    std::vector<Point3D> pointsB = {
+        {17624.99999998905, -8510.00401819123, 17830.90112057791},
+        {17374.99999998905, -8510.00401819123, 17830.90112057791},
+        {17374.99999998905, -8510.00401819123, 18408.25138976753},
+        {17624.99999998905, -8510.00401819123, 18408.25138976753},
+        {17624.99999998906, -3577.50000024273, 20678.68364311939},
+        {17374.99999998906, -3577.50000024273, 20678.68364311939},
+        {17374.99999998906, -3577.50000024273, 21256.03391230901},
+        {17624.99999998906, -3577.50000024273, 21256.03391230901}
+    };
+    std::vector<int>                 indicesB = {0, 1, 2, 3, 0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6, 3, 0, 4, 7, 4, 5, 6, 7};
+    std::vector<std::pair<int, int>> facesB   = {
+        {0,  4},
+        {4,  4},
+        {8,  4},
+        {12, 4},
+        {16, 4},
+        {20, 4}
+    };
+
+    const std::array<std::array<double, 3>, 3> rotationMatrix = {-0.6402, 0.6868, 0.3440,  0.6713, 0.7180,
+                                                                 -0.1840, 0.3734, -0.1131, 0.9208};
+    // for (auto& point : pointsA) { point.rotation(rotationMatrix); }
+
+    auto tag1 = loader.addMesh(pointsA, indicesA, objFaces[0]);
+    for (int i = 1; i < objFaces.size(); i++) {
+        auto tmpTag = loader.addMesh(pointsA, indicesA, objFaces[i]);
+        loader.unionNode(tag1, tmpTag);
+    }
+    auto tag2 = loader.addMesh(pointsB, indicesB, facesB);
+
     // loader.getAreaAndVolume(tag2);
 
     loader.differentNode(tag1, tag2);
+    timer.stop();
     loader.getAreaAndVolume(tag1);
 }
 
@@ -387,7 +415,8 @@ void loaderPolyThickPlane()
 {
     Loader loader;
 
-    auto tag1 = loader.addMesh(
+    Timer timer("loading and blob construction");
+    auto  tag1 = loader.addMesh(
         {
             {2629.82969004962, 6801.21455575771, 49030.00000000000},
             {2869.82968877607, 6801.18983125328, 49030.00000000000},
@@ -404,28 +433,34 @@ void loaderPolyThickPlane()
     },
         {0, 1, 2, 3, 4, 1, 0, 5, 6, 7, 6, 5, 8, 9, 5, 0, 3, 8, 8, 3, 2, 10, 11, 9, 1, 4, 10, 2, 4, 7, 11, 10, 7, 6, 9, 11},
         {{0, 4}, {4, 6}, {10, 4}, {14, 4}, {18, 6}, {24, 4}, {28, 4}, {32, 4}});
-    std::cout << "Des volume:" << std::endl;
-    loader.getAreaAndVolume(tag1);
 
-    // auto tag2 = loader.addMesh(
-    //     {
-    //         {2629.99453728099, 8340.01156816906, 52180.00000000000},
-    //         {2869.99453728099, 8340.01156816906, 52180.00000000000},
-    //         {2869.99453728099, 8401.36916293070, 52180.00000000000},
-    //         {2869.99453728099, 8580.01156816906, 52180.00000000000},
-    //         {2629.99453728099, 8580.01156816906, 52180.00000000000},
-    //         {2629.99453728099, 8340.02684557581, 52180.00000000000},
-    //         {2629.99453728099, 8340.01156816906, 49030.00000000000},
-    //         {2869.99453728099, 8340.01156816906, 49030.00000000000},
-    //         {2869.99453728099, 8401.36916293070, 49030.00000000000},
-    //         {2869.99453728099, 8580.01156816906, 49030.00000000000},
-    //         {2629.99453728099, 8580.01156816906, 49030.00000000000}
-    // },
-    //     {0, 1, 2, 3, 4, 5, 6, 7, 1, 0, 7, 8, 9, 3, 2, 1, 9, 10, 4, 3, 10, 6, 0, 5, 4, 6, 10, 9, 8, 7},
-    //     {{0, 6}, {6, 4}, {10, 6}, {16, 4}, {20, 5}, {25, 5}});
-    // std::cout << "Src volume:" << std::endl;
-    // loader.getAreaAndVolume(tag2);
 
-    // loader.differentNode(tag1, tag2);
+    /**** 处理凹多面体 */
+    auto& rep                           = loader.getVisible(tag1);
+    rep.subBlobTree.structure[5].nodeOp = algoim::organizer::OP_UNION;
+    // std::cout << "Des volume:" << std::endl;
     // loader.getAreaAndVolume(tag1);
+
+    auto tag2 = loader.addMesh(
+        {
+            {2629.99453728099, 8340.01156816906, 52180.00000000000},
+            {2869.99453728099, 8340.01156816906, 52180.00000000000},
+            {2869.99453728099, 8401.36916293070, 52180.00000000000},
+            {2869.99453728099, 8580.01156816906, 52180.00000000000},
+            {2629.99453728099, 8580.01156816906, 52180.00000000000},
+            {2629.99453728099, 8340.02684557581, 52180.00000000000},
+            {2629.99453728099, 8340.01156816906, 49030.00000000000},
+            {2869.99453728099, 8340.01156816906, 49030.00000000000},
+            {2869.99453728099, 8401.36916293070, 49030.00000000000},
+            {2869.99453728099, 8580.01156816906, 49030.00000000000},
+            {2629.99453728099, 8580.01156816906, 49030.00000000000}
+    },
+        {0, 1, 2, 3, 4, 5, 6, 7, 1, 0, 7, 8, 9, 3, 2, 1, 9, 10, 4, 3, 10, 6, 0, 5, 4, 6, 10, 9, 8, 7},
+        {{0, 6}, {6, 4}, {10, 6}, {16, 4}, {20, 5}, {25, 5}});
+
+    loader.differentNode(tag1, tag2);
+
+    timer.stop();
+
+    loader.getAreaAndVolume(tag1);
 }