From 1c4d06f102b0fac441071d404165307260fc9831 Mon Sep 17 00:00:00 2001 From: gjj Date: Wed, 18 Sep 2024 22:38:06 +0800 Subject: [PATCH] octree debug --- algoim/organizer/organizer.hpp | 147 +++++++++++++++++++++++------- algoim/organizer/primitive.hpp | 159 ++++++++++++++++++++++++++------- algoim/uvector.hpp | 11 ++- gjj/primitiveDebug.hpp | 93 ++++++++++++++++--- 4 files changed, 333 insertions(+), 77 deletions(-) diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 227a892..47e079b 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -37,7 +37,7 @@ void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) { assert(0 <= k && k < 3); assert(all(res.ext() == remove_component(phi.ext(), k))); - assert(0 < facePos && facePos < 1); + // assert(0 < facePos && facePos < 1); real* basisAlongK; int P = phi.ext(k); algoim_spark_alloc(real, &basisAlongK, P); @@ -93,13 +93,13 @@ 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) +std::pair, uvector> getOneEightCellAABB(const uvector3& min, + const uvector3& max, + const uvector side, + int negativeRep = -1) { - std::pair, uvector> res = {min, max}; - uvector3 mid = (min + max) / 2; + std::pair, uvector> res = {min, max}; + uvector3 mid = (min + max) / 2; for (int i = 0; i < 3; ++i) { if (side(i) == negativeRep) { res.second(i) = mid(i); @@ -138,15 +138,11 @@ 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]].originalPower, originPt); + primitiveInOuts[i] = isPointInside(scene.polys[polyIntersectIndices[i]].rawPower, originPt); } int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); assert(res == organizer::NODE_IN || res == organizer::NODE_OUT); - if (res == organizer::NODE_OUT) { - return false; - } else { - return true; - } + return res == organizer::NODE_IN; } // std::vector> primitives; @@ -167,15 +163,16 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, const int polyIntersectIndex, const tensor3& poly, - std::array& subNodes, - int startIdxToCheck, - uvector& mark) + // std::array& subNodes, + std::vector& subNodes, + int startIdxToCheck, + uvector& mark) { int zeroIdx = detail::findFirst(mark, 0, startIdxToCheck); if (zeroIdx == -1) { tensor3 halfCellPoly(nullptr, poly.ext()); algoim_spark_alloc(real, halfCellPoly); - int subIdx = detail::binaryToDecimal(mark); + int subIdx = detail::binaryToDecimal(mark, -1); auto& subNode = subNodes[subIdx]; bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 if (bernstein::uniformSign(halfCellPoly) != 1) { @@ -195,18 +192,19 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, // 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历 // 所以树结构是不需要的记录的 // void buildOcTree(const Scene& scene, std::vector& intermediateNodes, std::vector leaves, int nowNodeIdx) -void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector& leaves) +void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector& leaves) { const std::vector& polyIntersectIndices = node.polyIntersectIndices; if (polyIntersectIndices.size() <= 4) { leaves.emplace_back(node); return; } - const uvector3& nowNodeMin = node.min; - const uvector3& nowNodeMax = node.max; - std::array subNodes; + const uvector3& nowNodeMin = node.min; + const uvector3& nowNodeMax = node.max; + std::vector subNodes(CHILD_NUM); + // std::array subNodes; // intermediateNodes.resize(lastIdx + 8); - int subIdx = 0; + 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; @@ -274,7 +272,7 @@ void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector& leaves) +{ + const std::vector& polyIntersectIndices = node.polyIntersectIndices; + if (polyIntersectIndices.size() <= 4) { + 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].blobTree = node.blobTree; + subNodes[subIdx].polyIntersectIndices = std::vector(); + } + uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; + 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) { + tensor3 subcellPoly(nullptr, poly.compositedBernstein.ext()); + algoim_spark_alloc(real, subcellPoly); + bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].min, subNodes[subIdx].max, subcellPoly); + int sign = bernstein::uniformSign(subcellPoly); + if (sign == 1) { + // organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); + } else if (sign == -1) { + // organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_IN); + } else { + subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); + } + } + uvector3 testX(0.5); + real evalX = bernstein::evalBernsteinPoly(poly.compositedBernstein, testX); + // std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl; + } + for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { buildOcTreeV0(scene, subNodes[subIdx], leaves); } +} + /** single object quadrature */ void basicTask(const std::shared_ptr& p, int q = 20, real xmin = -1, real xmax = 1) { @@ -465,16 +508,15 @@ void quadratureScene(const std::vector>& primitiv const organizer::BlobTree& blobTree, const real q = 10) { - OcTreeNode rootNode(xmin, xmax, blobTree); std::vector leaves; std::vector*> tensorStacks; std::vector completeTensorReps; - real volume; + real volume = 0; auto integrand = [](const uvector& x) { return 1.0; }; uvector3 range = xmax - xmin; - uvector testX(0., 0.75, 0.2); + uvector testX(0.8, 0.8, 0.8); for (int i = 0; i < primitives.size(); i++) { if (auto pt = std::dynamic_pointer_cast(primitives[i])) { tensor3 originTensor(nullptr, 3); @@ -487,9 +529,10 @@ void quadratureScene(const std::vector>& primitiv algoim_spark_alloc(real, transformedTensor); makeSphere(*pt, originTensor); - detail::powerTransformation(range, xmin, originTensor, transformedTensor); + real testEvaOri = evalPower(originTensor, testX); + detail::power2BernsteinTensor(transformedTensor, phi); completeTensorReps.emplace_back(CompleteTensorRep{phi, {originTensor}}); } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { @@ -514,14 +557,54 @@ void quadratureScene(const std::vector>& primitiv real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); completeTensorReps.emplace_back(CompleteTensorRep{phi, planeTensors}); + } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { + std::vector rawTensors(3, tensor3(nullptr, 2)); + rawTensors[0].ext_ = 3; + rawTensors[0].ext_(pt->alignAxis) = 1; + algoimSparkAllocHeapVector(tensorStacks, rawTensors); + + uvector3i resExt = detail::getCompositeExt(rawTensors); + tensor3 phi(nullptr, resExt); + tensorStacks.emplace_back(algoim_spark_alloc_heap(real, phi)); + + tensor3 compositeTensor(nullptr, resExt); + algoim_spark_alloc(real, compositeTensor); + + makeCylinder(*pt, compositeTensor, rawTensors); + detail::powerTransformation(range, xmin, compositeTensor); + real testEvalPower = evalPower(compositeTensor, testX); + + detail::power2BernsteinTensor(compositeTensor, phi); + + completeTensorReps.emplace_back(CompleteTensorRep{phi, rawTensors}); } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { - tensor3 originTensor(nullptr, 3); - tensor3 transformedTensor(nullptr, 3); + std::vector rawTensors(3, tensor3(nullptr, 2)); + rawTensors[0].ext_ = 3; + algoimSparkAllocHeapVector(tensorStacks, rawTensors); + + uvector resExt = 2 + 1 + 1 + 1; + tensor3 phi(nullptr, resExt); + tensorStacks.emplace_back(algoim_spark_alloc_heap(real, phi)); + + tensor3 compositeTensor(nullptr, resExt); + algoim_spark_alloc(real, compositeTensor); + + makeCone(*pt, compositeTensor, rawTensors); + detail::powerTransformation(range, xmin, compositeTensor); + real testEvalPower = evalPower(compositeTensor, testX); + + detail::power2BernsteinTensor(compositeTensor, phi); + + completeTensorReps.emplace_back(CompleteTensorRep{phi, rawTensors}); + } else { + std::cerr << "unrecognized primitive type." << std::endl; } } - Scene scene{completeTensorReps, blobTree}; - buildOcTree(scene, rootNode, leaves); + Scene scene{completeTensorReps, blobTree}; + 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); diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index df32a49..aeb261a 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -112,19 +112,21 @@ void powerTransformation(const uvector& scale, const uvector& } } -static void compositePower(const std::vector>& powers, - int powerIdx, - uvector powerSum, - real factor, - xarray& res) +uvector3i getCompositeExt(const std::vector& tensors) +{ + uvector3i resExt = 1; + for (const auto& t : tensors) { resExt += t.ext() - 1; } + return resExt; +} +void compositePower(const std::vector>& powers, + int powerIdx, + uvector powerSum, + real factor, + xarray& res) { if (powerIdx == 0) { - { - uvector3 ext(1, 1, 1); - for (auto& t : powers) { ext += t.ext() - 1; } - assert(all(ext == res.ext())); - } + assert(all(res.ext() == getCompositeExt(powers))); xarrayInit(res); } if (powerIdx == powers.size()) { @@ -498,17 +500,19 @@ class CylinderDesc : virtual public PrimitiveDesc public: const static PrimitiveType type = Cylinder; uvector3 node1; - uvector3 node2; + // uvector3 node2; + real height; real radius; + int alignAxis; // CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} // 现在是只支持轴对齐x,y,z的圆柱 // TODO: 实现张量积的旋转后换成支持任意方向的 - CylinderDesc(const uvector3& n1, real r, real h, int alignAxis) : PrimitiveDesc(), node1(n1), radius(r) + CylinderDesc(const uvector3& n1, real r, real h, int ax) : PrimitiveDesc(), node1(n1), radius(r), height(h), alignAxis(ax) { assert(alignAxis >= 0 && alignAxis <= 2); - node2 = node1; - node2(alignAxis) += h; + // node2 = node1; + // node2(alignAxis) += h; } void print() override { std::cout << "Cylinder Description" << std::endl; } @@ -519,10 +523,15 @@ class ConeDesc : virtual public PrimitiveDesc public: const static PrimitiveType type = Cone; uvector3 node1; - uvector3 node2; + // uvector3 node2; real radius; + real height; + int alignAxis; - ConeDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} + ConeDesc(const uvector3& n1, real r, real h, int ax) : PrimitiveDesc(), node1(n1), radius(r), height(h), alignAxis(ax) + { + assert(alignAxis >= 0 && alignAxis <= 2); + } void print() override { std::cout << "Cone Description" << std::endl; } }; @@ -594,7 +603,36 @@ public: } }; -void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector>& planeTensors) +struct AABB { +public: + uvector3 min = uvector3(std::numeric_limits::max()); + uvector3 max = uvector3(std::numeric_limits::lowest()); + AABB() = default; + + AABB(const uvector3& min_, const uvector3& max_) : min(min_), max(max_) {} + + void extend(const uvector3& p) + { + for (int i = 0; i < 3; ++i) { + min(i) = std::min(min(i), p(i)); + max(i) = std::max(max(i), p(i)); + } + } + + void extend(const AABB& aabb) + { + for (int i = 0; i < 3; ++i) { + min(i) = std::min(min(i), aabb.min(i)); + max(i) = std::max(max(i), aabb.max(i)); + } + } + + uvector3 center() const { return (min + max) / 2; } + + uvector3 size() const { return max - min; } +}; + +void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector>& planeTensors, AABB& aabb) { uvector3 ext(1 + mesh.indexInclusiveScan.size()); @@ -625,10 +663,12 @@ void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector& tensor) +void makeSphere(const SphereDesc& sphereDesc, xarray& tensor, AABB& aabb) { uvector ext = 3; assert(all(ext == tensor.ext())); @@ -642,27 +682,84 @@ void makeSphere(const SphereDesc& sphereDesc, xarray& tensor) tensor.m(idx) = sphereDesc.amplitude(dim) * sphereDesc.amplitude(dim); } tensor.m(0) -= sphereDesc.radius * sphereDesc.radius; - // return PowerTensor(tensor); + // TODO consider the amplitude + aabb.min = sphereDesc.center - sphereDesc.radius; + aabb.max = sphereDesc.center + sphereDesc.radius; }; -void makeCylinder(const CylinderDesc& cylinder, xarray& tensor, std::vector& rawTensors) +void makeCylinder(const CylinderDesc& cylinderDesc, xarray& tensor, std::vector& rawTensors) { - assert(rawTensors.size() == 3 - && (all(rawTensors[0].ext() == uvector3(3, 3, 1)) || all(rawTensors[0].ext() == uvector3(3, 3, 1)) - || all(rawTensors[0].ext() == uvector3(3, 3, 1))) - && all(rawTensors[1].ext() == uvector3(2)) && all(rawTensors[1].ext() == uvector3(2))); - uvector ext = 5; - assert(all(ext == tensor.ext())); + auto & cylinderSrf = rawTensors[0], plane1 = rawTensors[1], plane2 = rawTensors[2]; + int alignAxis = cylinderDesc.alignAxis; + uvector cylinderExt = 3; + cylinderExt(alignAxis) = 1; + assert(rawTensors.size() == 3 && all(cylinderSrf.ext() == cylinderExt) && all(plane1.ext() == uvector3(2)) + && all(plane2.ext() == uvector3(2))); + assert(all(detail::getCompositeExt(rawTensors) == tensor.ext())); + + int dimA = (alignAxis + 1) % 3, dimB = (alignAxis + 2) % 3; + real a = cylinderDesc.node1(dimA), b = cylinderDesc.node1(dimB), c = cylinderDesc.node1(alignAxis), r = cylinderDesc.radius; + uvector idx = 0; + plane1.m(idx) = c; + plane2.m(idx) = -c - cylinderDesc.height; + cylinderSrf.m(idx) = a * a + b * b - r * r; + idx(alignAxis) = 1; + plane1.m(idx) = -1; + plane2.m(idx) = 1; + idx(alignAxis) = 0; + idx(dimA) = 1; + cylinderSrf.m(idx) = -2 * a; + idx(dimA) = 2; + cylinderSrf.m(idx) = 1; + + idx = 0; + idx(dimB) = 1; + cylinderSrf.m(idx) = -2 * b; + idx(dimB) = 2; + cylinderSrf.m(idx) = 1; + + detail::compositePower(rawTensors, 0, 0, 1, tensor); } -void makeCone(const CylinderDesc& cone, tensor3& tensor) +void makeCone(const ConeDesc& coneDesc, tensor3& tensor, std::vector& rawTensors) { - // TODO: + assert(rawTensors.size() == 3 && all(rawTensors[0].ext() == 3) && all(rawTensors[1].ext() == 2) + && all(rawTensors[2].ext() == 2)); + assert(all(tensor.ext() == detail::getCompositeExt(rawTensors))); + auto &coneSrf = rawTensors[0], &plane1 = rawTensors[1], &plane2 = rawTensors[2]; + int alignAxis = coneDesc.alignAxis; + + int dimA = (alignAxis + 1) % 3, dimB = (alignAxis + 2) % 3; + real a = coneDesc.node1(dimA), b = coneDesc.node1(dimB), c = coneDesc.node1(alignAxis); + real scale = coneDesc.height / coneDesc.radius, scaleSquare = scale * scale; + // f = (scale(x-a))^2 + (scale(y-b))^2 - (z-c)^2 + uvector idx = 0; + plane1.m(idx) = c; + plane2.m(idx) = -c - coneDesc.height; + coneSrf.m(idx) = scaleSquare * (a * a + b * b) - c * c; + idx(alignAxis) = 1; + plane1.m(idx) = -1; + plane2.m(idx) = 1; + coneSrf.m(idx) = 2 * c; + idx(alignAxis) = 2; + coneSrf.m(idx) = -1; + idx(alignAxis) = 0; + idx(dimA) = 1; + coneSrf.m(idx) = scaleSquare * (-2 * a); + idx(dimA) = 2; + coneSrf.m(idx) = scaleSquare; + idx = 0; + idx(dimB) = 1; + coneSrf.m(idx) = scaleSquare * (-2 * b); + idx(dimB) = 2; + coneSrf.m(idx) = scaleSquare; + + detail::compositePower(rawTensors, 0, 0, 1, tensor); } struct CompleteTensorRep { tensor3 compositedBernstein; // 合法输入在[0,1]^3 - std::vector originalPower; // 合法输入在场景 + std::vector rawPower; // 合法输入在场景 }; struct Scene { @@ -694,6 +791,8 @@ struct OcTreeNode { max(node.max), blobTree(node.blobTree) { + int a = 1; + int b = 2; } }; diff --git a/algoim/uvector.hpp b/algoim/uvector.hpp index d0e9024..dc6c832 100644 --- a/algoim/uvector.hpp +++ b/algoim/uvector.hpp @@ -71,8 +71,7 @@ using extent = std::integral_constant>::value>; // binary uvector expressions #define algoim_decl_binary_op(op) \ - template ::value > 0) || (detail::extent::value > 0), bool> = true> \ auto operator op(E1&& e1, E2&& e2) \ { \ @@ -162,8 +161,12 @@ public: #undef algoim_decl_assign_op }; -using uvector3 = uvector; -using uvector2 = uvector; +using uvector3f = uvector; +using uvector2f = uvector; +using uvector3i = uvector; +using uvector2i = uvector; +using uvector3 = uvector3f; +using uvector2 = uvector2f; // some methods, particularly those which remove components of a vector, benefit // from having a zero-length uvector; it is simply an empty class diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 5be5d4f..7e53327 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -124,21 +124,41 @@ void casePolyhedronSphere() void caseScene() { - std::vector> primitiveDescriptions(7); - primitiveDescriptions[0] = std::make_shared(CuboidDesc(0., 1.6)); - primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); - 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.)); - primitiveDescriptions[4] = std::make_shared(ConeDesc(uvector3(0., -0.2, 0.), uvector3(0., -0.9, 0.), 0.4)); + 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[3] = std::make_shared(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.)); std::vector pyramidBottomVertices = { uvector3{-1, -1, 0}, uvector3{1, -1, 0}, uvector3{1, 1, 0}, uvector3{-1, 1, 0} }; - primitiveDescriptions[5] = std::make_shared(pyramidBottomVertices, uvector3{0, 0, 1}); - primitiveDescriptions[4] = std::make_shared(CylinderDesc(uvector3(-0.3, 0.3, 2.3), 0.4, 3.6, 1)); - quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 2.3), organizer::BlobTree()); + primitiveDescriptions[4] = std::make_shared(pyramidBottomVertices, uvector3{0, 0, 1}); + primitiveDescriptions[5] = std::make_shared(ConeDesc(uvector3(0., -0.2, 0.), -0.7, 0.4, 1)); + // primitiveDescriptions[6] = std::make_shared(CylinderDesc(uvector3(-0.3, 0.3, 2.3), 0.4, 3.6, 1)); + organizer::BlobTree blobTree; + blobTree.structure.resize(PRIMITIVE_CNT * 2 - 1); + blobTree.primitiveNodeIdx.resize(PRIMITIVE_CNT); + blobTree.structure[0] = {1, 0, 0, 0, 1, 2}; // cube1 + blobTree.structure[1] = {1, 0, 0, 0, 0, 0}; // cube2 + blobTree.structure[2] = {0, 0, 0, 0, 1, 6}; // Union of cubes = opNode1 + blobTree.structure[3] = {1, 0, 0, 0, 1, 5}; // sphere1 + blobTree.structure[4] = {1, 0, 0, 0, 0, 0}; // sphere2 + blobTree.structure[5] = {0, 0, 0, 0, 0, 0}; // Union of spheres = opNode2 + blobTree.structure[6] = {0, 2, 0, 0, 1, 10}; // Difference of opNode1 and sphere2 = opNode3 + blobTree.structure[7] = {1, 0, 0, 0, 1, 9}; // Pyramid + blobTree.structure[8] = {1, 0, 0, 0, 0, 0}; // Cone + blobTree.structure[9] = {0, 1, 0, 0, 0, 0}; // Intersection of Pyramid and Cone = opNode4 + blobTree.structure[10] = {0, 0, 0, 0, 1, 12}; // Union of opNode3 and opNode4 = opNode5 + // blobTree.structure[11] = {1, 0, 0, 0, 0, 0}; // Cylinder + // blobTree.structure[12] = {0, 2, 0, 0, 1, 0}; // Difference of opNode5 and Cylinder + // blobTree.primitiveNodeIdx = {0, 1, 3, 4, 7, 8, 11}; + // quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 2.3), blobTree); + blobTree.primitiveNodeIdx = {0, 1, 3, 4, 7, 8}; + quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree); } void testDeCasteljau() @@ -168,6 +188,43 @@ void testDeCasteljau() std::cout << "sign = " << sign << std::endl; } +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); + 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); + xmin = 0, xmax = 1; + + // 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); + std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(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); + + std::cout << "evalX original rep = " << 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; +} + void testBlob() { organizer::Blob blob = {5, 3, 4, 5, 6}; @@ -181,10 +238,24 @@ void testBlob() std::cout << sizeof(int) << std::endl; } +void vectorDebug() +{ + struct A { + std::vector data; + A() = default; + }; + + A a; + int s = a.data.size(); + std::cout << "size: " << s << std::endl; +} + void testPrimitive() { // casePolyhedron1(); // casePolyhedronSphere(); - // testDeCasteljau(); - testBlob(); + // testSubDivideWithDeCasteljau(); + // testBlob(); + caseScene(); + // vectorDebug(); } \ No newline at end of file