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

blobtree

master
gjj 9 months ago
parent
b0ec3bd561
commit
0747ce10bd
5 changed files with 186 additions and 105 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. 21
      algoim/organizer/blobtree.hpp
  2. 115
      algoim/organizer/organizer.hpp
  3. 35
      algoim/organizer/primitive.hpp
  4. 14
      algoim/quadrature_multipoly.hpp
  5. 106
      gjj/primitiveDebug.hpp

21
algoim/organizer/blobtree.hpp
View File

@ -37,6 +37,21 @@ public:
std::vector<unsigned int> primitiveNodeIdx; // 由primitive数组下标指向node中对应leaf的下标 std::vector<unsigned int> primitiveNodeIdx; // 由primitive数组下标指向node中对应leaf的下标
BlobTree(const BlobTree& other)
{
structure = other.structure;
primitiveNodeIdx = other.primitiveNodeIdx;
}
BlobTree& operator=(const BlobTree& other)
{
if (this != &other) {
structure = other.structure;
primitiveNodeIdx = other.primitiveNodeIdx;
}
return *this;
}
BlobTree() {} BlobTree() {}
}; };
@ -44,7 +59,11 @@ 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); // if (node.inOut != NODE_IN_OUT_UNKNOWN) {
// int aaa = 1;
// int bb = 1;
// }
assert(node.inOut == NODE_IN_OUT_UNKNOWN);
node.inOut = in ? NODE_IN : NODE_OUT; node.inOut = in ? NODE_IN : NODE_OUT;
for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) { for (unsigned int nowIdx = nodeIdx; nowIdx != rootIdx;) {
const auto& nowNode = tree.structure[nowIdx]; const auto& nowNode = tree.structure[nowIdx];

115
algoim/organizer/organizer.hpp
View File

@ -93,41 +93,21 @@ int findFirst(const uvector<int, N>& v, int val, int startIdx = 0)
return -1; return -1;
} }
std::pair<uvector<real, 3>, uvector<real, 3>> getOneEightCellAABB(const uvector3& min, AABB getOneEightCellAABB(const AABB& fatherAABB, const uvector<int, 3> side, int negativeRep = -1)
const uvector3& max,
const uvector<int, 3> side,
int negativeRep = -1)
{ {
std::pair<uvector<real, 3>, uvector<real, 3>> res = {min, max}; AABB res = fatherAABB;
uvector3 mid = (min + max) / 2; uvector3 mid = res.center();
for (int i = 0; i < 3; ++i) { for (int i = 0; i < 3; ++i) {
if (side(i) == negativeRep) { if (side(i) == negativeRep) {
res.second(i) = mid(i); res.max(i) = mid(i);
} else { } else {
res.first(i) = mid(i); res.min(i) = mid(i);
} }
} }
return res; return res;
} }
} // namespace detail } // namespace detail
// TODO: delete it
bool isPointInside(const std::vector<tensor3>& originTensors, const uvector3& originPt)
{
for (auto& t : originTensors) {
if (evalPower(t, originPt) >= 0) { return false; }
}
return true;
}
// bool isPointInside(const std::vector<CompleteTensorRep>& completeTensorReps, const uvector3& originPt)
// {
// for (auto& completeTensorRep : completeTensorReps) {
// if (isPointInside(completeTensorRep.originalPower, originPt)) { return false; }
// }
// return true;
// }
// 这里blobTree是拷贝传参 // 这里blobTree是拷贝传参
bool keepQuadraturePoint(const Scene& scene, bool keepQuadraturePoint(const Scene& scene,
organizer::BlobTree blobTree, organizer::BlobTree blobTree,
@ -138,7 +118,7 @@ bool keepQuadraturePoint(const Scene& scene,
const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices;
std::vector<bool> primitiveInOuts(polyIntersectIndices.size()); std::vector<bool> primitiveInOuts(polyIntersectIndices.size());
for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) { for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) {
primitiveInOuts[i] = isPointInside(scene.polys[polyIntersectIndices[i]].rawPower, originPt); primitiveInOuts[i] = isInsidePowers(scene.polys[polyIntersectIndices[i]].rawPower, originPt);
} }
int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts);
assert(res == organizer::NODE_IN || res == organizer::NODE_OUT); assert(res == organizer::NODE_IN || res == organizer::NODE_OUT);
@ -149,7 +129,6 @@ bool keepQuadraturePoint(const Scene& scene,
// BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {}; // BasicTask(std::vector<std::shared_ptr<PrimitiveDesc>> ps) {};
const std::array<int, 2> sides = {-1, 1}; const std::array<int, 2> sides = {-1, 1};
struct BasicTaskRes { struct BasicTaskRes {
@ -174,11 +153,11 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid,
algoim_spark_alloc(real, halfCellPoly); algoim_spark_alloc(real, halfCellPoly);
int subIdx = detail::binaryToDecimal(mark, -1); int subIdx = detail::binaryToDecimal(mark, -1);
auto& subNode = subNodes[subIdx]; auto& subNode = subNodes[subIdx];
bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达 bernstein::deCasteljau(poly, subNode.aabb.min, subNode.aabb.max, halfCellPoly); // 求1/8空间下的表达
if (bernstein::uniformSign(halfCellPoly) != 1) { if (bernstein::uniformSign(halfCellPoly) != 1) {
subNode.polyIntersectIndices.emplace_back(polyIntersectIndex); subNode.polyIntersectIndices.emplace_back(polyIntersectIndex);
} else { } else {
organizer::traverse(subNode.blobTree, polyIntersectIndex, organizer::NODE_OUT); organizer::traverse(subNode.blobTree, polyIntersectIndex, false);
} }
} else { } else {
for (auto side : sides) { for (auto side : sides) {
@ -199,19 +178,15 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
leaves.emplace_back(node); leaves.emplace_back(node);
return; return;
} }
const uvector3& nowNodeMin = node.min;
const uvector3& nowNodeMax = node.max;
std::vector<OcTreeNode> subNodes(CHILD_NUM); std::vector<OcTreeNode> subNodes(CHILD_NUM);
// std::array<OcTreeNode, 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 = detail::getOneEightCellAABB(node.min, node.max, j(), 0); subNodes[subIdx].aabb = detail::getOneEightCellAABB(node.aabb, j(), 0);
subNodes[subIdx].min = nodeAABB.first;
subNodes[subIdx].max = nodeAABB.second;
subNodes[subIdx].blobTree = node.blobTree; subNodes[subIdx].blobTree = node.blobTree;
} }
uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; uvector3 nodeMid = node.aabb.center();
for (int i = 0; i < polyIntersectIndices.size(); ++i) { for (int i = 0; i < polyIntersectIndices.size(); ++i) {
const int polyIntersectIndex = polyIntersectIndices[i]; const int polyIntersectIndex = polyIntersectIndices[i];
const auto& poly = scene.polys[polyIntersectIndex]; const auto& poly = scene.polys[polyIntersectIndex];
@ -229,18 +204,18 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
} else if (signOnHalfPlane == 1) { } else if (signOnHalfPlane == 1) {
// primitive intersects either side or both sides of the centerface // primitive intersects either side or both sides of the centerface
// deCasteljau to transformation // deCasteljau to transformation
uvector3 halfCellMin = nowNodeMin, halfCellMax = nowNodeMax; AABB halfCell = node.aabb;
halfCellMax(faceAxis) = nodeMid(faceAxis); halfCell.max(faceAxis) = nodeMid(faceAxis);
tensor3 halfCellPoly(nullptr, poly.compositedBernstein.ext()); tensor3 halfCellPoly(nullptr, poly.compositedBernstein.ext());
algoim_spark_alloc(real, halfCellPoly); algoim_spark_alloc(real, halfCellPoly);
bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly); bernstein::deCasteljau(poly.compositedBernstein, halfCell.min, halfCell.max, halfCellPoly);
if (bernstein::uniformSign(halfCellPoly) != 1) { if (bernstein::uniformSign(halfCellPoly) != 1) {
// 负空间有 // 负空间有
mark(faceAxis) = -1; mark(faceAxis) = -1;
} }
halfCellMax(faceAxis) = nowNodeMax(faceAxis); halfCell.max(faceAxis) = node.aabb.max(faceAxis);
halfCellMin(faceAxis) = nodeMid(faceAxis); halfCell.min(faceAxis) = nodeMid(faceAxis);
bernstein::deCasteljau(poly.compositedBernstein, halfCellMin, halfCellMax, halfCellPoly); bernstein::deCasteljau(poly.compositedBernstein, halfCell.min, halfCell.max, halfCellPoly);
if (bernstein::uniformSign(halfCellPoly) != 1) { if (bernstein::uniformSign(halfCellPoly) != 1) {
// 正空间有 // 正空间有
mark(faceAxis) += 1; // 当正负空间都有，记0 mark(faceAxis) += 1; // 当正负空间都有，记0
@ -253,10 +228,10 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
// 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.compositedBernstein.ext()); tensor3 subPoly(nullptr, poly.compositedBernstein.ext());
bernstein::deCasteljau(poly.compositedBernstein, subNodes[i].min, subNodes[i].max, subPoly); bernstein::deCasteljau(poly.compositedBernstein, subNodes[i].aabb.min, subNodes[i].aabb.max, subPoly);
if (bernstein::uniformSign(subPoly) == -1) { if (bernstein::uniformSign(subPoly) == -1) {
// subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex); // subNodes[i].polyFullyContainedIndices.emplace_back(polyIntersectIndex);
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_IN); organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, true);
} else { } else {
subNodes[i].polyIntersectIndices.emplace_back(polyIntersectIndex); subNodes[i].polyIntersectIndices.emplace_back(polyIntersectIndex);
} }
@ -277,7 +252,7 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
if (subIdx == markIdx) if (subIdx == markIdx)
subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
else else
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
} }
} else if (zeroCount == 1) { } else if (zeroCount == 1) {
// poly related to 2 subcells // poly related to 2 subcells
@ -290,7 +265,7 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
subNodes[sideIdx2].polyIntersectIndices.emplace_back(polyIntersectIndex); subNodes[sideIdx2].polyIntersectIndices.emplace_back(polyIntersectIndex);
for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) {
if (subIdx != sideIdx1 && subIdx != sideIdx2) { if (subIdx != sideIdx1 && subIdx != sideIdx2) {
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
} }
} }
} else { } else {
@ -311,39 +286,44 @@ void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTre
void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves) void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves)
{ {
int a = 0;
const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices; const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices;
if (polyIntersectIndices.size() <= 4) { if (polyIntersectIndices.size() <= 3) {
leaves.emplace_back(node); leaves.emplace_back(node);
return; return;
} }
const uvector3& nowNodeMin = node.min;
const uvector3& nowNodeMax = node.max;
// std::array<OcTreeNode, CHILD_NUM> subNodes; // std::array<OcTreeNode, CHILD_NUM> subNodes;
std::vector<OcTreeNode> subNodes(CHILD_NUM); std::vector<OcTreeNode> subNodes(CHILD_NUM);
// 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 = detail::getOneEightCellAABB(node.min, node.max, j(), 0); subNodes[subIdx].aabb = detail::getOneEightCellAABB(node.aabb, j(), 0);
subNodes[subIdx].min = nodeAABB.first;
subNodes[subIdx].max = nodeAABB.second;
subNodes[subIdx].blobTree = node.blobTree; subNodes[subIdx].blobTree = node.blobTree;
subNodes[subIdx].polyIntersectIndices = std::vector<int>(); subNodes[subIdx].polyIntersectIndices = std::vector<int>();
} }
uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; uvector3 nodeMid = node.aabb.center();
for (int i = 0; i < polyIntersectIndices.size(); ++i) { for (int i = 0; i < polyIntersectIndices.size(); ++i) {
const int polyIntersectIndex = polyIntersectIndices[i]; const int polyIntersectIndex = polyIntersectIndices[i];
const auto& poly = scene.polys[polyIntersectIndex]; const auto& poly = scene.polys[polyIntersectIndex];
subIdx = 0; subIdx = 0;
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) {
if (!poly.aabb.intersect(subNodes[subIdx].aabb)) {
// out of the subcell
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
continue;
}
tensor3 subcellPoly(nullptr, poly.compositedBernstein.ext()); tensor3 subcellPoly(nullptr, poly.compositedBernstein.ext());
algoim_spark_alloc(real, subcellPoly); algoim_spark_alloc(real, subcellPoly);
bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].min, subNodes[subIdx].max, subcellPoly); bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].aabb.min, subNodes[subIdx].aabb.max, subcellPoly);
int sign = bernstein::uniformSign(subcellPoly); int sign = bernstein::uniformSign(subcellPoly);
if (sign == 1) { if (sign == 1) {
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT); organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
} else if (sign == -1) { } else if (sign == -1) {
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_IN); // 采样一个点，带入所有原始primitive看是否是冗余部分
uvector3 sampleX = subNodes[subIdx].aabb.center() * scene.boundary.size() + scene.boundary.min;
if (isInsidePowers(poly.rawPower, sampleX))
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, true);
else
organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, false);
} else { } else {
subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex);
} }
@ -479,7 +459,7 @@ void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, in
ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) { ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
auto realX = x * range + xmin; auto realX = x * range + xmin;
if (isPointInside(originTensors, realX)) volume += w * integrand(realX); if (isInsidePowers(originTensors, realX)) volume += w * integrand(realX);
}); });
volume *= pow(xmax - xmin, 3); volume *= pow(xmax - xmin, 3);
std::cout << "Volume xxx: " << volume << std::endl; std::cout << "Volume xxx: " << volume << std::endl;
@ -489,18 +469,18 @@ void basicTask(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitives, in
for (auto& p : originTensorStacks) delete p; for (auto& p : originTensorStacks) delete p;
}; };
BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 20) BasicTaskRes basicTask(const Scene& scene, const organizer::BlobTree& blobTree, const OcTreeNode& node, int q = 10)
{ {
auto integrand = [](const uvector<real, 3>& x) { return 1.0; }; auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
real volume = 0., surf = 0.; real volume = 0., surf = 0.;
auto range = node.max - node.min; auto range = node.aabb.size();
ImplicitPolyQuadrature<3> ipquad(scene.polys); ImplicitPolyQuadrature<3> ipquad(scene.polys, node.polyIntersectIndices);
ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) { ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
auto realX = x * range + node.min; auto realX = x * range + node.aabb.min;
if (keepQuadraturePoint(scene, blobTree, node, realX)) { volume += w * integrand(realX); } if (keepQuadraturePoint(scene, blobTree, node, realX)) { volume += w * integrand(realX); }
}); });
ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, 3>& x, real w, const uvector<real, 3>& wn) { ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, 3>& x, real w, const uvector<real, 3>& wn) {
surf += w * integrand(x * range + node.min); surf += w * integrand(x * range + node.aabb.min);
}); });
return {volume, surf}; return {volume, surf};
} }
@ -534,6 +514,7 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
AABB aabb; AABB aabb;
makeSphere(*pt, originTensor, aabb); makeSphere(*pt, originTensor, aabb);
detail::powerTransformation(range, xmin, originTensor, transformedTensor); detail::powerTransformation(range, xmin, originTensor, transformedTensor);
aabb.normalize(range, xmin);
real testEvaOri = evalPower(originTensor, testX); real testEvaOri = evalPower(originTensor, testX);
@ -553,7 +534,9 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
algoim_spark_alloc(real, compositeTensor); algoim_spark_alloc(real, compositeTensor);
AABB aabb; AABB aabb;
makeMesh(*pt, compositeTensor, planeTensors, aabb); makeMesh(*pt, compositeTensor, planeTensors, aabb);
detail::powerTransformation(range, xmin, compositeTensor); detail::powerTransformation(range, xmin, compositeTensor);
aabb.normalize(range, xmin);
real testEvalPower = evalPower(compositeTensor, testX); real testEvalPower = evalPower(compositeTensor, testX);
@ -576,6 +559,8 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
AABB aabb; AABB aabb;
makeCylinder(*pt, compositeTensor, rawTensors, aabb); makeCylinder(*pt, compositeTensor, rawTensors, aabb);
detail::powerTransformation(range, xmin, compositeTensor); detail::powerTransformation(range, xmin, compositeTensor);
aabb.normalize(range, xmin);
real testEvalPower = evalPower(compositeTensor, testX); real testEvalPower = evalPower(compositeTensor, testX);
detail::power2BernsteinTensor(compositeTensor, phi); detail::power2BernsteinTensor(compositeTensor, phi);
@ -595,6 +580,8 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
AABB aabb; AABB aabb;
makeCone(*pt, compositeTensor, rawTensors, aabb); makeCone(*pt, compositeTensor, rawTensors, aabb);
detail::powerTransformation(range, xmin, compositeTensor); detail::powerTransformation(range, xmin, compositeTensor);
aabb.normalize(range, xmin);
real testEvalPower = evalPower(compositeTensor, testX); real testEvalPower = evalPower(compositeTensor, testX);
detail::power2BernsteinTensor(compositeTensor, phi); detail::power2BernsteinTensor(compositeTensor, phi);
@ -605,14 +592,14 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv
} }
} }
Scene scene{completeTensorReps, blobTree}; Scene scene{completeTensorReps, AABB(xmin, xmax)};
OcTreeNode rootNode(0, 1, blobTree); OcTreeNode rootNode(0, 1, blobTree);
for (int i = 0; i < completeTensorReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } for (int i = 0; i < completeTensorReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); }
buildOcTreeV0(scene, rootNode, leaves); buildOcTreeV0(scene, rootNode, leaves);
for (const auto& leaf : leaves) { for (const auto& leaf : leaves) {
auto basicRes = basicTask(scene, blobTree, leaf, q); auto basicRes = basicTask(scene, blobTree, leaf, q);
volume += basicRes.volume * prod(leaf.max - leaf.min); volume += basicRes.volume * prod(leaf.aabb.max - leaf.aabb.min);
} }
volume *= prod(xmax - xmin); volume *= prod(xmax - xmin);

35
algoim/organizer/primitive.hpp
View File

@ -171,6 +171,11 @@ real evalBernstein(const xarray<real, N>& phi, const uvector<real, N>& x)
return bernstein::evalBernsteinPoly(phi, x); return bernstein::evalBernsteinPoly(phi, x);
} }
void inverseTensor(tensor3& tensor)
{
for (auto i = tensor.loop(); ~i; ++i) { tensor.l(i) = -tensor.l(i); }
}
// class PowerTensor : public Primitive // class PowerTensor : public Primitive
// { // {
// public: // public:
@ -603,7 +608,8 @@ public:
} }
}; };
struct AABB { class AABB
{
public: public:
uvector3 min = uvector3(std::numeric_limits<real>::max()); uvector3 min = uvector3(std::numeric_limits<real>::max());
uvector3 max = uvector3(std::numeric_limits<real>::lowest()); uvector3 max = uvector3(std::numeric_limits<real>::lowest());
@ -630,6 +636,20 @@ public:
uvector3 center() const { return (min + max) / 2; } uvector3 center() const { return (min + max) / 2; }
uvector3 size() const { return max - min; } uvector3 size() const { return max - min; }
bool intersect(const AABB& aabb) const
{
for (int i = 0; i < 3; ++i) {
if (min(i) > aabb.max(i) || max(i) < aabb.min(i)) { return false; }
}
return true;
}
void normalize(const uvector3& scale, const uvector3& boundaryMin)
{
min = (min - boundaryMin) / scale;
max = (max - boundaryMin) / scale;
}
}; };
void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<real, 3>>& planeTensors, AABB& aabb) void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<real, 3>>& planeTensors, AABB& aabb)
@ -663,7 +683,7 @@ void makeMesh(const MeshDesc& mesh, xarray<real, 3>& tensor, std::vector<xarray<
} }
} }
detail::compositePower(planeTensors, 0, 0, 1, tensor); detail::compositePower(planeTensors, 0, 0, 1, tensor);
if (mesh.indexInclusiveScan.size() % 2 == 0) { inverseTensor(tensor); }
// AABB // AABB
for (const auto& v : mesh.vertices) { aabb.extend(v); } for (const auto& v : mesh.vertices) { aabb.extend(v); }
}; };
@ -785,7 +805,8 @@ struct CompleteTensorRep {
struct Scene { struct Scene {
std::vector<CompleteTensorRep> polys; std::vector<CompleteTensorRep> polys;
organizer::BlobTree blobTree; // organizer::BlobTree blobTree;
AABB boundary;
// uvector3 min, max; // uvector3 min, max;
}; };
@ -795,21 +816,21 @@ 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; AABB aabb;
organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递 organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递
// Node() { children.fill(-1); } // Node() { children.fill(-1); }
OcTreeNode() = default; OcTreeNode() = default;
OcTreeNode(const uvector3& min_, const uvector3& max_, const organizer::BlobTree& blobTree_) OcTreeNode(const uvector3& min_, const uvector3& max_, const organizer::BlobTree& blobTree_)
: min(min_), max(max_), blobTree(blobTree_) {}; : aabb(min_, max_), blobTree(blobTree_) {};
OcTreeNode(const AABB& aabb_, const organizer::BlobTree& blobTree_) : aabb(aabb_), blobTree(blobTree_) {};
// 深拷贝 // 深拷贝
OcTreeNode(const OcTreeNode& node) OcTreeNode(const OcTreeNode& node)
: polyIntersectIndices(node.polyIntersectIndices), : polyIntersectIndices(node.polyIntersectIndices),
// polyFullyContainedIndices(node.polyFullyContainedIndices), // polyFullyContainedIndices(node.polyFullyContainedIndices),
min(node.min), aabb(node.aabb),
max(node.max),
blobTree(node.blobTree) blobTree(node.blobTree)
{ {
int a = 1; int a = 1;

14
algoim/quadrature_multipoly.hpp
View File

@ -569,15 +569,25 @@ struct ImplicitPolyQuadrature {
build(true, false); build(true, false);
} }
ImplicitPolyQuadrature(const std::vector<Organizer::CompleteTensorRep>& completeTensorRes) ImplicitPolyQuadrature(const std::vector<Organizer::CompleteTensorRep>& completeTensorReps)
{ {
for (const auto& ctr : completeTensorRes) { for (const auto& ctr : completeTensorReps) {
auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray<N, ALGOIM_M>(true)); auto mask = detail::nonzeroMask(ctr.compositedBernstein, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(ctr.compositedBernstein, mask); if (!detail::maskEmpty(mask)) phi.push_back(ctr.compositedBernstein, mask);
} }
build(true, false); build(true, false);
} }
ImplicitPolyQuadrature(const std::vector<Organizer::CompleteTensorRep>& completeTensorReps,
const std::vector<int>& polyIndices)
{
for (auto i : polyIndices) {
auto mask = detail::nonzeroMask(completeTensorReps[i].compositedBernstein, booluarray<N, ALGOIM_M>(true));
if (!detail::maskEmpty(mask)) phi.push_back(completeTensorReps[i].compositedBernstein, mask);
}
build(true, false);
}
// Build quadrature hierarchy for a given domain implicitly defined by two polynomials with user-defined masks // 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,

106
gjj/primitiveDebug.hpp
View File

@ -15,6 +15,7 @@
#include "binomial.hpp" #include "binomial.hpp"
#include "real.hpp" #include "real.hpp"
#include "sparkstack.hpp"
#include "uvector.hpp" #include "uvector.hpp"
#include "vector" #include "vector"
#include "xarray.hpp" #include "xarray.hpp"
@ -122,13 +123,28 @@ void casePolyhedronSphere()
basicTask(primitiveDescriptions); basicTask(primitiveDescriptions);
} }
void caseTwoCube()
{
auto phi0 = std::make_shared<CuboidDesc>(CuboidDesc(0.2, 0.6));
auto phi1 = std::make_shared<CuboidDesc>(CuboidDesc(-0.2, 0.6));
basicTask({phi0, phi1});
}
void caseCubeSphere()
{
auto phi0 = std::make_shared<CuboidDesc>(CuboidDesc(0.2, 0.6));
auto phi1 = std::make_shared<SphereDesc>(SphereDesc(0.6, -0.2, 1.));
basicTask({phi0, phi1});
}
void caseScene() void caseScene()
{ {
const int PRIMITIVE_CNT = 6; const int PRIMITIVE_CNT = 6;
std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions(PRIMITIVE_CNT); std::vector<std::shared_ptr<PrimitiveDesc>> primitiveDescriptions(PRIMITIVE_CNT);
primitiveDescriptions[2] = std::make_shared<CuboidDesc>(CuboidDesc(0., 1.6)); primitiveDescriptions[0] = std::make_shared<CuboidDesc>(CuboidDesc(0., 1.6));
primitiveDescriptions[1] = std::make_shared<CuboidDesc>(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); // primitiveDescriptions[1] = std::make_shared<CuboidDesc>(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.));
primitiveDescriptions[0] = std::make_shared<SphereDesc>(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.)); primitiveDescriptions[1] = std::make_shared<SphereDesc>(SphereDesc(1., uvector3(0.6, 0.6, -0.6), 1.));
primitiveDescriptions[2] = std::make_shared<SphereDesc>(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.));
primitiveDescriptions[3] = std::make_shared<SphereDesc>(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.)); primitiveDescriptions[3] = std::make_shared<SphereDesc>(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.));
std::vector<uvector3> pyramidBottomVertices = { std::vector<uvector3> pyramidBottomVertices = {
uvector3{-1, -1, 0}, uvector3{-1, -1, 0},
@ -166,7 +182,8 @@ void testDeCasteljau()
auto phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.8, uvector3(0), 1.)); auto phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.8, uvector3(0), 1.));
tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3); tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3);
algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein); algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein);
makeSphere(*phiDesc, tensor); AABB aabb;
makeSphere(*phiDesc, tensor, aabb);
uvector3 xmax = 1, xmin = -1, range = xmax - xmin; uvector3 xmax = 1, xmin = -1, range = xmax - xmin;
Organizer::detail::powerTransformation(range, xmin, tensor, tensor01); Organizer::detail::powerTransformation(range, xmin, tensor, tensor01);
Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein); Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein);
@ -193,7 +210,8 @@ void testSubDivideWithDeCasteljau()
auto phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.7, uvector3(0.8), 1.)); auto phiDesc = std::make_shared<SphereDesc>(SphereDesc(0.7, uvector3(0.8), 1.));
tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3); tensor3 tensor(nullptr, 3), tensor01(nullptr, 3), tensorBernstein(nullptr, 3), transformedTensorBernstein(nullptr, 3);
algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein); algoim_spark_alloc(real, tensor, tensor01, tensorBernstein, transformedTensorBernstein);
makeSphere(*phiDesc, tensor); AABB aabb;
makeSphere(*phiDesc, tensor, aabb);
uvector3 xmin = uvector3(-1., -1.3, -1.6), xmax = uvector3(1.6, 1.6, 2.3), range = xmax - xmin; uvector3 xmin = uvector3(-1., -1.3, -1.6), xmax = uvector3(1.6, 1.6, 2.3), range = xmax - xmin;
Organizer::detail::powerTransformation(range, xmin, tensor, tensor01); Organizer::detail::powerTransformation(range, xmin, tensor, tensor01);
Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein); Organizer::detail::power2BernsteinTensor(tensor01, tensorBernstein);
@ -202,27 +220,44 @@ void testSubDivideWithDeCasteljau()
// bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001), // bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001),
// transformedTensorBernstein); // transformedTensorBernstein);
for (MultiLoop<3> j(0, 2); ~j; ++j) { for (MultiLoop<3> j(0, 2); ~j; ++j) {
if (all(j() == uvector3i(1, 0, 0))) { auto aabbSub = Organizer::detail::getOneEightCellAABB(AABB{xmin, xmax}, j(), 0);
int aaa = 1; bernstein::deCasteljau(tensorBernstein, aabbSub.min, aabbSub.max, transformedTensorBernstein);
int bbb = 1; auto transformedTensorVec = xarray2StdVector(transformedTensorBernstein);
} int sign = bernstein::uniformSign(transformedTensorBernstein);
auto aabb = Organizer::detail::getOneEightCellAABB(xmin, xmax, j(), 0);
bernstein::deCasteljau(tensorBernstein, aabb.first, aabb.second, transformedTensorBernstein);
int sign = bernstein::uniformSign(transformedTensorBernstein);
std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(2) << "), sign: " << sign << std::endl; std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(2) << "), sign: " << sign << std::endl;
if (all(j() == uvector3i(0, 0, 1))) {
for (MultiLoop<3> k(0, 2); ~k; ++k) {
auto aabbSubSub = Organizer::detail::getOneEightCellAABB(aabbSub, k(), 0);
bernstein::deCasteljau(tensorBernstein, aabbSubSub.min, aabbSubSub.max, transformedTensorBernstein);
int sign = bernstein::uniformSign(transformedTensorBernstein);
std::cout << "k(): (" << k(0) << ", " << k(1) << ", " << k(2) << "), sign: " << sign << std::endl;
}
}
} }
uvector3 testX = uvector3(0.5);
bernstein::deCasteljau(tensorBernstein, uvector3(-1., -1.3, 0.35), uvector3(0.3, 0.15, 2.3), transformedTensorBernstein);
int sign = bernstein::uniformSign(transformedTensorBernstein);
real evalX = evalPower(tensor, testX);
uvector3 testX = uvector3(0.5);
real evalX = evalPower(tensor, testX);
std::cout << "evalX original rep = " << evalX << std::endl; std::cout << "evalX original rep = " << evalX << std::endl;
evalX = evalPower(tensor01, testX);
std::cout << "evalX power within 0-1 = " << evalX << std::endl;
evalX = bernstein::evalBernsteinPoly(tensorBernstein, testX); evalX = bernstein::evalBernsteinPoly(tensorBernstein, testX);
std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl; std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl;
evalX = bernstein::evalBernsteinPoly(transformedTensorBernstein, testX);
std::cout << "evalX bernstein after Decasteljau = " << evalX << std::endl; // 切比雪夫
std::cout << "sign = " << sign << std::endl; xmin = uvector3(-1., -1.3, -1.6), xmax = uvector3(1.6, 1.6, 2.3);
auto phi = [](const uvector<real, 3> &x) {
return (x(0) - 0.8) * (x(0) - 0.8) + (x(1) - 0.8) * (x(1) - 0.8) + (x(2) - 0.8) * (x(2) - 0.8) - 0.1 * 0.1;
};
xarray<real, 3> phipoly(nullptr, 3), subPhiPoly(nullptr, 3);
algoim_spark_alloc(real, phipoly, subPhiPoly);
bernstein::bernsteinInterpolate<3>([&](const uvector<real, 3> &x) { return phi(xmin + x * (xmax - xmin)); }, phipoly);
xmin = 0, xmax = 1;
for (MultiLoop<3> j(0, 2); ~j; ++j) {
auto aabbSub = Organizer::detail::getOneEightCellAABB(AABB{xmin, xmax}, j(), 0);
bernstein::deCasteljau(phipoly, aabbSub.min, aabbSub.max, subPhiPoly);
int sign = bernstein::uniformSign(subPhiPoly);
std::cout << "j(): (" << j(0) << ", " << j(1) << ", " << j(2) << "), sign: " << sign << std::endl;
}
} }
void testBlob() void testBlob()
@ -238,16 +273,23 @@ void testBlob()
std::cout << sizeof(int) << std::endl; std::cout << sizeof(int) << std::endl;
} }
void vectorDebug() void testTensorInverse()
{ {
struct A { auto phiDesc = std::make_shared<CuboidDesc>(uvector3(0.2, 0.3, 0.5), uvector3(0.8, 1.1, 0.9));
std::vector<real> data; const int faceCnt = 6;
A() = default; std::vector<tensor3> tensors(faceCnt, tensor3(nullptr, 2));
}; std::vector<SparkStack<real> *> sparkStackPtrs(faceCnt);
tensor3 tensor(nullptr, 1 + faceCnt);
A a; algoimSparkAllocHeapVector(sparkStackPtrs, tensors);
int s = a.data.size(); algoim_spark_alloc(real, tensor);
std::cout << "size: " << s << std::endl; AABB aabb;
makeMesh(*phiDesc, tensor, tensors, aabb);
uvector3 testX(0.5, 0.7, 0.2);
real evalX = evalPower(tensor, testX);
std::cout << "before inverse, evalX = " << evalX << std::endl;
inverseTensor(tensor);
evalX = evalPower(tensor, testX);
std::cout << "after inverse, evalX = " << evalX << std::endl;
} }
void testPrimitive() void testPrimitive()
@ -256,6 +298,8 @@ void testPrimitive()
// casePolyhedronSphere(); // casePolyhedronSphere();
// testSubDivideWithDeCasteljau(); // testSubDivideWithDeCasteljau();
// testBlob(); // testBlob();
caseScene(); // caseScene();
// vectorDebug(); caseCubeSphere();
// caseTwoCube();
// testTensorInverse();
} }
Write Preview
Loading…
Cancel
Save