From 5e8545d77083de3b8ac5e8c21a70fedd33a933fb Mon Sep 17 00:00:00 2001 From: gjj Date: Wed, 11 Sep 2024 19:48:16 +0800 Subject: [PATCH] blobtree octree integration --- algoim/kdtree.hpp | 437 +++++++++++++--------------- algoim/organizer/blobtree.hpp | 79 ++--- algoim/organizer/organizer.hpp | 493 ++++++++++++++++++++------------ algoim/organizer/primitive.hpp | 24 +- algoim/quadrature_multipoly.hpp | 12 +- 5 files changed, 579 insertions(+), 466 deletions(-) diff --git a/algoim/kdtree.hpp b/algoim/kdtree.hpp index ad8facd..741978e 100644 --- a/algoim/kdtree.hpp +++ b/algoim/kdtree.hpp @@ -2,11 +2,11 @@ #define ALGOIM_KDTREE_HPP /* algoim::KDTree constructs a k-d tree data structure for a given collection of - points of type uvector, where T is typically float or double, and N is - the dimension. This particular implementation of the k-d tree has been optimised for the + points of type uvector, where T is typically float or double, and N is + the dimension. This particular implementation of the k-d tree has been optimised for the case that the points are situated on a smooth codimension-one surface by using coordinate transformations that result in "tight" bounding boxes for some of the nodes in the tree. - + Refer to the paper R. I. Saye, High-order methods for computing distances to implicitly defined surfaces, Communications in Applied Mathematics and Computational Science, 9(1), 107-141 (2014), @@ -21,272 +21,243 @@ namespace algoim { - namespace detail - { - // Returns the squared-distance from a point x to the hyperrectangle [min,max] (if x is inside the rectangle, the distance is zero) - template - T sqrDistBox(const uvector& x, const uvector& min, const uvector& max) - { - T dsqr = T(0); - for (int i = 0; i < N; ++i) - if (x(i) < min(i)) - dsqr += util::sqr(x(i) - min(i)); - else if (x(i) > max(i)) - dsqr += util::sqr(x(i) - max(i)); - return dsqr; - } - } +namespace detail +{ +// Returns the squared-distance from a point x to the hyperrectangle [min,max] (if x is inside the rectangle, the distance is +// zero) +template +T sqrDistBox(const uvector& x, const uvector& min, const uvector& max) +{ + T dsqr = T(0); + for (int i = 0; i < N; ++i) + if (x(i) < min(i)) + dsqr += util::sqr(x(i) - min(i)); + else if (x(i) > max(i)) + dsqr += util::sqr(x(i) - max(i)); + return dsqr; +} +} // namespace detail - template - class KDTree - { - std::vector> points; - std::vector index; - static constexpr int leaf_size = 16; +template +class KDTree +{ + std::vector> points; + std::vector index; + static constexpr int leaf_size = 16; - // A Node of the tree is a leaf node iff type=-1. If type >= 0, then a coordinate transform - // is applied. If type < -1, the node is a standard splitting node with two children and no transform - struct Node - { - int type; - int i0, i1; - uvector xmin, xmax; - }; + // A Node of the tree is a leaf node iff type=-1. If type >= 0, then a coordinate transform + // is applied. If type < -1, the node is a standard splitting node with two children and no transform + struct Node { + int type; + int i0, i1; + uvector xmin, xmax; + }; - std::vector nodes; - std::vector,N>> transforms; + std::vector nodes; + std::vector, N>> transforms; - // Given a node and a range of points [lb,ub), build the tree - void build_tree(int nodeIndex, int lb, int ub, bool hasTransformed, int level) - { - assert(lb < ub); - Node& node = nodes[nodeIndex]; + // Given a node and a range of points [lb,ub), build the tree + void build_tree(int nodeIndex, int lb, int ub, bool hasTransformed, int level) + { + assert(lb < ub); + Node& node = nodes[nodeIndex]; - // Compute bounding box and mean - uvector mean = node.xmin = node.xmax = points[index[lb]]; - for (int i = lb + 1; i < ub; ++i) - { - const uvector& x = points[index[i]]; - mean += x; - for (int j = 0; j < N; ++j) - { - if (x(j) < node.xmin(j)) node.xmin(j) = x(j); - if (x(j) > node.xmax(j)) node.xmax(j) = x(j); - } + // Compute bounding box and mean + uvector mean = node.xmin = node.xmax = points[index[lb]]; + for (int i = lb + 1; i < ub; ++i) { + const uvector& x = points[index[i]]; + mean += x; + for (int j = 0; j < N; ++j) { + if (x(j) < node.xmin(j)) node.xmin(j) = x(j); + if (x(j) > node.xmax(j)) node.xmax(j) = x(j); } - mean /= static_cast(ub - lb); + } + mean /= static_cast(ub - lb); - // The node is a leaf iff point count is sufficiently small - if (ub - lb <= leaf_size) - { - node.type = -1; - node.i0 = lb; - node.i1 = ub; - return; - } + // The node is a leaf iff point count is sufficiently small + if (ub - lb <= leaf_size) { + node.type = -1; + node.i0 = lb; + node.i1 = ub; + return; + } - // Splitting node: default to splitting along greatest extent - node.type = -2; - int axis = argmax(node.xmax - node.xmin); + // Splitting node: default to splitting along greatest extent + node.type = -2; + int axis = argmax(node.xmax - node.xmin); - // Evaluate possibility for coordinate transformation - if (!hasTransformed && level > 5 && ub - lb >= leaf_size * (1 << 2)) - { - // Estimate normal - T holeRadiusSqr = util::sqr(0.05*max(node.xmax - node.xmin)); - uvector n = static_cast(0); + // Evaluate possibility for coordinate transformation + if (!hasTransformed && level > 5 && ub - lb >= leaf_size * (1 << 2)) { + // Estimate normal + T holeRadiusSqr = util::sqr(0.05 * max(node.xmax - node.xmin)); + uvector n = static_cast(0); + n(0) = 1.0; + for (int i = lb; i < ub; ++i) { + uvector x = points[index[i]] - mean; + T msqr = sqrnorm(x); + if (msqr > holeRadiusSqr) n -= x * (dot(x, n) / msqr); + } + T msqr = sqrnorm(n); + if (msqr == 0.0) n(0) = 1.0; - for (int i = lb; i < ub; ++i) - { - uvector x = points[index[i]] - mean; - T msqr = sqrnorm(x); - if (msqr > holeRadiusSqr) - n -= x * (dot(x,n)/msqr); - } - T msqr = sqrnorm(n); - if (msqr == 0.0) - n(0) = 1.0; - else - n /= sqrt(msqr); + else + n /= sqrt(msqr); - // Compute alpha range - T minAlpha = std::numeric_limits::max(); - T maxAlpha = -std::numeric_limits::max(); - for (int i = lb; i < ub; ++i) - { - T alpha = dot(points[index[i]], n); - if (alpha > maxAlpha) maxAlpha = alpha; - if (alpha < minAlpha) minAlpha = alpha; - } + // Compute alpha range + T minAlpha = std::numeric_limits::max(); + T maxAlpha = -std::numeric_limits::max(); + for (int i = lb; i < ub; ++i) { + T alpha = dot(points[index[i]], n); + if (alpha > maxAlpha) maxAlpha = alpha; + if (alpha < minAlpha) minAlpha = alpha; + } - if (maxAlpha - minAlpha < 0.1*max(node.xmax - node.xmin)) - { - // Perform transformation: calculate an orthonormal basis using the normal as first axis. - // A stable method for doing so is to compute the Householder matrix which maps n to ej, - // i.e., P = I - 2 uu^T where u = normalised(n - ej), where ej is the j-th basis vector, - // j chosen that that n != ej. - uvector,N> axes; - int j = argmin(abs(n)); - uvector u = n; u(j) -= 1.0; - u /= norm(u); - for (int dim = 0; dim < N; ++dim) - for (int i = 0; i < N; ++i) - axes(dim)(i) = (dim == i ? 1.0 : 0.0) - 2.0 * u(dim) * u(i); + if (maxAlpha - minAlpha < 0.1 * max(node.xmax - node.xmin)) { + // Perform transformation: calculate an orthonormal basis using the normal as first axis. + // A stable method for doing so is to compute the Householder matrix which maps n to ej, + // i.e., P = I - 2 uu^T where u = normalised(n - ej), where ej is the j-th basis vector, + // j chosen that that n != ej. + uvector, N> axes; + int j = argmin(abs(n)); + uvector u = n; + u(j) -= 1.0; + u /= norm(u); + for (int dim = 0; dim < N; ++dim) + for (int i = 0; i < N; ++i) axes(dim)(i) = (dim == i ? 1.0 : 0.0) - 2.0 * u(dim) * u(i); - // Swap the first row of axes with j, so that the normal is the first axis. This is likely - // unnecessary (but done in the name of consistency with old approach). - std::swap(axes(0), axes(j)); + // Swap the first row of axes with j, so that the normal is the first axis. This is likely + // unnecessary (but done in the name of consistency with old approach). + std::swap(axes(0), axes(j)); - // Apply coordinate transformation and calculate new bounding box in order to determine new split direction - uvector bmin = std::numeric_limits::max(); - uvector bmax = -std::numeric_limits::max(); - for (int i = lb; i < ub; ++i) - { - uvector x = points[index[i]]; - for (int dim = 0; dim < N; ++dim) - { - T alpha = dot(axes(dim), x); - points[index[i]](dim) = alpha; - if (alpha < bmin(dim)) bmin(dim) = alpha; - if (alpha > bmax(dim)) bmax(dim) = alpha; - } + // Apply coordinate transformation and calculate new bounding box in order to determine new split direction + uvector bmin = std::numeric_limits::max(); + uvector bmax = -std::numeric_limits::max(); + for (int i = lb; i < ub; ++i) { + uvector x = points[index[i]]; + for (int dim = 0; dim < N; ++dim) { + T alpha = dot(axes(dim), x); + points[index[i]](dim) = alpha; + if (alpha < bmin(dim)) bmin(dim) = alpha; + if (alpha > bmax(dim)) bmax(dim) = alpha; } - - node.type = static_cast(transforms.size()); - transforms.push_back(axes); - axis = argmax(bmax - bmin); - hasTransformed = true; } + + node.type = static_cast(transforms.size()); + transforms.push_back(axes); + axis = argmax(bmax - bmin); + hasTransformed = true; } + } - // Use median as the split - int m = (lb + ub)/2; + // Use median as the split + int m = (lb + ub) / 2; - // Rearrange points - std::nth_element(index.begin() + lb, index.begin() + m, index.begin() + ub, [&](int i0, int i1) { return points[i0](axis) < points[i1](axis); } ); + // Rearrange points + std::nth_element(index.begin() + lb, index.begin() + m, index.begin() + ub, [&](int i0, int i1) { + return points[i0](axis) < points[i1](axis); + }); - // Build child trees - int i0 = node.i0 = static_cast(nodes.size()); - int i1 = node.i1 = static_cast(nodes.size() + 1); - nodes.push_back(Node()); - nodes.push_back(Node()); - build_tree(i0, lb, m, hasTransformed, level + 1); - build_tree(i1, m, ub, hasTransformed, level + 1); - } + // Build child trees + int i0 = node.i0 = static_cast(nodes.size()); + int i1 = node.i1 = static_cast(nodes.size() + 1); + nodes.push_back(Node()); + nodes.push_back(Node()); + build_tree(i0, lb, m, hasTransformed, level + 1); + build_tree(i1, m, ub, hasTransformed, level + 1); + } - struct ClosestPoint - { - uvector x; - T distsqr; - int ind; - }; + struct ClosestPoint { + uvector x; + T distsqr; + int ind; + }; - // Recursive function for searching the tree for the closest point to a given point cp.x - void search(const Node& node, ClosestPoint& cp) const - { - if (node.type == -1) - { - // Leaf node - for (int j = node.i0; j < node.i1; ++j) - { - T dsqr = sqrnorm(points[j] - cp.x); - if (dsqr < cp.distsqr) - { - cp.distsqr = dsqr; - cp.ind = j; - } + // Recursive function for searching the tree for the closest point to a given point cp.x + void search(const Node& node, ClosestPoint& cp) const + { + if (node.type == -1) { + // Leaf node + for (int j = node.i0; j < node.i1; ++j) { + T dsqr = sqrnorm(points[j] - cp.x); + if (dsqr < cp.distsqr) { + cp.distsqr = dsqr; + cp.ind = j; } } - else - { - // Non-leaf node - if (node.type >= 0) - { - // Transform query point to new coordinate system - const uvector,N>& axes = transforms[node.type]; - uvector x = cp.x; - for (int dim = 0; dim < N; ++dim) - cp.x(dim) = dot(axes(dim), x); - } + } else { + // Non-leaf node + if (node.type >= 0) { + // Transform query point to new coordinate system + const uvector, N>& axes = transforms[node.type]; + uvector x = cp.x; + for (int dim = 0; dim < N; ++dim) cp.x(dim) = dot(axes(dim), x); + } - T dleft = detail::sqrDistBox(cp.x, nodes[node.i0].xmin, nodes[node.i0].xmax); - T dright = detail::sqrDistBox(cp.x, nodes[node.i1].xmin, nodes[node.i1].xmax); - if (dleft < dright) - { - if (dleft < cp.distsqr) - { - search(nodes[node.i0], cp); - if (dright < cp.distsqr) - search(nodes[node.i1], cp); - } + T dleft = detail::sqrDistBox(cp.x, nodes[node.i0].xmin, nodes[node.i0].xmax); + T dright = detail::sqrDistBox(cp.x, nodes[node.i1].xmin, nodes[node.i1].xmax); + if (dleft < dright) { + if (dleft < cp.distsqr) { + search(nodes[node.i0], cp); + if (dright < cp.distsqr) search(nodes[node.i1], cp); } - else - { - if (dright < cp.distsqr) - { - search(nodes[node.i1], cp); - if (dleft < cp.distsqr) - search(nodes[node.i0], cp); - } + } else { + if (dright < cp.distsqr) { + search(nodes[node.i1], cp); + if (dleft < cp.distsqr) search(nodes[node.i0], cp); } + } - if (node.type >= 0) - { - // Transform query point back to old coordinate system. This is about 5% faster than storing - // the old value of cp.x on the stack - const uvector,N>& axes = transforms[node.type]; - uvector x = cp.x; - cp.x = axes(0)*x(0); - for (int dim = 1; dim < N; ++dim) - cp.x += axes(dim)*x(dim); - } + if (node.type >= 0) { + // Transform query point back to old coordinate system. This is about 5% faster than storing + // the old value of cp.x on the stack + const uvector, N>& axes = transforms[node.type]; + uvector x = cp.x; + cp.x = axes(0) * x(0); + for (int dim = 1; dim < N; ++dim) cp.x += axes(dim) * x(dim); } } + } - public: - - // Construct a KDTree from a given collection of points p (the given points are not overwritten) - KDTree(const std::vector>& p) - { - assert(p.size() < std::numeric_limits::max()); // code currently uses int to index but could easily be adapted to size_t - int len = static_cast(p.size()); - if (len == 0) - return; +public: + // Construct a KDTree from a given collection of points p (the given points are not overwritten) + KDTree(const std::vector>& p) + { + assert(p.size() + < std::numeric_limits::max()); // code currently uses int to index but could easily be adapted to size_t + int len = static_cast(p.size()); + if (len == 0) return; - // Copy points - points = p; + // Copy points + points = p; - // Build initial index array, which shall soon be reordered - index.resize(len); - for (int i = 0; i < len; ++i) - index[i] = i; + // Build initial index array, which shall soon be reordered + index.resize(len); + for (int i = 0; i < len; ++i) index[i] = i; - // Recursively build tree starting from root node. This only manipulates index array - nodes.push_back(Node()); - build_tree(0, 0, len, false, 0); + // Recursively build tree starting from root node. This only manipulates index array + nodes.push_back(Node()); + build_tree(0, 0, len, false, 0); - // Rearrange so that points in leaf nodes are contiguous in memory. Based on tree - // construction, both nodes & point ranges will be organised depth-first. - std::vector> pointscopy(points); - for (size_t i = 0; i < len; ++i) - points[i] = pointscopy[index[i]]; - } + // Rearrange so that points in leaf nodes are contiguous in memory. Based on tree + // construction, both nodes & point ranges will be organised depth-first. + std::vector> pointscopy(points); + for (size_t i = 0; i < len; ++i) points[i] = pointscopy[index[i]]; + } - // Search the tree for the closest point to x, where the search is restricted to all points that are - // within a squared distance of rsqr to x. Returns -1 if no such point exists, otherwise the index - // of the closest point in the original array p passed to the KDTree constructor is returned. - int nearest(const uvector& x, T rsqr = std::numeric_limits::max()) const - { - if (nodes.empty()) - return -1; - ClosestPoint cp; - cp.x = x; - cp.distsqr = rsqr; - cp.ind = -1; - search(nodes[0], cp); - return cp.ind >= 0 ? index[cp.ind] : -1; - } - }; + // Search the tree for the closest point to x, where the search is restricted to all points that are + // within a squared distance of rsqr to x. Returns -1 if no such point exists, otherwise the index + // of the closest point in the original array p passed to the KDTree constructor is returned. + int nearest(const uvector& x, T rsqr = std::numeric_limits::max()) const + { + if (nodes.empty()) return -1; + ClosestPoint cp; + cp.x = x; + cp.distsqr = rsqr; + cp.ind = -1; + search(nodes[0], cp); + return cp.ind >= 0 ? index[cp.ind] : -1; + } +}; } // namespace algoim #endif diff --git a/algoim/organizer/blobtree.hpp b/algoim/organizer/blobtree.hpp index 822fbf5..33d56cd 100644 --- a/algoim/organizer/blobtree.hpp +++ b/algoim/organizer/blobtree.hpp @@ -5,6 +5,14 @@ namespace algoim::organizer { +const int NODE_IN_OUT_UNKNOWN = 0; +const int NODE_IN = 1; +const int NODE_OUT = 2; + +const int OP_UNION = 0; +const int OP_INTERSECTION = 1; +const int OP_DIFFERENCE = 2; + // bitfield struct Blob { unsigned int isPrimitive : 1; // 1 for leaf @@ -37,67 +45,66 @@ 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); - node.inOut = in ? 1 : 2; - for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx; ++nodeIdx) { + node.inOut = in ? NODE_IN : NODE_OUT; + for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) { const auto& nowNode = tree.structure[nowIdx]; unsigned int nextIdx = isLeft(nowNode) ? nowNode.ancestor : nodeIdx + 1; - auto& nextNode = tree.structure[nextIdx]; + auto& nextNode = tree.structure[nextIdx]; // father if (nextNode.inOut != 0) { return; } - if (nowNode.inOut == 2) { - // out - if (nextNode.nodeOp == 1) { - // intersection - } else if (nextNode.nodeOp == 2 && isLeft(nowNode)) { - // difference and left - nextNode.inOut = 2; - } else if (nextNode.oneChildInOut == 1) { + if (nowNode.inOut == NODE_OUT) { + if (nextNode.nodeOp == OP_INTERSECTION) { + nextNode.inOut = NODE_OUT; + } else if (nextNode.nodeOp == OP_DIFFERENCE && isLeft(nowNode)) { + nextNode.inOut = NODE_OUT; + } else if (nextNode.oneChildInOut == NODE_IN) { // 两种情况，Union: in, Difference: NowNode一定是right，in-out，还是in - nextNode.inOut = 1; - } else if (nextNode.oneChildInOut == 2) { + nextNode.inOut = NODE_IN; + } else if (nextNode.oneChildInOut == NODE_OUT) { // 两种情况，Union: out, Difference: out-out, 还是out - nextNode.inOut = 2; + nextNode.inOut = NODE_OUT; } else { - nextNode.oneChildInOut = 2; + nextNode.oneChildInOut = NODE_OUT; return; } - } else if (nowNode.inOut == 1) { - // in - if (nextNode.nodeOp == 0) { - // union - nextNode.inOut = 1; - } else if (nextNode.nodeOp == 2 && !isLeft(nowNode)) { + } else if (nowNode.inOut == NODE_IN) { + if (nextNode.nodeOp == OP_UNION) { + nextNode.inOut = NODE_IN; + } else if (nextNode.nodeOp == OP_DIFFERENCE && !isLeft(nowNode)) { // difference and right - nextNode.inOut = 2; - } else if (nextNode.oneChildInOut == 1) { + nextNode.inOut = NODE_OUT; + } else if (nextNode.oneChildInOut == NODE_IN) { // 两种情况，Intersection: in, Difference: in-in，out - nextNode.inOut = nextNode.nodeOp == 1 ? 1 : 2; - } else if (nextNode.oneChildInOut == 2) { + nextNode.inOut = nextNode.nodeOp == OP_INTERSECTION ? NODE_IN : NODE_OUT; + } else if (nextNode.oneChildInOut == NODE_OUT) { // 两种情况，Intersection: out, Difference: NowNode一定是left，in-out，in - nextNode.inOut = nextNode.nodeOp == 1 ? 2 : 1; + nextNode.inOut = nextNode.nodeOp == OP_INTERSECTION ? NODE_OUT : NODE_IN; } else { - nextNode.oneChildInOut = 1; + nextNode.oneChildInOut = NODE_IN; return; } } - nextIdx = nowIdx; + nowIdx = nextIdx; } return; } -// 这里传入的BlobTree是拷贝，因为要修改且不同traversal的修改不相互影响 -// TODO: 最好优化掉拷贝 -// TODO: std::vector 可以改成bitset之TT +// TODO: std::vector 可以改成bitset之类 /** - * relatedPrimitives 是当次traversal需要care的primitives的索引，例如OcTree子问题中管理的Primitive + * relatedPrimitives 是当次traversal需要care的primitives的索引，例如OcTree子问题中关联的Primitive */ -bool traverse(BlobTree tree, const std::vector& relatedPrimitives, const std::vector& primitiveInout) +int traverse(BlobTree& tree, const std::vector& relatedPrimitives, const std::vector& primitiveInout) { assert(relatedPrimitives.size() == primitiveInout.size()); for (int primitiveIdx : relatedPrimitives) { propagate(tree, tree.primitiveNodeIdx[primitiveIdx], primitiveInout[primitiveIdx]); - if (tree.structure.back().inOut != 0) { return tree.structure.back().inOut == 1; } + if (tree.structure.back().inOut != 0) { return tree.structure.back().inOut; } } - std::cerr << "should not be here" << std::endl; - return false; + return 0; +} + +int traverse(BlobTree& tree, const int relatedPrimitive, const bool primitiveInout) +{ + propagate(tree, tree.primitiveNodeIdx[relatedPrimitive], primitiveInout); + return tree.structure.back().inOut; } }; // namespace algoim::organizer diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index 1b6bc31..f889ab1 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -4,7 +4,6 @@ #include #include - #include #include #include @@ -12,6 +11,7 @@ #include #include "bernstein.hpp" #include "multiloop.hpp" +#include "organizer/blobtree.hpp" #include "quadrature_multipoly.hpp" #include "binomial.hpp" @@ -32,20 +32,7 @@ namespace algoim::Organizer namespace detail { void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {} -} // namespace detail - -bool keepQuadraturePoint(std::vector& originTensors, const uvector3& originPt) -{ - // TODO: using blobtree to filter tensors - for (auto& t : originTensors) { - if (evalPower(t, originPt) >= 0) { return false; } - } - return true; -} - -// std::vector> primitives; -// BasicTask(std::vector> ps) {}; void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) { assert(0 <= k && k < 3); @@ -63,143 +50,6 @@ void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) } } -void basicTask(const std::shared_ptr& p, int q = 20, real xmin = -1, real xmax = 1) -{ - 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); - - makeMesh(*pt, compositeTensor, planeTensors); - detail::powerTransformation(range, xmin, compositeTensor); - - auto planeStdVector1 = xarray2StdVector(planeTensors[0]); - auto planeStdVector2 = xarray2StdVector(planeTensors[1]); - auto compositeTensorStdVector = xarray2StdVector(compositeTensor); - - 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); - - 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); - 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)); - 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); - - makeMesh(*pt, compositeTensor, planeTensors); - detail::powerTransformation(range, xmin, compositeTensor); - - real testEvalPower = evalPower(compositeTensor, testX); - - originTensors.insert(originTensors.end(), planeTensors.begin(), planeTensors.end()); - - detail::power2BernsteinTensor(compositeTensor, phi); - real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); - phis.emplace_back(phi); - } - } - 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) {} - template int numOfZero(const uvector& v) { @@ -259,28 +109,60 @@ std::pair, uvector> getOneEightCellAABB(const uvector3& } return res; } +} // namespace detail + +// TODO: delete it +bool keepQuadraturePoint(std::vector& originTensors, const uvector3& originPt) +{ + for (auto& t : originTensors) { + if (evalPower(t, originPt) >= 0) { return false; } + } + return true; +} -std::array sides = {-1, 1}; +bool keepQuadraturePoint(const Scene& scene, const OcTreeNode& ocTreeNode, const uvector3& originPt) +{ + // 只需要考虑intersect polys，不用考虑fully contained polys + for (auto& polyIntersectIdx : ocTreeNode.polyIntersectIndices) { + // TODO: + } + return true; +} + +// std::vector> primitives; + +// BasicTask(std::vector> ps) {}; + + +const std::array sides = {-1, 1}; + +struct BasicTaskRes { + real volume; + real area; +}; // 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合 // TODO: 参数太多了，考虑换用std::function + lambda void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, - , - const int polyIntersectIndex, - const tensor3& poly, + const int polyIntersectIndex, + const tensor3& poly, - std::array& subNodes, - int startIdxToCheck, - uvector& mark) + std::array& subNodes, + int startIdxToCheck, + uvector& mark) { - int zeroIdx = findFirst(mark, 0, startIdxToCheck); + int zeroIdx = detail::findFirst(mark, 0, startIdxToCheck); if (zeroIdx == -1) { tensor3 halfCellPoly(nullptr, poly.ext()); algoim_spark_alloc(real, halfCellPoly); - int subIdx = binaryToDecimal(mark); + int subIdx = detail::binaryToDecimal(mark); auto& subNode = subNodes[subIdx]; bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 - if (bernstein::uniformSign(halfCellPoly) != 1) { subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); } + if (bernstein::uniformSign(halfCellPoly) != 1) { + subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); + } else { + organizer::traverse(subNode.blobTree, polyIntersectIndex, organizer::NODE_OUT); + } } else { for (auto side : sides) { mark(zeroIdx) = side; @@ -292,23 +174,24 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, // 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历 // 所以树结构是不需要的记录的 -// void build(const Scene& scene, std::vector& intermediateNodes, std::vector leaves, int nowNodeIdx) -void build(const Scene& scene, const Node& node, std::vector 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) { 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::array subNodes; // intermediateNodes.resize(lastIdx + 8); - int subIdx = 0; + int subIdx = 0; for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { - auto nodeAABB = getOneEightCellAABB(node.min, node.max, j(), 0); - subNodes[subIdx].min = nodeAABB.first; - subNodes[subIdx].max = nodeAABB.second; + 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; } uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; for (int i = 0; i < polyIntersectIndices.size(); ++i) { @@ -317,9 +200,10 @@ void build(const Scene& scene, const Node& node, std::vector leaves) uvector mark(0, 0, 0); for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { real centerPlane = nodeMid(faceAxis); - tensor2 restrictToCenterPlane(nullptr, remove_component(poly.ext(), faceAxis)); + tensor2 restrictToCenterPlane(nullptr, remove_component(poly.compositedBernstein.ext(), faceAxis)); algoim_spark_alloc(real, restrictToCenterPlane); - restrictToInnerFace(poly, faceAxis, centerPlane, restrictToCenterPlane); + // NOTE: 这里restrict To Face用的是合成的tensor，这可能导致产生很多实际不需要考虑的交线 + detail::restrictToInnerFace(poly.compositedBernstein, faceAxis, centerPlane, restrictToCenterPlane); int signOnHalfPlane = bernstein::uniformSign(restrictToCenterPlane); if (signOnHalfPlane == -1) { // primitive contain the centerface @@ -329,16 +213,16 @@ void build(const Scene& scene, const Node& node, std::vector leaves) // deCasteljau to transformation uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax; halfCellMax(faceAxis) = nodeMid(faceAxis); - tensor3 halfCellPoly(nullptr, poly.ext()); + tensor3 halfCellPoly(nullptr, poly.compositedBernstein.ext()); algoim_spark_alloc(real, halfCellPoly); - bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); + bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly); if (bernstein::uniformSign(halfCellPoly) != 1) { // 负空间有 mark(faceAxis) = -1; } halfCellMax(faceAxis) = nowNodeMax(faceAxis); halfCellMin(faceAxis) = nodeMid(faceAxis); - bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); + bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly); if (bernstein::uniformSign(halfCellPoly) != 1) { // 正空间有 mark(faceAxis) += 1; // 当正负空间都有，记0 @@ -350,8 +234,8 @@ void build(const Scene& scene, const Node& node, std::vector leaves) if (all(mark == uvector3(2, 2, 2))) { // fully containing cases for (int i = 0; i < CHILD_NUM; ++i) { - tensor3 subPoly(nullptr, poly.ext()); - bernstein::deCasteljau(poly, subNodes[i].min, subNodes[i].max, subPoly); + tensor3 subPoly(nullptr, poly.compositedBernstein.ext()); + bernstein::deCasteljau(poly.compositedBernstein, subNodes[i].min, subNodes[i].max, subPoly); if (bernstein::uniformSign(subPoly) == -1) { subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex); } else { @@ -368,22 +252,34 @@ void build(const Scene& scene, const Node& node, std::vector leaves) } continue; } - int zeroCount = numOfZero<3>(mark); + int zeroCount = detail::numOfZero<3>(mark); if (zeroCount == 0) { // mark has -1 or 1 only - real nodeMidVal = bernstein::evalBernsteinPoly(poly, nodeMid); - if (mark(0) == 2 && mark(1) == 1 && mark(2) == 1) {} - subNodes[binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); + real nodeMidVal = bernstein::evalBernsteinPoly(poly.compositedBernstein, nodeMid); + int markIdx = detail::binaryToDecimal(mark); + for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { + if (subIdx == markIdx) + subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); + else + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); + } } else if (zeroCount == 1) { // poly related to 2 subcells uvector sideMark = mark; - int zeroIdx = replaceFirst(sideMark, 0, 1); - subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); - sideMark(zeroIdx) = -1; - subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); + int zeroIdx = detail::replaceFirst(sideMark, 0, 1); + int sideIdx1 = detail::binaryToDecimal(sideMark, -1); + sideMark(zeroIdx) = -1; + int sideIdx2 = detail::binaryToDecimal(sideMark, -1); + subNodes[sideIdx1].polyIntersectIndices.emplace_back(polyIntersectIndex); + subNodes[sideIdx2].polyIntersectIndices.emplace_back(polyIntersectIndex); + for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { + if (subIdx != sideIdx1 && subIdx != sideIdx2) { + organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); + } + } } else { // zeroCount == 2 or 3 - visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, 0, mark); + visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly.compositedBernstein, subNodes, 0, mark); } } // launch subdivision in subcells @@ -393,14 +289,231 @@ void build(const Scene& scene, const Node& node, std::vector leaves) subNodes[subIdx].polyFullyContainedIndices.insert(subNodes[subIdx].polyFullyContainedIndices.end(), node.polyFullyContainedIndices.begin(), node.polyFullyContainedIndices.end()); - build(scene, subNodes[subIdx], leaves); + buildOcTree(scene, subNodes[subIdx], leaves); + } +} + +/** single object quadrature */ +void basicTask(const std::shared_ptr& p, int q = 20, real xmin = -1, real xmax = 1) +{ + 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); + 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); + + makeMesh(*pt, compositeTensor, planeTensors); + detail::powerTransformation(range, xmin, compositeTensor); + + auto planeStdVector1 = xarray2StdVector(planeTensors[0]); + auto planeStdVector2 = xarray2StdVector(planeTensors[1]); + auto compositeTensorStdVector = xarray2StdVector(compositeTensor); + + 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); + + 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; } - // TODO: 考虑fully contain 问题 + 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); + 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)); + 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); + + makeMesh(*pt, compositeTensor, planeTensors); + detail::powerTransformation(range, xmin, compositeTensor); + + real testEvalPower = evalPower(compositeTensor, testX); + + originTensors.insert(originTensors.end(), planeTensors.begin(), planeTensors.end()); + + detail::power2BernsteinTensor(compositeTensor, phi); + real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); + phis.emplace_back(phi); + } + } + 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, further deallocating memory of xarray + for (auto& p : phiStacks) delete p; + for (auto& p : originTensorStacks) delete p; +}; + +BasicTaskRes basicTask(const Scene& scene, const OcTreeNode& node, int q = 20) +{ + auto integrand = [](const uvector& x) { return 1.0; }; + real volume = 0., surf = 0.; + auto range = node.max - node.min; + ImplicitPolyQuadrature<3> ipquad(scene.polys); + ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { + auto realX = x * range + node.min; + if (keepQuadraturePoint(scene, node, realX)) volume += w * integrand(realX); + }); +} + +void quadratureScene(const std::vector>& primitives, + const uvector3& xmin, + const uvector3& xmax, + const organizer::BlobTree& blobTree) +{ + OcTreeNode rootNode(xmin, xmax, blobTree); + std::vector leaves; + + std::vector*> tensorStacks; + std::vector completeTensorReps; + + 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); + tensor3 transformedTensor(nullptr, 3); + tensorStacks.emplace_back(algoim_spark_alloc_heap(real, originTensor)); + + tensor3 phi(nullptr, 3); + tensorStacks.emplace_back( + algoim_spark_alloc_heap(real, phi)); // 必须先于algoim_spark_alloc，使transformedTensor的内存以栈形式释放 + + algoim_spark_alloc(real, transformedTensor); + makeSphere(*pt, originTensor); + + detail::powerTransformation(range, xmin, originTensor, transformedTensor); + + detail::power2BernsteinTensor(transformedTensor, phi); + completeTensorReps.emplace_back(CompleteTensorRep{phi, {originTensor}}); + } 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(tensorStacks, planeTensors); + + tensor3 phi(nullptr, 1 + faceCount); + tensorStacks.emplace_back( + algoim_spark_alloc_heap(real, phi)); // 必须先于algoim_spark_alloc，使compositeTensor的内存以栈形式释放 + + tensor3 compositeTensor(nullptr, 1 + faceCount); + algoim_spark_alloc(real, compositeTensor); + + makeMesh(*pt, compositeTensor, planeTensors); + detail::powerTransformation(range, xmin, compositeTensor); + + real testEvalPower = evalPower(compositeTensor, testX); + + detail::power2BernsteinTensor(compositeTensor, phi); + real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); + + completeTensorReps.emplace_back(CompleteTensorRep{phi, planeTensors}); + } + } + + Scene scene{completeTensorReps, blobTree}; + buildOcTree(scene, rootNode, leaves); + + for (const auto& leaf : leaves) { + auto basicRes = basicTask(scene, leaf, xmin, xmax, 10); + volume += basicRes.volume; + } + + volume *= prod(xmax - xmin); + // TODO: surface area + std::cout << "Volume xxx: " << volume << std::endl; + + + // free memory, further deallocating memory of xarray + for (auto& p : tensorStacks) delete p; } }; // namespace algoim::Organizer // clang-format off -/** +/** 最一开始的cell应当包含所有的primitive func(given cell) { 得到8个subcell @@ -408,7 +521,7 @@ void build(const Scene& scene, const Node& node, std::vector leaves) mark = [0,0,0] for k in 0,1,2 { 测试primitive和k轴向的centerface有无交 - if 无交 { + if 无交 {0 if primitive包裹centerface { mark[k] = 2 } else { diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 494eb68..3def130 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -12,11 +12,13 @@ #include "xarray.hpp" #include "binomial.hpp" #include "bernstein.hpp" +#include "blobtree.hpp" #include namespace algoim::Organizer { + namespace detail { @@ -589,29 +591,39 @@ void makeCylinder(xarray& tensor, uvector3 startPt, uvector3 endPt, rea // TODO: } +struct CompleteTensorRep { + tensor3 compositedBernstein; // 合法输入在[0,1]^3 + std::vector originalPower; // 合法输入在场景 +}; + struct Scene { - std::vector polys; - int* boolDescription; // blobtree + std::vector polys; + organizer::BlobTree blobTree; // uvector3 min, max; }; const int CHILD_NUM = 8; -struct Node { +struct OcTreeNode { std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 // std::array children; // octree uvector3 min, max; + organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递 // Node() { children.fill(-1); } - Node() = default; + OcTreeNode() = default; + + OcTreeNode(const uvector3& min_, const uvector3& max_, const organizer::BlobTree& blobTree_) + : min(min_), max(max_), blobTree(blobTree_) {}; // 深拷贝 - Node(const Node& node) + OcTreeNode(const OcTreeNode& node) : polyIntersectIndices(node.polyIntersectIndices), polyFullyContainedIndices(node.polyFullyContainedIndices), min(node.min), - max(node.max) + max(node.max), + blobTree(node.blobTree) { } }; diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index 7f333a2..befbb11 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -19,6 +19,7 @@ #include #include #include +#include "organizer/primitive.hpp" #define ALGOIM_M 8 #include @@ -568,6 +569,15 @@ struct ImplicitPolyQuadrature { build(true, false); } + ImplicitPolyQuadrature(const std::vector& completeTensorRes) + { + for (const auto& ctr : completeTensorRes) { + auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray(true)); + if (!detail::maskEmpty(mask)) phi.push_back(ctr.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, @@ -807,7 +817,7 @@ struct ImplicitPolyQuadrature { // Apply primary base integral k_active = k; base.integrate(strategy, q, integrand); - // situation the + // situation the // If in aggregate mode, apply to other dimensions as well if (type == OuterAggregate) { for (int i = 0; i < N - 1; ++i) {