polyhedron tests

algoim/organizer/organizer.hpp

@@ -15,6 +15,7 @@
 #include "binomial.hpp"
 
 #include "real.hpp"
+#include "sparkstack.hpp"
 #include "uvector.hpp"
 #include "vector"
 #include "xarray.hpp"
@@ -26,6 +27,16 @@
 
 namespace algoim::Organizer
 {
+
+bool keepQuadraturePoint(std::vector<tensor3>& originTensors, const uvector3& originPt)
+{
+    // TODO: using blobtree to filter tensors
+    for (auto& t : originTensors) {
+        if (evalPower(t, originPt) >= 0) { return false; }
+    }
+    return true;
+}
+
 class BasicTask
 {
 public:
@@ -33,6 +44,7 @@ public:
 
     // BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
 
+
     BasicTask(std::shared_ptr<PrimitiveDesc> p)
     {
         int      q         = 20;
@@ -74,36 +86,84 @@ public:
             makeMesh(*pt, compositeTensor, planeTensors);
             detail::powerTransformation(range, xmin, compositeTensor);
 
+            auto compositeTensorStdVector = xarray2StdVector(compositeTensor);
+
+            uvector<real, 3> testX(0.8, 0.8, 0.8);
+            real             textEvalPower = evalPower(compositeTensor, testX);
+
             tensor3 phi(nullptr, 1 + faceCount);
             algoim_spark_alloc(real, phi);
             detail::power2BernsteinTensor(compositeTensor, phi);
 
+            int                       quadraturePointCount = 0;
             ImplicitPolyQuadrature<3> ipquad(phi);
             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 << "quadraturePointCount: " << quadraturePointCount << std::endl;
         }
         volume *= pow(xmax - xmin, 3);
         std::cout << "Volume xxx: " << volume << std::endl;
     };
 
-    // BasicTask(std::shared_ptr<PowerTensorComplex> pc)
-    // {
-    //     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;
-
-    //     detail::powerTransformation(range, xmin, pc->compositeTensor);
-    //     auto                      primitive = BernsteinTensorComplex(*pc);
-    //     ImplicitPolyQuadrature<3> ipquad(primitive.compositeTensor);
-    //     ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
-    //         if (primitive.isInside(x)) volume += w * integrand(xmin + x * (xmax - xmin));
-    //     });
-    //     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;
+    };
 };
 }; // namespace algoim::Organizer

algoim/organizer/primitive.hpp

@@ -18,7 +18,6 @@ namespace algoim::Organizer
 
 namespace detail
 {
-// void compositePower(const std::vector<xarray<real, 3>>& powers) {}
 
 template <int N>
 void power2BernsteinTensor(const xarray<real, N>& phiPower, xarray<real, N>& phiBernsetin)
@@ -50,6 +49,38 @@ void power2BernsteinTensor(const xarray<real, N>& phiPower, xarray<real, N>& phi
     }
 }
 
+void powerTransformation(const uvector<real, 3>& scale,
+                         const uvector<real, 3>& bias,
+                         const xarray<real, 3>&  origin,
+                         xarray<real, 3>&        res)
+{
+    assert(all(origin.ext() == res.ext()));
+    std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
+    const auto&                                 ext = origin.ext();
+    for (int dim = 0; dim < 3; ++dim) {
+        dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
+        for (int degree = 0; degree < ext(dim); ++degree) {
+            const real* binomDegree = Binomial::row(degree);
+            dimOrderExpansion[dim][degree].reserve(degree + 1);
+            // 根据二项定理展开
+            for (int i = 0; i <= degree; ++i) {
+                dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
+            }
+        }
+    }
+    for (auto i = origin.loop(); ~i; ++i) {
+        // 迭代器必须按照坐标的升序进行访问，即，访问ijk时，(i-1)jk，i(j-1)k，ij(k-1)必须已经访问过
+        real base = origin.l(i);
+        res.l(i)  = 0;
+        auto exps = i();
+        for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
+            real item = base;
+            for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
+            res.m(j()) += item;
+        }
+    }
+}
+
 void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
 {
     std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
@@ -77,6 +108,35 @@ void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>&
         }
     }
 }
+
+static void compositePower(const std::vector<xarray<real, 3>>& powers,
+                           int                                 powerIdx,
+                           uvector<int, 3>                     powerSum,
+                           real                                factor,
+                           xarray<real, 3>&                    res)
+
+{
+    if (powerIdx == 0) {
+        {
+            uvector3 ext(1, 1, 1);
+            for (auto& t : powers) { ext += t.ext() - 1; }
+            assert(all(ext == res.ext()));
+        }
+        xarrayInit(res);
+    }
+    if (powerIdx == powers.size()) {
+        res.m(powerSum) += factor;
+        return;
+    }
+    auto& power = powers[powerIdx];
+    for (auto i = power.loop(); ~i; ++i) {
+        if (power.l(i) == 0) {
+            factor = 0;
+            continue;
+        }
+        compositePower(powers, powerIdx + 1, powerSum + i(), factor * power.l(i), res);
+    }
+}
 } // namespace detail
 
 class Primitive
@@ -501,7 +561,7 @@ void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<
             assert(dot(N, vertices[indices[j]]) + d < std::numeric_limits<real>::epsilon());
         }
     }
-    // compositePower(planeTensors, 0, 0, 1, tensor);
+    detail::compositePower(planeTensors, 0, 0, 1, tensor);
     // return PowerTensorComplex(planeTensors, tensor);
 };
 

algoim/quadrature_multipoly.hpp

@@ -41,6 +41,9 @@ inline namespace v1
 {
 namespace detail
 {
+/**
+ * g 和 gmask之所以用ptr，是因为其不一定存在，方便传入空
+ */
 template <int N>
 booluarray<N, ALGOIM_M> mask_driver(const xarray<real, N>&         f,
                                     const booluarray<N, ALGOIM_M>& fmask,
@@ -54,7 +57,7 @@ booluarray<N, ALGOIM_M> mask_driver(const xarray<real, N>&         f,
             if (fmask(i()) && (!gmask || (*gmask)(i()))) overlap = true;
         if (!overlap) return;
 
-        real             eps = 0.015625 / ALGOIM_M;
+        real             eps = 0.05 / ALGOIM_M;
         uvector<real, N> xa, xb;
         for (int dim = 0; dim < N; ++dim) {
             xa(dim) = real(a(dim)) / ALGOIM_M - eps;
@@ -530,7 +533,11 @@ struct ImplicitPolyQuadrature {
     ImplicitPolyQuadrature(const xarray<real, N>& p)
     {
         auto mask = detail::nonzeroMask(p, booluarray<N, ALGOIM_M>(true));
-        if (!detail::maskEmpty(mask)) phi.push_back(p, mask);
+
+        if (!detail::maskEmpty(mask)) {
+            phi.push_back(p, mask);
+            int a = 1;
+        }
         build(true, false);
     }
 

algoim/sparkstack.hpp

@@ -107,7 +107,7 @@ public:
     ~SparkStack()
     {
         pos() -= len_;
-        // std::cout << "Here!" << std::endl;
+        std::cout << "Here!" << std::endl;
     }
 };
 
@@ -120,6 +120,10 @@ public:
 #define algoim_spark_alloc_heap(T, ...)          new SparkStack<T>(__VA_ARGS__)
 #define algoim_spark_release_heap(SparkStackPtr) delete SparkStackPtr
 
+void algoimSparkAllocHeapVector(std::vector<SparkStack<real>*>& sparkStackPtrs, std::vector<xarray<real, 3>>& tensors)
+{
+    for (int i = 0; i < tensors.size(); ++i) { sparkStackPtrs.push_back(algoim_spark_alloc_heap(real, tensors[i])); }
+}
 } // namespace algoim
 
 #endif

gjj/primitiveDebug.hpp

@@ -28,7 +28,7 @@
 using namespace algoim::Organizer;
 using namespace algoim;
 
-void case0()
+void casePolyhedron1()
 {
     // std::vector<std::shared_ptr<PrimitiveDesc>> ps;
 
@@ -40,21 +40,52 @@ void case0()
                                       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
+                                      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
+    std::vector<int>      scan     = {4, 8, 12, 16, 20, 24};
+
+    // ps.emplace_back(std::make_shared<MeshDesc>(vertices, indices, scan));
+    // 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)));
+}
 
+void casePolyhedron2()
+{
+    // std::vector<std::shared_ptr<PrimitiveDesc>> ps;
+
+    // mesh
+    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
     };
-    std::vector<int>      indices  = {2, 3, 1, 0, // force break here
-                                      4, 0, 1, 5, //
-                                      4, 6, 2, 0, //
-                                      6, 7, 3, 2, //
-                                      7, 6, 4, 5, //
-                                      3, 7, 5, 1};
-    std::vector<int>      scan     = {4, 8, 12, 16, 20, 24};
+    std::vector<int> scan = {4, 8};
 
-    // ps.emplace_back(std::make_shared<PowerTensorComplex>(makeMesh(vertices, indixes)));
+    // ps.emplace_back(std::make_shared<MeshDesc>(vertices, indices, scan));
+    // 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(ps);
+    auto basicTask = BasicTask(std::make_shared<MeshDesc>(MeshDesc(vertices, indices, scan)));
 }
 
 void case1()
@@ -64,4 +95,26 @@ void case1()
     auto basicTask = BasicTask(phi0);
 }
 
-void testPrimitive() { case1(); }
+void case2()
+{
+    std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions;
+    std::vector<uvector3>                       vertices = {uvector3(-1.6, 0, 0),
+                                                            uvector3(-1.6, 0, 1.6),
+                                                            uvector3(-1.6, 1.6, 0),
+                                                            uvector3(-1.6, 1.6, 1.6),
+                                                            uvector3(1.6, 0, 0),
+                                                            uvector3(1.6, 0, 1.6),
+                                                            uvector3(1.6, 1.6, 0),
+                                                            uvector3(1.6, 1.6, 1.6)
+
+    };
+    std::vector<int>                            indices  = {2, 3, 1, 0, // force break here
+                                                            4, 0, 1, 5, //
+                                                            4, 6, 2, 0, //
+                                                            6, 7, 3, 2, //
+                                                            7, 6, 4, 5, //
+                                                            3, 7, 5, 1};
+    std::vector<int>                            scan     = {4, 8, 12, 16, 20, 24};
+}
+
+void testPrimitive() { casePolyhedron2(); }