fix debug of wrong SparkStack allocation order

algoim/organizer/organizer.hpp

@@ -37,48 +37,113 @@ bool keepQuadraturePoint(std::vector<tensor3>& originTensors, const uvector3& or
     return true;
 }
 
-class BasicTask
+// std::vector<std::shared_ptr<PrimitiveDesc>> primitives;
+
+// BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
+
+
+void basicTask(const std::shared_ptr<PrimitiveDesc>& p, int q = 20, real xmin = -1, real xmax = 1)
 {
-public:
-    // std::vector<std::shared_ptr<PrimitiveDesc>> primitives;
+    real     volume    = 0;
+    auto     integrand = [](const uvector<real, 3>& x) { return 1.0; };
+    uvector3 range     = xmax - xmin;
+
+
+    if (auto pt = std::dynamic_pointer_cast<SphereDesc>(p)) {
+        tensor3 tensor(nullptr, 3);
+        algoim_spark_alloc(real, tensor);
+        makeSphere(*pt, tensor);
+        detail::powerTransformation(range, xmin, tensor);
+
+        tensor3 phi(nullptr, 3);
+        algoim_spark_alloc(real, phi);
+        detail::power2BernsteinTensor(tensor, phi);
+
+        uvector<real, 3> testX(0., 0., 0.25);
+        real             testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
+        // auto             vec1              = xarray2StdVector(phi);
+        std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;
+
+        ImplicitPolyQuadrature<3> ipquad(phi);
+        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
+            if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin));
+        });
+    } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
+        const int faceCount = pt->indexInclusiveScan.size();
+        assert(faceCount > 1);
+        std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 2));
+        algoim_spark_alloc(real, planeTensors);
+
+        tensor3 compositeTensor(nullptr, 1 + faceCount);
+        algoim_spark_alloc(real, compositeTensor);
 
-    // BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
+        makeMesh(*pt, compositeTensor, planeTensors);
+        detail::powerTransformation(range, xmin, compositeTensor);
 
+        auto planeStdVector1          = xarray2StdVector(planeTensors[0]);
+        auto planeStdVector2          = xarray2StdVector(planeTensors[1]);
+        auto compositeTensorStdVector = xarray2StdVector(compositeTensor);
 
-    BasicTask(std::shared_ptr<PrimitiveDesc> p)
-    {
-        int      q         = 20;
-        real     volume    = 0;
-        real     xmin      = -1;
-        real     xmax      = 1;
-        auto     integrand = [](const uvector<real, 3>& x) { return 1.0; };
-        uvector3 range     = xmax - xmin;
+        uvector<real, 3> testX(0., 0.75, 0.2);
+        real             textEvalPower = evalPower(compositeTensor, testX);
 
+        tensor3 phi(nullptr, 1 + faceCount);
+        algoim_spark_alloc(real, phi);
+        detail::power2BernsteinTensor(compositeTensor, phi);
 
-        if (auto pt = std::dynamic_pointer_cast<SphereDesc>(p)) {
-            tensor3 tensor(nullptr, 3);
-            algoim_spark_alloc(real, tensor);
-            makeSphere(*pt, tensor);
-            detail::powerTransformation(range, xmin, tensor);
+        int                       quadraturePointCount = 0;
+        ImplicitPolyQuadrature<3> ipquad(phi);
+
+        real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
+        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
+            quadraturePointCount++;
+            auto realX = x * range + xmin;
+            if (isInsidePowers(planeTensors, realX)) volume += w * integrand(realX);
+        });
+        std::cout << "textEvalPower: " << textEvalPower << std::endl;
+        std::cout << "quadraturePointCount: " << quadraturePointCount << std::endl;
+    }
+    volume *= pow(xmax - xmin, 3);
+    std::cout << "Volume xxx: " << volume << std::endl;
+};
+
+void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, int q = 20, real xmin = -1, real xmax = 1)
+{
+    std::vector<SparkStack<real>*> phiStacks;
+    std::vector<tensor3>           phis;
+
+    std::vector<SparkStack<real>*> originTensorStacks;
+    std::vector<tensor3>           originTensors;
+    real                           volume;
+    auto                           integrand = [](const uvector<real, 3>& x) { return 1.0; };
+    uvector3                       range     = xmax - xmin;
+    uvector<real, 3>               testX(0., 0.75, 0.2);
+    for (int i = 0; i < primitives.size(); i++) {
+        if (auto pt = std::dynamic_pointer_cast<SphereDesc>(primitives[i])) {
+            tensor3 originTensor(nullptr, 3), transformedTensor(nullptr, 3);
+            originTensorStacks.emplace_back(algoim_spark_alloc_heap(real, originTensor)); // 记录，用以最后bool
 
             tensor3 phi(nullptr, 3);
-            algoim_spark_alloc(real, phi);
-            detail::power2BernsteinTensor(tensor, phi);
-
-            uvector<real, 3> testX(0., 0., 0.25);
-            real             testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
-            // auto             vec1              = xarray2StdVector(phi);
-            std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;
-
-            ImplicitPolyQuadrature<3> ipquad(phi);
-            ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
-                if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin));
-            });
-        } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
+            phiStacks.emplace_back(
+                algoim_spark_alloc_heap(real, phi)); // 必须先于algoim_spark_alloc，使transformedTensor的内存以栈形式释放
+
+            algoim_spark_alloc(real, transformedTensor);
+            makeSphere(*pt, originTensor);
+            originTensors.emplace_back(originTensor);
+            detail::powerTransformation(range, xmin, originTensor, transformedTensor);
+
+
+            detail::power2BernsteinTensor(transformedTensor, phi);
+            phis.emplace_back(phi);
+        } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
             const int faceCount = pt->indexInclusiveScan.size();
             assert(faceCount > 1);
             std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 2));
-            algoim_spark_alloc(real, planeTensors);
+            algoimSparkAllocHeapVector(originTensorStacks, planeTensors);
+
+            tensor3 phi(nullptr, 1 + faceCount);
+            phiStacks.emplace_back(
+                algoim_spark_alloc_heap(real, phi)); // 必须先于algoim_spark_alloc，使compositeTensor的内存以栈形式释放
 
             tensor3 compositeTensor(nullptr, 1 + faceCount);
             algoim_spark_alloc(real, compositeTensor);
@@ -86,89 +151,31 @@ public:
             makeMesh(*pt, compositeTensor, planeTensors);
             detail::powerTransformation(range, xmin, compositeTensor);
 
-            auto planeStdVector1          = xarray2StdVector(planeTensors[0]);
-            auto planeStdVector2          = xarray2StdVector(planeTensors[1]);
-            auto compositeTensorStdVector = xarray2StdVector(compositeTensor);
+            real testEvalPower = evalPower(compositeTensor, testX);
 
-            uvector<real, 3> testX(0., 0.75, 0.2);
-            real             textEvalPower = evalPower(compositeTensor, testX);
+            originTensors.insert(originTensors.end(), planeTensors.begin(), planeTensors.end());
 
-            tensor3 phi(nullptr, 1 + faceCount);
-            algoim_spark_alloc(real, phi);
             detail::power2BernsteinTensor(compositeTensor, phi);
-
-            int                       quadraturePointCount = 0;
-            ImplicitPolyQuadrature<3> ipquad(phi);
-
             real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
-            ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
-                quadraturePointCount++;
-                auto realX = x * range + xmin;
-                if (isInsidePowers(planeTensors, realX)) volume += w * integrand(realX);
-            });
-            std::cout << "textEvalPower: " << textEvalPower << std::endl;
-            std::cout << "quadraturePointCount: " << quadraturePointCount << std::endl;
+            phis.emplace_back(phi);
         }
-        volume *= pow(xmax - xmin, 3);
-        std::cout << "Volume xxx: " << volume << std::endl;
-    };
-
-    BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> primitives)
-    {
-        std::vector<SparkStack<real>*> phiStacks;
-        std::vector<tensor3>           phis;
-
-        std::vector<SparkStack<real>*> originTensorStacks;
-        std::vector<tensor3>           originTensors;
-        int                            q = 10;
-        real                           volume;
-        uvector3                       xmin      = 0;
-        uvector3                       xmax      = 1;
-        auto                           integrand = [](const uvector<real, 3>& x) { return 1.0; };
-        uvector3                       range     = xmax - xmin;
-
-        for (int i = 0; i < primitives.size(); i++) {
-            if (auto pt = std::dynamic_pointer_cast<SphereDesc>(primitives[i])) {
-                tensor3 originTensor(nullptr, 3), transformedTensor(nullptr, 3);
-                originTensorStacks.emplace_back(algoim_spark_alloc_heap(real, originTensor)); // 记录，用以最后bool
-                algoim_spark_alloc(real, transformedTensor);
-                makeSphere(*pt, originTensor);
-                originTensors.emplace_back(originTensor);
-                detail::powerTransformation(range, xmin, originTensor, transformedTensor);
-
-                tensor3 phi(nullptr, 3);
-                phiStacks.emplace_back(algoim_spark_alloc_heap(real, phi));
-                detail::power2BernsteinTensor(transformedTensor, phi);
-                phis.emplace_back(phi);
-            } else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) {
-                const int faceCount = pt->indexInclusiveScan.size();
-                assert(faceCount > 1);
-                std::vector<tensor3> planeTensors(faceCount, tensor3(nullptr, 2));
-                algoimSparkAllocHeapVector(originTensorStacks, planeTensors);
-
-                tensor3 compositeTensor(nullptr, 1 + faceCount);
-                algoim_spark_alloc(real, compositeTensor);
-
-                makeMesh(*pt, compositeTensor, planeTensors);
-                detail::powerTransformation(range, xmin, compositeTensor);
-                originTensors.insert(originTensors.end(), planeTensors.begin(), planeTensors.end());
-
-                tensor3 phi(nullptr, 1 + faceCount);
-                phiStacks.emplace_back(algoim_spark_alloc_heap(real, phi));
-                detail::power2BernsteinTensor(compositeTensor, phi);
-                phis.emplace_back(phi);
-            }
-        }
-        ImplicitPolyQuadrature<3> ipquad(phis);
-        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
-            auto realX = x * range + xmin;
-            if (keepQuadraturePoint(originTensors, realX)) volume += w * integrand(realX);
-        });
-        std::cout << "Volume xxx: " << volume << std::endl;
-
-        // free memory
-        for (auto& p : phiStacks) delete p;
-        for (auto& p : originTensorStacks) delete p;
-    };
+    }
+    real                      testEvalBernstein = bernstein::evalBernsteinPoly(phis[0], testX);
+    ImplicitPolyQuadrature<3> ipquad(phis);
+
+    ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
+        auto realX = x * range + xmin;
+        if (keepQuadraturePoint(originTensors, realX)) volume += w * integrand(realX);
+    });
+    volume *= pow(xmax - xmin, 3);
+    std::cout << "Volume xxx: " << volume << std::endl;
+
+    // free memory, thus deallocate memory of xarray
+    for (auto& p : phiStacks) delete p;
+    for (auto& p : originTensorStacks) delete p;
 };
+
+void quadratureTask(const Scene& scene) {}
+
+void buildOctree(const Scene& scene, std::vector<Node>& nodes, const uvector3& min, const uvector3& max) {}
 }; // namespace algoim::Organizer

algoim/organizer/primitive.hpp

@@ -12,6 +12,7 @@
 #include "xarray.hpp"
 #include "binomial.hpp"
 #include "bernstein.hpp"
+#include <memory>
 
 namespace algoim::Organizer
 {
@@ -588,5 +589,17 @@ void makeCylinder(xarray<real, 3>& tensor, uvector3 startPt, uvector3 endPt, rea
     // TODO:
 }
 
+struct Scene {
+    std::vector<tensor3> polys;
+    real*                boolDescription; // blobtree
+    uvector3             min, max;
+};
+
+struct Node {
+    std::vector<int> polyIntersectIndices;      // 同一poly会出现在不同node的polyIndices中
+    std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
+    int              children[8];               // octree
+    uvector3         min, max;
+};
 
 } // namespace algoim::Organizer

algoim/sparkstack.hpp

@@ -107,7 +107,7 @@ public:
     ~SparkStack()
     {
         pos() -= len_;
-        // std::cout << "Here!" << std::endl;
+        // std::cout << "Here! len = " << len_ << std::endl;
     }
 };
 

algoim/xarray.hpp

@@ -182,12 +182,15 @@ public:
 
     xarray(const xarray& x)
     {
-        this->data_ = nullptr;
-        this->ext_  = x.ext();
-        if (x.data_) {
-            algoim_spark_alloc(T, *this);
-            *this = x;
-        };
+        // this->data_ = nullptr;
+        // this->ext_  = x.ext();
+        // if (x.data_) {
+        //     algoim_spark_alloc(T, *this);
+        //     *this = x;
+        // };
+        // 这里只能浅拷贝，因为不能在这里用algoim_spark_alloc分配新的内存
+        this->data_ = x.data_;
+        this->ext_  = x.ext_;
     };
 
     xarray()

gjj/primitiveDebug.hpp

@@ -53,7 +53,7 @@ void casePolyhedron1()
     // ps.emplace_back(std::make_shared<SphereDesc>(0.8, 0., 1.));
     // ps.emplace_back(std::make_shared<PowerTensor>(makeSphere(0.2)));
     // ps.emplace_back(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
-    auto basicTask = BasicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
+    basicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
 }
 
 void casePolyhedron2()
@@ -85,7 +85,7 @@ void casePolyhedron2()
     // ps.emplace_back(std::make_shared<SphereDesc>(0.8, 0., 1.));
     // ps.emplace_back(std::make_shared<PowerTensor>(makeSphere(0.2)));
     // ps.emplace_back(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
-    auto basicTask = BasicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
+    basicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
 }
 
 void casePolyhedron3()
@@ -115,23 +115,51 @@ void casePolyhedron3()
                                                             2,
                                                             3,
                                                             7, // top
-                                                            7,
-                                                            3,
-                                                            1,
-                                                            5, // front
+                                                            // 7,
+                                                            // 3,
+                                                            // 1,
+                                                            // 5, // front
                                                             2,
                                                             6,
                                                             4,
                                                             0}; // back
-    std::vector<int>                            scan     = {4, 8, 12, 16};
-    auto basicTask = BasicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
+    std::vector<int>                            scan     = {4, 8, 12};
+    basicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
 }
 
-void case1()
+void casePolyhedronSphere()
 {
-    auto phi0      = std::make_shared<SphereDesc>(SphereDesc(0.8, 0., 1.));
-    // SphereDesc sphere(0.8, 0., 1.);
-    auto basicTask = BasicTask(phi0);
+    // PI * r^3 / 6
+    auto                  phi0     = std::make_shared<SphereDesc>(SphereDesc(0.8, uvector3(-0.8, 0.8, -0.8), 1.));
+    std::vector<uvector3> vertices = {uvector3(-0.8, -0.8, -0.8),
+                                      uvector3(-0.8, -0.8, 0.8),
+                                      uvector3(-0.8, 0.8, -0.8),
+                                      uvector3(-0.8, 0.8, 0.8),
+                                      uvector3(0.8, -0.8, -0.8),
+                                      uvector3(0.8, -0.8, 0.8),
+                                      uvector3(0.8, 0.8, -0.8),
+                                      uvector3(0.8, 0.8, 0.8)};
+    std::vector<int>      indices  = {3,
+                                      2,
+                                      0,
+                                      1, // left
+                                      6,
+                                      2,
+                                      3,
+                                      7, // top
+                                      2,
+                                      6,
+                                      4,
+                                      0}; // back
+    std::vector<int>      scan     = {4, 8, 12};
+    std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions(2);
+    primitiveDescriptions[0] = phi0;
+    primitiveDescriptions[1] = std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan));
+    basicTask(primitiveDescriptions);
 }
 
-void testPrimitive() { casePolyhedron3(); }
+void testPrimitive()
+{
+    casePolyhedron1();
+    // casePolyhedronSphere();
+}