diff --git a/.clang-format b/.clang-format index 7b54908..25cde52 100644 --- a/.clang-format +++ b/.clang-format @@ -1,4 +1,5 @@ BasedOnStyle: Google + AccessModifierOffset: -4 AlignAfterOpenBracket: Align AlignArrayOfStructures: Left @@ -51,7 +52,7 @@ AlwaysBreakBeforeMultilineStrings: false AlwaysBreakTemplateDeclarations: Yes AttributeMacros: - __capability -BinPackArguments: false +BinPackArguments: true BinPackParameters: false BitFieldColonSpacing: Both BraceWrapping: diff --git a/algoim/organizer/organizer.hpp b/algoim/organizer/organizer.hpp index b91d419..227a892 100644 --- a/algoim/organizer/organizer.hpp +++ b/algoim/organizer/organizer.hpp @@ -158,7 +158,7 @@ const std::array sides = {-1, 1}; struct BasicTaskRes { real volume; - real area; + real surf; }; // 对于mark含2个及以上0的情况，尝试对每个0填1或-1的所有组合 @@ -257,7 +257,8 @@ void buildOcTree(const Scene& scene, const OcTreeNode& node, std::vector ipquad(scene.polys); ipquad.integrate(AutoMixed, q, [&](const uvector& x, real w) { auto realX = x * range + node.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& x, real w, const uvector& wn) { + surf += w * integrand(x * range + node.min); + }); + return {volume, surf}; } void quadratureScene(const std::vector>& primitives, const uvector3& xmin, const uvector3& xmax, - const organizer::BlobTree& blobTree) + const organizer::BlobTree& blobTree, + const real q = 10) { OcTreeNode rootNode(xmin, xmax, blobTree); std::vector leaves; @@ -511,6 +514,9 @@ void quadratureScene(const std::vector>& primitiv real testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX); completeTensorReps.emplace_back(CompleteTensorRep{phi, planeTensors}); + } else if (auto pt = std::dynamic_pointer_cast(primitives[i])) { + tensor3 originTensor(nullptr, 3); + tensor3 transformedTensor(nullptr, 3); } } @@ -518,15 +524,13 @@ void quadratureScene(const std::vector>& primitiv buildOcTree(scene, rootNode, leaves); for (const auto& leaf : leaves) { - auto basicRes = basicTask(scene, blobTree, leaf, 10); - volume += basicRes.volume; + auto basicRes = basicTask(scene, blobTree, leaf, q); + volume += basicRes.volume * prod(leaf.max - leaf.min); } 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; } diff --git a/algoim/organizer/primitive.hpp b/algoim/organizer/primitive.hpp index 3def130..df32a49 100644 --- a/algoim/organizer/primitive.hpp +++ b/algoim/organizer/primitive.hpp @@ -451,12 +451,12 @@ private: class BRepDesc : virtual public PrimitiveDesc { +public: const static PrimitiveType type = BRep; std::vector vertices; std::vector curves; std::vector surfaces; -public: void print() override { std::cout << "BRep Description" << std::endl; } real eval(const uvector3& p) @@ -495,18 +495,28 @@ public: class CylinderDesc : virtual public PrimitiveDesc { +public: const static PrimitiveType type = Cylinder; uvector3 node1; uvector3 node2; real radius; - CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} + // CylinderDesc(const uvector3& n1_, const uvector3& n2_, real r_) : PrimitiveDesc(), node1(n1_), node2(n2_), radius(r_) {} + // 现在是只支持轴对齐x,y,z的圆柱 + // TODO: 实现张量积的旋转后换成支持任意方向的 + CylinderDesc(const uvector3& n1, real r, real h, int alignAxis) : PrimitiveDesc(), node1(n1), radius(r) + { + assert(alignAxis >= 0 && alignAxis <= 2); + node2 = node1; + node2(alignAxis) += h; + } void print() override { std::cout << "Cylinder Description" << std::endl; } }; class ConeDesc : virtual public PrimitiveDesc { +public: const static PrimitiveType type = Cone; uvector3 node1; uvector3 node2; @@ -532,10 +542,60 @@ public: { } + MeshDesc() = default; + void print() override { std::cout << "Mesh Description" << std::endl; } }; +class PyramidDesc : public MeshDesc +{ +public: + PyramidDesc(const std::vector& bottomVertices, const uvector3& topVertex) : PrimitiveDesc() + { + vertices = bottomVertices; + int bottomSize = bottomVertices.size(); + vertices.emplace_back(topVertex); + indices.reserve(bottomSize * 3 + bottomSize); + indexInclusiveScan.reserve(bottomSize + 1); + for (int i = 0; i < bottomSize; ++i) { + indices.emplace_back(i); + indices.emplace_back(i == bottomSize - 1 ? 0 : indices.emplace_back(i + 1)); + indices.emplace_back(bottomSize); + indexInclusiveScan.emplace_back(indexInclusiveScan.empty() ? 3 : indexInclusiveScan.back() + 3); + } + for (int i = 0; i < bottomSize; ++i) { + indices.emplace_back(i); + indices.emplace_back(i); + } + indexInclusiveScan.emplace_back(indexInclusiveScan.back() + bottomSize); + } +}; + +class CuboidDesc : public MeshDesc +{ +public: + CuboidDesc(const uvector3& center, const uvector3& size) : MeshDesc() + { + vertices.resize(8); + uvector3 halfSize = size / 2; + for (int i = 0; i < 8; ++i) { + vertices[i] = center; + vertices[i](0) += (i & 1) ? halfSize(0) : -halfSize(0); + vertices[i](1) += (i & 2) ? halfSize(1) : -halfSize(1); + vertices[i](2) += (i & 4) ? halfSize(2) : -halfSize(2); + } + indices = {3, 2, 0, 1, // left + 4, 6, 7, 5, // right + 6, 2, 3, 7, // top + 1, 0, 4, 5, // bottom + 7, 3, 1, 5, // front + 2, 6, 4, 0}; // back + indexInclusiveScan = {4, 8, 12, 16, 20, 24}; + } +}; + void makeMesh(const MeshDesc& mesh, xarray& tensor, std::vector>& planeTensors) + { uvector3 ext(1 + mesh.indexInclusiveScan.size()); assert(all(ext == tensor.ext())); @@ -585,9 +645,18 @@ void makeSphere(const SphereDesc& sphereDesc, xarray& tensor) // return PowerTensor(tensor); }; -void makeCylinder(xarray& tensor, uvector3 startPt, uvector3 endPt, real r) +void makeCylinder(const CylinderDesc& cylinder, xarray& tensor, std::vector& rawTensors) +{ + assert(rawTensors.size() == 3 + && (all(rawTensors[0].ext() == uvector3(3, 3, 1)) || all(rawTensors[0].ext() == uvector3(3, 3, 1)) + || all(rawTensors[0].ext() == uvector3(3, 3, 1))) + && all(rawTensors[1].ext() == uvector3(2)) && all(rawTensors[1].ext() == uvector3(2))); + uvector ext = 5; + assert(all(ext == tensor.ext())); +} + +void makeCone(const CylinderDesc& cone, tensor3& tensor) { - // PowerTensor pt; // TODO: } @@ -605,8 +674,8 @@ struct Scene { const int CHILD_NUM = 8; struct OcTreeNode { - std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 - std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 + std::vector polyIntersectIndices; // 同一poly会出现在不同node的polyIndices中 + // std::vector polyFullyContainedIndices; // 完全包含该node的poly, 同样会出现在不同node的polyIndices中 // std::array children; // octree uvector3 min, max; organizer::BlobTree blobTree; // 部分遍历过的octree，那些不属于该node的primitive的out信息已经在该blobtree种尽可能地向上传递 @@ -620,7 +689,7 @@ struct OcTreeNode { // 深拷贝 OcTreeNode(const OcTreeNode& node) : polyIntersectIndices(node.polyIntersectIndices), - polyFullyContainedIndices(node.polyFullyContainedIndices), + // polyFullyContainedIndices(node.polyFullyContainedIndices), min(node.min), max(node.max), blobTree(node.blobTree) diff --git a/algoim/quadrature_multipoly.hpp b/algoim/quadrature_multipoly.hpp index befbb11..2eccac2 100644 --- a/algoim/quadrature_multipoly.hpp +++ b/algoim/quadrature_multipoly.hpp @@ -623,7 +623,7 @@ struct ImplicitPolyQuadrature { assert(max(abs(score)) > 0); score /= 2 * max(abs(score)); for (int i = 0; i < N; ++i) - if (!has_disc(i)) score(i) += 1.0; + if (!has_disc(i)) score(i) += 1.0; // 1个与k相切的多项式都没有，才+1 // Choose height direction and form base polynomials; if tanh-sinh is being used at this // level, suggest the same all the way down; moreover, suggest tanh-sinh if a non-empty @@ -689,7 +689,7 @@ struct ImplicitPolyQuadrature { // For every phi(i) ... for (size_t i = 0; i < phi.count(); ++i) { const auto& p = phi.poly(i); - const auto& mask = phi.mask(i); // 这个mask干什么的？ + const auto& mask = phi.mask(i); int P = p.ext(k); // Ignore phi if its mask is void everywhere above the base point @@ -702,6 +702,7 @@ struct ImplicitPolyQuadrature { int rcount = bernstein::bernsteinUnitIntervalRealRoots(pline, P, roots); // Add all real roots in [0,1] which are also within masked region of phi + // 为什么[0,1]的real roots还必须要在mask为true的subregion？ for (int j = 0; j < rcount; ++j) { auto x = add_component(xbase, k, roots[j]); if (detail::pointWithinMask(mask, x)) nodes[count++] = roots[j]; diff --git a/gjj/primitiveDebug.hpp b/gjj/primitiveDebug.hpp index 030e3e4..5be5d4f 100644 --- a/gjj/primitiveDebug.hpp +++ b/gjj/primitiveDebug.hpp @@ -35,14 +35,9 @@ void casePolyhedron1() // std::vector> ps; // mesh - std::vector vertices = {uvector3(-0.8, -0.8, -0.8), - uvector3(-0.8, -0.8, 0.8), - uvector3(-0.8, 0.8, -0.8), - uvector3(-0.8, 0.8, 0.8), - uvector3(0.8, -0.8, -0.8), - uvector3(0.8, -0.8, 0.8), - uvector3(0.8, 0.8, -0.8), - uvector3(0.8, 0.8, 0.8)}; + std::vector vertices = {uvector3(-0.8, -0.8, -0.8), uvector3(-0.8, -0.8, 0.8), uvector3(-0.8, 0.8, -0.8), + uvector3(-0.8, 0.8, 0.8), uvector3(0.8, -0.8, -0.8), uvector3(0.8, -0.8, 0.8), + uvector3(0.8, 0.8, -0.8), uvector3(0.8, 0.8, 0.8)}; std::vector indices = {3, 2, 0, 1, // left 4, 6, 7, 5, // right 6, 2, 3, 7, // top @@ -63,23 +58,14 @@ void casePolyhedron2() // std::vector> ps; // mesh - std::vector vertices = {uvector3(-0.8, -0.8, -0.8), - uvector3(-0.8, -0.8, 0.8), - uvector3(-0.8, 0.8, -0.8), - uvector3(-0.8, 0.8, 0.8), - uvector3(0.8, -0.8, -0.8), - uvector3(0.8, -0.8, 0.8), - uvector3(0.8, 0.8, -0.8), - uvector3(0.8, 0.8, 0.8)}; + std::vector vertices = {uvector3(-0.8, -0.8, -0.8), uvector3(-0.8, -0.8, 0.8), uvector3(-0.8, 0.8, -0.8), + uvector3(-0.8, 0.8, 0.8), uvector3(0.8, -0.8, -0.8), uvector3(0.8, -0.8, 0.8), + uvector3(0.8, 0.8, -0.8), uvector3(0.8, 0.8, 0.8)}; std::vector indices = { - 3, - 2, - 0, + 3, 2, 0, 1, // left - 6, - 2, - 3, - 7 // top + 6, 2, 3, + 7 // top }; std::vector scan = {4, 8}; @@ -101,31 +87,19 @@ void casePolyhedron3() // uvector3(1.6, 0, 1.6), // uvector3(1.6, 1.6, 0), // uvector3(1.6, 1.6, 1.6)}; - std::vector vertices = {uvector3(-0.8, -0.8, -0.8), - uvector3(-0.8, -0.8, 0.8), - uvector3(-0.8, 0.8, -0.8), - uvector3(-0.8, 0.8, 0.8), - uvector3(0.8, -0.8, -0.8), - uvector3(0.8, -0.8, 0.8), - uvector3(0.8, 0.8, -0.8), - uvector3(0.8, 0.8, 0.8)}; - std::vector indices = {3, - 2, - 0, - 1, // left - 6, - 2, - 3, - 7, // top - // 7, - // 3, - // 1, - // 5, // front - 2, - 6, - 4, - 0}; // back - std::vector scan = {4, 8, 12}; + std::vector vertices = {uvector3(-0.8, -0.8, -0.8), uvector3(-0.8, -0.8, 0.8), uvector3(-0.8, 0.8, -0.8), + uvector3(-0.8, 0.8, 0.8), uvector3(0.8, -0.8, -0.8), uvector3(0.8, -0.8, 0.8), + uvector3(0.8, 0.8, -0.8), uvector3(0.8, 0.8, 0.8)}; + std::vector indices = {3, 2, 0, + 1, // left + 6, 2, 3, + 7, // top + // 7, + // 3, + // 1, + // 5, // front + 2, 6, 4, 0}; // back + std::vector scan = {4, 8, 12}; basicTask(std::make_shared(MeshDesc(vertices, indices, scan))); } @@ -133,26 +107,14 @@ void casePolyhedronSphere() { // PI * r^3 / 6 auto phi0 = std::make_shared(SphereDesc(0.8, uvector3(-0.8, 0.8, -0.8), 1.)); - std::vector vertices = {uvector3(-0.8, -0.8, -0.8), - uvector3(-0.8, -0.8, 0.8), - uvector3(-0.8, 0.8, -0.8), - uvector3(-0.8, 0.8, 0.8), - uvector3(0.8, -0.8, -0.8), - uvector3(0.8, -0.8, 0.8), - uvector3(0.8, 0.8, -0.8), - uvector3(0.8, 0.8, 0.8)}; - std::vector indices = {3, - 2, - 0, - 1, // left - 6, - 2, - 3, - 7, // top - 2, - 6, - 4, - 0}; // back + std::vector vertices = {uvector3(-0.8, -0.8, -0.8), uvector3(-0.8, -0.8, 0.8), uvector3(-0.8, 0.8, -0.8), + uvector3(-0.8, 0.8, 0.8), uvector3(0.8, -0.8, -0.8), uvector3(0.8, -0.8, 0.8), + uvector3(0.8, 0.8, -0.8), uvector3(0.8, 0.8, 0.8)}; + std::vector indices = {3, 2, 0, + 1, // left + 6, 2, 3, + 7, // top + 2, 6, 4, 0}; // back std::vector scan = {4, 8, 12}; std::vector> primitiveDescriptions(2); primitiveDescriptions[0] = phi0; @@ -160,6 +122,25 @@ void casePolyhedronSphere() basicTask(primitiveDescriptions); } +void caseScene() +{ + std::vector> primitiveDescriptions(7); + primitiveDescriptions[0] = std::make_shared(CuboidDesc(0., 1.6)); + primitiveDescriptions[1] = std::make_shared(CuboidDesc(uvector3(0.6, 0.6, -0.6), 2.)); + primitiveDescriptions[2] = std::make_shared(SphereDesc(0.7, uvector3(0.8, 0.8, 0.8), 1.)); + primitiveDescriptions[3] = std::make_shared(SphereDesc(0.5, uvector3(-0.3, -0.8, 0.8), 1.)); + primitiveDescriptions[4] = std::make_shared(ConeDesc(uvector3(0., -0.2, 0.), uvector3(0., -0.9, 0.), 0.4)); + std::vector pyramidBottomVertices = { + uvector3{-1, -1, 0}, + uvector3{1, -1, 0}, + uvector3{1, 1, 0}, + uvector3{-1, 1, 0} + }; + primitiveDescriptions[5] = std::make_shared(pyramidBottomVertices, uvector3{0, 0, 1}); + primitiveDescriptions[4] = std::make_shared(CylinderDesc(uvector3(-0.3, 0.3, 2.3), 0.4, 3.6, 1)); + quadratureScene(primitiveDescriptions, uvector3(-1., -1.3, -1.6), uvector3(1.6, 1.6, 2.3), organizer::BlobTree()); +} + void testDeCasteljau() { auto phiDesc = std::make_shared(SphereDesc(0.8, uvector3(0), 1.)); @@ -172,10 +153,8 @@ void testDeCasteljau() // bernstein::deCasteljau(tensorBernstein, uvector3(0), uvector3(0.2690598923241497 + 0.0000001), // transformedTensorBernstein); - bernstein::deCasteljau(tensorBernstein, - uvector3(0.2690598923241497 + 0.0000001), - uvector3(1 - 0.2690598923241497 - 0.0000001), - transformedTensorBernstein); + bernstein::deCasteljau(tensorBernstein, uvector3(0.2690598923241497 + 0.0000001), + uvector3(1 - 0.2690598923241497 - 0.0000001), transformedTensorBernstein); uvector3 testX = uvector3(0.5); int sign = bernstein::uniformSign(transformedTensorBernstein); @@ -191,7 +170,7 @@ void testDeCasteljau() void testBlob() { - orgzer::Blob blob = {5, 3, 4, 5, 6}; + organizer::Blob blob = {5, 3, 4, 5, 6}; // std::cout << blob.isPrimitive << std::endl; // std::cout << blob.nodeOp << std::endl;