first add

CMakeLists.txt

@@ -0,0 +1,42 @@
+cmake_minimum_required (VERSION 3.2.0)
+
+# Project name
+project (OCCTMeshless CXX)
+
+list(APPEND CMAKE_PREFIX_PATH "F:/OCC/demo/opencascade-install" )
+
+# Enable C++17
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+# OpenCascade
+find_package(OpenCASCADE REQUIRED)
+
+# Eigen3
+set(EIGEN3_INCLUDE_DIR "F:/eigen-3.4.0")
+
+# Configure C++ compiler's includes dir
+include_directories (SYSTEM ${OpenCASCADE_INCLUDE_DIR} )
+include_directories(${EIGEN3_INCLUDE_DIR})
+include_directories(${PROJECT_SOURCE_DIR}/include)
+
+
+file(GLOB HEADER_FILES "include/*.h*")
+file(GLOB SRC_FILES "src/*.cpp")
+source_group("include" FILES HEADER_FILES)
+
+# Add executable
+add_executable (${PROJECT_NAME}
+   ${SRC_FILES}
+  ${HEADER_FILES}
+)
+
+# Add linker options
+foreach (LIB ${OpenCASCADE_LIBRARIES})
+  target_link_libraries(${PROJECT_NAME} PUBLIC ${OpenCASCADE_LIBRARY_DIR}/${LIB}.lib)
+  target_link_libraries(${PROJECT_NAME} PUBLIC ${OpenCASCADE_LIBRARY_DIR}d/${LIB}.lib)
+endforeach()
+
+# Adjust runtime environment
+set_property(TARGET ${PROJECT_NAME} PROPERTY VS_DEBUGGER_ENVIRONMENT "PATH=$<$<CONFIG:DEBUG>:${OpenCASCADE_BINARY_DIR}d>$<$<NOT:$<CONFIG:DEBUG>>:${OpenCASCADE_BINARY_DIR}>;%PATH%")

include/DomainDiscretization.hpp

@@ -0,0 +1,96 @@
+#pragma once
+#include "assert.hpp"
+#include "DomainDiscretization_fwd.hpp"
+
+namespace meshless {
+	template<typename vec>
+	const vec& DomainDiscretization<vec>::normal(int i) const {
+		assert_msg(0 <= i && i <= size(), "Index %d out of range [0, %d)", i, size());
+		assert_msg(types_[i] < 0, "Node %d must be a boundary node, got type %d.", i, types_[i]);
+		assert_msg(boundaryMap_[i] != -1, "Node %d does not have a normal. Maybe you manually set"
+			" supports and positions instead of using addInternalNode* methods?",
+			i);
+		return normals_[boundaryMap_[i]];
+
+	}
+	template<typename vec>
+	vec& DomainDiscretization<vec>::normal(int i) {
+		assert_msg(0 <= i && i <= size(), "Index %d out of range [0, %d)", i, size());
+		assert_msg(types_[i] < 0, "Node %d must be a boundary node, got type %d.", i, types_[i]);
+		assert_msg(boundaryMap_[i] != -1, "Node %d does not have a normal. Maybe you manually set"
+			" supports and positions instead of using addInternalNode* methods?",
+			i);
+		return normals_[boundaryMap_[i]];
+
+	}
+
+	template<typename vec>
+	int DomainDiscretization<vec>::addInternalNode(const vec& point, int type) {
+
+		assert_msg(type > 0, "This function is for adding internal points, but got type %d, which is "
+			"not positive. Use addBoundaryNode to add boundary nodes.",
+			type);
+		return addNode(point, type);
+	}
+
+	template<typename vec>
+	int DomainDiscretization<vec>::addInternalNodeWithT(const vec& point, double t, int type) {
+
+		assert_msg(type > 0, "This function is for adding internal points, but got type %d, which is "
+			"not positive. Use addBoundaryNode to add boundary nodes.",
+			type);
+		return addNodeWithT(point, t, type);
+	}
+
+	/**
+	 * Adds a boundary node with given type and normal to the domain.
+	 * @param point Coordinates of the node to add.
+	 * @param type Type of the point, must be negative.
+	 * @param normal Outside unit normal to the boundary at point `point`.
+	 * @return The index of the new node.
+	 * @sa addInternalNode
+	 */
+	template<typename vec>
+	int DomainDiscretization<vec>::addBoundaryNode(const vec& point, int type, const vec& normal) {
+		assert_msg(type < 0, "Type of boundary points must be negative, got %d.", type);
+		int idx = addNode(point, type);
+		boundaryMap_[idx] = normals_.size();
+		normals_.push_back(normal);
+		return idx;
+
+	}
+
+	template<typename vec>
+	int DomainDiscretization<vec>::addBoundaryNodeWithT(const vec& point, double t, int type, const vec& normal) {
+		assert_msg(type < 0, "Type of boundary points must be negative, got %d.", type);
+		int idx = addNodeWithT(point, t, type);
+		boundaryMap_[idx] = normals_.size();
+		normals_.push_back(normal);
+		return idx;
+
+	}
+
+	template<typename vec>
+	int DomainDiscretization<vec>::addNode(const vec& point, int type) {
+		positions_.push_back(point);
+		types_.push_back(type);
+		support_.emplace_back();
+		boundaryMap_.push_back(-1);
+		return positions_.size() - 1;
+	}
+
+	template<typename vec>
+	int DomainDiscretization<vec>::addNodeWithT(const vec& point, double t, int type) {
+		positions_.push_back(point);
+		inCurve_.push_back(t);
+		types_.push_back(type);
+		support_.emplace_back();
+		boundaryMap_.push_back(-1);
+		return positions_.size() - 1;
+	}
+	//bool DomainDiscretization::contains(const vec& point)const;
+
+
+
+
+};

include/DomainDiscretization_fwd.hpp

@@ -0,0 +1,197 @@
+#pragma once
+#include <vector>
+#include <Eigen/Core>
+#include "OccShape.hpp"
+
+namespace  meshless {
+
+
+
+	template<typename vec>
+	class DomainDiscretization {
+	protected:
+		std::vector<vec> positions_;
+		std::vector<int> types_;
+		std::vector<std::vector<int>> support_;
+		std::vector<int> boundaryMap_;
+		std::vector<vec> normals_;
+		OccShape shape_;
+		std::vector<double> inCurve_;
+
+	public:
+		DomainDiscretization() {}
+		DomainDiscretization(const OccShape& shape) {
+			shape_.max_points = shape.max_points;
+			shape_.seed_ = shape.seed_;
+			shape_.n_samples = shape.n_samples;
+			shape_.zeta = shape.zeta;
+			shape_.myShape = shape.myShape;
+
+		}
+		const double getT(int i)const {
+			return inCurve_[i];
+		}
+		const std::vector<vec>& positions()const {
+			return positions_;
+		}
+
+		const vec& pos(int i)const {
+			return positions_[i];
+		}
+
+		vec& pos(int i) {
+			return positions_[i];
+		}
+		double pos(int i, int j) const {
+			return positions_[i][j];
+		}
+
+		const std::vector<std::vector<int> >& supports()const {
+			return support_;
+		}
+
+		const std::vector<int>& support(int i)const {
+			return support_[i];
+		}
+
+		std::vector<int>& support(int i) {
+			return support_[i];
+		}
+
+		int support(int i, int j) const {
+			return support_[i][j];
+		}
+
+		std::vector<vec> supportNodes(int i) const {
+			auto sp = support_[i];
+			std::vector<vec> retsP;
+			for(auto it : sp) {
+				retsP.push_back(positions_[it]);
+
+			}
+			return retsP;
+		}
+		/// Returns position of `j`-th support node of `i`-th node.
+		vec supportNode(int i, int j) const {
+			return positions_[support_[i][j]];
+		}
+		/// Returns Euclidean distance to the second support node.
+		double dr(int i) const {
+			return (positions_[i] - positions_[support_[i][1]]).norm();
+		}
+		/// Returns size of `i`-th node support.
+		int supportSize(int i) const {
+			return support_[i].size();
+		}
+		/// Returns a vector of support sizes for each node.
+		//std::vector<int> supportSizes() const;
+
+		/// Returns types of all nodes.
+		const std::vector<int>& types() const {
+			return types_;
+		}
+		/// Returns mutable types of all nodes.
+		std::vector<int>& types() {
+			return types_;
+		}
+		/// Returns type of `i`-th node.
+		int type(int i) const {
+			return types_[i];
+		}
+		/// Returns writeable type of `i`-th node.
+		int& type(int i) {
+			return types_[i];
+		}
+
+
+		/// Returns boundary map. @sa boundary_map_
+		const std::vector<int>& bmap() const {
+			return boundaryMap_;
+		}
+		/**
+		 * Returns index of node `node` among only boundary nodes. The returned index is
+		 * in range `[0, boundary().size())` if `node` is a boundary node, and `-1` otherwise.
+		 */
+		int bmap(int node) const {
+			return boundaryMap_[node];
+		}
+		/// Returns normals of all boundary nodes.
+		const std::vector<vec>& normals() const {
+			return normals_;
+		}
+		/**
+		 * Returns outside unit normal of `i`-th node. The node must be a boundary node.
+		 * @throw Assertion fails if the noe is not a boundary node, i.e.\ `type(i) < 0` must hold.
+		 */
+		const vec& normal(int i) const;
+		/// Returns writable outside unit normal of `i`-th node. @sa normal
+		vec& normal(int i);
+
+		/// Returns indexes of all boundary nodes.
+		std::vector<int> boundary() const {
+			std::vector<int> ret;
+			for(int i = 0; i < types_.size(); i++) {
+				if(types_[i] < 0) {
+					ret.push_back(i);
+				}
+			}
+			return ret;
+		}
+		/// Returns indexes of all internal nodes.
+		std::vector<int> interior() const {
+			std::vector<int> ret;
+			for(int i = 0; i < types_.size(); i++) {
+				if(types_[i] > 0) {
+					ret.push_back(i);
+				}
+			}
+			return ret;
+		}
+		/// Returns indexes of all nodes, i.e.\ `{0, 1, ..., N-1}`.
+		std::vector<int> all() const {
+			std::vector<int> ret;
+			for(int i = 0; i < types_.size(); i++) {
+				if(types_[i] != 0) {
+					ret.push_back(i);
+				}
+			}
+			return ret;
+		}
+
+		/// Returns `N`, the number of nodes in this discretization.
+		int size() const {
+			return positions_.size();
+		}
+
+		/**
+		 * Adds a single interior node with specified type to this discretization.
+		 * @param point Coordinates of the node to add.
+		 * @param type Type of the node to add. Must be positive.
+		 * @return The index of the new node.
+		 * @sa addBoundaryNode
+		 */
+		int addInternalNode(const vec& point, int type);
+
+		/**
+		 * Adds a boundary node with given type and normal to the domain.
+		 * @param point Coordinates of the node to add.
+		 * @param type Type of the point, must be negative.
+		 * @param normal Outside unit normal to the boundary at point `point`.
+		 * @return The index of the new node.
+		 * @sa addInternalNode
+		 */
+		int addBoundaryNode(const vec& point, int type, const vec& normal);
+
+		int addBoundaryNodeWithT(const vec& point, double t, int type, const vec& normal);
+		int addInternalNodeWithT(const vec& point, double t, int type);
+	private:
+		int addNode(const vec& point, int type);
+		int addNodeWithT(const vec& point, double t, int type);
+
+	public:
+		//bool discreteContains(const vec& point, )const;
+	};
+
+
+
+};

include/FillBoundary.hpp

@@ -0,0 +1,193 @@
+#pragma once
+#include "FillBoundary_fwd.hpp"
+#include "KDTree.hpp"
+#include "KDTreeMutable.hpp"
+#include <iomanip>
+#include <random>
+#include <map>
+
+namespace meshless {
+
+	DomainDiscretization<Eigen::Vector3d> discretizeBoundaryWithDensity(const OccShape& shape, const std::function<double(Eigen::Vector3d)>& h, int type) {
+		if(type == 0)type = -1;
+		std::mt19937 gen(shape.seed_);
+		DomainDiscretization<Eigen::Vector3d> domainGlobal(shape);
+		KDTreeMutable treeGlobal;
+		TopExp_Explorer expFace(shape.myShape, TopAbs_FACE);
+		int numFace = 0;
+
+		std::cout.precision(10);
+
+		for(; expFace.More(); expFace.Next()) {
+			std::cout << "正在离散化第" << ++numFace << "张曲面" << std::endl;
+			//if(numFace != 5)continue;
+			bool reverseNormal = false;
+			const TopoDS_Face& face = TopoDS::Face(expFace.Current());
+			Handle(Geom_Surface) geomSurface = BRep_Tool::Surface(face);
+			if(face.Orientable() == TopAbs_REVERSED) {
+				std::cout << "该面法线需反向\n";
+				reverseNormal = true;
+			}
+
+			//开始离散化每个wire
+			auto outerWire = BRepTools::OuterWire(face);
+			int numWire = 0;
+			TopExp_Explorer expWire(face, TopAbs_WIRE);
+			for(; expWire.More(); expWire.Next()) {
+				std::cout << "正在离散化第" << ++numWire << "个Wire" << std::endl;
+				const TopoDS_Wire& wire = TopoDS::Wire(expWire.Current());
+				if(wire == outerWire) {
+					std::cout << "此wire为外环\n";
+				}
+
+
+				DomainDiscretization<Eigen::Vector2d> domainParamLocalPatch(shape);
+
+				TopExp_Explorer expEdge(wire, TopAbs_EDGE);
+				int numEdge = 0;
+				int genPoints = 0;
+				for(; expEdge.More(); expEdge.Next()) {
+
+					std::cout << "正在离散化第" << ++numEdge << "条Edge" << std::endl;
+					//if(numEdge != 5)continue;
+					const TopoDS_Edge& edge = TopoDS::Edge(expEdge.Current());
+
+					KDTreeMutable treeLocal;
+					Standard_Real firstParam, lastParam;
+					Handle(Geom_Curve) geomCurve = BRep_Tool::Curve(edge, firstParam, lastParam);
+					//std::cout << "param: " << firstParam << " -> " << lastParam << '\n';
+
+					//曲线参数点在曲面参数域上的坐标
+					Standard_Real pCurveFirstParam, pCurveLastParam;
+					Handle(Geom2d_Curve) geomPCurve = BRep_Tool::CurveOnSurface(edge, face, pCurveFirstParam, pCurveLastParam);
+					Handle(Geom2d_TrimmedCurve) geomTrimmedCurve = new Geom2d_TrimmedCurve(geomPCurve, pCurveFirstParam, pCurveLastParam);
+
+					//GeomTools::Dump(geomCurve, std::cout);
+					//std::cout << "In param Pcurve:\n";
+					//GeomTools::Dump(geomPCurve, std::cout);
+					std::cout << "GeomPcurve P:" << pCurveFirstParam << ',' << pCurveLastParam << '\n';
+
+					std::uniform_real_distribution<> dis(pCurveFirstParam, pCurveLastParam);
+					//曲线参数点
+					double tSeedInParam = dis(gen);
+					gp_Pnt2d tSeedInPCurve;
+					geomPCurve->D0(tSeedInParam, tSeedInPCurve);
+
+
+					//std::cout << "t in param: " << tSeedInParam << "\n";
+					//std::cout << "t in PCurve: " << tSeedInPCurve.Coord().X() << "," << tSeedInPCurve.Coord().Y() << "\n";
+
+					//!这里两种计算方法得到的结果不一致，第一种更准确
+					gp_Pnt tSeedInCurve;
+					geomCurve->D0(tSeedInParam, tSeedInCurve);
+					gp_Pnt tSeedInSurf;
+					geomSurface->D0(tSeedInPCurve.X(), tSeedInPCurve.Y(), tSeedInSurf);
+					//! Tip: use GeomLib::NormEstim() to calculate surface normal at specified (U, V) point.
+
+					std::cout << "t in Surface:#1 " << tSeedInCurve.Coord().X() << "," << tSeedInCurve.Coord().Y() << "," << tSeedInCurve.Coord().Z() << "\n";
+					std::cout << "t in Surface:#2 " << tSeedInSurf.Coord().X() << "," << tSeedInSurf.Coord().Y() << "," << tSeedInSurf.Coord().Z() << "\n";
+
+					Eigen::Vector3d eigenSeedPnt(tSeedInCurve.X(), tSeedInCurve.Y(), tSeedInCurve.Z());
+					Eigen::Vector2d eigenSeedPntParam(tSeedInPCurve.X(), tSeedInPCurve.Y());
+					domainParamLocalPatch.addInternalNodeWithT(transPnt2d(tSeedInPCurve), tSeedInParam, 1);
+					treeLocal.insert(transPnt(tSeedInCurve));
+
+
+					//Insert into global structures.
+					double checkRadiusSeed = h(eigenSeedPnt);
+
+					double d_sq = treeGlobal.size() == 0 ? 10 * checkRadiusSeed * checkRadiusSeed : treeGlobal.query(eigenSeedPnt).second[0];
+					if(d_sq >= (shape.zeta * checkRadiusSeed) * (shape.zeta * checkRadiusSeed)) {
+						//domainGlobal.addBoundaryNode(eigenSeedPnt, -1);
+						gp_Dir tSeedNormal;
+						auto retStatus = GeomLib::NormEstim(geomSurface, tSeedInPCurve, 1e-6, tSeedNormal);
+						if(retStatus >= 2) {
+							std::cout << "calculate wrong!\n";
+							exit(-1);
+						}
+						if(reverseNormal) tSeedNormal = -tSeedNormal;
+						Eigen::Vector3d eigenNormal(tSeedNormal.X(), tSeedNormal.Y(), tSeedNormal.Z());
+						domainGlobal.addBoundaryNode(eigenSeedPnt, type, eigenNormal);
+						treeGlobal.insert(eigenSeedPnt);
+					}
+
+					int curNode = domainParamLocalPatch.size() - 1;
+					int endNode = domainParamLocalPatch.size();
+					while(curNode < endNode && endNode < shape.max_points) {
+						genPoints++;
+						//std::cout << "curNode / endNode : " << curNode << '/' << endNode << '\n';
+						auto param = domainParamLocalPatch.pos(curNode);
+						double t = domainParamLocalPatch.getT(curNode);
+						std::vector<double> candidates{ 1.0, -1.0 };
+
+						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';
+						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";
+								continue;
+							}
+							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';
+							double checkradius = pt.SquareDistance(ptNew);
+							std::cout << "CheckRadius: " << checkradius << '\n';
+							double d_sq_loc = treeLocal.query(transPnt(ptNew)).second[0];
+							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));
+								endNode++;
+								double d_sq_global = treeGlobal.query(transPnt(ptNew)).second[0];
+								if(d_sq_global >= (shape.zeta * shape.zeta * checkradius)) {
+									gp_Dir tSeedNormalNew;
+									auto retStatus = GeomLib::NormEstim(geomSurface, tNewSeedInPCurve, 1e-6, tSeedNormalNew);
+									if(retStatus >= 2) {
+										std::cout << "calculate wrong!\n";
+										exit(-1);
+									}
+									if(reverseNormal) {
+										tSeedNormalNew = -tSeedNormalNew;
+									}
+									domainGlobal.addBoundaryNode(transPnt(ptNew), type, transVec(tSeedNormalNew));
+									treeGlobal.insert(transPnt(ptNew));
+								} else {
+									std::cout << "全局不添加点\n";
+								}
+							} else {
+								std::cout << "局部不添加点\n";
+							}
+						}
+						curNode++;
+					}
+					std::cout << "该边已经生成：" << genPoints << "个点\n";
+				}
+			}
+
+			//fill boundary surface!
+			KDTree paramBoundarySearch;
+
+
+
+
+
+
+		}
+		return domainGlobal;
+	}
+
+};

include/FillBoundary_fwd.hpp

@@ -0,0 +1,23 @@
+#pragma once
+#include "OccHelper.hpp"
+#include "OccShape.hpp"
+#include "DomainDiscretization.hpp"
+#include <vector>
+#include <Eigen/Core>
+
+namespace meshless {
+
+
+	DomainDiscretization<Eigen::Vector3d> discretizeBoundaryWithDensity(const OccShape& shape, const std::function<double(Eigen::Vector3d)>& h, int type);
+
+	DomainDiscretization<Eigen::Vector3d> discretizeBoundaryWithStep(const OccShape& shape, double step, int type) {
+		auto f = [=](Eigen::Vector3d) {return step; };
+		return discretizeBoundaryWithDensity(shape, f, type);
+
+	}
+
+
+
+
+
+};

include/KDTree.hpp

@@ -0,0 +1,30 @@
+#pragma once
+#include <KDTree_fwd.hpp>
+#include <assert.hpp>
+
+namespace meshless {
+	std::pair<std::vector<int>, std::vector<double>> KDTree::query(const Eigen::Vector3d& point, int k)const {
+
+		assert_msg(point.array().isFinite().prod() == 1, "Invalid point.");
+		std::vector<int> ret_index(k);
+		std::vector<double> out_dist_sqr(k);
+		int actual_k = tree.knnSearch(point.data(), k, &ret_index[0], &out_dist_sqr[0]);
+		assert_msg(actual_k == k, "There were not enough points in the tree, you requested %d "
+			"points, the tree only contains %d points.", k, actual_k);
+		return { ret_index, out_dist_sqr };
+	}
+
+	std::pair<std::vector<int>, std::vector<double>> KDTree::query(const Eigen::Vector3d& point, const double& radius_squared)const {
+		assert_msg(point.array().isFinite().prod() == 1, "Invalid point.");
+		std::vector<std::pair<int, double>> idx_dist;
+		int k = tree.radiusSearch(point.data(), radius_squared, idx_dist, nanoflann::SearchParams());
+		std::vector<int> idx(k); std::vector<double> dists(k);
+		for(int i = 0; i < k; ++i) {
+			std::tie(idx[i], dists[i]) = idx_dist[i];
+		}
+		return { idx, dists };
+	}
+
+
+
+};

include/KDTreeMutable.hpp

@@ -0,0 +1,40 @@
+#pragma once
+#include "assert.hpp"
+#include "KDTreeMutable_fwd.hpp"
+
+namespace meshless {
+
+	void KDTreeMutable::insert(const Eigen::Vector3d& point) {
+		assert_msg(point.array().isFinite().prod() == 1, "Invalid point.");
+		auto n = points_.kdtree_get_point_count();
+		points_.add(point);
+		tree.addPoints(n, n);
+		++size_;
+
+	}
+
+	void KDTreeMutable::insert(const std::vector<Eigen::Vector3d>& points) {
+		auto n = points_.kdtree_get_point_count();
+		for(const auto& p : points) {
+			assert_msg(p.array().isFinite().prod() == 1, "One of the points is invalid.");
+			points_.add(p);
+		}
+		size_ += points.size();
+		tree.addPoints(n, n + points.size() - 1);
+
+	}
+
+	std::pair<std::vector<int>, std::vector<double>> KDTreeMutable::query(const Eigen::Vector3d& point, int k) {
+		assert_msg(point.array().isFinite().prod() == 1, "Invalid query point.");
+		nanoflann::KNNResultSet<double, int> resultSet(k);
+		std::vector<int> ret_index(k);
+		std::vector<double> out_dist_sqr(k);
+		resultSet.init(&ret_index[0], &out_dist_sqr[0]);
+		tree.findNeighbors(resultSet, point.data(), nanoflann::SearchParams(k));
+		assert_msg(resultSet.full(), "Not enough points in the tree, you requested %d points, "
+			"but the tree contains only %d points.",
+			k, size());
+		return { ret_index, out_dist_sqr };
+	}
+
+};

include/KDTreeMutable_fwd.hpp

@@ -0,0 +1,69 @@
+#pragma once
+#include "nanoflann.hpp"
+#include <array>
+#include <iosfwd>
+#include "PointCloud.hpp"
+
+namespace meshless {
+	class KDTreeMutable {
+
+	private:
+		typedef nanoflann::KDTreeSingleIndexDynamicAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloud>, PointCloud, 3, int> kd_tree_t;
+		int size_;
+		PointCloud points_;
+		kd_tree_t tree;
+
+	public:
+		explicit KDTreeMutable(const std::vector<Eigen::Vector3d>& points) :points_(points), size_(points.size()), tree(3, points_, nanoflann::KDTreeSingleIndexAdaptorParams(20)) {}
+
+		KDTreeMutable() :points_(), size_(0), tree(3, points_, nanoflann::KDTreeSingleIndexAdaptorParams(20)) {}
+
+		void reset(const std::vector<Eigen::Vector3d>& points) {
+			points_.setPts(points);
+			size_ = points.size();
+			tree.reset();
+
+		}
+
+		void insert(const Eigen::Vector3d& point);
+		void insert(const std::vector<Eigen::Vector3d>& points);
+
+		/// Check if any point exists in sphere centered at `p` with radius `r`.
+		bool existsPointInSphere(const Eigen::Vector3d& p, double r) {
+			if(size_ == 0)
+				return false;
+			return query(p).second[0] <= r * r;
+		}
+
+		/**
+		 * Removes a point with given index from the tree. The indexes of the points are given
+		 * sequentially at insertion and do not change. The removal is lazy and point still takes
+		 * up memory.
+		 */
+		void remove(int index) {
+			size_ -= tree.removePoint(index);
+		}
+
+		/// Returns number of points in the tree.
+		int size() const {
+			return size_;
+		}
+
+		/**
+		 * Find `k` nearest neighbors to given point.
+		 * @param point Find closest points to this point.
+		 * @param k How many nearest points to find.
+		 *
+		 * @return A pair of two vectors of size `k` containing
+		 * indices of nearest neighbors and squared distances to said neighbors.
+		 * @throw Assertion fails if there are not enough points in the tree.
+		 */
+		std::pair<std::vector<int>, std::vector<double>> query(const Eigen::Vector3d& point, int k = 1);
+
+
+
+	};
+
+
+
+};

include/KDTree_fwd.hpp

@@ -0,0 +1,49 @@
+#pragma once
+#include <PointCloud.hpp>
+#include <nanoflann.hpp>
+#include <Eigen/Core>
+
+namespace meshless {
+	class KDTree {
+	private:
+		typedef nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloud>, PointCloud, 3, int> kd_tree_t;
+		PointCloud points_;
+		kd_tree_t tree;
+	public:
+		explicit KDTree(const std::vector<Eigen::Vector3d>& points) :points_(points), tree(3, points_, nanoflann::KDTreeSingleIndexAdaptorParams(20)) {
+			tree.buildIndex();
+		}
+
+		KDTree() :points_(), tree(3, points_, nanoflann::KDTreeSingleIndexAdaptorParams(20)) {}
+
+		void reset(const std::vector<Eigen::Vector3d>& points) {
+			points_.setPts(points);
+			tree.buildIndex();
+		}
+		std::pair<std::vector<int>, std::vector<double>> query(const Eigen::Vector3d& point, int k = 1)const;
+
+		std::pair<std::vector<int>, std::vector<double>> query(const Eigen::Vector3d& point, const double& radius_squared) const;
+
+		Eigen::Vector3d get(int index) const {
+			return points_.get(index);
+		}
+
+		std::vector<Eigen::Vector3d> get(const std::vector<int>& indexes) const {
+			const int n = indexes.size();
+			std::vector<Eigen::Vector3d> result(n);
+			for(int i = 0; i < n; i++) {
+				result[i] = points_.get(indexes[i]);
+			}
+			return result;
+		}
+
+		int size() const {
+			return points_.kdtree_get_point_count();
+		}
+
+
+	};
+
+
+
+};

include/OccHelper.hpp

@@ -0,0 +1,65 @@
+#pragma once
+
+// OpenCascade includes
+#include <Geom2d_TrimmedCurve.hxx>
+#include <BRepClass_FaceClassifier.hxx>
+#include <Interface_Static.hxx>
+#include <STEPCAFControl_Writer.hxx>
+#include <STEPControl_Reader.hxx>
+#include <TopExp.hxx>
+#include <TopoDS.hxx>
+#include <TopExp_Explorer.hxx>
+#include <iostream>
+#include <TopoDS_Shape.hxx>
+#include <Geom_Surface.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBuilderAPI_NurbsConvert.hxx>
+#include <Brep_Tool.hxx>
+#include <GeomConvert.hxx>
+#include <BRepLib_FindSurface.hxx>
+#include <Geom_BSplineSurface.hxx>
+#include <BRepAdaptor_CompCurve.hxx>
+#include <BRepTools.hxx>
+#include <ShapeFix_Wire.hxx>
+#include <Adaptor3d_Curve.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <BOPTools_AlgoTools3D.hxx>
+#include <BRepBuilderAPI_Sewing.hxx>
+#include <ShapeFix_Shape.hxx>
+#include <ShapeFix_Shell.hxx>
+#include <ShapeAnalysis.hxx>
+#include <TopoDS_TShape.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_CylindricalSurface.hxx>
+#include <Geom_ConicalSurface.hxx>
+#include <Geom_SphericalSurface.hxx>
+#include <Geom_ToroidalSurface.hxx>
+#include <Geom_SurfaceOfLinearExtrusion.hxx>
+#include <Geom_SurfaceOfRevolution.hxx>
+#include <Geom_BezierSurface.hxx>
+#include <Geom_BSplineSurface.hxx>
+#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_OffsetSurface.hxx>
+#include <TopoDS_Compound.hxx>
+#include <BRepBuilderAPI_Copy.hxx>
+#include <TColgp_Array2OfPnt.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <TColStd_Array2OfReal.hxx>
+#include <TColStd_Array1OfInteger.hxx>
+#include <BRepBuilderAPI_MakeSolid.hxx>
+#include <GeomTools.hxx>
+#include <Geom_Circle.hxx>
+#include <BRep_Builder.hxx>
+#include <GeomLib.hxx>
+
+Eigen::Vector3d transVec(const gp_Vec& vec) {
+	return Eigen::Vector3d(vec.X(), vec.Y(), vec.Z());
+}
+
+Eigen::Vector3d transPnt(const gp_Pnt& Pnt) {
+	return Eigen::Vector3d(Pnt.X(), Pnt.Y(), Pnt.Z());
+}
+
+Eigen::Vector2d transPnt2d(const gp_Pnt2d& Pnt) {
+	return Eigen::Vector2d(Pnt.X(), Pnt.Y());
+}

include/OccShape.hpp

@@ -0,0 +1,57 @@
+#pragma once
+#include "OccHelper.hpp"
+#include "assert.hpp"
+
+namespace meshless {
+	class OccShape {
+	public:
+		int max_points = 500000;
+		int seed_;
+		int n_samples = 15;
+		double zeta = 1 - 1e-10;
+		double epsilon = 0;
+		TopoDS_Shape myShape;
+	public:
+		OccShape() {}
+		OccShape(const TopoDS_Shape& myShape_) {
+			max_points = 500000;
+			seed_;
+			n_samples = 15;
+			zeta = 1 - 1e-10;
+			epsilon = 0;
+			myShape = myShape_;
+		}
+		OccShape& maxPoints(int max_points) {
+			this->max_points = max_points;
+			return *this;
+		}
+
+		OccShape& seed(int seed) {
+			seed_ = seed;
+			return *this;
+		}
+
+		OccShape& numSamples(int n_samples) {
+			this->n_samples = n_samples;
+			return *this;
+		}
+
+
+		OccShape& proximityTolerance(double zeta) {
+			assert_msg((0 < zeta && zeta < 1), "Zeta must be between 0 and 1, got %f.", zeta);
+			this->zeta = zeta;
+			return *this;
+		}
+
+		OccShape& boundaryProximity(double epsilon) {
+			assert_msg((0 < epsilon && epsilon < 1), "Epsilon must be between 0 and 1, got %f.", epsilon);
+			this->epsilon = epsilon;
+			return *this;
+		}
+
+
+	};
+
+
+
+};

include/PointCloud.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+/**
+ * @file
+ * Implementation of KDTree storage class.
+ */
+
+#include <vector>
+#include <Eigen/Core>
+
+namespace meshless {
+
+	/// Helper class for KDTree with appropriate accessors containing a set of points. @ingroup utils
+	struct PointCloud {
+		std::vector<Eigen::Vector3d> pts; ///< Points, contained in the tree.
+
+		/// Construct an empty point set.
+		PointCloud() = default;
+
+		/// Construct from an array of points.
+		PointCloud(const std::vector<Eigen::Vector3d>& pts) : pts(pts) {}
+
+		/// Reset contained points.
+		void setPts(const std::vector<Eigen::Vector3d>& pts) {
+			PointCloud::pts = pts;
+		}
+
+		/// Interface requirement: returns number of data points.
+		inline int kdtree_get_point_count() const {
+			return pts.size();
+		}
+
+		/// Interface requirement: returns `dim`-th coordinate of `idx`-th point.
+		inline typename double kdtree_get_pt(const size_t idx, int dim) const {
+			return pts[idx][dim];
+		}
+
+		/// Access the points.
+		inline Eigen::Vector3d get(const size_t idx) const {
+			return pts[idx];
+		}
+
+		/// Add a point to the cloud.
+		inline void add(const Eigen::Vector3d& p) {
+			pts.push_back(p);
+		}
+
+		/// Comply with the interface.
+		template <class BBOX>
+		bool kdtree_get_bbox(BBOX& /* bb */) const {
+			return false;
+		}
+	};
+
+} // namespace mm
+

include/assert.hpp

@@ -0,0 +1,110 @@
+#pragma once
+
+/**
+ * @file
+ * Implementation of custom assert and debug utilities.
+ */
+
+#include <tinyformat.h>
+#include "print.hpp"
+
+namespace meshless {
+
+#ifdef _MSC_VER // Check if using MSVC
+#define FUNCTION_NAME __FUNCSIG__ // Use __FUNCSIG__ for function name
+#else
+#define FUNCTION_NAME __PRETTY_FUNCTION__ // Use __PRETTY_FUNCTION__ for other compilers
+#endif
+
+
+	// print macro
+	/// @cond
+#define VA_NUM_ARGS(...) VA_NUM_ARGS_IMPL(__VA_ARGS__, 5, 4, 3, 2, 1, 0)
+#define VA_NUM_ARGS_IMPL(_1, _2, _3, _4, _5, N, ...) N
+#define macro_dispatcher(func, ...)     macro_dispatcher_(func, VA_NUM_ARGS(__VA_ARGS__))
+#define macro_dispatcher_(func, nargs)  macro_dispatcher__(func, nargs)
+#define macro_dispatcher__(func, nargs) func ## nargs
+#define addflag(a) {std::cerr << "flags=[flags, " << (a) << "];" << std::endl;}
+#define prnv2(a, b) {std::cerr << a << " = " << (b) << ";" << std::endl;}
+#define prnv1(a)   {std::cerr << #a << " = " << (a) << ";" << std::endl;}
+/// @endcond
+/**
+ * Prints a variable name and value to standard output. Can take one or two parameters.
+ * Example:
+ * @code
+ * int a = 6;
+ * prn(a) // prints 'a = 6;'
+ * prn("value", a) // prints 'value = 6;'
+ * @endcode
+ */
+#define prn(...) macro_dispatcher(prnv, __VA_ARGS__)(__VA_ARGS__)
+
+	using tinyformat::printf;
+	using tinyformat::format;
+
+	/// Namespace holding custom assert implementation.
+	namespace assert_internal {
+		/**
+		 * Actual assert implementation.
+		 * @param condition Condition to test, e.g.\ `n > 0`.
+		 * @param file File where the assertion failed.
+		 * @param func_name Function name where the assertion failed.
+		 * @param line Line on which the assertion failed.
+		 * @param message Message as specified in the `assert_msg` macro.
+		 * @param format_list List of format field values to pass to `tinyformat::format` function.
+		 */
+		bool assert_handler_implementation(const char* condition, const char* file, const char* func_name,
+			int line, const char* message, tfm::FormatListRef format_list);
+		/// Assert handler that unpacks varargs.
+		template<typename... Args>
+		bool assert_handler(const char* condition, const char* file, const char* func_name, int line,
+			const char* message, const Args&... args) {  // unpacks first argument
+			tfm::FormatListRef arg_list = tfm::makeFormatList(args...);
+			return assert_handler_implementation(condition, file, func_name, line, message, arg_list);
+		}
+	}  // namespace assert_internal
+
+#ifdef NDEBUG
+#define assert_msg(cond, ...) ((void)sizeof(cond))
+#else
+/**
+ * @brief Assert with better error reporting.
+ * @param cond Conditions to test.
+ * @param ... The second parameter is also required and represents the message to print on failure.
+ * For every %* field in message one additional parameter must be present.
+ *
+ * Example:
+ * @code
+ * assert_msg(n > 0, "n must be positive, got %d.", n);
+ * @endcode
+ */
+#define assert_msg(cond, ...) ((void)(!(cond) && \
+    meshless::assert_internal::assert_handler( \
+            #cond, __FILE__, FUNCTION_NAME, __LINE__, __VA_ARGS__) && (assert(0), 1)))
+ //            #cond, __FILE__, __PRETTY_FUNCTION__, __LINE__, __VA_ARGS__) && (exit(1), 1)))
+#endif
+
+/**
+ * Prints given text in bold red.
+ * @param s text to print.
+ */
+	inline void print_red(const std::string& s) {
+		std::cout << "\x1b[31;1m" << s << "\x1b[37;0m";
+	}
+	/**
+	 * Prints given text in bold white.
+	 * @param s text to print.
+	 */
+	inline void print_white(const std::string& s) {
+		std::cout << "\x1b[37;1m" << s << "\x1b[37;0m";
+	}
+	/**
+	 * Prints given text in bold green.
+	 * @param s text to print.
+	 */
+	inline void print_green(const std::string& s) {
+		std::cout << "\x1b[32;1m" << s << "\x1b[37;0m";
+	}
+
+}  // namespace meshless
+

include/print.hpp

@@ -0,0 +1,103 @@
+#ifndef MEDUSA_BITS_UTILS_PRINT_HPP_
+#define MEDUSA_BITS_UTILS_PRINT_HPP_
+
+/**
+ * @file
+ * Printing helpers for std types.
+ */
+
+#include <iostream>
+#include <vector>
+#include <array>
+#include <utility>
+#include <tuple>
+
+ // additional ostream operators
+namespace std {
+
+	/// Output pairs as `(1, 2)`. @ingroup utils
+	template <class T, class U>
+	std::ostream& operator<<(std::ostream& xx, const std::pair<T, U>& par) {
+		return xx << "(" << par.first << "," << par.second << ")";
+	}
+
+	/// Output arrays as `[1, 2, 3]`. @ingroup utils
+	template <class T, size_t N>
+	std::ostream& operator<<(std::ostream& xx, const std::array<T, N>& arr) {
+		xx << "[";
+		for(size_t i = 0; i < N; ++i) {
+			xx << arr[i];
+			if(i < N - 1)
+				xx << ", ";
+		}
+		xx << "]";
+		return xx;
+	}
+
+	/// Output vectors as `[1, 2, 3]`. @ingroup utils
+	template <class T, class A>
+	std::ostream& operator<<(std::ostream& xx, const std::vector<T, A>& arr) {
+		// do it like the matlab does it.
+		xx << "[";
+		for(size_t i = 0; i < arr.size(); ++i) {
+			xx << arr[i];
+			if(i < arr.size() - 1)
+				xx << ", ";
+		}
+		xx << "]";
+		return xx;
+	}
+
+	/// Output nested vectors as `[[1, 2]; [3, 4]]`. @ingroup utils
+	template <class T, class A>
+	std::ostream& operator<<(std::ostream& xx, const std::vector<std::vector<T, A>>& arr) {
+		xx << "[";
+		for(size_t i = 0; i < arr.size(); ++i) {
+			for(size_t j = 0; j < arr[i].size(); ++j) {
+				xx << arr[i][j];
+				if(j < arr[i].size() - 1)
+					xx << ", ";
+			}
+			if(i < arr.size() - 1)
+				xx << "; ";
+		}
+		xx << "]";
+		return xx;
+	}
+
+	/// @cond
+	namespace tuple_print_internal {
+		template <class Tuple, std::size_t N>
+		struct TuplePrinter {
+			static void print(std::ostream& os, const Tuple& t) { // recursive
+				TuplePrinter<Tuple, N - 1>::print(os, t);
+				os << ", " << std::get<N - 1>(t);
+			}
+		};
+
+		template <class Tuple>
+		struct TuplePrinter<Tuple, 1> {
+			static void print(std::ostream& os, const Tuple& t) { // one element
+				os << std::get<0>(t);
+			}
+		};
+
+		template <class Tuple>
+		struct TuplePrinter<Tuple, 0> {
+			static void print(std::ostream&, const Tuple&) {}
+		}; // zero elt
+
+	} // namespace tuple_print_internal
+	/// @endcond
+
+	/// Print a tuple as (1, 4.5, abc). @ingroup utils
+	template <class... Args>
+	std::ostream& operator<<(std::ostream& os, const std::tuple<Args...>& t) {
+		os << "(";
+		tuple_print_internal::TuplePrinter<decltype(t), sizeof...(Args)>::print(os, t);
+		return os << ")";
+	}
+
+} // namespace std
+
+#endif // MEDUSA_BITS_UTILS_PRINT_HPP_

include/writeVTK.hpp

@@ -0,0 +1,54 @@
+#pragma once
+#include <Eigen/Core>
+#include <cstring>
+namespace meshless {
+	template <int dim>
+	void writePntVTK(const std::string& path, const Eigen::MatrixXd& pnts, const std::vector<double>& attri) {
+		std::ofstream out(path);
+		out << "# vtk DataFile Version 3.0\n"
+			"Volume Mesh\n"
+			"ASCII\n"
+			"DATASET UNSTRUCTURED_GRID"
+			<< std::endl;
+		out << "POINTS " << pnts.rows() << " float" << std::endl;
+		for(int i = 0; i < pnts.rows(); ++i) {
+			for(int j = 0; j < dim; ++j) {
+				out << std::setprecision(4) << pnts(i, j) << " ";
+			}
+			for(int j = dim; j < 3; ++j) {
+				out << "0 ";
+			}
+			out << std::endl;
+		}
+		int innersize = 0;
+		int interSize = 0;
+		int bounSize = 0;
+
+		out << "CELLS " << pnts.rows() << " " << pnts.rows() * (1 + 1) << std::endl;
+		for(int i = 0; i < pnts.rows(); ++i) {
+			out << "1 " << i << std::endl;
+		}
+		out << "CELL_TYPES " << pnts.rows() << std::endl;
+		for(int i = 0; i < pnts.rows(); ++i) {
+			out << 1 << std::endl;
+		}
+
+		if(!attri.empty()) {
+			out << "POINT_DATA " << attri.size() << "\n"
+				<< "SCALARS point_scalars double 1\n"
+				<< "LOOKUP_TABLE default" << std::endl;
+			for(auto& d : attri) {
+				out << d << std::endl;
+				if(d == 1)
+					innersize++;
+				else
+					bounSize++;
+			}
+		}
+		std::cout << "ÄÚ²¿µã£º " << innersize << "\n";
+		std::cout << "±ß½çµã£º " << bounSize << "\n";
+	}
+
+
+
+};

src/assert.cpp

@@ -0,0 +1,26 @@
+#include "assert.hpp"
+
+/**
+ * @file
+ * Implementation of custom assert utilities.
+ */
+
+namespace meshless {
+
+	namespace assert_internal {
+
+		bool assert_handler_implementation(const char* condition, const char* file, const char* func_name,
+			int line, const char* message, tfm::FormatListRef format_list) {
+			std::cerr << "\x1b[37;1m";  // white bold
+			tfm::format(std::cerr, "%s:%d: %s: Assertion `%s' failed with message:\n",
+				file, line, func_name, condition);
+			std::cerr << "\x1b[31;1m";  // red bold
+			tfm::vformat(std::cerr, message, format_list);
+			std::cerr << "\x1b[37;0m\n";  // no color
+			std::cerr.flush();
+			return true;
+		}
+
+	}  // namespace assert_internal
+
+}  // namespace meshless

src/main.cpp

@@ -0,0 +1,38 @@
+#include <iostream>
+#include <fstream>
+#include <Eigen/Core>
+#include "FillBoundary.hpp"
+#include "writeVTK.hpp"
+
+using namespace std;
+using namespace meshless;
+
+int main() {
+	Eigen::MatrixXd mat;
+	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);
+	OccShape occShape(myShape);
+	auto domain = discretizeBoundaryWithStep(occShape, 1, -1);
+	Eigen::MatrixXd quaPoints;
+
+	quaPoints.resize(domain.positions().size(), 3);
+	std::vector<double> attri;
+	for(int i = 0; i < domain.positions().size(); i++) {
+		if(domain.types()[i] == 1) {
+			attri.push_back(1.0);
+		} else {
+			attri.push_back(0);
+		}
+		quaPoints(i, 0) = domain.positions()[i].x();
+		quaPoints(i, 1) = domain.positions()[i].y();
+		quaPoints(i, 2) = domain.positions()[i].z();
+	}
+	writePntVTK<3>("Exhaust_Manifold.vtk", quaPoints, attri);
+
+
+
+
+	return 0;
+}