GJJ
/
HeteQuadrature
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

blobtree octree integration

master
gjj 1 year ago
parent
c490eea4bf
commit
5e8545d770
5 changed files with 579 additions and 466 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 437
      algoim/kdtree.hpp
  2. 79
      algoim/organizer/blobtree.hpp
  3. 493
      algoim/organizer/organizer.hpp
  4. 24
      algoim/organizer/primitive.hpp
  5. 12
      algoim/quadrature_multipoly.hpp

437
algoim/kdtree.hpp
View File

@ -2,11 +2,11 @@
#define ALGOIM_KDTREE_HPP #define ALGOIM_KDTREE_HPP
/* algoim::KDTree<T,N> constructs a k-d tree data structure for a given collection of /* algoim::KDTree<T,N> constructs a k-d tree data structure for a given collection of
points of type uvector<T,N>, where T is typically float or double, and N is points of type uvector<T,N>, 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 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 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. transformations that result in "tight" bounding boxes for some of the nodes in the tree.
Refer to the paper Refer to the paper
R. I. Saye, High-order methods for computing distances to implicitly defined surfaces, 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), Communications in Applied Mathematics and Computational Science, 9(1), 107-141 (2014),
@ -21,272 +21,243 @@
namespace algoim namespace algoim
{ {
namespace detail 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) // Returns the squared-distance from a point x to the hyperrectangle [min,max] (if x is inside the rectangle, the distance is
template<typename T, int N> // zero)
T sqrDistBox(const uvector<T,N>& x, const uvector<T,N>& min, const uvector<T,N>& max) template <typename T, int N>
{ T sqrDistBox(const uvector<T, N>& x, const uvector<T, N>& min, const uvector<T, N>& max)
T dsqr = T(0); {
for (int i = 0; i < N; ++i) T dsqr = T(0);
if (x(i) < min(i)) for (int i = 0; i < N; ++i)
dsqr += util::sqr(x(i) - min(i)); if (x(i) < min(i))
else if (x(i) > max(i)) dsqr += util::sqr(x(i) - min(i));
dsqr += util::sqr(x(i) - max(i)); else if (x(i) > max(i))
return dsqr; dsqr += util::sqr(x(i) - max(i));
} return dsqr;
} }
} // namespace detail
template<typename T, int N> template <typename T, int N>
class KDTree class KDTree
{ {
std::vector<uvector<T,N>> points; std::vector<uvector<T, N>> points;
std::vector<int> index; std::vector<int> index;
static constexpr int leaf_size = 16; 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 // 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 // is applied. If type < -1, the node is a standard splitting node with two children and no transform
struct Node struct Node {
{ int type;
int type; int i0, i1;
int i0, i1; uvector<T, N> xmin, xmax;
uvector<T,N> xmin, xmax; };
};
std::vector<Node> nodes; std::vector<Node> nodes;
std::vector<uvector<uvector<T,N>,N>> transforms; std::vector<uvector<uvector<T, N>, N>> transforms;
// Given a node and a range of points [lb,ub), build the tree // 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) void build_tree(int nodeIndex, int lb, int ub, bool hasTransformed, int level)
{ {
assert(lb < ub); assert(lb < ub);
Node& node = nodes[nodeIndex]; Node& node = nodes[nodeIndex];
// Compute bounding box and mean // Compute bounding box and mean
uvector<T,N> mean = node.xmin = node.xmax = points[index[lb]]; uvector<T, N> mean = node.xmin = node.xmax = points[index[lb]];
for (int i = lb + 1; i < ub; ++i) for (int i = lb + 1; i < ub; ++i) {
{ const uvector<T, N>& x = points[index[i]];
const uvector<T,N>& x = points[index[i]]; mean += x;
mean += x; for (int j = 0; j < N; ++j) {
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);
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<T>(ub - lb); }
mean /= static_cast<T>(ub - lb);
// The node is a leaf iff point count is sufficiently small // The node is a leaf iff point count is sufficiently small
if (ub - lb <= leaf_size) if (ub - lb <= leaf_size) {
{ node.type = -1;
node.type = -1; node.i0 = lb;
node.i0 = lb; node.i1 = ub;
node.i1 = ub; return;
return; }
}
// Splitting node: default to splitting along greatest extent // Splitting node: default to splitting along greatest extent
node.type = -2; node.type = -2;
int axis = argmax(node.xmax - node.xmin); int axis = argmax(node.xmax - node.xmin);
// Evaluate possibility for coordinate transformation // Evaluate possibility for coordinate transformation
if (!hasTransformed && level > 5 && ub - lb >= leaf_size * (1 << 2)) if (!hasTransformed && level > 5 && ub - lb >= leaf_size * (1 << 2)) {
{ // Estimate normal
// Estimate normal T holeRadiusSqr = util::sqr(0.05 * max(node.xmax - node.xmin));
T holeRadiusSqr = util::sqr(0.05*max(node.xmax - node.xmin)); uvector<T, N> n = static_cast<T>(0);
uvector<T,N> n = static_cast<T>(0); n(0) = 1.0;
for (int i = lb; i < ub; ++i) {
uvector<T, N> 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; n(0) = 1.0;
for (int i = lb; i < ub; ++i) else
{ n /= sqrt(msqr);
uvector<T,N> 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);
// Compute alpha range // Compute alpha range
T minAlpha = std::numeric_limits<T>::max(); T minAlpha = std::numeric_limits<T>::max();
T maxAlpha = -std::numeric_limits<T>::max(); T maxAlpha = -std::numeric_limits<T>::max();
for (int i = lb; i < ub; ++i) for (int i = lb; i < ub; ++i) {
{ T alpha = dot(points[index[i]], n);
T alpha = dot(points[index[i]], n); if (alpha > maxAlpha) maxAlpha = alpha;
if (alpha > maxAlpha) maxAlpha = alpha; if (alpha < minAlpha) minAlpha = alpha;
if (alpha < minAlpha) minAlpha = alpha; }
}
if (maxAlpha - minAlpha < 0.1*max(node.xmax - node.xmin)) if (maxAlpha - minAlpha < 0.1 * max(node.xmax - node.xmin)) {
{ // Perform transformation: calculate an orthonormal basis using the normal as first axis.
// 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,
// 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,
// 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.
// j chosen that that n != ej. uvector<uvector<T, N>, N> axes;
uvector<uvector<T,N>,N> axes; int j = argmin(abs(n));
int j = argmin(abs(n)); uvector<T, N> u = n;
uvector<T,N> u = n; u(j) -= 1.0; u(j) -= 1.0;
u /= norm(u); u /= norm(u);
for (int dim = 0; dim < N; ++dim) for (int dim = 0; dim < N; ++dim)
for (int i = 0; i < N; ++i) for (int i = 0; i < N; ++i) axes(dim)(i) = (dim == i ? 1.0 : 0.0) - 2.0 * u(dim) * u(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 // 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). // unnecessary (but done in the name of consistency with old approach).
std::swap(axes(0), axes(j)); std::swap(axes(0), axes(j));
// Apply coordinate transformation and calculate new bounding box in order to determine new split direction // Apply coordinate transformation and calculate new bounding box in order to determine new split direction
uvector<T,N> bmin = std::numeric_limits<T>::max(); uvector<T, N> bmin = std::numeric_limits<T>::max();
uvector<T,N> bmax = -std::numeric_limits<T>::max(); uvector<T, N> bmax = -std::numeric_limits<T>::max();
for (int i = lb; i < ub; ++i) for (int i = lb; i < ub; ++i) {
{ uvector<T, N> x = points[index[i]];
uvector<T,N> x = points[index[i]]; for (int dim = 0; dim < N; ++dim) {
for (int dim = 0; dim < N; ++dim) T alpha = dot(axes(dim), x);
{ points[index[i]](dim) = alpha;
T alpha = dot(axes(dim), x); if (alpha < bmin(dim)) bmin(dim) = alpha;
points[index[i]](dim) = alpha; if (alpha > bmax(dim)) bmax(dim) = alpha;
if (alpha < bmin(dim)) bmin(dim) = alpha;
if (alpha > bmax(dim)) bmax(dim) = alpha;
}
} }
node.type = static_cast<int>(transforms.size());
transforms.push_back(axes);
axis = argmax(bmax - bmin);
hasTransformed = true;
} }
node.type = static_cast<int>(transforms.size());
transforms.push_back(axes);
axis = argmax(bmax - bmin);
hasTransformed = true;
} }
}
// Use median as the split // Use median as the split
int m = (lb + ub)/2; int m = (lb + ub) / 2;
// Rearrange points // 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); } ); 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 // Build child trees
int i0 = node.i0 = static_cast<int>(nodes.size()); int i0 = node.i0 = static_cast<int>(nodes.size());
int i1 = node.i1 = static_cast<int>(nodes.size() + 1); int i1 = node.i1 = static_cast<int>(nodes.size() + 1);
nodes.push_back(Node()); nodes.push_back(Node());
nodes.push_back(Node()); nodes.push_back(Node());
build_tree(i0, lb, m, hasTransformed, level + 1); build_tree(i0, lb, m, hasTransformed, level + 1);
build_tree(i1, m, ub, hasTransformed, level + 1); build_tree(i1, m, ub, hasTransformed, level + 1);
} }
struct ClosestPoint struct ClosestPoint {
{ uvector<T, N> x;
uvector<T,N> x; T distsqr;
T distsqr; int ind;
int ind; };
};
// Recursive function for searching the tree for the closest point to a given point cp.x // Recursive function for searching the tree for the closest point to a given point cp.x
void search(const Node& node, ClosestPoint& cp) const void search(const Node& node, ClosestPoint& cp) const
{ {
if (node.type == -1) if (node.type == -1) {
{ // Leaf node
// Leaf node for (int j = node.i0; j < node.i1; ++j) {
for (int j = node.i0; j < node.i1; ++j) T dsqr = sqrnorm(points[j] - cp.x);
{ if (dsqr < cp.distsqr) {
T dsqr = sqrnorm(points[j] - cp.x); cp.distsqr = dsqr;
if (dsqr < cp.distsqr) cp.ind = j;
{
cp.distsqr = dsqr;
cp.ind = j;
}
} }
} }
else } else {
{ // Non-leaf node
// Non-leaf node if (node.type >= 0) {
if (node.type >= 0) // Transform query point to new coordinate system
{ const uvector<uvector<T, N>, N>& axes = transforms[node.type];
// Transform query point to new coordinate system uvector<T, N> x = cp.x;
const uvector<uvector<T,N>,N>& axes = transforms[node.type]; for (int dim = 0; dim < N; ++dim) cp.x(dim) = dot(axes(dim), x);
uvector<T,N> 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 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); T dright = detail::sqrDistBox(cp.x, nodes[node.i1].xmin, nodes[node.i1].xmax);
if (dleft < dright) if (dleft < dright) {
{ if (dleft < cp.distsqr) {
if (dleft < cp.distsqr) search(nodes[node.i0], cp);
{ if (dright < cp.distsqr) search(nodes[node.i1], cp);
search(nodes[node.i0], cp);
if (dright < cp.distsqr)
search(nodes[node.i1], cp);
}
} }
else } else {
{ if (dright < cp.distsqr) {
if (dright < cp.distsqr) search(nodes[node.i1], cp);
{ if (dleft < cp.distsqr) search(nodes[node.i0], cp);
search(nodes[node.i1], cp);
if (dleft < cp.distsqr)
search(nodes[node.i0], cp);
}
} }
}
if (node.type >= 0) if (node.type >= 0) {
{ // Transform query point back to old coordinate system. This is about 5% faster than storing
// Transform query point back to old coordinate system. This is about 5% faster than storing // the old value of cp.x on the stack
// the old value of cp.x on the stack const uvector<uvector<T, N>, N>& axes = transforms[node.type];
const uvector<uvector<T,N>,N>& axes = transforms[node.type]; uvector<T, N> x = cp.x;
uvector<T,N> x = cp.x; cp.x = axes(0) * x(0);
cp.x = axes(0)*x(0); for (int dim = 1; dim < N; ++dim) cp.x += axes(dim) * x(dim);
for (int dim = 1; dim < N; ++dim)
cp.x += axes(dim)*x(dim);
}
} }
} }
}
public: public:
// Construct a KDTree from a given collection of points p (the given points are not overwritten)
// Construct a KDTree from a given collection of points p (the given points are not overwritten) KDTree(const std::vector<uvector<T, N>>& p)
KDTree(const std::vector<uvector<T,N>>& p) {
{ assert(p.size()
assert(p.size() < std::numeric_limits<int>::max()); // code currently uses int to index but could easily be adapted to size_t < std::numeric_limits<int>::max()); // code currently uses int to index but could easily be adapted to size_t
int len = static_cast<int>(p.size()); int len = static_cast<int>(p.size());
if (len == 0) if (len == 0) return;
return;
// Copy points // Copy points
points = p; points = p;
// Build initial index array, which shall soon be reordered // Build initial index array, which shall soon be reordered
index.resize(len); index.resize(len);
for (int i = 0; i < len; ++i) for (int i = 0; i < len; ++i) index[i] = i;
index[i] = i;
// Recursively build tree starting from root node. This only manipulates index array // Recursively build tree starting from root node. This only manipulates index array
nodes.push_back(Node()); nodes.push_back(Node());
build_tree(0, 0, len, false, 0); build_tree(0, 0, len, false, 0);
// Rearrange so that points in leaf nodes are contiguous in memory. Based on tree // Rearrange so that points in leaf nodes are contiguous in memory. Based on tree
// construction, both nodes & point ranges will be organised depth-first. // construction, both nodes & point ranges will be organised depth-first.
std::vector<uvector<T,N>> pointscopy(points); std::vector<uvector<T, N>> pointscopy(points);
for (size_t i = 0; i < len; ++i) for (size_t i = 0; i < len; ++i) points[i] = pointscopy[index[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 // 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 // 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. // of the closest point in the original array p passed to the KDTree constructor is returned.
int nearest(const uvector<T,N>& x, T rsqr = std::numeric_limits<T>::max()) const int nearest(const uvector<T, N>& x, T rsqr = std::numeric_limits<T>::max()) const
{ {
if (nodes.empty()) if (nodes.empty()) return -1;
return -1; ClosestPoint cp;
ClosestPoint cp; cp.x = x;
cp.x = x; cp.distsqr = rsqr;
cp.distsqr = rsqr; cp.ind = -1;
cp.ind = -1; search(nodes[0], cp);
search(nodes[0], cp); return cp.ind >= 0 ? index[cp.ind] : -1;
return cp.ind >= 0 ? index[cp.ind] : -1; }
} };
};
} // namespace algoim } // namespace algoim
#endif #endif

79
algoim/organizer/blobtree.hpp
View File

@ -5,6 +5,14 @@
namespace algoim::organizer 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 // bitfield
struct Blob { struct Blob {
unsigned int isPrimitive : 1; // 1 for leaf 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; const std::size_t rootIdx = tree.structure.size() - 1;
auto& node = tree.structure[nodeIdx]; auto& node = tree.structure[nodeIdx];
assert(node.inOut == 0); assert(node.inOut == 0);
node.inOut = in ? 1 : 2; node.inOut = in ? NODE_IN : NODE_OUT;
for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx; ++nodeIdx) { for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) {
const auto& nowNode = tree.structure[nowIdx]; const auto& nowNode = tree.structure[nowIdx];
unsigned int nextIdx = isLeft(nowNode) ? nowNode.ancestor : nodeIdx + 1; 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 (nextNode.inOut != 0) { return; }
if (nowNode.inOut == 2) { if (nowNode.inOut == NODE_OUT) {
// out if (nextNode.nodeOp == OP_INTERSECTION) {
if (nextNode.nodeOp == 1) { nextNode.inOut = NODE_OUT;
// intersection } else if (nextNode.nodeOp == OP_DIFFERENCE && isLeft(nowNode)) {
} else if (nextNode.nodeOp == 2 && isLeft(nowNode)) { nextNode.inOut = NODE_OUT;
// difference and left } else if (nextNode.oneChildInOut == NODE_IN) {
nextNode.inOut = 2;
} else if (nextNode.oneChildInOut == 1) {
// 两种情况，Union: in, Difference: NowNode一定是right，in-out，还是in // 两种情况，Union: in, Difference: NowNode一定是right，in-out，还是in
nextNode.inOut = 1; nextNode.inOut = NODE_IN;
} else if (nextNode.oneChildInOut == 2) { } else if (nextNode.oneChildInOut == NODE_OUT) {
// 两种情况，Union: out, Difference: out-out, 还是out // 两种情况，Union: out, Difference: out-out, 还是out
nextNode.inOut = 2; nextNode.inOut = NODE_OUT;
} else { } else {
nextNode.oneChildInOut = 2; nextNode.oneChildInOut = NODE_OUT;
return; return;
} }
} else if (nowNode.inOut == 1) { } else if (nowNode.inOut == NODE_IN) {
// in if (nextNode.nodeOp == OP_UNION) {
if (nextNode.nodeOp == 0) { nextNode.inOut = NODE_IN;
// union } else if (nextNode.nodeOp == OP_DIFFERENCE && !isLeft(nowNode)) {
nextNode.inOut = 1;
} else if (nextNode.nodeOp == 2 && !isLeft(nowNode)) {
// difference and right // difference and right
nextNode.inOut = 2; nextNode.inOut = NODE_OUT;
} else if (nextNode.oneChildInOut == 1) { } else if (nextNode.oneChildInOut == NODE_IN) {
// 两种情况，Intersection: in, Difference: in-in，out // 两种情况，Intersection: in, Difference: in-in，out
nextNode.inOut = nextNode.nodeOp == 1 ? 1 : 2; nextNode.inOut = nextNode.nodeOp == OP_INTERSECTION ? NODE_IN : NODE_OUT;
} else if (nextNode.oneChildInOut == 2) { } else if (nextNode.oneChildInOut == NODE_OUT) {
// 两种情况，Intersection: out, Difference: NowNode一定是left，in-out，in // 两种情况，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 { } else {
nextNode.oneChildInOut = 1; nextNode.oneChildInOut = NODE_IN;
return; return;
} }
} }
nextIdx = nowIdx; nowIdx = nextIdx;
} }
return; return;
} }
// 这里传入的BlobTree是拷贝，因为要修改且不同traversal的修改不相互影响 // TODO: std::vector<bool> 可以改成bitset之类
// TODO: 最好优化掉拷贝
// TODO: std::vector<bool> 可以改成bitset之TT
/** /**
* relatedPrimitives 是当次traversal需要care的primitives的索引，例如OcTree子问题中管理的Primitive * relatedPrimitives 是当次traversal需要care的primitives的索引，例如OcTree子问题中关联的Primitive
*/ */
bool traverse(BlobTree tree, const std::vector<int>& relatedPrimitives, const std::vector<bool>& primitiveInout) int traverse(BlobTree& tree, const std::vector<int>& relatedPrimitives, const std::vector<bool>& primitiveInout)
{ {
assert(relatedPrimitives.size() == primitiveInout.size()); assert(relatedPrimitives.size() == primitiveInout.size());
for (int primitiveIdx : relatedPrimitives) { for (int primitiveIdx : relatedPrimitives) {
propagate(tree, tree.primitiveNodeIdx[primitiveIdx], primitiveInout[primitiveIdx]); 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 0;
return false; }
int traverse(BlobTree& tree, const int relatedPrimitive, const bool primitiveInout)
{
propagate(tree, tree.primitiveNodeIdx[relatedPrimitive], primitiveInout);
return tree.structure.back().inOut;
} }
}; // namespace algoim::organizer }; // namespace algoim::organizer

493
algoim/organizer/organizer.hpp
View File

@ -4,7 +4,6 @@
#include <iostream> #include <iostream>
#include <booluarray.hpp> #include <booluarray.hpp>
#include <cstddef> #include <cstddef>
#include <iostream> #include <iostream>
#include <iomanip> #include <iomanip>
@ -12,6 +11,7 @@
#include <vector> #include <vector>
#include "bernstein.hpp" #include "bernstein.hpp"
#include "multiloop.hpp" #include "multiloop.hpp"
#include "organizer/blobtree.hpp"
#include "quadrature_multipoly.hpp" #include "quadrature_multipoly.hpp"
#include "binomial.hpp" #include "binomial.hpp"
@ -32,20 +32,7 @@ namespace algoim::Organizer
namespace detail namespace detail
{ {
void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {} void bernstein2PowerTensor(const tensor3& phiBernstein, tensor3& phiPower) {}
} // namespace detail
bool keepQuadraturePoint(std::vector<tensor3>& 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<std::shared_ptr<PrimitiveDesc>> primitives;
// BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res)
{ {
assert(0 <= k && k < 3); 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<PrimitiveDesc>& p, int q = 20, real xmin = -1, real xmax = 1)
{
real volume = 0;
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
uvector3 range = xmax - xmin;
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(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<real, 3> 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<real, 3>& x, real w) {
if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin));
});
} else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
const int faceCount = pt->indexInclusiveScan.size();
assert(faceCount > 1);
std::vector<tensor3> 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<real, 3> 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<real, 3>& 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<std::shared_ptr<PrimitiveDesc>>& primitives, int q = 20, real xmin = -1, real xmax = 1)
{
std::vector<SparkStack<real>*> phiStacks;
std::vector<tensor3> phis;
std::vector<SparkStack<real>*> originTensorStacks;
std::vector<tensor3> originTensors;
real volume;
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
uvector3 range = xmax - xmin;
uvector<real, 3> testX(0., 0.75, 0.2);
for (int i = 0; i < primitives.size(); i++) {
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(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<MeshDesc>(primitives[i])) {
const int faceCount = pt->indexInclusiveScan.size();
assert(faceCount > 1);
std::vector<tensor3> 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<real, 3>& 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<Node>& nodes, const uvector3& min, const uvector3& max) {}
template <int N> template <int N>
int numOfZero(const uvector<int, N>& v) int numOfZero(const uvector<int, N>& v)
{ {
@ -259,28 +109,60 @@ std::pair<uvector<int, 3>, uvector<int, 3>> getOneEightCellAABB(const uvector3&
} }
return res; return res;
} }
} // namespace detail
// TODO: delete it
bool keepQuadraturePoint(std::vector<tensor3>& originTensors, const uvector3& originPt)
{
for (auto& t : originTensors) {
if (evalPower(t, originPt) >= 0) { return false; }
}
return true;
}
std::array<int, 2> 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<std::shared_ptr<PrimitiveDesc>> primitives;
// BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
const std::array<int, 2> sides = {-1, 1};
struct BasicTaskRes {
real volume;
real area;
};
// 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合 // 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合
// TODO: 参数太多了，考虑换用std::function + lambda // TODO: 参数太多了，考虑换用std::function + lambda
void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid,
, const int polyIntersectIndex,
const int polyIntersectIndex, const tensor3& poly,
const tensor3& poly,
std::array<Node, CHILD_NUM>& subNodes, std::array<OcTreeNode, CHILD_NUM>& subNodes,
int startIdxToCheck, int startIdxToCheck,
uvector<int, 3>& mark) uvector<int, 3>& mark)
{ {
int zeroIdx = findFirst(mark, 0, startIdxToCheck); int zeroIdx = detail::findFirst(mark, 0, startIdxToCheck);
if (zeroIdx == -1) { if (zeroIdx == -1) {
tensor3 halfCellPoly(nullptr, poly.ext()); tensor3 halfCellPoly(nullptr, poly.ext());
algoim_spark_alloc(real, halfCellPoly); algoim_spark_alloc(real, halfCellPoly);
int subIdx = binaryToDecimal(mark); int subIdx = detail::binaryToDecimal(mark);
auto& subNode = subNodes[subIdx]; auto& subNode = subNodes[subIdx];
bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 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 { } else {
for (auto side : sides) { for (auto side : sides) {
mark(zeroIdx) = side; mark(zeroIdx) = side;
@ -292,23 +174,24 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid,
// 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历 // 叶节点单独存在leaves中，是因为最后只需要叶节点信息，没有自顶向下/自底向上遍历
// 所以树结构是不需要的记录的 // 所以树结构是不需要的记录的
// void build(const Scene& scene, std::vector<Node>& intermediateNodes, std::vector<Node> leaves, int nowNodeIdx) // void buildOcTree(const Scene& scene, std::vector<Node>& intermediateNodes, std::vector<Node> leaves, int nowNodeIdx)
void build(const Scene& scene, const Node& node, std::vector<Node> leaves) void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves)
{ {
const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices; const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
if (polyIntersectIndices.size() <= 4) { if (polyIntersectIndices.size() <= 4) {
leaves.emplace_back(node); leaves.emplace_back(node);
return; return;
} }
const uvector3& nowNodeMin = node.min; const uvector3& nowNodeMin = node.min;
const uvector3& nowNodeMax = node.max; const uvector3& nowNodeMax = node.max;
std::array<Node, CHILD_NUM> subNodes; std::array<OcTreeNode, CHILD_NUM> subNodes;
// intermediateNodes.resize(lastIdx + 8); // intermediateNodes.resize(lastIdx + 8);
int subIdx = 0; int subIdx = 0;
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
auto nodeAABB = getOneEightCellAABB(node.min, node.max, j(), 0); auto nodeAABB = detail::getOneEightCellAABB(node.min, node.max, j(), 0);
subNodes[subIdx].min = nodeAABB.first; subNodes[subIdx].min = nodeAABB.first;
subNodes[subIdx].max = nodeAABB.second; subNodes[subIdx].max = nodeAABB.second;
subNodes[subIdx].blobTree = node.blobTree;
} }
uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2;
for (int i = 0; i < polyIntersectIndices.size(); ++i) { for (int i = 0; i < polyIntersectIndices.size(); ++i) {
@ -317,9 +200,10 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
uvector<int, 3> mark(0, 0, 0); uvector<int, 3> mark(0, 0, 0);
for (int faceAxis = 0; faceAxis < 3; ++faceAxis) { for (int faceAxis = 0; faceAxis < 3; ++faceAxis) {
real centerPlane = nodeMid(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); 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); int signOnHalfPlane = bernstein::uniformSign(restrictToCenterPlane);
if (signOnHalfPlane == -1) { if (signOnHalfPlane == -1) {
// primitive contain the centerface // primitive contain the centerface
@ -329,16 +213,16 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
// deCasteljau to transformation // deCasteljau to transformation
uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax; uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax;
halfCellMax(faceAxis) = nodeMid(faceAxis); halfCellMax(faceAxis) = nodeMid(faceAxis);
tensor3 halfCellPoly(nullptr, poly.ext()); tensor3 halfCellPoly(nullptr, poly.compositedBernstein.ext());
algoim_spark_alloc(real, halfCellPoly); algoim_spark_alloc(real, halfCellPoly);
bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly);
if (bernstein::uniformSign(halfCellPoly) != 1) { if (bernstein::uniformSign(halfCellPoly) != 1) {
// 负空间有 // 负空间有
mark(faceAxis) = -1; mark(faceAxis) = -1;
} }
halfCellMax(faceAxis) = nowNodeMax(faceAxis); halfCellMax(faceAxis) = nowNodeMax(faceAxis);
halfCellMin(faceAxis) = nodeMid(faceAxis); halfCellMin(faceAxis) = nodeMid(faceAxis);
bernstein::deCasteljau(poly, halfCellMin, halfCellMax, halfCellPoly); bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly);
if (bernstein::uniformSign(halfCellPoly) != 1) { if (bernstein::uniformSign(halfCellPoly) != 1) {
// 正空间有 // 正空间有
mark(faceAxis) += 1; // 当正负空间都有，记0 mark(faceAxis) += 1; // 当正负空间都有，记0
@ -350,8 +234,8 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
if (all(mark == uvector3(2, 2, 2))) { if (all(mark == uvector3(2, 2, 2))) {
// fully containing cases // fully containing cases
for (int i = 0; i < CHILD_NUM; ++i) { for (int i = 0; i < CHILD_NUM; ++i) {
tensor3 subPoly(nullptr, poly.ext()); tensor3 subPoly(nullptr, poly.compositedBernstein.ext());
bernstein::deCasteljau(poly, subNodes[i].min, subNodes[i].max, subPoly); bernstein::deCasteljau(poly.compositedBernstein, subNodes[i].min, subNodes[i].max, subPoly);
if (bernstein::uniformSign(subPoly) == -1) { if (bernstein::uniformSign(subPoly) == -1) {
subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex); subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex);
} else { } else {
@ -368,22 +252,34 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
} }
continue; continue;
} }
int zeroCount = numOfZero<3>(mark); int zeroCount = detail::numOfZero<3>(mark);
if (zeroCount == 0) { if (zeroCount == 0) {
// mark has -1 or 1 only // mark has -1 or 1 only
real nodeMidVal = bernstein::evalBernsteinPoly(poly, nodeMid); real nodeMidVal = bernstein::evalBernsteinPoly(poly.compositedBernstein, nodeMid);
if (mark(0) == 2 && mark(1) == 1 && mark(2) == 1) {} int markIdx = detail::binaryToDecimal(mark);
subNodes[binaryToDecimal(mark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); 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) { } else if (zeroCount == 1) {
// poly related to 2 subcells // poly related to 2 subcells
uvector<int, 3> sideMark = mark; uvector<int, 3> sideMark = mark;
int zeroIdx = replaceFirst(sideMark, 0, 1); int zeroIdx = detail::replaceFirst(sideMark, 0, 1);
subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); int sideIdx1 = detail::binaryToDecimal(sideMark, -1);
sideMark(zeroIdx) = -1; sideMark(zeroIdx) = -1;
subNodes[binaryToDecimal(sideMark, -1)].polyIntersectIndices.emplace_back(polyIntersectIndex); 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 { } else {
// zeroCount == 2 or 3 // zeroCount == 2 or 3
visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly, subNodes, 0, mark); visitSubcellOnBothSidesOfDir(nodeMid, polyIntersectIndex, poly.compositedBernstein, subNodes, 0, mark);
} }
} }
// launch subdivision in subcells // launch subdivision in subcells
@ -393,14 +289,231 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
subNodes[subIdx].polyFullyContainedIndices.insert(subNodes[subIdx].polyFullyContainedIndices.end(), subNodes[subIdx].polyFullyContainedIndices.insert(subNodes[subIdx].polyFullyContainedIndices.end(),
node.polyFullyContainedIndices.begin(), node.polyFullyContainedIndices.begin(),
node.polyFullyContainedIndices.end()); node.polyFullyContainedIndices.end());
build(scene, subNodes[subIdx], leaves); buildOcTree(scene, subNodes[subIdx], leaves);
}
}
/** single object quadrature */
void basicTask(const std::shared_ptr<PrimitiveDesc>& p, int q = 20, real xmin = -1, real xmax = 1)
{
real volume = 0;
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
uvector3 range = xmax - xmin;
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(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<real, 3> 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<real, 3>& x, real w) {
if (isInsideBernstein(phi, x)) volume += w * integrand(xmin + x * (xmax - xmin));
});
} else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(p)) {
const int faceCount = pt->indexInclusiveScan.size();
assert(faceCount > 1);
std::vector<tensor3> 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<real, 3> 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<real, 3>& 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<std::shared_ptr<PrimitiveDesc>>& primitives, int q = 20, real xmin = -1, real xmax = 1)
{
std::vector<SparkStack<real>*> phiStacks;
std::vector<tensor3> phis;
std::vector<SparkStack<real>*> originTensorStacks;
std::vector<tensor3> originTensors;
real volume;
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
uvector3 range = xmax - xmin;
uvector<real, 3> testX(0., 0.75, 0.2);
for (int i = 0; i < primitives.size(); i++) {
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(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<MeshDesc>(primitives[i])) {
const int faceCount = pt->indexInclusiveScan.size();
assert(faceCount > 1);
std::vector<tensor3> 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<real, 3>& 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<real, 3>& 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<real, 3>& x, real w) {
auto realX = x * range + node.min;
if (keepQuadraturePoint(scene, node, realX)) volume += w * integrand(realX);
});
}
void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives,
const uvector3& xmin,
const uvector3& xmax,
const organizer::BlobTree& blobTree)
{
OcTreeNode rootNode(xmin, xmax, blobTree);
std::vector<OcTreeNode> leaves;
std::vector<SparkStack<real>*> tensorStacks;
std::vector<CompleteTensorRep> completeTensorReps;
real volume;
auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
uvector3 range = xmax - xmin;
uvector<real, 3> testX(0., 0.75, 0.2);
for (int i = 0; i < primitives.size(); i++) {
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(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<MeshDesc>(primitives[i])) {
const int faceCount = pt->indexInclusiveScan.size();
assert(faceCount > 1);
std::vector<tensor3> 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 }; // namespace algoim::Organizer
// clang-format off // clang-format off
/** /**
最一开始的cell应当包含所有的primitive 最一开始的cell应当包含所有的primitive
func(given cell) { func(given cell) {
得到8个subcell 得到8个subcell
@ -408,7 +521,7 @@ void build(const Scene& scene, const Node& node, std::vector<Node> leaves)
mark = [0,0,0] mark = [0,0,0]
for k in 0,1,2 { for k in 0,1,2 {
测试primitive和k轴向的centerface有无交 测试primitive和k轴向的centerface有无交
if 无交 { if 无交 {0
if primitive包裹centerface { if primitive包裹centerface {
mark[k] = 2 mark[k] = 2
} else { } else {

24
algoim/organizer/primitive.hpp
View File

@ -12,11 +12,13 @@
#include "xarray.hpp" #include "xarray.hpp"
#include "binomial.hpp" #include "binomial.hpp"
#include "bernstein.hpp" #include "bernstein.hpp"
#include "blobtree.hpp"
#include <memory> #include <memory>
namespace algoim::Organizer namespace algoim::Organizer
{ {
namespace detail namespace detail
{ {
@ -589,29 +591,39 @@ void makeCylinder(xarray<real, 3>& tensor, uvector3 startPt, uvector3 endPt, rea
// TODO: // TODO:
} }
struct CompleteTensorRep {
tensor3 compositedBernstein; // 合法输入在[0,1]^3
std::vector<tensor3> originalPower; // 合法输入在场景<min, max>
};
struct Scene { struct Scene {
std::vector<tensor3> polys; std::vector<CompleteTensorRep> polys;
int* boolDescription; // blobtree organizer::BlobTree blobTree;
// uvector3 min, max; // uvector3 min, max;
}; };
const int CHILD_NUM = 8; const int CHILD_NUM = 8;
struct Node { struct OcTreeNode {
std::vector<int> polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 std::vector<int> polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中
std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 std::vector<int> polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中
// std::array<int, CHILD_NUM> children; // octree // std::array<int, CHILD_NUM> children; // octree
uvector3 min, max; uvector3 min, max;
organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递
// Node() { children.fill(-1); } // 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), : polyIntersectIndices(node.polyIntersectIndices),
polyFullyContainedIndices(node.polyFullyContainedIndices), polyFullyContainedIndices(node.polyFullyContainedIndices),
min(node.min), min(node.min),
max(node.max) max(node.max),
blobTree(node.blobTree)
{ {
} }
}; };

12
algoim/quadrature_multipoly.hpp
View File

@ -19,6 +19,7 @@
#include <cassert> #include <cassert>
#include <cmath> #include <cmath>
#include <vector> #include <vector>
#include "organizer/primitive.hpp"
#define ALGOIM_M 8 #define ALGOIM_M 8
#include <algorithm> #include <algorithm>
@ -568,6 +569,15 @@ struct ImplicitPolyQuadrature {
build(true, false); build(true, false);
} }
ImplicitPolyQuadrature(const std::vector<Organizer::CompleteTensorRep>& completeTensorRes)
{
for (const auto& ctr : completeTensorRes) {
auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray<N, ALGOIM_M>(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 // Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks
ImplicitPolyQuadrature(const xarray<real, N>& p, ImplicitPolyQuadrature(const xarray<real, N>& p,
const booluarray<N, ALGOIM_M>& pmask, const booluarray<N, ALGOIM_M>& pmask,
@ -807,7 +817,7 @@ struct ImplicitPolyQuadrature {
// Apply primary base integral // Apply primary base integral
k_active = k; k_active = k;
base.integrate(strategy, q, integrand); base.integrate(strategy, q, integrand);
// situation the // situation the
// If in aggregate mode, apply to other dimensions as well // If in aggregate mode, apply to other dimensions as well
if (type == OuterAggregate) { if (type == OuterAggregate) {
for (int i = 0; i < N - 1; ++i) { for (int i = 0; i < N - 1; ++i) {

Write Preview
Loading…
Cancel
Save