diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index f70912c..09f6656 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -108,14 +108,19 @@ AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector side, int } } // namespace detail -bool keepQuadraturePoint(const Scene& scene, const OcTreeNode& ocTreeNode, const uvector3& originPt) +bool keepQuadraturePoint(const std::vector& tensors, const OcTreeNode& ocTreeNode, const uvector3& point) { // 只需要考虑intersect polys，不用考虑fully contained polys - const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; + const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; + assert(tensors.size() == polyIntersectIndices.size()); + if (polyIntersectIndices.size() == 0) { + assert(ocTreeNode.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN); + return ocTreeNode.blobTree.structure.back().inOut; + } std::vector primitiveInOuts(polyIntersectIndices.size()); for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) { // primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt); - primitiveInOuts[i] = isInsideBernstein(scene.tensors[polyIntersectIndices[i]], originPt); + primitiveInOuts[i] = isInsideBernstein(tensors[i], point); } // 这里blobTree是拷贝传参 auto blobTree = ocTreeNode.blobTree; @@ -188,7 +193,7 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector mark(0, 0, 0); for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { real centerPlane = nodeMid(faceAxis); @@ -296,7 +301,7 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector& polyIntersectIndices = node.polyIntersectIndices; - if (polyIntersectIndices.size() <= 5) { + if (polyIntersectIndices.size() <= 3) { leaves.emplace_back(node); return; } @@ -326,21 +331,21 @@ void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector 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; - // } + if (!minimalRep.aabb.intersect(subNodes[subIdx].aabb)) { + // out of the subcell + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false); + continue; + } if (subIdx == 1) { int aaa = 1; int bbb = 1; } - tensor3 subcellPoly(nullptr, poly.ext()); + tensor3 subcellPoly(nullptr, minimalRep.tensor.ext()); algoim_spark_alloc(real, subcellPoly); - bernstein::deCasteljau(poly, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly); + bernstein::deCasteljau(minimalRep.tensor, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly); int sign = bernstein::uniformSign(subcellPoly); if (sign == 1) { organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false); @@ -458,18 +463,36 @@ void basicTask(const std::vector>& primitives, in BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 10) { - auto integrand = [](const uvector& x) { return 1.0; }; - real volume = 0., surf = 0.; - auto range = node.aabb.size(); - ImplicitPolyQuadrature<3> ipquad(scene.tensors, node.polyIntersectIndices); + if (node.polyIntersectIndices.size() == 0) { + assert(node.blobTree.structure.back().inOut != NODE_IN_OUT_UNKNOWN); + if (node.blobTree.structure.back().inOut == NODE_IN) { + return {node.aabb.volume(), 0}; + } else { + return {0, 0}; + } + } + + auto integrand = [](const uvector& x) { return 1.0; }; + real volume = 0., surf = 0.; + auto range = node.aabb.size(); + std::vector phis; // phis using subcell as [0,1]^3 + std::vector*> phiStacks; + for (int i = 0; i < node.polyIntersectIndices.size(); i++) { + const auto& polyWithinScene = scene.minimalReps[node.polyIntersectIndices[i]].tensor; + phis.emplace_back(tensor3(nullptr, polyWithinScene.ext())); + phiStacks.emplace_back(algoim_spark_alloc_heap(real, phis.back())); + bernstein::deCasteljau(polyWithinScene, node.aabb.min, node.aabb.max, phis.back()); + } + // ImplicitPolyQuadrature<3> ipquad(scene.minimalReps, node.polyIntersectIndices); + ImplicitPolyQuadrature<3> ipquad(phis); + ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { - auto realX = x * range + node.aabb.min; - if (keepQuadraturePoint(scene, node, realX)) { volume += w * integrand(realX); } + auto realX = x * range + node.aabb.min; // 这里realX应该是最原始空间下的点？不过因为算体积，所以不影响 + if (keepQuadraturePoint(phis, node, x)) { volume += w * integrand(realX); } }); - ipquad.integrate_surf(AutoMixed, q, [&](const uvector& x, real w, const uvector& wn) { - volume += w * integrand(x * range + node.aabb.min); - }); - return {volume, surf}; + + for (auto& p : phiStacks) delete p; + return {volume * node.aabb.volume(), surf}; } void quadratureScene(const std::vector>& primitives, @@ -513,7 +536,7 @@ void quadratureScene(const std::vector>& primitiv } visiblePrimitiveReps[i].aabb.normalize(range, xmin); } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { - visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 2), tensor3(nullptr, 2)}; + visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)}; auto& tensors = visiblePrimitiveReps[i].tensors; tensors[0].ext_(pt->alignAxis) = 1; algoimSparkAllocHeapVector(tensorStacks, tensors); @@ -526,7 +549,7 @@ void quadratureScene(const std::vector>& primitiv } visiblePrimitiveReps[i].aabb.normalize(range, xmin); } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { - visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 2), tensor3(nullptr, 2)}; + visiblePrimitiveReps[i].tensors = {tensor3(nullptr, 3), tensor3(nullptr, 3), tensor3(nullptr, 3)}; auto& tensors = visiblePrimitiveReps[i].tensors; algoimSparkAllocHeapVector(tensorStacks, tensors); @@ -542,9 +565,9 @@ void quadratureScene(const std::vector>& primitiv } } - std::vector realPrimitives; + std::vector minimalReps; /*** merge subtrees to main tree ***/ - std::vector realLeafIndices; + std::vector realLeafIndices; for (int i = 0; i < visiblePrimitiveReps.size() - 1; ++i) { // 必须以自左向右的顺序访问所有叶节点 assert(blobTree.primitiveNodeIdx[i] < blobTree.primitiveNodeIdx[i + 1]); @@ -574,26 +597,33 @@ void quadratureScene(const std::vector>& primitiv for (auto primitiveIdx : visiblePrimitiveReps[i].subBlobTree.primitiveNodeIdx) { realLeafIndices.push_back(primitiveIdx + originLeafIdx); } - - realPrimitives.insert(realPrimitives.end(), visiblePrimitiveReps[i].tensors.begin(), - visiblePrimitiveReps[i].tensors.end()); + minimalReps.reserve(minimalReps.size() + visiblePrimitiveReps[i].tensors.size()); + const auto& aabb = visiblePrimitiveReps[i].aabb; + for (const auto& tensor : visiblePrimitiveReps[i].tensors) { + minimalReps.emplace_back(MinimalPrimitiveRep{tensor, aabb}); + } } else { blobTree.structure[originLeafIdx].isPrimitive = true; realLeafIndices.push_back(originLeafIdx); - realPrimitives.emplace_back(visiblePrimitiveReps[i].tensors[0]); + minimalReps.emplace_back(MinimalPrimitiveRep{visiblePrimitiveReps[i].tensors[0], visiblePrimitiveReps[i].aabb}); } } blobTree.primitiveNodeIdx = realLeafIndices; /*** merge subtrees to main tree ***/ - Scene scene{realPrimitives, AABB(xmin, xmax)}; + Scene scene{minimalReps, AABB(xmin, xmax)}; OcTreeNode rootNode(0, 1, blobTree); - for (int i = 0; i < realPrimitives.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } + for (int i = 0; i < minimalReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } int cnt = 1; buildOcTreeV0(scene, rootNode, leaves, 1, cnt); + basicTask(scene, leaves[14], q); + + int i = 0; for (const auto& leaf : leaves) { - auto basicRes = basicTask(scene, leaf, q); - volume += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min); + auto basicRes = basicTask(scene, leaf, q); + if (std::isinf(basicRes.volume)) { std::cout << "inf volume when solving leaf: " << i << std::endl; } + volume += basicRes.volume; + std::cout << "Solved leaves: " << ++i << "/" << leaves.size() << std::endl; } volume *= prod(xmax - xmin); diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 8fb0f5d..e03a87c 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -635,9 +635,9 @@ public: uvector3 halfSize = size / 2; for (int i = 0; i < 8; ++i) { vertices[i] = center; - vertices[i](0) += (i & 1) ? halfSize(0) : -halfSize(0); + vertices[i](0) += (i & 4) ? halfSize(0) : -halfSize(0); vertices[i](1) += (i & 2) ? halfSize(1) : -halfSize(1); - vertices[i](2) += (i & 4) ? halfSize(2) : -halfSize(2); + vertices[i](2) += (i & 1) ? halfSize(2) : -halfSize(2); } indices = {3, 2, 0, 1, // left 4, 6, 7, 5, // right @@ -678,6 +678,8 @@ public: uvector3 size() const { return max - min; } + real volume() const { return prod(size()); } + bool intersect(const AABB& aabb) const { for (int i = 0; i < 3; ++i) { @@ -708,7 +710,7 @@ struct Scene { // std::vector polys; // std::vector visiblePrimitives; // organizer::BlobTree blobTree; - std::vector tensors; + std::vector minimalReps; AABB boundary; }; diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index 0979334..2c38d79 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -18,6 +18,7 @@ // a good choice, and this is specified by the following macro def. #include #include +#include #include #include "organizer/primitive.hpp" #define ALGOIM_M 8 @@ -614,6 +615,15 @@ struct ImplicitPolyQuadrature { build(true, false); } + ImplicitPolyQuadrature(const std::vector& minimalReps, const std::vector& polyIndices) + { + for (auto i : polyIndices) { + auto mask = detail::nonzeroMask(minimalReps[i].tensor, booluarray(true)); + if (!detail::maskEmpty(mask)) phi.push_back(minimalReps[i].tensor, 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, @@ -654,12 +664,13 @@ struct ImplicitPolyQuadrature { // Compute score; penalise any directions which likely contain vertical tangents uvector has_disc; uvector score = detail::score_estimate(phi, has_disc); - // if (max(abs(score)) == 0) - // score(0) = 1.0; + /**gjj */ + if (max(abs(score)) == 0) score(0) = 1.0; + /**gjj */ assert(max(abs(score)) > 0); score /= 2 * max(abs(score)); for (int i = 0; i < N; ++i) - if (!has_disc(i)) score(i) += 1.0; // 1个与k相切的多项式都没有，才+1 + if (!has_disc(i)) score(i) += 1.0; // 一个与k相切的多项式都没有时，才+1 // Choose height direction and form base polynomials; if tanh-sinh is being used at this // level, suggest the same all the way down; moreover, suggest tanh-sinh if a non-empty @@ -748,9 +759,20 @@ struct ImplicitPolyQuadrature { if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j]; } }; - + /**gjj */ + auto newEnd = std::remove_if(nodes, nodes + count, [](real x) { return std::isnan(x); }); + count = std::distance(nodes, newEnd); + // for (int i = 0; i < count; ++i) { + // if (std::isnan(nodes[i])) { nanCount++; } + // } + // count -= nanCount; + /**gjj */ // Sort the nodes std::sort(nodes, nodes + count); + if (!(nodes[0] == real(0) && nodes[count - 1] == real(1))) { + int aaa = 1; + int bbb = 1; + } assert(nodes[0] == real(0) && nodes[count - 1] == real(1)); // Force nearly-degenerate sub-intervals to be exactly degenerate diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 2a3723b..fb1e213 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -151,8 +151,8 @@ void caseScene() uvector3{-1, 1, 0} }; // primitiveDescriptions[4] = std::make_shared(pyramidBottomVertices, uvector3{0, 0, 1}); - primitiveDescriptions[4] = std::make_shared(SphereDesc(0.2, uvector3(0., 0., -1), 1.)); - primitiveDescriptions[5] = std::make_shared(ConeDesc(uvector3(0., -0.2, 0.), -0.7, 0.4, 1)); + primitiveDescriptions[4] = std::make_shared(SphereDesc(0.2, uvector3(0.2, -0.7, 0.), 1.)); + primitiveDescriptions[5] = std::make_shared(ConeDesc(uvector3(0., -0.2, 0.), 0.4, -0.7, 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); @@ -164,16 +164,16 @@ void caseScene() 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 - 6}; // Difference of opNode1 and opNode2 = opNode3 - blobTree.structure[7] = {1, 0, 0, 0, 1, 9 - 7}; // Pyramid + blobTree.structure[7] = {1, 0, 0, 0, 1, 9 - 7}; // Pyramid (sphere3) 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 - 10}; // Union of opNode3 and opNode4 = opNode5 + blobTree.structure[9] = {0, 0, 0, 0, 0, 0}; // UNION of Pyramid (sphere3) and Cone = opNode4 + blobTree.structure[10] = {0, 2, 0, 0, 1, 12 - 10}; // Difference 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); + quadratureScene(primitiveDescriptions, uvector3(-0.8, -1.3, -1.6), uvector3(1.6, 1.6, 1.5), blobTree); } void testDeCasteljau() @@ -279,6 +279,30 @@ void testPlaneUniformSign() std::cout << "sign = " << sign << std::endl; } +void testPlaneWithinSubcell() +{ + AABB sceneAABB{uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 1.5)}; + AABB subCell{ + uvector3{0.076923076923076983, 0.17241379310344829, 0.77419354838709675}, + uvector3{0.076923076923077094, 0.1724137931034484, 0.77419354838709686} + }; + tensor3 planeTensor(nullptr, 3), subPlaneTensor(nullptr, 3); + algoim_spark_alloc(real, planeTensor, subPlaneTensor); + xarrayInit(planeTensor); + planeTensor.m(0) = -0.8; + planeTensor.m(uvector3{0, 0, 1}) = 1; + organizer::detail::powerTransformation(sceneAABB.size(), sceneAABB.min, planeTensor); + organizer::detail::power2BernsteinTensor(planeTensor); + bernstein::deCasteljau(planeTensor, subCell.min, subCell.max, subPlaneTensor); + int sign = bernstein::uniformSign(subPlaneTensor); + std::cout << sign << std::endl; + ImplicitPolyQuadrature<3> ipquad(subPlaneTensor); + real volume = 0; + ipquad.integrate(AutoMixed, 10, [&](const uvector &x, real w) { + if (evalBernstein(subPlaneTensor, x)) volume += w * 1; + }); +} + void testBlob() { organizer::Blob blob = {5, 3, 4, 5, 6}; @@ -320,6 +344,7 @@ void testPrimitive() // testSubDivideWithDeCasteljau(); // testBlob(); caseScene(); + // testPlaneWithinSubcell(); // caseCubeSphere(); // caseTwoCube(); // testTensorInverse();