diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4019e82 --- /dev/null +++ b/.gitignore @@ -0,0 +1,13 @@ +.DS_store +.idea + +build +cmake-build-debug +cmake-build-release + +input +output +result +tools + +test diff --git a/include/DomainDiscretization.hpp b/include/DomainDiscretization.hpp index c6f2e90..424cd63 100644 --- a/include/DomainDiscretization.hpp +++ b/include/DomainDiscretization.hpp @@ -88,9 +88,48 @@ namespace meshless { boundaryMap_.push_back(-1); return positions_.size() - 1; } - //bool DomainDiscretization::contains(const vec& point)const; - - + template + template + void DomainDiscretization::makeDiscreteContainsStructure(search_struc& search) const { + int size = boundary().size(); + std::vector boundary_points(size); + for(int i = 0; i < boundaryMap_.size(); i++) { + if(boundaryMap_[i] == -1) + continue; + boundary_points[boundaryMap_[i]] = positions_[i]; + } + search.reset(boundary_points); + } + //bool DomainDiscretization::contains(const vec& point)const; + template + template + bool DomainDiscretization::discreteContains(const Eigen::Vector3d& point, search_struc& search, int num)const { + assert_msg((search.size() != 0), "Search structure should not be empty."); + auto closet = search.query(point, num); + auto closet_indexs = closet.first; + for(int i = 0; i < closet_indexs.size(); i++) { + auto idx = closet_indexs[i]; + auto n = normals_[idx]; + auto closet_point = search.get(idx); + if(n.dot(point - closet_point) >= 0) { + return false; + } + } + // for (auto idx : closet_indexs) + // { + // auto n = normals_[idx]; + // auto closet_point = search.get(idx); + // if (n.dot(point - closet_point) >= 0) + // { + // return false; + // } + // } + return true; + // int closest_index = search.query(point, 1).first[0]; + // auto closest_point = search.get(closest_index); + // auto n = normals_[closest_index]; + // return n.dot(point - closest_point) < 0; + } }; \ No newline at end of file diff --git a/include/DomainDiscretization_fwd.hpp b/include/DomainDiscretization_fwd.hpp index 0b897a0..4a9a2bd 100644 --- a/include/DomainDiscretization_fwd.hpp +++ b/include/DomainDiscretization_fwd.hpp @@ -2,6 +2,7 @@ #include #include #include "OccShape.hpp" +#include "KDTreeMutable.hpp" namespace meshless { @@ -21,7 +22,8 @@ namespace meshless { public: DomainDiscretization() {} DomainDiscretization(const OccShape& shape) { - shape_.max_points = shape.max_points; + shape_.max_points_boundary = shape.max_points_boundary; + shape_.max_points_interior = shape.max_points_interior; shape_.seed_ = shape.seed_; shape_.n_samples = shape.n_samples; shape_.zeta = shape.zeta; @@ -184,6 +186,13 @@ namespace meshless { int addBoundaryNodeWithT(const vec& point, double t, int type, const vec& normal); int addInternalNodeWithT(const vec& point, double t, int type); + + template + bool discreteContains(const Eigen::Vector3d& point, search_struc& search, int num = 1)const; + + template + void makeDiscreteContainsStructure(search_struc& search) const; + private: int addNode(const vec& point, int type); int addNodeWithT(const vec& point, double t, int type); diff --git a/include/FillBoundary.hpp b/include/FillBoundary.hpp index 59b4c8c..5c09c1b 100644 --- a/include/FillBoundary.hpp +++ b/include/FillBoundary.hpp @@ -1,4 +1,5 @@ #pragma once +#include "discretization_helpers.hpp" #include "FillBoundary_fwd.hpp" #include "KDTree.hpp" #include "KDTreeMutable.hpp" @@ -28,6 +29,7 @@ namespace meshless { std::cout << "该面法线需反向\n"; reverseNormal = true; } + DomainDiscretization domainParamLocalPatch(shape); //开始离散化每个wire auto outerWire = BRepTools::OuterWire(face); @@ -39,10 +41,6 @@ namespace meshless { if(wire == outerWire) { std::cout << "此wire为外环\n"; } - - - DomainDiscretization domainParamLocalPatch(shape); - TopExp_Explorer expEdge(wire, TopAbs_EDGE); int numEdge = 0; int genPoints = 0; @@ -113,7 +111,7 @@ namespace meshless { int curNode = domainParamLocalPatch.size() - 1; int endNode = domainParamLocalPatch.size(); - while(curNode < endNode && endNode < shape.max_points) { + while(curNode < endNode && endNode < shape.max_points_boundary) { genPoints++; //std::cout << "curNode / endNode : " << curNode << '/' << endNode << '\n'; auto param = domainParamLocalPatch.pos(curNode); @@ -123,31 +121,33 @@ namespace meshless { gp_Pnt pt; gp_Vec der; geomCurve->D1(t, pt, der); - std::cout << "pt: " << pt.X() << "," << pt.Y() << "," << pt.Z() << '\n'; - std::cout << "der: " << der.X() << "," << der.Y() << "," << der.Z() << '\n'; - - - std::cout << "der.Mag:" << der.Magnitude() << '\n'; + //std::cout << "pt: " << pt.X() << "," << pt.Y() << "," << pt.Z() << '\n'; + //std::cout << "der: " << der.X() << "," << der.Y() << "," << der.Z() << '\n'; + //std::cout << "der.Mag:" << der.Magnitude() << '\n'; double alpha = h(transPnt(pt)) / der.Magnitude(); for(const auto& uCan : candidates) { + //参数点必须在规定的参数域内 double tNew = t + alpha * uCan; if(tNew > pCurveLastParam || tNew < pCurveFirstParam) { - std::cout << "排除\n"; + //std::cout << "排除\n"; continue; } - std::cout << "tNew:" << tNew << "\n"; + //std::cout << "tNew:" << tNew << "\n"; gp_Pnt2d tNewSeedInPCurve; geomPCurve->D0(tNew, tNewSeedInPCurve); gp_Pnt ptNew; geomCurve->D0(tNew, ptNew); gp_Pnt ptNew2; geomSurface->D0(tNewSeedInPCurve.X(), tNewSeedInPCurve.Y(), ptNew2); - std::cout << "ptNew: " << ptNew.X() << "," << ptNew.Y() << "," << ptNew.Z() << '\n'; - std::cout << "ptNew2: " << ptNew2.X() << "," << ptNew2.Y() << "," << ptNew2.Z() << '\n'; + //std::cout << "ptNew: " << ptNew.X() << "," << ptNew.Y() << "," << ptNew.Z() << '\n'; + //std::cout << "ptNew2: " << ptNew2.X() << "," << ptNew2.Y() << "," << ptNew2.Z() << '\n'; double checkradius = pt.SquareDistance(ptNew); - std::cout << "CheckRadius: " << checkradius << '\n'; + //std::cout << "CheckRadius: " << checkradius << '\n'; + + + //距离已有的点不能太近 double d_sq_loc = treeLocal.query(transPnt(ptNew)).second[0]; - std::cout << "d_sq_loc: " << d_sq_loc << '\n'; + //std::cout << "d_sq_loc: " << d_sq_loc << '\n'; if(d_sq_loc >= (shape.zeta * shape.zeta * checkradius)) { domainParamLocalPatch.addInternalNodeWithT(transPnt2d(tNewSeedInPCurve), tNew, 1); treeLocal.insert(transPnt(ptNew)); @@ -166,10 +166,10 @@ namespace meshless { domainGlobal.addBoundaryNode(transPnt(ptNew), type, transVec(tSeedNormalNew)); treeGlobal.insert(transPnt(ptNew)); } else { - std::cout << "全局不添加点\n"; + //std::cout << "全局不添加点\n"; } } else { - std::cout << "局部不添加点\n"; + //std::cout << "局部不添加点\n"; } } curNode++; @@ -177,17 +177,79 @@ namespace meshless { std::cout << "该边已经生成：" << genPoints << "个点\n"; } } - //fill boundary surface! KDTree paramBoundarySearch; + int curNode = 0, endNode = domainParamLocalPatch.size(); + while(curNode < endNode && endNode < shape.max_points_boundary) { + auto paramPoint = domainParamLocalPatch.pos(curNode); + gp_Pnt initPt; + gp_Vec derU, derV; + //geomSurface->D0(paramPoint.x(), paramPoint.y(), initPt); + geomSurface->D1(paramPoint.x(), paramPoint.y(), initPt, derU, derV); + Eigen::MatrixXd jacobMat; + jacobMat.resize(3, 2); + jacobMat(0, 0) = derU.X(); + jacobMat(0, 1) = derV.X(); + jacobMat(1, 0) = derU.Y(); + jacobMat(1, 1) = derV.Y(); + jacobMat(2, 0) = derU.Z(); + jacobMat(2, 1) = derV.Z(); + + auto unitCandidates = SphereDiscretization::construct(1, shape.n_samples, gen); + + for(auto& uCand : unitCandidates) { + double alpha = h(transPnt(initPt)) / (jacobMat * uCand).norm(); + auto node = paramPoint + alpha * uCand; + gp_Pnt2d nodePnt(node.x(), node.y()); + + //该参数点必须在参数域内 + BRepClass_FaceClassifier brepF; + brepF.Perform(face, nodePnt, 1e-6); + if(brepF.State() == TopAbs_State::TopAbs_IN) { + //std::cout << "TopAbs_IN\n"; + } + if(brepF.State() == TopAbs_State::TopAbs_ON) { + //std::cout << "TopAbs_ON\n"; + continue; + } + if(brepF.State() == TopAbs_State::TopAbs_OUT) { + //std::cout << "TopAbs_OUT\n"; + continue; + } + if(brepF.State() == TopAbs_State::TopAbs_UNKNOWN) { + //std::cout << "TopAbs_UNKNOWN\n"; + continue; + } + gp_Pnt newPt; + gp_Vec derU, derV; + geomSurface->D1(node.x(), node.y(), newPt, derU, derV); - - - - + double d_sq = treeGlobal.query(transPnt(newPt)).second[0]; + double checkradius = newPt.SquareDistance(initPt); + if(d_sq < (shape.zeta * shape.zeta * checkradius)) { + continue; + } + domainParamLocalPatch.addInternalNode(node, 1); + treeGlobal.insert(transPnt(newPt)); + gp_Dir normal; + auto retStatus = GeomLib::NormEstim(geomSurface, nodePnt, 1e-6, normal); + if(retStatus >= 2) { + std::cout << "calculate wrong!\n"; + exit(-1); + } + if(reverseNormal) { + normal = -normal; + } + domainGlobal.addBoundaryNode(transPnt(newPt), type, transVec(normal)); + endNode++; + } + curNode++; + } + if(endNode >= shape.max_points_boundary) { + std::cerr << "Maximum number of points reached, fill may be incomplete." << std::endl; + } } return domainGlobal; } - }; \ No newline at end of file diff --git a/include/FillInterior.hpp b/include/FillInterior.hpp new file mode 100644 index 0000000..ea321ea --- /dev/null +++ b/include/FillInterior.hpp @@ -0,0 +1,46 @@ +#pragma once + +#include "assert.hpp" +#include "discretization_helpers.hpp" +#include "FillInterior_fwd.hpp" + +namespace meshless { + template + void fillInteriorWithDensity(const OccShape& shape, DomainDiscretization& domain, const std::function& h, KDTreeMutable& search, contains_func& containsFunc, int type) { + if(type == 0) + type = 1; + std::mt19937 gen(shape.seed_); + + std::cout << "Filling the domain interior...\n"; + int curNode = 0; + int endNode = domain.size(); + + search.insert(domain.positions()); + while((curNode < endNode) && (endNode < shape.max_points_interior)) { + std::cout << "CurNode/EndNode: " << curNode << "/" << endNode << '\n'; + auto p = domain.pos(curNode); + double r = h(p); + assert_msg(r > 0, "Nodal spacing radius should be > 0, got %g.", r); + + auto candidates = SphereDiscretization::construct(r, shape.n_samples, gen); + for(const auto& f : candidates) { + auto node = p + f; + if(!containsFunc(node))continue; + + if(search.existsPointInSphere(node, r * shape.zeta)) { + continue; + } + + domain.addInternalNode(node, type); + search.insert(node); + endNode++; + } + curNode++; + } + if(endNode >= shape.max_points_interior) { + std::cerr << "Maximum number of points reached, fill may be incomplete." << std::endl; + } + std::cout << "Filling the domain interior done!\n"; + } + +}; \ No newline at end of file diff --git a/include/FillInterior_fwd.hpp b/include/FillInterior_fwd.hpp new file mode 100644 index 0000000..edd6c21 --- /dev/null +++ b/include/FillInterior_fwd.hpp @@ -0,0 +1,24 @@ + +#pragma once +#include "OccHelper.hpp" +#include "OccShape.hpp" +#include "DomainDiscretization.hpp" +#include "KDTree.hpp" +#include "KDTreeMutable.hpp" +#include +#include + +namespace meshless { + + + template + void fillInteriorWithDensity(const OccShape& shape, DomainDiscretization& domain, const std::function& h, KDTreeMutable& search, contains_func& containsFunc, int type); + + template + void fillInteriorWithStep(const OccShape& shape, DomainDiscretization& domain, double step, KDTreeMutable& search, contains_func& containsFunc, int type) { + auto f = [=](Eigen::Vector3d) {return step; }; + fillInteriorWithDensity(shape, domain, f, search, containsFunc, type); + } + +}; + diff --git a/include/OccHelper.hpp b/include/OccHelper.hpp index bf40075..19ed383 100644 --- a/include/OccHelper.hpp +++ b/include/OccHelper.hpp @@ -2,6 +2,7 @@ // OpenCascade includes #include +#include #include #include #include diff --git a/include/OccShape.hpp b/include/OccShape.hpp index 83b9137..54ead7b 100644 --- a/include/OccShape.hpp +++ b/include/OccShape.hpp @@ -5,7 +5,9 @@ namespace meshless { class OccShape { public: - int max_points = 500000; + int max_points_boundary = 500000; + int max_points_interior = 500000; + int seed_; int n_samples = 15; double zeta = 1 - 1e-10; @@ -14,15 +16,19 @@ namespace meshless { public: OccShape() {} OccShape(const TopoDS_Shape& myShape_) { - max_points = 500000; + max_points_boundary = 500000; seed_; n_samples = 15; zeta = 1 - 1e-10; epsilon = 0; myShape = myShape_; } - OccShape& maxPoints(int max_points) { - this->max_points = max_points; + OccShape& maxPointsBoundary(int max_points_boundary) { + this->max_points_boundary = max_points_boundary; + return *this; + } + OccShape& maxPointsInterior(int max_points_interior) { + this->max_points_interior = max_points_interior; return *this; } diff --git a/include/discretization_helpers.hpp b/include/discretization_helpers.hpp new file mode 100644 index 0000000..1227296 --- /dev/null +++ b/include/discretization_helpers.hpp @@ -0,0 +1,128 @@ +#pragma once + +#include +#include + +namespace meshless { + + template + struct Pi { + static constexpr T value = T(3.14159265358979323846264338327950L); ///< Value of @f$\pi@f$. + }; + + static const double PI = Pi::value; ///< Mathematical constant pi in double precision. + static const double INF = 1.0 / 0.0; ///< Infinite floating point value. + static const double NaN = 0.0 / 0.0; ///< Not-a-number floating point value. + + template + struct SphereDiscretization { + /** + * Construct a randomized discretization. + * @param radius Radius of the sphere. + * @param num_samples Number of points on the equator, implies nodal spacing `dp = 2*pi*r/n`. + * @param generator A random number generator. + * @return A vector of discretization points. + */ + template + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int num_samples, generator_t& generator) { + scalar_t dphi = 2 * Pi::value / num_samples; + std::uniform_real_distribution distribution(0, Pi::value); + scalar_t offset = distribution(generator); + std::vector, + Eigen::aligned_allocator>> result; + for(int i = 0; i < num_samples / 2; ++i) { + scalar_t phi = i * dphi + offset; + if(phi > Pi::value) phi -= Pi::value; + int slice_n = static_cast(std::ceil(num_samples * std::sin(phi))); + if(slice_n == 0) continue; + auto slice = SphereDiscretization::construct( + radius * std::sin(phi), slice_n, generator); + Eigen::Matrix v; + for(const auto& p : slice) { + v[0] = radius * std::cos(phi); + v.template tail() = p; + result.push_back(v); + } + } + return result; + } + /// Construct the discretization. + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int num_samples) { + scalar_t dphi = 2 * Pi::value / num_samples; + std::vector, + Eigen::aligned_allocator>> result; + for(int i = 0; i < num_samples / 2; ++i) { + scalar_t phi = i * dphi; + if(phi > Pi::value) phi -= Pi::value; + int slice_n = static_cast(std::ceil(num_samples * std::sin(phi))); + if(slice_n == 0) continue; + auto slice = SphereDiscretization::construct( + radius * std::sin(phi), slice_n); + Eigen::Matrix v; + for(const auto& p : slice) { + v[0] = radius * std::cos(phi); + v.template tail() = p; + result.push_back(v); + } + } + return result; + } + }; + + /// Two-dimensional base case of the discretization. + template + struct SphereDiscretization { + /// Construct a randomized discretization. + template + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int num_samples, generator_t& generator) { + scalar_t dphi = 2 * Pi::value / num_samples; + std::uniform_real_distribution distribution(0, 2 * Pi::value); + scalar_t offset = distribution(generator); + std::vector, + Eigen::aligned_allocator>> result; + for(int i = 0; i < num_samples; ++i) { + scalar_t phi = i * dphi + offset; + result.emplace_back(radius * std::cos(phi), radius * std::sin(phi)); + } + return result; + } + /// Construct the discretization. + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int num_samples) { + scalar_t dphi = 2 * Pi::value / num_samples; + std::vector, + Eigen::aligned_allocator>> result; + for(int i = 0; i < num_samples; ++i) { + scalar_t phi = i * dphi; + result.emplace_back(radius * std::cos(phi), radius * std::sin(phi)); + } + return result; + } + }; + + /// One-dimensional base case of the discretization. + template + struct SphereDiscretization { + template + /// Construct a randomized discretization. + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int, generator_t&) { + return { Eigen::Matrix(-radius), Eigen::Matrix(radius) }; + } + /// Construct the discretization. + static std::vector, + Eigen::aligned_allocator>> construct( + scalar_t radius, int) { + return { Eigen::Matrix(-radius), Eigen::Matrix(radius) }; + } + }; + +}; \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 8073c76..c348269 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,7 +1,9 @@ #include #include #include +#include "KDTree.hpp" #include "FillBoundary.hpp" +#include "FillInterior.hpp" #include "writeVTK.hpp" using namespace std; @@ -12,11 +14,24 @@ int main() { BRep_Builder brepBuilder; TopoDS_Shape myShape; //auto status = BRepTools::Read(myShape, "F:/CADMeshless/StepAndBrepModel/DWE_ADM102_ASSIGN.brep", brepBuilder); - auto status = BRepTools::Read(myShape, "F:/CADMeshless/StepAndBrepModel/Exhaust_Manifold.brep", brepBuilder); + auto status = BRepTools::Read(myShape, "F:/CADMeshless/StepAndBrepModel/DWE_ADM102_ASSIGN.brep", brepBuilder); OccShape occShape(myShape); - auto domain = discretizeBoundaryWithStep(occShape, 1, -1); - Eigen::MatrixXd quaPoints; + double h1 = 1.0; + double h2 = 4.0; + DomainDiscretization domain = discretizeBoundaryWithStep(occShape, h1, -1); + DomainDiscretization domainOverSampled = discretizeBoundaryWithStep(occShape, 0.2, -1); + + + KDTree containsTree; + domainOverSampled.makeDiscreteContainsStructure(containsTree); + auto containsFunc = [&](const Eigen::Vector3d p) { + return domainOverSampled.discreteContains(p, containsTree, 15); + }; + KDTreeMutable tree; + fillInteriorWithStep(myShape, domain, h2, tree, containsFunc, 1); + + Eigen::MatrixXd quaPoints; quaPoints.resize(domain.positions().size(), 3); std::vector attri; for(int i = 0; i < domain.positions().size(); i++) { @@ -29,7 +44,7 @@ int main() { quaPoints(i, 1) = domain.positions()[i].y(); quaPoints(i, 2) = domain.positions()[i].z(); } - writePntVTK<3>("Exhaust_Manifold.vtk", quaPoints, attri); + writePntVTK<3>("DWE_ADM102_ASSIGN.vtk", quaPoints, attri);