|
|
@ -37,7 +37,7 @@ void restrictToInnerFace(const tensor3& phi, int k, real facePos, tensor2& res) |
|
|
|
{ |
|
|
|
assert(0 <= k && k < 3); |
|
|
|
assert(all(res.ext() == remove_component(phi.ext(), k))); |
|
|
|
assert(0 < facePos && facePos < 1); |
|
|
|
// assert(0 < facePos && facePos < 1);
|
|
|
|
real* basisAlongK; |
|
|
|
int P = phi.ext(k); |
|
|
|
algoim_spark_alloc(real, &basisAlongK, P); |
|
|
@ -93,13 +93,13 @@ int findFirst(const uvector<int, N>& v, int val, int startIdx = 0) |
|
|
|
return -1; |
|
|
|
} |
|
|
|
|
|
|
|
std::pair<uvector<int, 3>, uvector<int, 3>> getOneEightCellAABB(const uvector3& min, |
|
|
|
const uvector3& max, |
|
|
|
const uvector<int, 3> side, |
|
|
|
int negativeRep = -1) |
|
|
|
std::pair<uvector<real, 3>, uvector<real, 3>> getOneEightCellAABB(const uvector3& min, |
|
|
|
const uvector3& max, |
|
|
|
const uvector<int, 3> side, |
|
|
|
int negativeRep = -1) |
|
|
|
{ |
|
|
|
std::pair<uvector<int, 3>, uvector<int, 3>> res = {min, max}; |
|
|
|
uvector3 mid = (min + max) / 2; |
|
|
|
std::pair<uvector<real, 3>, uvector<real, 3>> res = {min, max}; |
|
|
|
uvector3 mid = (min + max) / 2; |
|
|
|
for (int i = 0; i < 3; ++i) { |
|
|
|
if (side(i) == negativeRep) { |
|
|
|
res.second(i) = mid(i); |
|
|
@ -138,15 +138,11 @@ bool keepQuadraturePoint(const Scene& scene, |
|
|
|
const auto& polyIntersectIndices = ocTreeNode.polyIntersectIndices; |
|
|
|
std::vector<bool> primitiveInOuts(polyIntersectIndices.size()); |
|
|
|
for (int i = 0; i < ocTreeNode.polyIntersectIndices.size(); ++i) { |
|
|
|
primitiveInOuts[i] = isPointInside(scene.polys[polyIntersectIndices[i]].originalPower, originPt); |
|
|
|
primitiveInOuts[i] = isPointInside(scene.polys[polyIntersectIndices[i]].rawPower, originPt); |
|
|
|
} |
|
|
|
int res = organizer::traverse(blobTree, ocTreeNode.polyIntersectIndices, primitiveInOuts); |
|
|
|
assert(res == organizer::NODE_IN || res == organizer::NODE_OUT); |
|
|
|
if (res == organizer::NODE_OUT) { |
|
|
|
return false; |
|
|
|
} else { |
|
|
|
return true; |
|
|
|
} |
|
|
|
return res == organizer::NODE_IN; |
|
|
|
} |
|
|
|
|
|
|
|
// std::vector<std::shared_ptr<PrimitiveDesc>> primitives;
|
|
|
@ -167,15 +163,16 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, |
|
|
|
const int polyIntersectIndex, |
|
|
|
const tensor3& poly, |
|
|
|
|
|
|
|
std::array<OcTreeNode, CHILD_NUM>& subNodes, |
|
|
|
int startIdxToCheck, |
|
|
|
uvector<int, 3>& mark) |
|
|
|
// std::array<OcTreeNode, CHILD_NUM>& subNodes,
|
|
|
|
std::vector<OcTreeNode>& subNodes, |
|
|
|
int startIdxToCheck, |
|
|
|
uvector<int, 3>& mark) |
|
|
|
{ |
|
|
|
int zeroIdx = detail::findFirst(mark, 0, startIdxToCheck); |
|
|
|
if (zeroIdx == -1) { |
|
|
|
tensor3 halfCellPoly(nullptr, poly.ext()); |
|
|
|
algoim_spark_alloc(real, halfCellPoly); |
|
|
|
int subIdx = detail::binaryToDecimal(mark); |
|
|
|
int subIdx = detail::binaryToDecimal(mark, -1); |
|
|
|
auto& subNode = subNodes[subIdx]; |
|
|
|
bernstein::deCasteljau(poly, subNode.min, subNode.max, halfCellPoly); // 求1/8空间下的表达
|
|
|
|
if (bernstein::uniformSign(halfCellPoly) != 1) { |
|
|
@ -195,18 +192,19 @@ void visitSubcellOnBothSidesOfDir(const uvector3& nodeMid, |
|
|
|
// 叶节点单独存在leaves中,是因为最后只需要叶节点信息,没有自顶向下/自底向上遍历
|
|
|
|
// 所以树结构是不需要的记录的
|
|
|
|
// void buildOcTree(const Scene& scene, std::vector<Node>& intermediateNodes, std::vector<Node> leaves, int nowNodeIdx)
|
|
|
|
void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves) |
|
|
|
void buildOcTreeV1(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves) |
|
|
|
{ |
|
|
|
const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices; |
|
|
|
if (polyIntersectIndices.size() <= 4) { |
|
|
|
leaves.emplace_back(node); |
|
|
|
return; |
|
|
|
} |
|
|
|
const uvector3& nowNodeMin = node.min; |
|
|
|
const uvector3& nowNodeMax = node.max; |
|
|
|
std::array<OcTreeNode, CHILD_NUM> subNodes; |
|
|
|
const uvector3& nowNodeMin = node.min; |
|
|
|
const uvector3& nowNodeMax = node.max; |
|
|
|
std::vector<OcTreeNode> subNodes(CHILD_NUM); |
|
|
|
// std::array<OcTreeNode, CHILD_NUM> subNodes;
|
|
|
|
// intermediateNodes.resize(lastIdx + 8);
|
|
|
|
int subIdx = 0; |
|
|
|
int subIdx = 0; |
|
|
|
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { |
|
|
|
auto nodeAABB = detail::getOneEightCellAABB(node.min, node.max, j(), 0); |
|
|
|
subNodes[subIdx].min = nodeAABB.first; |
|
|
@ -274,7 +272,7 @@ void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeN |
|
|
|
if (zeroCount == 0) { |
|
|
|
// mark has -1 or 1 only
|
|
|
|
real nodeMidVal = bernstein::evalBernsteinPoly(poly.compositedBernstein, nodeMid); |
|
|
|
int markIdx = detail::binaryToDecimal(mark); |
|
|
|
int markIdx = detail::binaryToDecimal(mark, -1); |
|
|
|
for (int subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { |
|
|
|
if (subIdx == markIdx) |
|
|
|
subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); |
|
|
@ -307,10 +305,55 @@ void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeN |
|
|
|
// subNodes[subIdx].polyFullyContainedIndices.insert(subNodes[subIdx].polyFullyContainedIndices.end(),
|
|
|
|
// node.polyFullyContainedIndices.begin(),
|
|
|
|
// node.polyFullyContainedIndices.end());
|
|
|
|
buildOcTree(scene, subNodes[subIdx], leaves); |
|
|
|
buildOcTreeV1(scene, subNodes[subIdx], leaves); |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
void buildOcTreeV0(const Scene& scene, const OcTreeNode& node, std::vector<OcTreeNode>& leaves) |
|
|
|
{ |
|
|
|
const std::vector<int>& polyIntersectIndices = node.polyIntersectIndices; |
|
|
|
if (polyIntersectIndices.size() <= 4) { |
|
|
|
leaves.emplace_back(node); |
|
|
|
return; |
|
|
|
} |
|
|
|
const uvector3& nowNodeMin = node.min; |
|
|
|
const uvector3& nowNodeMax = node.max; |
|
|
|
// std::array<OcTreeNode, CHILD_NUM> subNodes;
|
|
|
|
std::vector<OcTreeNode> subNodes(CHILD_NUM); |
|
|
|
// intermediateNodes.resize(lastIdx + 8);
|
|
|
|
int subIdx = 0; |
|
|
|
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { |
|
|
|
auto nodeAABB = detail::getOneEightCellAABB(node.min, node.max, j(), 0); |
|
|
|
subNodes[subIdx].min = nodeAABB.first; |
|
|
|
subNodes[subIdx].max = nodeAABB.second; |
|
|
|
subNodes[subIdx].blobTree = node.blobTree; |
|
|
|
subNodes[subIdx].polyIntersectIndices = std::vector<int>(); |
|
|
|
} |
|
|
|
uvector3 nodeMid = (nowNodeMin + nowNodeMax) / 2; |
|
|
|
for (int i = 0; i < polyIntersectIndices.size(); ++i) { |
|
|
|
const int polyIntersectIndex = polyIntersectIndices[i]; |
|
|
|
const auto& poly = scene.polys[polyIntersectIndex]; |
|
|
|
subIdx = 0; |
|
|
|
for (MultiLoop<3> j(0, 2); ~j; ++j, ++subIdx) { |
|
|
|
tensor3 subcellPoly(nullptr, poly.compositedBernstein.ext()); |
|
|
|
algoim_spark_alloc(real, subcellPoly); |
|
|
|
bernstein::deCasteljau(poly.compositedBernstein, subNodes[subIdx].min, subNodes[subIdx].max, subcellPoly); |
|
|
|
int sign = bernstein::uniformSign(subcellPoly); |
|
|
|
if (sign == 1) { |
|
|
|
// organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_OUT);
|
|
|
|
} else if (sign == -1) { |
|
|
|
// organizer::traverse(subNodes[subIdx].blobTree, polyIntersectIndex, organizer::NODE_IN);
|
|
|
|
} else { |
|
|
|
subNodes[subIdx].polyIntersectIndices.emplace_back(polyIntersectIndex); |
|
|
|
} |
|
|
|
} |
|
|
|
uvector3 testX(0.5); |
|
|
|
real evalX = bernstein::evalBernsteinPoly(poly.compositedBernstein, testX); |
|
|
|
// std::cout << "evalX bernstein within 0-1 = " << evalX << std::endl;
|
|
|
|
} |
|
|
|
for (subIdx = 0; subIdx < CHILD_NUM; ++subIdx) { buildOcTreeV0(scene, subNodes[subIdx], leaves); } |
|
|
|
} |
|
|
|
|
|
|
|
/** single object quadrature */ |
|
|
|
void basicTask(const std::shared_ptr<PrimitiveDesc>& p, int q = 20, real xmin = -1, real xmax = 1) |
|
|
|
{ |
|
|
@ -465,16 +508,15 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv |
|
|
|
const organizer::BlobTree& blobTree, |
|
|
|
const real q = 10) |
|
|
|
{ |
|
|
|
OcTreeNode rootNode(xmin, xmax, blobTree); |
|
|
|
std::vector<OcTreeNode> leaves; |
|
|
|
|
|
|
|
std::vector<SparkStack<real>*> tensorStacks; |
|
|
|
std::vector<CompleteTensorRep> completeTensorReps; |
|
|
|
|
|
|
|
real volume; |
|
|
|
real volume = 0; |
|
|
|
auto integrand = [](const uvector<real, 3>& x) { return 1.0; }; |
|
|
|
uvector3 range = xmax - xmin; |
|
|
|
uvector<real, 3> testX(0., 0.75, 0.2); |
|
|
|
uvector<real, 3> testX(0.8, 0.8, 0.8); |
|
|
|
for (int i = 0; i < primitives.size(); i++) { |
|
|
|
if (auto pt = std::dynamic_pointer_cast<SphereDesc>(primitives[i])) { |
|
|
|
tensor3 originTensor(nullptr, 3); |
|
|
@ -487,9 +529,10 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv |
|
|
|
|
|
|
|
algoim_spark_alloc(real, transformedTensor); |
|
|
|
makeSphere(*pt, originTensor); |
|
|
|
|
|
|
|
detail::powerTransformation(range, xmin, originTensor, transformedTensor); |
|
|
|
|
|
|
|
real testEvaOri = evalPower(originTensor, testX); |
|
|
|
|
|
|
|
detail::power2BernsteinTensor(transformedTensor, phi); |
|
|
|
completeTensorReps.emplace_back(CompleteTensorRep{phi, {originTensor}}); |
|
|
|
} else if (auto pt = std::dynamic_pointer_cast<MeshDesc>(primitives[i])) { |
|
|
@ -514,14 +557,54 @@ void quadratureScene(const std::vector<std::shared_ptr<PrimitiveDesc>>& primitiv |
|
|
|
real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); |
|
|
|
|
|
|
|
completeTensorReps.emplace_back(CompleteTensorRep{phi, planeTensors}); |
|
|
|
} else if (auto pt = std::dynamic_pointer_cast<CylinderDesc>(primitives[i])) { |
|
|
|
std::vector<tensor3> rawTensors(3, tensor3(nullptr, 2)); |
|
|
|
rawTensors[0].ext_ = 3; |
|
|
|
rawTensors[0].ext_(pt->alignAxis) = 1; |
|
|
|
algoimSparkAllocHeapVector(tensorStacks, rawTensors); |
|
|
|
|
|
|
|
uvector3i resExt = detail::getCompositeExt(rawTensors); |
|
|
|
tensor3 phi(nullptr, resExt); |
|
|
|
tensorStacks.emplace_back(algoim_spark_alloc_heap(real, phi)); |
|
|
|
|
|
|
|
tensor3 compositeTensor(nullptr, resExt); |
|
|
|
algoim_spark_alloc(real, compositeTensor); |
|
|
|
|
|
|
|
makeCylinder(*pt, compositeTensor, rawTensors); |
|
|
|
detail::powerTransformation(range, xmin, compositeTensor); |
|
|
|
real testEvalPower = evalPower(compositeTensor, testX); |
|
|
|
|
|
|
|
detail::power2BernsteinTensor(compositeTensor, phi); |
|
|
|
|
|
|
|
completeTensorReps.emplace_back(CompleteTensorRep{phi, rawTensors}); |
|
|
|
} else if (auto pt = std::dynamic_pointer_cast<ConeDesc>(primitives[i])) { |
|
|
|
tensor3 originTensor(nullptr, 3); |
|
|
|
tensor3 transformedTensor(nullptr, 3); |
|
|
|
std::vector<tensor3> rawTensors(3, tensor3(nullptr, 2)); |
|
|
|
rawTensors[0].ext_ = 3; |
|
|
|
algoimSparkAllocHeapVector(tensorStacks, rawTensors); |
|
|
|
|
|
|
|
uvector<int, 3> resExt = 2 + 1 + 1 + 1; |
|
|
|
tensor3 phi(nullptr, resExt); |
|
|
|
tensorStacks.emplace_back(algoim_spark_alloc_heap(real, phi)); |
|
|
|
|
|
|
|
tensor3 compositeTensor(nullptr, resExt); |
|
|
|
algoim_spark_alloc(real, compositeTensor); |
|
|
|
|
|
|
|
makeCone(*pt, compositeTensor, rawTensors); |
|
|
|
detail::powerTransformation(range, xmin, compositeTensor); |
|
|
|
real testEvalPower = evalPower(compositeTensor, testX); |
|
|
|
|
|
|
|
detail::power2BernsteinTensor(compositeTensor, phi); |
|
|
|
|
|
|
|
completeTensorReps.emplace_back(CompleteTensorRep{phi, rawTensors}); |
|
|
|
} else { |
|
|
|
std::cerr << "unrecognized primitive type." << std::endl; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
Scene scene{completeTensorReps, blobTree}; |
|
|
|
buildOcTree(scene, rootNode, leaves); |
|
|
|
Scene scene{completeTensorReps, blobTree}; |
|
|
|
OcTreeNode rootNode(0, 1, blobTree); |
|
|
|
for (int i = 0; i < completeTensorReps.size(); ++i) { rootNode.polyIntersectIndices.emplace_back(i); } |
|
|
|
buildOcTreeV0(scene, rootNode, leaves); |
|
|
|
|
|
|
|
for (const auto& leaf : leaves) { |
|
|
|
auto basicRes = basicTask(scene, blobTree, leaf, q); |
|
|
|