diff --git a/algoim/organizer/blobtree.hpp b/algoim/organizer/blobtree.hpp index 73e1e6e..b185743 100644 --- a/algoim/organizer/blobtree.hpp +++ b/algoim/organizer/blobtree.hpp @@ -37,6 +37,21 @@ public: std::vector primitiveNodeIdx; // 由primitive数组下标指向node中对应leaf的下标 + BlobTree(const BlobTree& other) + { + structure = other.structure; + primitiveNodeIdx = other.primitiveNodeIdx; + } + + BlobTree& operator=(const BlobTree& other) + { + if (this != &other) { + structure = other.structure; + primitiveNodeIdx = other.primitiveNodeIdx; + } + return *this; + } + BlobTree() {} }; @@ -44,7 +59,11 @@ void propagate(BlobTree& tree, unsigned int nodeIdx, bool in) { const std::size_t rootIdx = tree.structure.size() - 1; auto& node = tree.structure[nodeIdx]; - assert(node.inOut == 0); + // if (node.inOut != NODE_IN_OUT_UNKNOWN) { + // int aaa = 1; + // int bb = 1; + // } + assert(node.inOut == NODE_IN_OUT_UNKNOWN); node.inOut = in ? NODE_IN : NODE_OUT; for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) { const auto& nowNode = tree.structure[nowIdx]; diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 0092912..98800d5 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -93,41 +93,21 @@ int findFirst(const uvector& v, int val, int startIdx = 0) return -1; } -std::pair, uvector> getOneEightCellAABB(const uvector3& min, - const uvector3& max, - const uvector side, - int negativeRep = -1) +AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector side, int negativeRep = -1) { - std::pair, uvector> res = {min, max}; - uvector3 mid = (min + max) / 2; + AABB res = fatherAABB; + uvector3 mid = res.center(); for (int i = 0; i < 3; ++i) { if (side(i) == negativeRep) { - res.second(i) = mid(i); + res.max(i) = mid(i); } else { - res.first(i) = mid(i); + res.min(i) = mid(i); } } return res; } } // namespace detail -// TODO: delete it -bool isPointInside(const std::vector& originTensors, const uvector3& originPt) -{ - for (auto& t : originTensors) { - if (evalPower(t, originPt) >= 0) { return false; } - } - return true; -} - -// bool isPointInside(const std::vector& completeTensorReps, const uvector3& originPt) -// { -// for (auto& completeTensorRep : completeTensorReps) { -// if (isPointInside(completeTensorRep.originalPower, originPt)) { return false; } -// } -// return true; -// } - // 这里blobTree是拷贝传参 bool keepQuadraturePoint(const Scene& scene, organizer::BlobTree blobTree, @@ -138,7 +118,7 @@ bool keepQuadraturePoint(const Scene& scene, const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; std::vector primitiveInOuts(polyIntersectIndices.size()); for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) { - primitiveInOuts[i] = isPointInside(scene.polys[polyIntersectIndices[i]].rawPower, originPt); + primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt); } int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); assert(res == organizer::NODE_IN || res == organizer::NODE_OUT); @@ -149,7 +129,6 @@ bool keepQuadraturePoint(const Scene& scene, // BasicTask(std::vector> ps) {}; - const std::array sides = {-1, 1}; struct BasicTaskRes { @@ -174,11 +153,11 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, algoim_spark_alloc(real, halfCellPoly); int subIdx = detail::binaryToDecimal(mark, -1); auto& subNode = subNodes[subIdx]; - bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 + bernstein::deCasteljau(poly, subNode.aabb.min, subNode.aabb.max, halfCellPoly); // 求1/8空间下的表达 if (bernstein::uniformSign(halfCellPoly) != 1) { subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); } else { - organizer::traverse(subNode.blobTree, polyIntersectIndex, organizer::NODE_OUT); + organizer::traverse(subNode.blobTree, polyIntersectIndex, false); } } else { for (auto side : sides) { @@ -199,19 +178,15 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector subNodes(CHILD_NUM); // std::array subNodes; // intermediateNodes.resize(lastIdx + 8); int subIdx = 0; for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { - auto nodeAABB = detail::getOneEightCellAABB(node.min, node.max, j(), 0); - subNodes[subIdx].min = nodeAABB.first; - subNodes[subIdx].max = nodeAABB.second; + subNodes[subIdx].aabb = detail::getOneEightCellAABB(node.aabb, j(), 0); subNodes[subIdx].blobTree = node.blobTree; } - uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; + uvector3 nodeMid = node.aabb.center(); for (int i = 0; i < polyIntersectIndices.size(); ++i) { const int polyIntersectIndex = polyIntersectIndices[i]; const auto& poly = scene.polys[polyIntersectIndex]; @@ -229,18 +204,18 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector& leaves) { - int a = 0; const std::vector& polyIntersectIndices = node.polyIntersectIndices; - if (polyIntersectIndices.size() <= 4) { + if (polyIntersectIndices.size() <= 3) { leaves.emplace_back(node); return; } - const uvector3& nowNodeMin = node.min; - const uvector3& nowNodeMax = node.max; // std::array subNodes; std::vector subNodes(CHILD_NUM); // intermediateNodes.resize(lastIdx + 8); int subIdx = 0; for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { - auto nodeAABB = detail::getOneEightCellAABB(node.min, node.max, j(), 0); - subNodes[subIdx].min = nodeAABB.first; - subNodes[subIdx].max = nodeAABB.second; + subNodes[subIdx].aabb = detail::getOneEightCellAABB(node.aabb, j(), 0); subNodes[subIdx].blobTree = node.blobTree; subNodes[subIdx].polyIntersectIndices = std::vector(); } - uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; + uvector3 nodeMid = node.aabb.center(); for (int i = 0; i < polyIntersectIndices.size(); ++i) { const int polyIntersectIndex = polyIntersectIndices[i]; const auto& poly = scene.polys[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; + } tensor3 subcellPoly(nullptr, poly.compositedBernstein.ext()); algoim_spark_alloc(real, subcellPoly); - bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].min, subNodes[subIdx].max, subcellPoly); + bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly); int sign = bernstein::uniformSign(subcellPoly); if (sign == 1) { - organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false); } else if (sign == -1) { - organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_IN); + // 采样一个点，带入所有原始primitive看是否是冗余部分 + uvector3 sampleX = subNodes[subIdx].aabb.center() * scene.boundary.size() + scene.boundary.min; + if (isInsidePowers(poly.rawPower, sampleX)) + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, true); + else + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false); } else { subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); } @@ -479,7 +459,7 @@ void basicTask(const std::vector>& primitives, in ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { auto realX = x * range + xmin; - if (isPointInside(originTensors, realX)) volume += w * integrand(realX); + if (isInsidePowers(originTensors, realX)) volume += w * integrand(realX); }); volume *= pow(xmax - xmin, 3); std::cout << "Volume xxx: " << volume << std::endl; @@ -489,18 +469,18 @@ void basicTask(const std::vector>& primitives, in for (auto& p : originTensorStacks) delete p; }; -BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 20) +BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 10) { auto integrand = [](const uvector& x) { return 1.0; }; real volume = 0., surf = 0.; - auto range = node.max - node.min; - ImplicitPolyQuadrature<3> ipquad(scene.polys); + auto range = node.aabb.size(); + ImplicitPolyQuadrature<3> ipquad(scene.polys, node.polyIntersectIndices); ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { - auto realX = x * range + node.min; + auto realX = x * range + node.aabb.min; if (keepQuadraturePoint(scene, blobTree, node, realX)) { volume += w * integrand(realX); } }); ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { - surf += w * integrand(x * range + node.min); + surf += w * integrand(x * range + node.aabb.min); }); return {volume, surf}; } @@ -534,6 +514,7 @@ void quadratureScene(const std::vector>& primitiv AABB aabb; makeSphere(*pt, originTensor, aabb); detail::powerTransformation(range, xmin, originTensor, transformedTensor); + aabb.normalize(range, xmin); real testEvaOri = evalPower(originTensor, testX); @@ -553,7 +534,9 @@ void quadratureScene(const std::vector>& primitiv algoim_spark_alloc(real, compositeTensor); AABB aabb; makeMesh(*pt, compositeTensor, planeTensors, aabb); + detail::powerTransformation(range, xmin, compositeTensor); + aabb.normalize(range, xmin); real testEvalPower = evalPower(compositeTensor, testX); @@ -576,6 +559,8 @@ void quadratureScene(const std::vector>& primitiv AABB aabb; makeCylinder(*pt, compositeTensor, rawTensors, aabb); detail::powerTransformation(range, xmin, compositeTensor); + aabb.normalize(range, xmin); + real testEvalPower = evalPower(compositeTensor, testX); detail::power2BernsteinTensor(compositeTensor, phi); @@ -595,6 +580,8 @@ void quadratureScene(const std::vector>& primitiv AABB aabb; makeCone(*pt, compositeTensor, rawTensors, aabb); detail::powerTransformation(range, xmin, compositeTensor); + aabb.normalize(range, xmin); + real testEvalPower = evalPower(compositeTensor, testX); detail::power2BernsteinTensor(compositeTensor, phi); @@ -605,14 +592,14 @@ void quadratureScene(const std::vector>& primitiv } } - Scene scene{completeTensorReps, blobTree}; + Scene scene{completeTensorReps, AABB(xmin, xmax)}; OcTreeNode rootNode(0, 1, blobTree); for (int i = 0; i < completeTensorReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } buildOcTreeV0(scene, rootNode, leaves); for (const auto& leaf : leaves) { auto basicRes = basicTask(scene, blobTree, leaf, q); - volume += basicRes.volume * prod(leaf.max - leaf.min); + volume += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min); } volume *= prod(xmax - xmin); diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index b5933d4..fd271e1 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -171,6 +171,11 @@ real evalBernstein(const xarray& phi, const uvector& x) return bernstein::evalBernsteinPoly(phi, x); } +void inverseTensor(tensor3& tensor) +{ + for (auto i = tensor.loop(); ~i; ++i) { tensor.l(i) = -tensor.l(i); } +} + // class PowerTensor : public Primitive // { // public: @@ -603,7 +608,8 @@ public: } }; -struct AABB { +class AABB +{ public: uvector3 min = uvector3(std::numeric_limits::max()); uvector3 max = uvector3(std::numeric_limits::lowest()); @@ -630,6 +636,20 @@ public: uvector3 center() const { return (min + max) / 2; } uvector3 size() const { return max - min; } + + bool intersect(const AABB& aabb) const + { + for (int i = 0; i < 3; ++i) { + if (min(i) > aabb.max(i) || max(i) < aabb.min(i)) { return false; } + } + return true; + } + + void normalize(const uvector3& scale, const uvector3& boundaryMin) + { + min = (min - boundaryMin) / scale; + max = (max - boundaryMin) / scale; + } }; void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector>& planeTensors, AABB& aabb) @@ -663,7 +683,7 @@ void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector polys; - organizer::BlobTree blobTree; + // organizer::BlobTree blobTree; + AABB boundary; // uvector3 min, max; }; @@ -795,21 +816,21 @@ struct OcTreeNode { std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 // std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 // std::array children; // octree - uvector3 min, max; + AABB aabb; organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递 // Node() { children.fill(-1); } OcTreeNode() = default; OcTreeNode(const uvector3& min_, const uvector3& max_, const organizer::BlobTree& blobTree_) - : min(min_), max(max_), blobTree(blobTree_) {}; + : aabb(min_, max_), blobTree(blobTree_) {}; + OcTreeNode(const AABB& aabb_, const organizer::BlobTree& blobTree_) : aabb(aabb_), blobTree(blobTree_) {}; // 深拷贝 OcTreeNode(const OcTreeNode& node) : polyIntersectIndices(node.polyIntersectIndices), // polyFullyContainedIndices(node.polyFullyContainedIndices), - min(node.min), - max(node.max), + aabb(node.aabb), blobTree(node.blobTree) { int a = 1; diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index 2eccac2..4e94891 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -569,15 +569,25 @@ struct ImplicitPolyQuadrature { build(true, false); } - ImplicitPolyQuadrature(const std::vector& completeTensorRes) + ImplicitPolyQuadrature(const std::vector& completeTensorReps) { - for (const auto& ctr : completeTensorRes) { + for (const auto& ctr : completeTensorReps) { auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray(true)); if (!detail::maskEmpty(mask)) phi.push_back(ctr.compositedBernstein, mask); } build(true, false); } + ImplicitPolyQuadrature(const std::vector& completeTensorReps, + const std::vector& polyIndices) + { + for (auto i : polyIndices) { + auto mask = detail::nonzeroMask(completeTensorReps[i].compositedBernstein, booluarray(true)); + if (!detail::maskEmpty(mask)) phi.push_back(completeTensorReps[i].compositedBernstein, mask); + } + build(true, false); + } + // Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks ImplicitPolyQuadrature(const xarray& p, const booluarray& pmask, diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 7e53327..e85e941 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -15,6 +15,7 @@ #include "binomial.hpp" #include "real.hpp" +#include "sparkstack.hpp" #include "uvector.hpp" #include "vector" #include "xarray.hpp" @@ -122,13 +123,28 @@ void casePolyhedronSphere() basicTask(primitiveDescriptions); } +void caseTwoCube() +{ + auto phi0 = std::make_shared(CuboidDesc(0.2, 0.6)); + auto phi1 = std::make_shared(CuboidDesc(-0.2, 0.6)); + basicTask({phi0, phi1}); +} + +void caseCubeSphere() +{ + auto phi0 = std::make_shared(CuboidDesc(0.2, 0.6)); + auto phi1 = std::make_shared(SphereDesc(0.6, -0.2, 1.)); + basicTask({phi0, phi1}); +} + void caseScene() { const int PRIMITIVE_CNT = 6; std::vector> primitiveDescriptions(PRIMITIVE_CNT); - primitiveDescriptions[2] = std::make_shared(CuboidDesc(0., 1.6)); - primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); - primitiveDescriptions[0] = std::make_shared(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.)); + primitiveDescriptions[0] = std::make_shared(CuboidDesc(0., 1.6)); + // primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); + primitiveDescriptions[1] = std::make_shared(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.)); + primitiveDescriptions[2] = std::make_shared(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.)); primitiveDescriptions[3] = std::make_shared(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.)); std::vector pyramidBottomVertices = { uvector3{-1, -1, 0}, @@ -166,7 +182,8 @@ void testDeCasteljau() auto phiDesc = std::make_shared(SphereDesc(0.8, uvector3(0), 1.)); tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3); algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein); - makeSphere(*phiDesc, tensor); + AABB aabb; + makeSphere(*phiDesc, tensor, aabb); uvector3 xmax = 1, xmin = -1, range = xmax - xmin; Organizer::detail::powerTransformation(range, xmin, tensor, tensor01); Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein); @@ -193,7 +210,8 @@ void testSubDivideWithDeCasteljau() auto phiDesc = std::make_shared(SphereDesc(0.7, uvector3(0.8), 1.)); tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3); algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein); - makeSphere(*phiDesc, tensor); + AABB aabb; + makeSphere(*phiDesc, tensor, aabb); uvector3 xmin = uvector3(-1., -1.3, -1.6), xmax = uvector3(1.6, 1.6, 2.3), range = xmax - xmin; Organizer::detail::powerTransformation(range, xmin, tensor, tensor01); Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein); @@ -202,27 +220,44 @@ void testSubDivideWithDeCasteljau() // bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001), // transformedTensorBernstein); for (MultiLoop<3> j(0, 2); ~j; ++j) { - if (all(j() == uvector3i(1, 0, 0))) { - int aaa = 1; - int bbb = 1; - } - auto aabb = Organizer::detail::getOneEightCellAABB(xmin, xmax, j(), 0); - bernstein::deCasteljau(tensorBernstein, aabb.first, aabb.second, transformedTensorBernstein); - int sign = bernstein::uniformSign(transformedTensorBernstein); + auto aabbSub = Organizer::detail::getOneEightCellAABB(AABB{xmin, xmax}, j(), 0); + bernstein::deCasteljau(tensorBernstein, aabbSub.min, aabbSub.max, transformedTensorBernstein); + auto transformedTensorVec = xarray2StdVector(transformedTensorBernstein); + int sign = bernstein::uniformSign(transformedTensorBernstein); std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(2) << "), sign: " << sign << std::endl; + if (all(j() == uvector3i(0, 0, 1))) { + for (MultiLoop<3> k(0, 2); ~k; ++k) { + auto aabbSubSub = Organizer::detail::getOneEightCellAABB(aabbSub, k(), 0); + bernstein::deCasteljau(tensorBernstein, aabbSubSub.min, aabbSubSub.max, transformedTensorBernstein); + int sign = bernstein::uniformSign(transformedTensorBernstein); + std::cout << "k(): (" << k(0) << ", " << k(1) << ", " << k(2) << "), sign: " << sign << std::endl; + } + } } - uvector3 testX = uvector3(0.5); - - bernstein::deCasteljau(tensorBernstein, uvector3(-1., -1.3, 0.35), uvector3(0.3, 0.15, 2.3), transformedTensorBernstein); - int sign = bernstein::uniformSign(transformedTensorBernstein); - real evalX = evalPower(tensor, testX); + uvector3 testX = uvector3(0.5); + real evalX = evalPower(tensor, testX); std::cout << "evalX original rep = " << evalX << std::endl; + evalX = evalPower(tensor01, testX); + std::cout << "evalX power within 0-1 = " << evalX << std::endl; evalX = bernstein::evalBernsteinPoly(tensorBernstein, testX); std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl; - evalX = bernstein::evalBernsteinPoly(transformedTensorBernstein, testX); - std::cout << "evalX bernstein after Decasteljau = " << evalX << std::endl; - std::cout << "sign = " << sign << std::endl; + + // 切比雪夫 + xmin = uvector3(-1., -1.3, -1.6), xmax = uvector3(1.6, 1.6, 2.3); + auto phi = [](const uvector &x) { + return (x(0) - 0.8) * (x(0) - 0.8) + (x(1) - 0.8) * (x(1) - 0.8) + (x(2) - 0.8) * (x(2) - 0.8) - 0.1 * 0.1; + }; + xarray phipoly(nullptr, 3), subPhiPoly(nullptr, 3); + algoim_spark_alloc(real, phipoly, subPhiPoly); + bernstein::bernsteinInterpolate<3>([&](const uvector &x) { return phi(xmin + x * (xmax - xmin)); }, phipoly); + xmin = 0, xmax = 1; + for (MultiLoop<3> j(0, 2); ~j; ++j) { + auto aabbSub = Organizer::detail::getOneEightCellAABB(AABB{xmin, xmax}, j(), 0); + bernstein::deCasteljau(phipoly, aabbSub.min, aabbSub.max, subPhiPoly); + int sign = bernstein::uniformSign(subPhiPoly); + std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(2) << "), sign: " << sign << std::endl; + } } void testBlob() @@ -238,16 +273,23 @@ void testBlob() std::cout << sizeof(int) << std::endl; } -void vectorDebug() +void testTensorInverse() { - struct A { - std::vector data; - A() = default; - }; - - A a; - int s = a.data.size(); - std::cout << "size: " << s << std::endl; + auto phiDesc = std::make_shared(uvector3(0.2, 0.3, 0.5), uvector3(0.8, 1.1, 0.9)); + const int faceCnt = 6; + std::vector tensors(faceCnt, tensor3(nullptr, 2)); + std::vector *> sparkStackPtrs(faceCnt); + tensor3 tensor(nullptr, 1 + faceCnt); + algoimSparkAllocHeapVector(sparkStackPtrs, tensors); + algoim_spark_alloc(real, tensor); + AABB aabb; + makeMesh(*phiDesc, tensor, tensors, aabb); + uvector3 testX(0.5, 0.7, 0.2); + real evalX = evalPower(tensor, testX); + std::cout << "before inverse, evalX = " << evalX << std::endl; + inverseTensor(tensor); + evalX = evalPower(tensor, testX); + std::cout << "after inverse, evalX = " << evalX << std::endl; } void testPrimitive() @@ -256,6 +298,8 @@ void testPrimitive() // casePolyhedronSphere(); // testSubDivideWithDeCasteljau(); // testBlob(); - caseScene(); - // vectorDebug(); + // caseScene(); + caseCubeSphere(); + // caseTwoCube(); + // testTensorInverse(); } \ No newline at end of file