diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 9c7dc02..caecf50 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -37,48 +37,113 @@ bool keepQuadraturePoint(std::vector& originTensors, const uvector3& or return true; } -class BasicTask +// std::vector> primitives; + +// BasicTask(std::vector> ps) {}; + + +void basicTask(const std::shared_ptr& p, int q = 20, real xmin = -1, real xmax = 1) { -public: - // std::vector> primitives; + real volume = 0; + auto integrand = [](const uvector& x) { return 1.0; }; + uvector3 range = xmax - xmin; + + + if (auto pt = std::dynamic_pointer_cast(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 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& x, real w) { + if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin)); + }); + } else if (auto pt = std::dynamic_pointer_cast(p)) { + const int faceCount = pt->indexInclusiveScan.size(); + assert(faceCount > 1); + std::vector planeTensors(faceCount, tensor3(nullptr, 2)); + algoim_spark_alloc(real, planeTensors); + + tensor3 compositeTensor(nullptr, 1 + faceCount); + algoim_spark_alloc(real, compositeTensor); - // BasicTask(std::vector> 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 p) - { - int q = 20; - real volume = 0; - real xmin = -1; - real xmax = 1; - auto integrand = [](const uvector& x) { return 1.0; }; - uvector3 range = xmax - xmin; + uvector 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(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& 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>& primitives, int q = 20, real xmin = -1, real xmax = 1) +{ + std::vector*> phiStacks; + std::vector phis; + + std::vector*> originTensorStacks; + std::vector originTensors; + real volume; + auto integrand = [](const uvector& x) { return 1.0; }; + uvector3 range = xmax - xmin; + uvector testX(0., 0.75, 0.2); + for (int i = 0; i < primitives.size(); i++) { + if (auto pt = std::dynamic_pointer_cast(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 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& x, real w) { - if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin)); - }); - } else if (auto pt = std::dynamic_pointer_cast(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(primitives[i])) { const int faceCount = pt->indexInclusiveScan.size(); assert(faceCount > 1); std::vector 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 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& 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> primitives) - { - std::vector*> phiStacks; - std::vector phis; - - std::vector*> originTensorStacks; - std::vector originTensors; - int q = 10; - real volume; - uvector3 xmin = 0; - uvector3 xmax = 1; - auto integrand = [](const uvector& x) { return 1.0; }; - uvector3 range = xmax - xmin; - - for (int i = 0; i < primitives.size(); i++) { - if (auto pt = std::dynamic_pointer_cast(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(primitives[i])) { - const int faceCount = pt->indexInclusiveScan.size(); - assert(faceCount > 1); - std::vector 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& 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& 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& nodes, const uvector3& min, const uvector3& max) {} }; // namespace algoim::Organizer diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 0e38e40..642ef74 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -12,6 +12,7 @@ #include "xarray.hpp" #include "binomial.hpp" #include "bernstein.hpp" +#include namespace algoim::Organizer { @@ -588,5 +589,17 @@ void makeCylinder(xarray& tensor, uvector3 startPt, uvector3 endPt, rea // TODO: } +struct Scene { + std::vector polys; + real* boolDescription; // blobtree + uvector3 min, max; +}; + +struct Node { + std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 + std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 + int children[8]; // octree + uvector3 min, max; +}; } // namespace algoim::Organizer \ No newline at end of file diff --git a/algoim/sparkstack.hpp b/algoim/sparkstack.hpp index 8cb4573..7016af5 100644 --- a/algoim/sparkstack.hpp +++ b/algoim/sparkstack.hpp @@ -107,7 +107,7 @@ public: ~SparkStack() { pos() -= len_; - // std::cout << "Here!" << std::endl; + // std::cout << "Here! len = " << len_ << std::endl; } }; diff --git a/algoim/xarray.hpp b/algoim/xarray.hpp index c586bab..78861ac 100644 --- a/algoim/xarray.hpp +++ b/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() diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 0e6a0e0..cc30cdb 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -53,7 +53,7 @@ void casePolyhedron1() // ps.emplace_back(std::make_shared(0.8, 0., 1.)); // ps.emplace_back(std::make_shared(makeSphere(0.2))); // ps.emplace_back(std::make_shared(MeshDesc(vertices, indices, scan))); - auto basicTask = BasicTask(std::make_shared(MeshDesc(vertices, indices, scan))); + basicTask(std::make_shared(MeshDesc(vertices, indices, scan))); } void casePolyhedron2() @@ -85,7 +85,7 @@ void casePolyhedron2() // ps.emplace_back(std::make_shared(0.8, 0., 1.)); // ps.emplace_back(std::make_shared(makeSphere(0.2))); // ps.emplace_back(std::make_shared(MeshDesc(vertices, indices, scan))); - auto basicTask = BasicTask(std::make_shared(MeshDesc(vertices, indices, scan))); + basicTask(std::make_shared(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 scan = {4, 8, 12, 16}; - auto basicTask = BasicTask(std::make_shared(MeshDesc(vertices, indices, scan))); + std::vector scan = {4, 8, 12}; + basicTask(std::make_shared(MeshDesc(vertices, indices, scan))); } -void case1() +void casePolyhedronSphere() { - auto phi0 = std::make_shared(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(0.8, uvector3(-0.8, 0.8, -0.8), 1.)); + std::vector 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 indices = {3, + 2, + 0, + 1, // left + 6, + 2, + 3, + 7, // top + 2, + 6, + 4, + 0}; // back + std::vector scan = {4, 8, 12}; + std::vector> primitiveDescriptions(2); + primitiveDescriptions[0] = phi0; + primitiveDescriptions[1] = std::make_shared(MeshDesc(vertices, indices, scan)); + basicTask(primitiveDescriptions); } -void testPrimitive() { casePolyhedron3(); } \ No newline at end of file +void testPrimitive() +{ + casePolyhedron1(); + // casePolyhedronSphere(); +} \ No newline at end of file