diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp index 602a67e1498b4bb591801ca61fccb969a08c3c0e..44624dd8727bfce03f1a1f3015fbe433fae5cc0f 100644 --- a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp +++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp @@ -28,63 +28,63 @@ namespace MeshGeoToolsLib std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines( const MeshLib::Mesh& mesh, const GeoLib::PolylineVec& ply_vec) { - // copy existing nodes and elements - std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(mesh.getNodes()); - std::vector<MeshLib::Element*> vec_new_eles = MeshLib::copyElementVector(mesh.getElements(), vec_new_nodes); + // copy existing nodes and elements + std::vector<MeshLib::Node*> vec_new_nodes = MeshLib::copyNodeVector(mesh.getNodes()); + std::vector<MeshLib::Element*> vec_new_eles = MeshLib::copyElementVector(mesh.getElements(), vec_new_nodes); - std::vector<int> new_mat_ids; - { - auto mat_ids = mesh.getProperties().getPropertyVector<int>("MaterialIDs"); - if (mat_ids) { - new_mat_ids.reserve((*mat_ids).size()); - std::copy((*mat_ids).cbegin(), (*mat_ids).cend(), - std::back_inserter(new_mat_ids)); - } - } - int max_matID(0); - if (!new_mat_ids.empty()) - max_matID = *(std::max_element(new_mat_ids.cbegin(), new_mat_ids.cend())); + std::vector<int> new_mat_ids; + { + auto mat_ids = mesh.getProperties().getPropertyVector<int>("MaterialIDs"); + if (mat_ids) { + new_mat_ids.reserve((*mat_ids).size()); + std::copy((*mat_ids).cbegin(), (*mat_ids).cend(), + std::back_inserter(new_mat_ids)); + } + } + int max_matID(0); + if (!new_mat_ids.empty()) + max_matID = *(std::max_element(new_mat_ids.cbegin(), new_mat_ids.cend())); - const std::size_t n_ply (ply_vec.size()); - // for each polyline - for (std::size_t k(0); k < n_ply; k++) - { - const GeoLib::Polyline* ply = (*ply_vec.getVector())[k]; + const std::size_t n_ply (ply_vec.size()); + // for each polyline + for (std::size_t k(0); k < n_ply; k++) + { + const GeoLib::Polyline* ply = (*ply_vec.getVector())[k]; - // search nodes on the polyline - MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly(mesh, *ply, mesh.getMinEdgeLength()*0.5); - auto &vec_nodes_on_ply = mshNodesAlongPoly.getNodeIDs(); - if (vec_nodes_on_ply.empty()) { - std::string ply_name; - ply_vec.getNameOfElementByID(k, ply_name); - INFO("No nodes found on polyline %s", ply_name.c_str()); - continue; - } + // search nodes on the polyline + MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly(mesh, *ply, mesh.getMinEdgeLength()*0.5); + auto &vec_nodes_on_ply = mshNodesAlongPoly.getNodeIDs(); + if (vec_nodes_on_ply.empty()) { + std::string ply_name; + ply_vec.getNameOfElementByID(k, ply_name); + INFO("No nodes found on polyline %s", ply_name.c_str()); + continue; + } - // add line elements - for (std::size_t i=0; i<vec_nodes_on_ply.size()-1; i++) { - std::array<MeshLib::Node*, 2> element_nodes; - element_nodes[0] = vec_new_nodes[vec_nodes_on_ply[i]]; - element_nodes[1] = vec_new_nodes[vec_nodes_on_ply[i+1]]; - vec_new_eles.push_back( - new MeshLib::Line(element_nodes, vec_new_eles.size())); - new_mat_ids.push_back(max_matID+k+1); - } - } + // add line elements + for (std::size_t i=0; i<vec_nodes_on_ply.size()-1; i++) { + std::array<MeshLib::Node*, 2> element_nodes; + element_nodes[0] = vec_new_nodes[vec_nodes_on_ply[i]]; + element_nodes[1] = vec_new_nodes[vec_nodes_on_ply[i+1]]; + vec_new_eles.push_back( + new MeshLib::Line(element_nodes, vec_new_eles.size())); + new_mat_ids.push_back(max_matID+k+1); + } + } - // generate a mesh - const std::string name = mesh.getName() + "_with_lines"; - std::unique_ptr<MeshLib::Mesh> new_mesh( - new MeshLib::Mesh(name, vec_new_nodes, vec_new_eles)); - auto opt_mat_pv = new_mesh->getProperties().createNewPropertyVector<int>( - "MaterialIDs", MeshLib::MeshItemType::Cell); - if (opt_mat_pv) { - auto & mat_pv = *opt_mat_pv; - mat_pv.reserve(new_mat_ids.size()); - std::copy(new_mat_ids.cbegin(), new_mat_ids.cend(), - std::back_inserter(mat_pv)); - } - return new_mesh; + // generate a mesh + const std::string name = mesh.getName() + "_with_lines"; + std::unique_ptr<MeshLib::Mesh> new_mesh( + new MeshLib::Mesh(name, vec_new_nodes, vec_new_eles)); + auto opt_mat_pv = new_mesh->getProperties().createNewPropertyVector<int>( + "MaterialIDs", MeshLib::MeshItemType::Cell); + if (opt_mat_pv) { + auto & mat_pv = *opt_mat_pv; + mat_pv.reserve(new_mat_ids.size()); + std::copy(new_mat_ids.cbegin(), new_mat_ids.cend(), + std::back_inserter(mat_pv)); + } + return new_mesh; } } // MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp index 5d24cc43445c6f02ccb9eb74f167ffafa9cb28c8..3d76ba00bd409cc35d0aa349fe3a5550e0522c75 100644 --- a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp +++ b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp @@ -26,77 +26,77 @@ namespace MeshGeoToolsLib BoundaryElementsAlongPolyline::BoundaryElementsAlongPolyline(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Polyline const& ply) : _mesh(mesh), _ply(ply) { - // search nodes and elements located along the polyline - auto node_ids_on_poly = mshNodeSearcher.getMeshNodeIDsAlongPolyline(ply); - MeshLib::ElementSearch es(_mesh); - es.searchByNodeIDs(node_ids_on_poly); - auto &ele_ids_near_ply = es.getSearchedElementIDs(); + // search nodes and elements located along the polyline + auto node_ids_on_poly = mshNodeSearcher.getMeshNodeIDsAlongPolyline(ply); + MeshLib::ElementSearch es(_mesh); + es.searchByNodeIDs(node_ids_on_poly); + auto &ele_ids_near_ply = es.getSearchedElementIDs(); - // check all edges of the elements near the polyline - for (auto ele_id : ele_ids_near_ply) { - auto* e = _mesh.getElement(ele_id); - // skip internal elements - if (!e->isBoundaryElement()) - continue; - // find edges on the polyline - for (unsigned i=0; i<e->getNEdges(); i++) { - auto* edge = e->getEdge(i); - // check if all edge nodes are along the polyline (if yes, store a distance) - std::vector<std::size_t> edge_node_distances_along_ply; - if (includesAllEdgeNodeIDs(node_ids_on_poly, *edge, edge_node_distances_along_ply)) { - auto* new_edge = modifyEdgeNodeOrdering(*edge, ply, edge_node_distances_along_ply, node_ids_on_poly); - if (edge != new_edge) - delete edge; - _boundary_elements.push_back(new_edge); - } else { - delete edge; - } - } - } + // check all edges of the elements near the polyline + for (auto ele_id : ele_ids_near_ply) { + auto* e = _mesh.getElement(ele_id); + // skip internal elements + if (!e->isBoundaryElement()) + continue; + // find edges on the polyline + for (unsigned i=0; i<e->getNEdges(); i++) { + auto* edge = e->getEdge(i); + // check if all edge nodes are along the polyline (if yes, store a distance) + std::vector<std::size_t> edge_node_distances_along_ply; + if (includesAllEdgeNodeIDs(node_ids_on_poly, *edge, edge_node_distances_along_ply)) { + auto* new_edge = modifyEdgeNodeOrdering(*edge, ply, edge_node_distances_along_ply, node_ids_on_poly); + if (edge != new_edge) + delete edge; + _boundary_elements.push_back(new_edge); + } else { + delete edge; + } + } + } - // sort picked edges according to a distance of their first node along the polyline - std::sort(_boundary_elements.begin(), _boundary_elements.end(), - [&](MeshLib::Element*e1, MeshLib::Element*e2) - { - std::size_t dist1 = std::distance(node_ids_on_poly.begin(), - std::find(node_ids_on_poly.begin(), node_ids_on_poly.end(), e1->getNodeIndex(0))); - std::size_t dist2 = std::distance(node_ids_on_poly.begin(), - std::find(node_ids_on_poly.begin(), node_ids_on_poly.end(), e2->getNodeIndex(0))); - return (dist1 < dist2); - }); + // sort picked edges according to a distance of their first node along the polyline + std::sort(_boundary_elements.begin(), _boundary_elements.end(), + [&](MeshLib::Element*e1, MeshLib::Element*e2) + { + std::size_t dist1 = std::distance(node_ids_on_poly.begin(), + std::find(node_ids_on_poly.begin(), node_ids_on_poly.end(), e1->getNodeIndex(0))); + std::size_t dist2 = std::distance(node_ids_on_poly.begin(), + std::find(node_ids_on_poly.begin(), node_ids_on_poly.end(), e2->getNodeIndex(0))); + return (dist1 < dist2); + }); } BoundaryElementsAlongPolyline::~BoundaryElementsAlongPolyline() { - for (auto p : _boundary_elements) - delete p; + for (auto p : _boundary_elements) + delete p; } bool BoundaryElementsAlongPolyline::includesAllEdgeNodeIDs(const std::vector<std::size_t> &vec_node_ids, const MeshLib::Element &edge, std::vector<std::size_t> &edge_node_distances) const { - unsigned j=0; - for (; j<edge.getNBaseNodes(); j++) { - auto itr = std::find(vec_node_ids.begin(), vec_node_ids.end(), edge.getNodeIndex(j)); - if (itr != vec_node_ids.end()) - edge_node_distances.push_back(std::distance(vec_node_ids.begin(), itr)); - else - break; - } - return (j==edge.getNBaseNodes()); + unsigned j=0; + for (; j<edge.getNBaseNodes(); j++) { + auto itr = std::find(vec_node_ids.begin(), vec_node_ids.end(), edge.getNodeIndex(j)); + if (itr != vec_node_ids.end()) + edge_node_distances.push_back(std::distance(vec_node_ids.begin(), itr)); + else + break; + } + return (j==edge.getNBaseNodes()); } MeshLib::Element* BoundaryElementsAlongPolyline::modifyEdgeNodeOrdering(const MeshLib::Element &edge, const GeoLib::Polyline &ply, const std::vector<std::size_t> &edge_node_distances_along_ply, const std::vector<std::size_t> &node_ids_on_poly) const { - MeshLib::Element* new_edge = const_cast<MeshLib::Element*>(&edge); - // The first node of the edge should be always closer to the beginning of the polyline than other nodes. - // Otherwise, create a new element with reversed local node index - if (edge_node_distances_along_ply.front() > edge_node_distances_along_ply.back() - || (ply.isClosed() && edge_node_distances_along_ply.back() == node_ids_on_poly.size()-1)) { - MeshLib::Node** new_nodes = new MeshLib::Node*[edge.getNBaseNodes()]; - std::reverse_copy(edge.getNodes(), edge.getNodes()+edge.getNBaseNodes(), new_nodes); - new_edge = new MeshLib::Line(new_nodes); - } - return new_edge; + MeshLib::Element* new_edge = const_cast<MeshLib::Element*>(&edge); + // The first node of the edge should be always closer to the beginning of the polyline than other nodes. + // Otherwise, create a new element with reversed local node index + if (edge_node_distances_along_ply.front() > edge_node_distances_along_ply.back() + || (ply.isClosed() && edge_node_distances_along_ply.back() == node_ids_on_poly.size()-1)) { + MeshLib::Node** new_nodes = new MeshLib::Node*[edge.getNBaseNodes()]; + std::reverse_copy(edge.getNodes(), edge.getNodes()+edge.getNBaseNodes(), new_nodes); + new_edge = new MeshLib::Line(new_nodes); + } + return new_edge; } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h index 713139ba9fdee714853cf1c13af695bb916347cd..574b6ee9923f01701293fc4150a4bfd349fcd9d5 100644 --- a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h +++ b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h @@ -32,56 +32,56 @@ class MeshNodeSearcher; class BoundaryElementsAlongPolyline { public: - /** - * Constructor - * @param mesh a mesh object - * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes - * @param ply a polyline object where edges are searched - */ - BoundaryElementsAlongPolyline(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Polyline const& ply); + /** + * Constructor + * @param mesh a mesh object + * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes + * @param ply a polyline object where edges are searched + */ + BoundaryElementsAlongPolyline(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Polyline const& ply); - /// destructor - virtual ~BoundaryElementsAlongPolyline(); + /// destructor + virtual ~BoundaryElementsAlongPolyline(); - /// return the mesh object - MeshLib::Mesh const& getMesh() const {return _mesh;} + /// return the mesh object + MeshLib::Mesh const& getMesh() const {return _mesh;} - /** - * Deploying this method the user can get access to the underlying - * GeoLib::Polyline. - * @return the underlying GeoLib::Polyline - */ - GeoLib::Polyline const& getPolyline () const {return _ply;} + /** + * Deploying this method the user can get access to the underlying + * GeoLib::Polyline. + * @return the underlying GeoLib::Polyline + */ + GeoLib::Polyline const& getPolyline () const {return _ply;} - /** - * Return the vector of boundary elements (i.e. edges). The elements are sorted - * according to their distance to the starting point of the given polyline. - */ - std::vector<MeshLib::Element*> const& getBoundaryElements() const {return _boundary_elements; } + /** + * Return the vector of boundary elements (i.e. edges). The elements are sorted + * according to their distance to the starting point of the given polyline. + */ + std::vector<MeshLib::Element*> const& getBoundaryElements() const {return _boundary_elements; } private: - /** - * Check if a vector of node IDs includes all nodes of a given element - * @param vec_node_ids a vector of Node IDs - * @param edge Edge object whose node IDs are checked - * @param edge_node_distances a vector of distances of the edge nodes from the beginning of the given node ID vector - * @return true if all element nodes are included in the vector - */ - bool includesAllEdgeNodeIDs(const std::vector<std::size_t> &vec_node_ids, const MeshLib::Element &edge, std::vector<std::size_t> &edge_node_distances) const; + /** + * Check if a vector of node IDs includes all nodes of a given element + * @param vec_node_ids a vector of Node IDs + * @param edge Edge object whose node IDs are checked + * @param edge_node_distances a vector of distances of the edge nodes from the beginning of the given node ID vector + * @return true if all element nodes are included in the vector + */ + bool includesAllEdgeNodeIDs(const std::vector<std::size_t> &vec_node_ids, const MeshLib::Element &edge, std::vector<std::size_t> &edge_node_distances) const; - /** - * Modify node ordering of an edge so that its first node is closer to the beginning of a polyline than others - * @param edge Element object - * @param ply Polyline object - * @param edge_node_distances_along_ply A vector of current edge node distances along poly - * @param node_ids_on_poly A vector of node IDs along the polyine - * @return A pointer to the new modified edge object. A pointer to the original edge is returned if the modification is unnecessary. - */ - MeshLib::Element* modifyEdgeNodeOrdering(const MeshLib::Element &edge, const GeoLib::Polyline &ply, const std::vector<std::size_t> &edge_node_distances_along_ply, const std::vector<std::size_t> &node_ids_on_poly) const; + /** + * Modify node ordering of an edge so that its first node is closer to the beginning of a polyline than others + * @param edge Element object + * @param ply Polyline object + * @param edge_node_distances_along_ply A vector of current edge node distances along poly + * @param node_ids_on_poly A vector of node IDs along the polyine + * @return A pointer to the new modified edge object. A pointer to the original edge is returned if the modification is unnecessary. + */ + MeshLib::Element* modifyEdgeNodeOrdering(const MeshLib::Element &edge, const GeoLib::Polyline &ply, const std::vector<std::size_t> &edge_node_distances_along_ply, const std::vector<std::size_t> &node_ids_on_poly) const; - MeshLib::Mesh const& _mesh; - GeoLib::Polyline const& _ply; - std::vector<MeshLib::Element*> _boundary_elements; + MeshLib::Mesh const& _mesh; + GeoLib::Polyline const& _ply; + std::vector<MeshLib::Element*> _boundary_elements; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp index 6430b5fa817a3f2bbbbc19d09b89fa1334cb2dc9..a875960f356ccc42567796c0e93328aba91b6955 100644 --- a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp +++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp @@ -24,18 +24,18 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint( GeoLib::Point const &point) : _mesh(mesh), _point(point) { - auto const node_ids = mshNodeSearcher.getMeshNodeIDs(_point); - assert(node_ids.size() == 1); - std::array<MeshLib::Node*, 1> const nodes = {{ - const_cast<MeshLib::Node*>(_mesh.getNode(node_ids[0]))}}; + auto const node_ids = mshNodeSearcher.getMeshNodeIDs(_point); + assert(node_ids.size() == 1); + std::array<MeshLib::Node*, 1> const nodes = {{ + const_cast<MeshLib::Node*>(_mesh.getNode(node_ids[0]))}}; - _boundary_elements.push_back(new MeshLib::Point{nodes, node_ids[0]}); + _boundary_elements.push_back(new MeshLib::Point{nodes, node_ids[0]}); } BoundaryElementsAtPoint::~BoundaryElementsAtPoint() { - for (auto p : _boundary_elements) - delete p; + for (auto p : _boundary_elements) + delete p; } } // MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.h b/MeshGeoToolsLib/BoundaryElementsAtPoint.h index fd1a6d2eeb293e91fc42e314b073eba3e5d8f265..7f7f88de22e1891434d3b92f8dbbc3b9dde11391 100644 --- a/MeshGeoToolsLib/BoundaryElementsAtPoint.h +++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.h @@ -30,37 +30,37 @@ class MeshNodeSearcher; class BoundaryElementsAtPoint final { public: - /// Constructor - /// \param mesh a mesh object - /// \param mshNodeSearcher a MeshNodeSearcher object which is internally - /// used to search mesh nodes - /// \param point a point object where edges are searched - BoundaryElementsAtPoint(MeshLib::Mesh const& mesh, - MeshNodeSearcher& mshNodeSearcher, - GeoLib::Point const& point); + /// Constructor + /// \param mesh a mesh object + /// \param mshNodeSearcher a MeshNodeSearcher object which is internally + /// used to search mesh nodes + /// \param point a point object where edges are searched + BoundaryElementsAtPoint(MeshLib::Mesh const& mesh, + MeshNodeSearcher& mshNodeSearcher, + GeoLib::Point const& point); - ~BoundaryElementsAtPoint(); + ~BoundaryElementsAtPoint(); - MeshLib::Mesh const& getMesh() const - { - return _mesh; - } + MeshLib::Mesh const& getMesh() const + { + return _mesh; + } - GeoLib::Point const& getPoint() const - { - return _point; - } + GeoLib::Point const& getPoint() const + { + return _point; + } - /// Return the vector of boundary elements (i.e. points). - std::vector<MeshLib::Element*> const& getBoundaryElements() const - { - return _boundary_elements; - } + /// Return the vector of boundary elements (i.e. points). + std::vector<MeshLib::Element*> const& getBoundaryElements() const + { + return _boundary_elements; + } private: - MeshLib::Mesh const& _mesh; - GeoLib::Point const& _point; - std::vector<MeshLib::Element*> _boundary_elements; + MeshLib::Mesh const& _mesh; + GeoLib::Point const& _point; + std::vector<MeshLib::Element*> _boundary_elements; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp b/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp index 92a9a1f1263d99bf4147b4c7d63de0bcc0505a91..5ef22526e4e9e03d04e0210789ef5a92c93d5e09 100644 --- a/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp +++ b/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp @@ -22,42 +22,42 @@ namespace MeshGeoToolsLib BoundaryElementsOnSurface::BoundaryElementsOnSurface(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Surface const& sfc) : _mesh(mesh), _sfc(sfc) { - // search elements near the polyline - auto node_ids_on_sfc = mshNodeSearcher.getMeshNodeIDsAlongSurface(sfc); - MeshLib::ElementSearch es(_mesh); - es.searchByNodeIDs(node_ids_on_sfc); - auto &ele_ids_near_sfc = es.getSearchedElementIDs(); + // search elements near the polyline + auto node_ids_on_sfc = mshNodeSearcher.getMeshNodeIDsAlongSurface(sfc); + MeshLib::ElementSearch es(_mesh); + es.searchByNodeIDs(node_ids_on_sfc); + auto &ele_ids_near_sfc = es.getSearchedElementIDs(); - // get a list of faces made of the nodes - for (auto ele_id : ele_ids_near_sfc) { - auto* e = _mesh.getElement(ele_id); - // skip internal elements - if (!e->isBoundaryElement()) - continue; - // find faces on surface - for (unsigned i=0; i<e->getNFaces(); i++) { - auto* face = e->getFace(i); - // check - std::size_t cnt_match = 0; - for (std::size_t j=0; j<face->getNBaseNodes(); j++) { - if (std::find(node_ids_on_sfc.begin(), node_ids_on_sfc.end(), face->getNodeIndex(j)) != node_ids_on_sfc.end()) - cnt_match++; - else - break; - } - // update the list - if (cnt_match==face->getNBaseNodes()) - _boundary_elements.push_back(const_cast<MeshLib::Element*>(face)); - else - delete face; - } - } + // get a list of faces made of the nodes + for (auto ele_id : ele_ids_near_sfc) { + auto* e = _mesh.getElement(ele_id); + // skip internal elements + if (!e->isBoundaryElement()) + continue; + // find faces on surface + for (unsigned i=0; i<e->getNFaces(); i++) { + auto* face = e->getFace(i); + // check + std::size_t cnt_match = 0; + for (std::size_t j=0; j<face->getNBaseNodes(); j++) { + if (std::find(node_ids_on_sfc.begin(), node_ids_on_sfc.end(), face->getNodeIndex(j)) != node_ids_on_sfc.end()) + cnt_match++; + else + break; + } + // update the list + if (cnt_match==face->getNBaseNodes()) + _boundary_elements.push_back(const_cast<MeshLib::Element*>(face)); + else + delete face; + } + } } BoundaryElementsOnSurface::~BoundaryElementsOnSurface() { - for (auto p : _boundary_elements) - delete p; + for (auto p : _boundary_elements) + delete p; } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsOnSurface.h b/MeshGeoToolsLib/BoundaryElementsOnSurface.h index 64f07994dcb46949a3112b96533cc296e017ddf9..6dde46912661b612bd5984d7af9520b0ce6226bf 100644 --- a/MeshGeoToolsLib/BoundaryElementsOnSurface.h +++ b/MeshGeoToolsLib/BoundaryElementsOnSurface.h @@ -32,36 +32,36 @@ class MeshNodeSearcher; class BoundaryElementsOnSurface { public: - /** - * Constructor - * @param mesh a mesh object - * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes - * @param sfc a surface object where face elements are searched for - */ - BoundaryElementsOnSurface(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Surface const& sfc); + /** + * Constructor + * @param mesh a mesh object + * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes + * @param sfc a surface object where face elements are searched for + */ + BoundaryElementsOnSurface(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher, GeoLib::Surface const& sfc); - /// destructor - virtual ~BoundaryElementsOnSurface(); + /// destructor + virtual ~BoundaryElementsOnSurface(); - /// return the mesh object - MeshLib::Mesh const& getMesh() const {return _mesh;} + /// return the mesh object + MeshLib::Mesh const& getMesh() const {return _mesh;} - /** - * Deploying this method the user can get access to the underlying - * GeoLib::Surface. - * @return the underlying GeoLib::Surface - */ - GeoLib::Surface const& getSurface() const {return _sfc;} + /** + * Deploying this method the user can get access to the underlying + * GeoLib::Surface. + * @return the underlying GeoLib::Surface + */ + GeoLib::Surface const& getSurface() const {return _sfc;} - /** - * Return the vector of boundary elements (i.e. faces). The elements are unsorted. - */ - std::vector<MeshLib::Element*> const& getBoundaryElements() const {return _boundary_elements;} + /** + * Return the vector of boundary elements (i.e. faces). The elements are unsorted. + */ + std::vector<MeshLib::Element*> const& getBoundaryElements() const {return _boundary_elements;} private: - MeshLib::Mesh const& _mesh; - GeoLib::Surface const& _sfc; - std::vector<MeshLib::Element*> _boundary_elements; + MeshLib::Mesh const& _mesh; + GeoLib::Surface const& _sfc; + std::vector<MeshLib::Element*> _boundary_elements; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp index 1062d06be49f8c2404cbf6589d28d95ae4fa9217..5685729ba02c7be357f6465924e8745539f1a232 100644 --- a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp +++ b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp @@ -31,76 +31,76 @@ BoundaryElementsSearcher::BoundaryElementsSearcher(MeshLib::Mesh const& mesh, Me BoundaryElementsSearcher::~BoundaryElementsSearcher() { - for (auto p : _boundary_elements_at_point) - delete p; - for (auto p : _boundary_elements_along_polylines) - delete p; - for (auto p : _boundary_elements_along_surfaces) - delete p; + for (auto p : _boundary_elements_at_point) + delete p; + for (auto p : _boundary_elements_along_polylines) + delete p; + for (auto p : _boundary_elements_along_surfaces) + delete p; } std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElements(GeoLib::GeoObject const& geoObj) { - switch (geoObj.getGeoType()) { - case GeoLib::GEOTYPE::POINT: - return this->getBoundaryElementsAtPoint(*dynamic_cast<const GeoLib::Point*>(&geoObj)); - break; - case GeoLib::GEOTYPE::POLYLINE: - return this->getBoundaryElementsAlongPolyline(*dynamic_cast<const GeoLib::Polyline*>(&geoObj)); - break; - case GeoLib::GEOTYPE::SURFACE: - return this->getBoundaryElementsOnSurface(*dynamic_cast<const GeoLib::Surface*>(&geoObj)); - break; - default: - const static std::vector<MeshLib::Element*> dummy(0); - return dummy; - } + switch (geoObj.getGeoType()) { + case GeoLib::GEOTYPE::POINT: + return this->getBoundaryElementsAtPoint(*dynamic_cast<const GeoLib::Point*>(&geoObj)); + break; + case GeoLib::GEOTYPE::POLYLINE: + return this->getBoundaryElementsAlongPolyline(*dynamic_cast<const GeoLib::Polyline*>(&geoObj)); + break; + case GeoLib::GEOTYPE::SURFACE: + return this->getBoundaryElementsOnSurface(*dynamic_cast<const GeoLib::Surface*>(&geoObj)); + break; + default: + const static std::vector<MeshLib::Element*> dummy(0); + return dummy; + } } std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElementsAtPoint(GeoLib::Point const& point) { - // look for already saved points and return if found. - for (auto const& boundaryElements : _boundary_elements_at_point) - { - if (boundaryElements->getPoint() == point) - return boundaryElements->getBoundaryElements(); - } - - // create new boundary elements at points. - _boundary_elements_at_point.push_back( - new BoundaryElementsAtPoint(_mesh, _mshNodeSearcher, point)); - return _boundary_elements_at_point.back()->getBoundaryElements(); + // look for already saved points and return if found. + for (auto const& boundaryElements : _boundary_elements_at_point) + { + if (boundaryElements->getPoint() == point) + return boundaryElements->getBoundaryElements(); + } + + // create new boundary elements at points. + _boundary_elements_at_point.push_back( + new BoundaryElementsAtPoint(_mesh, _mshNodeSearcher, point)); + return _boundary_elements_at_point.back()->getBoundaryElements(); } std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElementsAlongPolyline(GeoLib::Polyline const& ply) { - std::vector<BoundaryElementsAlongPolyline*>::const_iterator it(_boundary_elements_along_polylines.begin()); - for (; it != _boundary_elements_along_polylines.end(); ++it) { - if (&(*it)->getPolyline() == &ply) { - // we calculated mesh nodes for this polyline already - return (*it)->getBoundaryElements(); - } - } - - _boundary_elements_along_polylines.push_back( - new BoundaryElementsAlongPolyline(_mesh, _mshNodeSearcher, ply)); - return _boundary_elements_along_polylines.back()->getBoundaryElements(); + std::vector<BoundaryElementsAlongPolyline*>::const_iterator it(_boundary_elements_along_polylines.begin()); + for (; it != _boundary_elements_along_polylines.end(); ++it) { + if (&(*it)->getPolyline() == &ply) { + // we calculated mesh nodes for this polyline already + return (*it)->getBoundaryElements(); + } + } + + _boundary_elements_along_polylines.push_back( + new BoundaryElementsAlongPolyline(_mesh, _mshNodeSearcher, ply)); + return _boundary_elements_along_polylines.back()->getBoundaryElements(); } std::vector<MeshLib::Element*> const& BoundaryElementsSearcher::getBoundaryElementsOnSurface(GeoLib::Surface const& sfc) { - std::vector<BoundaryElementsOnSurface*>::const_iterator it(_boundary_elements_along_surfaces.begin()); - for (; it != _boundary_elements_along_surfaces.end(); ++it) { - if (&(*it)->getSurface() == &sfc) { - // we calculated mesh nodes for this surface already - return (*it)->getBoundaryElements(); - } - } - - _boundary_elements_along_surfaces.push_back( - new BoundaryElementsOnSurface(_mesh, _mshNodeSearcher, sfc)); - return _boundary_elements_along_surfaces.back()->getBoundaryElements(); + std::vector<BoundaryElementsOnSurface*>::const_iterator it(_boundary_elements_along_surfaces.begin()); + for (; it != _boundary_elements_along_surfaces.end(); ++it) { + if (&(*it)->getSurface() == &sfc) { + // we calculated mesh nodes for this surface already + return (*it)->getBoundaryElements(); + } + } + + _boundary_elements_along_surfaces.push_back( + new BoundaryElementsOnSurface(_mesh, _mshNodeSearcher, sfc)); + return _boundary_elements_along_surfaces.back()->getBoundaryElements(); } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsSearcher.h b/MeshGeoToolsLib/BoundaryElementsSearcher.h index 5c2127d29583396741630ca916c95672a8b94dfb..34d86c45721913543abbcee707af8d6ba1f09241 100644 --- a/MeshGeoToolsLib/BoundaryElementsSearcher.h +++ b/MeshGeoToolsLib/BoundaryElementsSearcher.h @@ -38,52 +38,52 @@ class BoundaryElementsOnSurface; class BoundaryElementsSearcher { public: - /** - * Constructor - * @param mesh a mesh object - * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes - */ - BoundaryElementsSearcher(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher); + /** + * Constructor + * @param mesh a mesh object + * @param mshNodeSearcher a MeshNodeSearcher object which is internally used to search mesh nodes + */ + BoundaryElementsSearcher(MeshLib::Mesh const& mesh, MeshNodeSearcher &mshNodeSearcher); - /// destructor - virtual ~BoundaryElementsSearcher(); + /// destructor + virtual ~BoundaryElementsSearcher(); - /** - * generate boundary elements on the given geometric object (point, polyline, surface). - * - * @param geoObj a GeoLib::GeoObject where the nearest mesh node is searched for - * @return a vector of boundary element objects - */ - std::vector<MeshLib::Element*> const& getBoundaryElements(GeoLib::GeoObject const& geoObj); + /** + * generate boundary elements on the given geometric object (point, polyline, surface). + * + * @param geoObj a GeoLib::GeoObject where the nearest mesh node is searched for + * @return a vector of boundary element objects + */ + std::vector<MeshLib::Element*> const& getBoundaryElements(GeoLib::GeoObject const& geoObj); - /** - * generate boundary elements at the given point. - * @param point Search the mesh for given point - * @return a vector of boundary elements - */ - std::vector<MeshLib::Element*> const& getBoundaryElementsAtPoint( - GeoLib::Point const& point); - /** - * generate boundary elements on the given polyline. - * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for - * @return a vector of boundary element objects - */ - std::vector<MeshLib::Element*> const& getBoundaryElementsAlongPolyline(GeoLib::Polyline const& ply); + /** + * generate boundary elements at the given point. + * @param point Search the mesh for given point + * @return a vector of boundary elements + */ + std::vector<MeshLib::Element*> const& getBoundaryElementsAtPoint( + GeoLib::Point const& point); + /** + * generate boundary elements on the given polyline. + * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for + * @return a vector of boundary element objects + */ + std::vector<MeshLib::Element*> const& getBoundaryElementsAlongPolyline(GeoLib::Polyline const& ply); - /** - * generate boundary elements on the given surface. - * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for - * @return a vector of boundary element objects - */ - std::vector<MeshLib::Element*> const& getBoundaryElementsOnSurface(GeoLib::Surface const& sfc); + /** + * generate boundary elements on the given surface. + * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for + * @return a vector of boundary element objects + */ + std::vector<MeshLib::Element*> const& getBoundaryElementsOnSurface(GeoLib::Surface const& sfc); private: - MeshLib::Mesh const& _mesh; - MeshNodeSearcher &_mshNodeSearcher; - std::vector<BoundaryElementsAtPoint*> _boundary_elements_at_point; - std::vector<BoundaryElementsAlongPolyline*> _boundary_elements_along_polylines; - std::vector<BoundaryElementsOnSurface*> _boundary_elements_along_surfaces; + MeshLib::Mesh const& _mesh; + MeshNodeSearcher &_mshNodeSearcher; + std::vector<BoundaryElementsAtPoint*> _boundary_elements_at_point; + std::vector<BoundaryElementsAlongPolyline*> _boundary_elements_along_polylines; + std::vector<BoundaryElementsOnSurface*> _boundary_elements_along_surfaces; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/CMakeLists.txt b/MeshGeoToolsLib/CMakeLists.txt index 728c484ac92aae3ad9bb17cc1554e493907abb74..0a91b15d96425a185401907d3379964e41684c7c 100644 --- a/MeshGeoToolsLib/CMakeLists.txt +++ b/MeshGeoToolsLib/CMakeLists.txt @@ -5,13 +5,13 @@ GET_SOURCE_FILES(SOURCES_MeshGeoToolsLib) add_library(MeshGeoToolsLib ${SOURCES_MeshGeoToolsLib}) target_link_libraries(MeshGeoToolsLib - BaseLib - MathLib - MeshLib - GeoLib - FileIO + BaseLib + MathLib + MeshLib + GeoLib + FileIO ) if(TARGET Boost) - add_dependencies(MeshGeoToolsLib Boost) + add_dependencies(MeshGeoToolsLib Boost) endif() diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp index e74f7fba4d1176b60349cb60bbe02713162cb16c..8401b02499315c4d52a9a608a434319bae7c31c2 100644 --- a/MeshGeoToolsLib/GeoMapper.cpp +++ b/MeshGeoToolsLib/GeoMapper.cpp @@ -37,395 +37,395 @@ namespace MeshGeoToolsLib { GeoMapper::GeoMapper(GeoLib::GEOObjects &geo_objects, const std::string &geo_name) - : _geo_objects(geo_objects), _geo_name(const_cast<std::string&>(geo_name)), - _surface_mesh(nullptr), _grid(nullptr), _raster(nullptr) + : _geo_objects(geo_objects), _geo_name(const_cast<std::string&>(geo_name)), + _surface_mesh(nullptr), _grid(nullptr), _raster(nullptr) { } GeoMapper::~GeoMapper() { - delete _surface_mesh; - delete _raster; + delete _surface_mesh; + delete _raster; } void GeoMapper::mapOnDEM(const std::string &file_name) { - _raster = GeoLib::IO::AsciiRasterInterface::getRasterFromASCFile(file_name); - if (! _raster) { - ERR("GeoMapper::mapOnDEM(): failed to load %s", file_name.c_str()); - return; - } - - std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); - if (! pnts) { - ERR("Geometry \"%s\" does not exist.", _geo_name.c_str()); - return; - } - if (GeoLib::isStation((*pnts)[0])) { - mapStationData(*pnts); - } else { - mapPointDataToDEM(*pnts); - } + _raster = GeoLib::IO::AsciiRasterInterface::getRasterFromASCFile(file_name); + if (! _raster) { + ERR("GeoMapper::mapOnDEM(): failed to load %s", file_name.c_str()); + return; + } + + std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); + if (! pnts) { + ERR("Geometry \"%s\" does not exist.", _geo_name.c_str()); + return; + } + if (GeoLib::isStation((*pnts)[0])) { + mapStationData(*pnts); + } else { + mapPointDataToDEM(*pnts); + } } void GeoMapper::mapOnMesh(const std::string &file_name) { - MeshLib::Mesh *mesh (MeshLib::IO::readMeshFromFile(file_name)); - mapOnMesh(mesh); - delete mesh; + MeshLib::Mesh *mesh (MeshLib::IO::readMeshFromFile(file_name)); + mapOnMesh(mesh); + delete mesh; } void GeoMapper::mapOnMesh(const MeshLib::Mesh* mesh) { - std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); - if (! pnts) { - ERR("Geometry \"%s\" does not exist.", _geo_name.c_str()); - return; - } - - // the variable _surface_mesh is reused below, so first the existing - // _surface_mesh has to be cleaned up - if (_surface_mesh) - delete _surface_mesh; - - if (mesh->getDimension()<3) - _surface_mesh = new MeshLib::Mesh(*mesh); - else - { - const MathLib::Vector3 dir(0,0,-1); - _surface_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90); - } - - // init grid - MathLib::Point3d origin(std::array<double,3>{{0,0,0}}); - MathLib::Vector3 normal(0,0,-1); - std::vector<MeshLib::Node> flat_nodes; - flat_nodes.reserve(_surface_mesh->getNNodes()); - // copy nodes and project the copied nodes to the x-y-plane, i.e. set - // z-coordinate to zero - for (auto n_ptr : _surface_mesh->getNodes()) { - flat_nodes.emplace_back(*n_ptr); - flat_nodes.back()[2] = 0.0; - } - _grid = new GeoLib::Grid<MeshLib::Node>(flat_nodes.cbegin(), flat_nodes.cend()); - - if (GeoLib::isStation((*pnts)[0])) { - mapStationData(*pnts); - } else { - mapPointDataToMeshSurface(*pnts); - } - - delete _grid; + std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); + if (! pnts) { + ERR("Geometry \"%s\" does not exist.", _geo_name.c_str()); + return; + } + + // the variable _surface_mesh is reused below, so first the existing + // _surface_mesh has to be cleaned up + if (_surface_mesh) + delete _surface_mesh; + + if (mesh->getDimension()<3) + _surface_mesh = new MeshLib::Mesh(*mesh); + else + { + const MathLib::Vector3 dir(0,0,-1); + _surface_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface(*mesh, dir, 90); + } + + // init grid + MathLib::Point3d origin(std::array<double,3>{{0,0,0}}); + MathLib::Vector3 normal(0,0,-1); + std::vector<MeshLib::Node> flat_nodes; + flat_nodes.reserve(_surface_mesh->getNNodes()); + // copy nodes and project the copied nodes to the x-y-plane, i.e. set + // z-coordinate to zero + for (auto n_ptr : _surface_mesh->getNodes()) { + flat_nodes.emplace_back(*n_ptr); + flat_nodes.back()[2] = 0.0; + } + _grid = new GeoLib::Grid<MeshLib::Node>(flat_nodes.cbegin(), flat_nodes.cend()); + + if (GeoLib::isStation((*pnts)[0])) { + mapStationData(*pnts); + } else { + mapPointDataToMeshSurface(*pnts); + } + + delete _grid; } void GeoMapper::mapToConstantValue(double value) { - std::vector<GeoLib::Point*> const* points (this->_geo_objects.getPointVec(this->_geo_name)); - if (points == nullptr) - { - ERR ("Geometry \"%s\" not found.", this->_geo_name.c_str()); - return; - } - std::for_each(points->begin(), points->end(), [value](GeoLib::Point* pnt){ (*pnt)[2] = value; }); + std::vector<GeoLib::Point*> const* points (this->_geo_objects.getPointVec(this->_geo_name)); + if (points == nullptr) + { + ERR ("Geometry \"%s\" not found.", this->_geo_name.c_str()); + return; + } + std::for_each(points->begin(), points->end(), [value](GeoLib::Point* pnt){ (*pnt)[2] = value; }); } void GeoMapper::mapStationData(std::vector<GeoLib::Point*> const& points) { - double min_val(0), max_val(0); - if (_surface_mesh) - { - GeoLib::AABB bounding_box( - _surface_mesh->getNodes().begin(), _surface_mesh->getNodes().end()); - min_val = bounding_box.getMinPoint()[2]; - max_val = bounding_box.getMaxPoint()[2]; - } - - for (auto * pnt : points) - { - double offset = - (_grid) - ? (getMeshElevation((*pnt)[0], (*pnt)[1], min_val, max_val) - - (*pnt)[2]) - : getDemElevation(*pnt); - - if (!GeoLib::isBorehole(pnt)) - continue; - auto const& layers = static_cast<GeoLib::StationBorehole*>(pnt)->getProfile(); - for (auto * layer_pnt : layers) { - (*layer_pnt)[2] = (*layer_pnt)[2] + offset; - } - } + double min_val(0), max_val(0); + if (_surface_mesh) + { + GeoLib::AABB bounding_box( + _surface_mesh->getNodes().begin(), _surface_mesh->getNodes().end()); + min_val = bounding_box.getMinPoint()[2]; + max_val = bounding_box.getMaxPoint()[2]; + } + + for (auto * pnt : points) + { + double offset = + (_grid) + ? (getMeshElevation((*pnt)[0], (*pnt)[1], min_val, max_val) - + (*pnt)[2]) + : getDemElevation(*pnt); + + if (!GeoLib::isBorehole(pnt)) + continue; + auto const& layers = static_cast<GeoLib::StationBorehole*>(pnt)->getProfile(); + for (auto * layer_pnt : layers) { + (*layer_pnt)[2] = (*layer_pnt)[2] + offset; + } + } } void GeoMapper::mapPointDataToDEM(std::vector<GeoLib::Point*> const& points) { - for (auto * pnt : points) - { - GeoLib::Point &p(*pnt); - p[2] = getDemElevation(p); - } + for (auto * pnt : points) + { + GeoLib::Point &p(*pnt); + p[2] = getDemElevation(p); + } } void GeoMapper::mapPointDataToMeshSurface(std::vector<GeoLib::Point*> const& pnts) { - GeoLib::AABB const aabb( - _surface_mesh->getNodes().cbegin(), _surface_mesh->getNodes().cend()); - double const min_val(aabb.getMinPoint()[2]); - double const max_val(aabb.getMaxPoint()[2]); - - for (auto * pnt : pnts) { - // check if pnt is inside of the bounding box of the _surface_mesh - // projected onto the y-x plane - GeoLib::Point &p(*pnt); - if (p[0] < aabb.getMinPoint()[0] || aabb.getMaxPoint()[0] < p[0]) - continue; - if (p[1] < aabb.getMinPoint()[1] || aabb.getMaxPoint()[1] < p[1]) - continue; - - p[2] = getMeshElevation(p[0], p[1], min_val, max_val); - } + GeoLib::AABB const aabb( + _surface_mesh->getNodes().cbegin(), _surface_mesh->getNodes().cend()); + double const min_val(aabb.getMinPoint()[2]); + double const max_val(aabb.getMaxPoint()[2]); + + for (auto * pnt : pnts) { + // check if pnt is inside of the bounding box of the _surface_mesh + // projected onto the y-x plane + GeoLib::Point &p(*pnt); + if (p[0] < aabb.getMinPoint()[0] || aabb.getMaxPoint()[0] < p[0]) + continue; + if (p[1] < aabb.getMinPoint()[1] || aabb.getMaxPoint()[1] < p[1]) + continue; + + p[2] = getMeshElevation(p[0], p[1], min_val, max_val); + } } float GeoMapper::getDemElevation(GeoLib::Point const& pnt) const { - double const elevation (_raster->getValueAtPoint(pnt)); - if (std::abs(elevation-_raster->getHeader().no_data) < std::numeric_limits<double>::epsilon()) - return 0.0; - return static_cast<float>(elevation); + double const elevation (_raster->getValueAtPoint(pnt)); + if (std::abs(elevation-_raster->getHeader().no_data) < std::numeric_limits<double>::epsilon()) + return 0.0; + return static_cast<float>(elevation); } double GeoMapper::getMeshElevation( - double x, double y, double min_val, double max_val) const + double x, double y, double min_val, double max_val) const { - const MeshLib::Node* pnt = - _grid->getNearestPoint(MathLib::Point3d{{{x, y, 0}}}); - const std::vector<MeshLib::Element*> elements( - _surface_mesh->getNode(pnt->getID())->getElements()); - std::unique_ptr<GeoLib::Point> intersection; - - for (auto const & element : elements) - { - if (intersection == nullptr && - element->getGeomType() != MeshLib::MeshElemType::LINE) - intersection = GeoLib::triangleLineIntersection( - *element->getNode(0), *element->getNode(1), - *element->getNode(2), GeoLib::Point(x, y, max_val), - GeoLib::Point(x, y, min_val)); - - if (intersection == nullptr && - element->getGeomType() == MeshLib::MeshElemType::QUAD) - intersection = GeoLib::triangleLineIntersection( - *element->getNode(0), *element->getNode(2), - *element->getNode(3), GeoLib::Point(x, y, max_val), - GeoLib::Point(x, y, min_val)); - } - if (intersection) - return (*intersection)[2]; - // if something goes wrong, simply take the elevation of the nearest mesh node - return (*(_surface_mesh->getNode(pnt->getID())))[2]; + const MeshLib::Node* pnt = + _grid->getNearestPoint(MathLib::Point3d{{{x, y, 0}}}); + const std::vector<MeshLib::Element*> elements( + _surface_mesh->getNode(pnt->getID())->getElements()); + std::unique_ptr<GeoLib::Point> intersection; + + for (auto const & element : elements) + { + if (intersection == nullptr && + element->getGeomType() != MeshLib::MeshElemType::LINE) + intersection = GeoLib::triangleLineIntersection( + *element->getNode(0), *element->getNode(1), + *element->getNode(2), GeoLib::Point(x, y, max_val), + GeoLib::Point(x, y, min_val)); + + if (intersection == nullptr && + element->getGeomType() == MeshLib::MeshElemType::QUAD) + intersection = GeoLib::triangleLineIntersection( + *element->getNode(0), *element->getNode(2), + *element->getNode(3), GeoLib::Point(x, y, max_val), + GeoLib::Point(x, y, min_val)); + } + if (intersection) + return (*intersection)[2]; + // if something goes wrong, simply take the elevation of the nearest mesh node + return (*(_surface_mesh->getNode(pnt->getID())))[2]; } std::unique_ptr<std::vector<GeoLib::Polyline*>> copyPolylinesVector( - std::vector<GeoLib::Polyline*> const& polylines, - std::vector<GeoLib::Point*> const& points) + std::vector<GeoLib::Polyline*> const& polylines, + std::vector<GeoLib::Point*> const& points) { - std::size_t nLines = polylines.size(); - auto new_lines = std::unique_ptr<std::vector<GeoLib::Polyline*>>( - new std::vector<GeoLib::Polyline*>(nLines)); - - for (std::size_t i=0; i<nLines; ++i) - { - (*new_lines)[i] = new GeoLib::Polyline(points); - std::size_t nLinePnts (polylines[i]->getNumberOfPoints()); - for (std::size_t j=0; j<nLinePnts; ++j) - (*new_lines)[i]->addPoint(polylines[i]->getPointID(j)); - } - return new_lines; + std::size_t nLines = polylines.size(); + auto new_lines = std::unique_ptr<std::vector<GeoLib::Polyline*>>( + new std::vector<GeoLib::Polyline*>(nLines)); + + for (std::size_t i=0; i<nLines; ++i) + { + (*new_lines)[i] = new GeoLib::Polyline(points); + std::size_t nLinePnts (polylines[i]->getNumberOfPoints()); + for (std::size_t j=0; j<nLinePnts; ++j) + (*new_lines)[i]->addPoint(polylines[i]->getPointID(j)); + } + return new_lines; } void GeoMapper::advancedMapOnMesh( - MeshLib::Mesh const* mesh, std::string const& new_geo_name) + MeshLib::Mesh const* mesh, std::string const& new_geo_name) { - const std::vector<GeoLib::Point*> *points(_geo_objects.getPointVec(_geo_name)); - const std::vector<GeoLib::Polyline*> *org_lines(_geo_objects.getPolylineVec(_geo_name)); - - const GeoLib::AABB aabb(points->begin(), points->end()); - const double eps = sqrt(std::numeric_limits<float>::epsilon()) * - sqrt(MathLib::sqrDist(aabb.getMinPoint(),aabb.getMaxPoint())) ; - - // copy geometry (and set z=0 for all points) - auto new_points = std::unique_ptr<std::vector<GeoLib::Point*>>( - new std::vector<GeoLib::Point*>); - new_points->reserve(points->size()); - std::transform(points->cbegin(), points->cend(), std::back_inserter(*new_points), - [](GeoLib::Point* p) { return new GeoLib::Point((*p)[0],(*p)[1],0.0); }); - - auto new_lines = copyPolylinesVector(*_geo_objects.getPolylineVec(_geo_name), - *new_points); - - GeoLib::Grid<GeoLib::Point> grid(new_points->begin(), new_points->end()); - double max_segment_length(getMaxSegmentLength(*new_lines)); - // squared so it can be compared to the squared distances calculated later - max_segment_length *= max_segment_length; - - const unsigned nMeshNodes ( mesh->getNNodes() ); - // index of closest geo point for each mesh node in (x,y)-plane - std::vector<int> closest_geo_point(nMeshNodes, -1); - // distance between geo points and mesh nodes in (x,y)-plane - std::vector<double> dist(nMeshNodes); - auto zero_coords = GeoLib::Point{}; // All coordinates zero. - for (std::size_t i=0; i<nMeshNodes; ++i) { - zero_coords[0] = (*mesh->getNode(i))[0]; - zero_coords[1] = (*mesh->getNode(i))[1]; - GeoLib::Point* pnt = grid.getNearestPoint(zero_coords); - dist[i] = MathLib::sqrDist(*pnt, zero_coords); - closest_geo_point[i] = (dist[i]<=max_segment_length) ? pnt->getID() : -1; - } - - // store for each point the line segment to which it was added. - const std::size_t nLines(new_lines->size()); - std::vector< std::vector<unsigned> > line_segment_map(nLines); - for (std::size_t i=0; i<nLines; ++i) - { - line_segment_map[i] = std::vector<unsigned>((*new_lines)[i]->getNumberOfPoints(),0); - std::iota(line_segment_map[i].begin(), line_segment_map[i].end(), 0); - } - - for (std::size_t i=0; i<nMeshNodes; ++i) - { - // if mesh node too far away or exactly at point position - if (closest_geo_point[i] == -1 || dist[i] < eps) continue; - - const MeshLib::Node* node (mesh->getNode(i)); - for (std::size_t l=0; l<nLines; ++l) - { - // find relevant polylines - if (!(*org_lines)[l]->isPointIDInPolyline(closest_geo_point[i])) continue; - - // find point position of closest geo point in original polyline - GeoLib::Polyline* ply ((*org_lines)[l]); - std::size_t nLinePnts ( ply->getNumberOfPoints() ); - std::size_t node_index_in_ply (0); - for (node_index_in_ply=0; node_index_in_ply<nLinePnts; ++node_index_in_ply) - if (ply->getPoint(node_index_in_ply) == (*points)[closest_geo_point[i]]) - break; - const GeoLib::Point* geo_point (ply->getPoint(node_index_in_ply)); - - // check if line segments connected to closest geo point intersect connected elements of current node - const std::vector<MeshLib::Element*> elements (node->getElements()); - const std::size_t nElems = elements.size(); - for (std::size_t e=0; e<nElems; ++e) - { - const unsigned nEdges (elements[e]->getNEdges()); - unsigned intersection_count (0); - - for (unsigned n=0; n<nEdges; ++n) - { - if (intersection_count>1) break; //already two intersections - - const MeshLib::Element* line = elements[e]->getEdge(n); - unsigned index_offset(0); // default: add to first line segment - GeoLib::Point* intersection (nullptr); - if (node_index_in_ply>0) // test line segment before closest point - intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply-1)); - if (intersection == nullptr && node_index_in_ply<(nLinePnts-1)) // test line segment after closest point - { - intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply+1)); - index_offset = 1; // add to second segment - } - if (intersection) - { - intersection_count++; - unsigned start_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset-1;}))); - unsigned end_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset;}))); - std::size_t pos = getPointPosInLine((*new_lines)[l], start_point_idx, end_point_idx, intersection, eps); - - if (pos) - { - const std::size_t pnt_pos (new_points->size()); - new_points->push_back(intersection); - (*new_lines)[l]->insertPoint(pos, pnt_pos); - line_segment_map[l].insert(line_segment_map[l].begin()+pos, node_index_in_ply+index_offset-1); - } - } - } - } - } - } - - _geo_objects.addPointVec(std::move(new_points), const_cast<std::string&>(new_geo_name)); - std::vector<std::size_t> pnt_id_map = this->_geo_objects.getPointVecObj(new_geo_name)->getIDMap(); - for (auto & new_line : *new_lines) - new_line->updatePointIDs(pnt_id_map); - _geo_objects.addPolylineVec(std::move(new_lines), new_geo_name); - - // map new geometry incl. additional point using the normal mapping method - this->_geo_name = new_geo_name; - this->mapOnMesh(mesh); + const std::vector<GeoLib::Point*> *points(_geo_objects.getPointVec(_geo_name)); + const std::vector<GeoLib::Polyline*> *org_lines(_geo_objects.getPolylineVec(_geo_name)); + + const GeoLib::AABB aabb(points->begin(), points->end()); + const double eps = sqrt(std::numeric_limits<float>::epsilon()) * + sqrt(MathLib::sqrDist(aabb.getMinPoint(),aabb.getMaxPoint())) ; + + // copy geometry (and set z=0 for all points) + auto new_points = std::unique_ptr<std::vector<GeoLib::Point*>>( + new std::vector<GeoLib::Point*>); + new_points->reserve(points->size()); + std::transform(points->cbegin(), points->cend(), std::back_inserter(*new_points), + [](GeoLib::Point* p) { return new GeoLib::Point((*p)[0],(*p)[1],0.0); }); + + auto new_lines = copyPolylinesVector(*_geo_objects.getPolylineVec(_geo_name), + *new_points); + + GeoLib::Grid<GeoLib::Point> grid(new_points->begin(), new_points->end()); + double max_segment_length(getMaxSegmentLength(*new_lines)); + // squared so it can be compared to the squared distances calculated later + max_segment_length *= max_segment_length; + + const unsigned nMeshNodes ( mesh->getNNodes() ); + // index of closest geo point for each mesh node in (x,y)-plane + std::vector<int> closest_geo_point(nMeshNodes, -1); + // distance between geo points and mesh nodes in (x,y)-plane + std::vector<double> dist(nMeshNodes); + auto zero_coords = GeoLib::Point{}; // All coordinates zero. + for (std::size_t i=0; i<nMeshNodes; ++i) { + zero_coords[0] = (*mesh->getNode(i))[0]; + zero_coords[1] = (*mesh->getNode(i))[1]; + GeoLib::Point* pnt = grid.getNearestPoint(zero_coords); + dist[i] = MathLib::sqrDist(*pnt, zero_coords); + closest_geo_point[i] = (dist[i]<=max_segment_length) ? pnt->getID() : -1; + } + + // store for each point the line segment to which it was added. + const std::size_t nLines(new_lines->size()); + std::vector< std::vector<unsigned> > line_segment_map(nLines); + for (std::size_t i=0; i<nLines; ++i) + { + line_segment_map[i] = std::vector<unsigned>((*new_lines)[i]->getNumberOfPoints(),0); + std::iota(line_segment_map[i].begin(), line_segment_map[i].end(), 0); + } + + for (std::size_t i=0; i<nMeshNodes; ++i) + { + // if mesh node too far away or exactly at point position + if (closest_geo_point[i] == -1 || dist[i] < eps) continue; + + const MeshLib::Node* node (mesh->getNode(i)); + for (std::size_t l=0; l<nLines; ++l) + { + // find relevant polylines + if (!(*org_lines)[l]->isPointIDInPolyline(closest_geo_point[i])) continue; + + // find point position of closest geo point in original polyline + GeoLib::Polyline* ply ((*org_lines)[l]); + std::size_t nLinePnts ( ply->getNumberOfPoints() ); + std::size_t node_index_in_ply (0); + for (node_index_in_ply=0; node_index_in_ply<nLinePnts; ++node_index_in_ply) + if (ply->getPoint(node_index_in_ply) == (*points)[closest_geo_point[i]]) + break; + const GeoLib::Point* geo_point (ply->getPoint(node_index_in_ply)); + + // check if line segments connected to closest geo point intersect connected elements of current node + const std::vector<MeshLib::Element*> elements (node->getElements()); + const std::size_t nElems = elements.size(); + for (std::size_t e=0; e<nElems; ++e) + { + const unsigned nEdges (elements[e]->getNEdges()); + unsigned intersection_count (0); + + for (unsigned n=0; n<nEdges; ++n) + { + if (intersection_count>1) break; //already two intersections + + const MeshLib::Element* line = elements[e]->getEdge(n); + unsigned index_offset(0); // default: add to first line segment + GeoLib::Point* intersection (nullptr); + if (node_index_in_ply>0) // test line segment before closest point + intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply-1)); + if (intersection == nullptr && node_index_in_ply<(nLinePnts-1)) // test line segment after closest point + { + intersection = calcIntersection(line->getNode(0), line->getNode(1), geo_point, ply->getPoint(node_index_in_ply+1)); + index_offset = 1; // add to second segment + } + if (intersection) + { + intersection_count++; + unsigned start_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset-1;}))); + unsigned end_point_idx = static_cast<unsigned>(std::distance(line_segment_map[l].begin(), std::find_if(line_segment_map[l].begin(), line_segment_map[l].end(), [&node_index_in_ply, &index_offset](unsigned a){return a==node_index_in_ply+index_offset;}))); + std::size_t pos = getPointPosInLine((*new_lines)[l], start_point_idx, end_point_idx, intersection, eps); + + if (pos) + { + const std::size_t pnt_pos (new_points->size()); + new_points->push_back(intersection); + (*new_lines)[l]->insertPoint(pos, pnt_pos); + line_segment_map[l].insert(line_segment_map[l].begin()+pos, node_index_in_ply+index_offset-1); + } + } + } + } + } + } + + _geo_objects.addPointVec(std::move(new_points), const_cast<std::string&>(new_geo_name)); + std::vector<std::size_t> pnt_id_map = this->_geo_objects.getPointVecObj(new_geo_name)->getIDMap(); + for (auto & new_line : *new_lines) + new_line->updatePointIDs(pnt_id_map); + _geo_objects.addPolylineVec(std::move(new_lines), new_geo_name); + + // map new geometry incl. additional point using the normal mapping method + this->_geo_name = new_geo_name; + this->mapOnMesh(mesh); } GeoLib::Point* GeoMapper::calcIntersection(MathLib::Point3d const*const p1, MathLib::Point3d const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const { - const double x1 = (*p1)[0], x2 = (*p2)[0], x3 = (*q1)[0], x4 = (*q2)[0]; - const double y1 = (*p1)[1], y2 = (*p2)[1], y3 = (*q1)[1], y4 = (*q2)[1]; + const double x1 = (*p1)[0], x2 = (*p2)[0], x3 = (*q1)[0], x4 = (*q2)[0]; + const double y1 = (*p1)[1], y2 = (*p2)[1], y3 = (*q1)[1], y4 = (*q2)[1]; - const double det = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); - if (std::abs(det) < std::numeric_limits<double>::epsilon()) return nullptr; + const double det = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); + if (std::abs(det) < std::numeric_limits<double>::epsilon()) return nullptr; - const double pre = (x1*y2 - y1*x2); - const double post = (x3*y4 - y3*x4); - const double x = ( pre * (x3 - x4) - (x1 - x2) * post ) / det; - const double y = ( pre * (y3 - y4) - (y1 - y2) * post ) / det; + const double pre = (x1*y2 - y1*x2); + const double post = (x3*y4 - y3*x4); + const double x = ( pre * (x3 - x4) - (x1 - x2) * post ) / det; + const double y = ( pre * (y3 - y4) - (y1 - y2) * post ) / det; - // Check if the x and y coordinates are within both line segments - if (isPntInBoundingBox(x1,y1,x2,y2,x,y) && isPntInBoundingBox(x3,y3,x4,y4,x,y)) - return new GeoLib::Point(x, y, 0); - return nullptr; + // Check if the x and y coordinates are within both line segments + if (isPntInBoundingBox(x1,y1,x2,y2,x,y) && isPntInBoundingBox(x3,y3,x4,y4,x,y)) + return new GeoLib::Point(x, y, 0); + return nullptr; } unsigned GeoMapper::getPointPosInLine(GeoLib::Polyline const*const line, unsigned start, unsigned end, GeoLib::Point const*const point, double eps) const { - const double* first = line->getPoint(start)->getCoords(); - const double* pnt = point->getCoords(); - const double max_dist = MathLib::sqrDist(first, pnt); - - // if point is at start or end of line segment - if (max_dist<eps && MathLib::sqrDist(pnt, line->getPoint(end)->getCoords())) return 0; - - for (std::size_t i=start+1; i<end; ++i) - { - const double* current = (*line->getPoint(i)).getCoords(); - if (MathLib::sqrDist(pnt, current) < eps) return 0; - if (MathLib::sqrDist(first, current) > max_dist) return i; - } - return end; // last point of segment + const double* first = line->getPoint(start)->getCoords(); + const double* pnt = point->getCoords(); + const double max_dist = MathLib::sqrDist(first, pnt); + + // if point is at start or end of line segment + if (max_dist<eps && MathLib::sqrDist(pnt, line->getPoint(end)->getCoords())) return 0; + + for (std::size_t i=start+1; i<end; ++i) + { + const double* current = (*line->getPoint(i)).getCoords(); + if (MathLib::sqrDist(pnt, current) < eps) return 0; + if (MathLib::sqrDist(first, current) > max_dist) return i; + } + return end; // last point of segment } bool GeoMapper::isPntInBoundingBox(double ax, double ay, double bx, double by, double px, double py) const { - if ( px < (std::min(ax, bx)-std::numeric_limits<double>::epsilon()) || px > (std::max(ax, bx)+std::numeric_limits<double>::epsilon()) || - py < (std::min(ay, by)-std::numeric_limits<double>::epsilon()) || py > (std::max(ay, by)+std::numeric_limits<double>::epsilon()) ) - return false; - return true; + if ( px < (std::min(ax, bx)-std::numeric_limits<double>::epsilon()) || px > (std::max(ax, bx)+std::numeric_limits<double>::epsilon()) || + py < (std::min(ay, by)-std::numeric_limits<double>::epsilon()) || py > (std::max(ay, by)+std::numeric_limits<double>::epsilon()) ) + return false; + return true; } double GeoMapper::getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &lines) const { - double max_segment_length (0); - const std::size_t nPlys ( lines.size() ); - for (std::size_t i=0; i<nPlys; ++i) - { - const GeoLib::Polyline* line = lines[i]; - const std::size_t nPlyPoints = line->getNumberOfPoints(); - for (std::size_t j=1; j<nPlyPoints; ++j) - { - const double dist (line->getLength(j)-line->getLength(j-1)); - if (dist>max_segment_length) - max_segment_length=dist; - } - } - return max_segment_length; + double max_segment_length (0); + const std::size_t nPlys ( lines.size() ); + for (std::size_t i=0; i<nPlys; ++i) + { + const GeoLib::Polyline* line = lines[i]; + const std::size_t nPlyPoints = line->getNumberOfPoints(); + for (std::size_t j=1; j<nPlyPoints; ++j) + { + const double dist (line->getLength(j)-line->getLength(j-1)); + if (dist>max_segment_length) + max_segment_length=dist; + } + } + return max_segment_length; } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/GeoMapper.h b/MeshGeoToolsLib/GeoMapper.h index f1e93fecc078b18fac6df440985f9d7075296a02..601350b33c5b9f3f97ce3dab36380571afaf3c98 100644 --- a/MeshGeoToolsLib/GeoMapper.h +++ b/MeshGeoToolsLib/GeoMapper.h @@ -24,12 +24,12 @@ #include "MathLib/Point3d.h" namespace MeshLib { - class Mesh; - class Node; + class Mesh; + class Node; } namespace GeoLib { - class Raster; + class Raster; } namespace MeshGeoToolsLib { @@ -40,73 +40,73 @@ namespace MeshGeoToolsLib { class GeoMapper { public: - GeoMapper(GeoLib::GEOObjects &geo_objects, const std::string &geo_name); - ~GeoMapper(); - - /// Maps geometry based on a raster file - void mapOnDEM(const std::string &file_name); - - /// Maps geometry based on a mesh file - void mapOnMesh(const std::string &file_name); - - /** - * Maps the geometry based on the given mesh file. The elevation value of all geometric - * points are modified such that they are located on the mesh surface. - */ - void mapOnMesh(const MeshLib::Mesh* mesh); - - /// Maps geometry to a constant elevation value - void mapToConstantValue(double value); - - /** - * Maps the geometry based on the given mesh file. A new geometry is created where all - * geometric points are assigned an elevation value on the mesh surface. Additional points are - * inserted whenever a polyline from the original geometry intersects a mesh node or the edge - * of a mesh element. - * \param mesh Mesh the geometry is mapped on - * \param new_geo_name Name of the newly created geometry - * \result A new geometry with the given name is inserted into _geo_objects - */ - void advancedMapOnMesh(const MeshLib::Mesh* mesh, const std::string &new_geo_name); + GeoMapper(GeoLib::GEOObjects &geo_objects, const std::string &geo_name); + ~GeoMapper(); + + /// Maps geometry based on a raster file + void mapOnDEM(const std::string &file_name); + + /// Maps geometry based on a mesh file + void mapOnMesh(const std::string &file_name); + + /** + * Maps the geometry based on the given mesh file. The elevation value of all geometric + * points are modified such that they are located on the mesh surface. + */ + void mapOnMesh(const MeshLib::Mesh* mesh); + + /// Maps geometry to a constant elevation value + void mapToConstantValue(double value); + + /** + * Maps the geometry based on the given mesh file. A new geometry is created where all + * geometric points are assigned an elevation value on the mesh surface. Additional points are + * inserted whenever a polyline from the original geometry intersects a mesh node or the edge + * of a mesh element. + * \param mesh Mesh the geometry is mapped on + * \param new_geo_name Name of the newly created geometry + * \result A new geometry with the given name is inserted into _geo_objects + */ + void advancedMapOnMesh(const MeshLib::Mesh* mesh, const std::string &new_geo_name); private: - /// Mapping stations, boreholes on a raster or mesh. - void mapStationData(std::vector<GeoLib::Point*> const& points); + /// Mapping stations, boreholes on a raster or mesh. + void mapStationData(std::vector<GeoLib::Point*> const& points); - /// Mapping points on a raster. - void mapPointDataToDEM(std::vector<GeoLib::Point*> const& points); + /// Mapping points on a raster. + void mapPointDataToDEM(std::vector<GeoLib::Point*> const& points); - /// Mapping points on mesh. - void mapPointDataToMeshSurface(std::vector<GeoLib::Point*> const& points); + /// Mapping points on mesh. + void mapPointDataToMeshSurface(std::vector<GeoLib::Point*> const& points); - /// Returns the elevation at Point (x,y) based on a mesh. This uses collision detection for triangles and nearest neighbor for quads. - /// NOTE: This medhod only returns correct values if the node numbering of the elements is correct! - double getMeshElevation(double x, double y, double min_val, double max_val) const; + /// Returns the elevation at Point (x,y) based on a mesh. This uses collision detection for triangles and nearest neighbor for quads. + /// NOTE: This medhod only returns correct values if the node numbering of the elements is correct! + double getMeshElevation(double x, double y, double min_val, double max_val) const; - /// Returns the elevation at Point (x,y) based on a raster - float getDemElevation(GeoLib::Point const& pnt) const; + /// Returns the elevation at Point (x,y) based on a raster + float getDemElevation(GeoLib::Point const& pnt) const; - /// Calculates the intersection of two lines embedded in the xy-plane - GeoLib::Point* calcIntersection(MathLib::Point3d const*const p1, MathLib::Point3d const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const; + /// Calculates the intersection of two lines embedded in the xy-plane + GeoLib::Point* calcIntersection(MathLib::Point3d const*const p1, MathLib::Point3d const*const p2, GeoLib::Point const*const q1, GeoLib::Point const*const q2) const; - /// Returns the position of a point within a line-segment - unsigned getPointPosInLine(GeoLib::Polyline const*const line, unsigned start, unsigned end, GeoLib::Point const*const point, double eps) const; + /// Returns the position of a point within a line-segment + unsigned getPointPosInLine(GeoLib::Polyline const*const line, unsigned start, unsigned end, GeoLib::Point const*const point, double eps) const; - /// Returns the maximum segment length in a polyline vector - double getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &lines) const; + /// Returns the maximum segment length in a polyline vector + double getMaxSegmentLength(const std::vector<GeoLib::Polyline*> &lines) const; - /// Returns if a point p is within a bounding box defined by a and b - bool isPntInBoundingBox(double ax, double ay, double bx, double by, double px, double py) const; + /// Returns if a point p is within a bounding box defined by a and b + bool isPntInBoundingBox(double ax, double ay, double bx, double by, double px, double py) const; - GeoLib::GEOObjects& _geo_objects; - std::string& _geo_name; + GeoLib::GEOObjects& _geo_objects; + std::string& _geo_name; - /// only necessary for mapping on mesh - MeshLib::Mesh* _surface_mesh; - GeoLib::Grid<MeshLib::Node>* _grid; + /// only necessary for mapping on mesh + MeshLib::Mesh* _surface_mesh; + GeoLib::Grid<MeshLib::Node>* _grid; - /// only necessary for mapping on DEM - GeoLib::Raster *_raster; + /// only necessary for mapping on DEM + GeoLib::Raster *_raster; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/HeuristicSearchLength.cpp b/MeshGeoToolsLib/HeuristicSearchLength.cpp index 12398609681255a08da623052fd0a324acd5d21b..915f96d29b08a5231a0e6f8b37ee8c213b3ec4cd 100644 --- a/MeshGeoToolsLib/HeuristicSearchLength.cpp +++ b/MeshGeoToolsLib/HeuristicSearchLength.cpp @@ -22,51 +22,51 @@ namespace MeshGeoToolsLib HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthType length_type) : _mesh(mesh) { - double sum (0.0); - double sum_of_sqr (0.0); // total length of edges + double sum (0.0); + double sum_of_sqr (0.0); // total length of edges - std::size_t n_sampling(0); // total length of edges squared - std::vector<MeshLib::Element*> const& elements(_mesh.getElements()); + std::size_t n_sampling(0); // total length of edges squared + std::vector<MeshLib::Element*> const& elements(_mesh.getElements()); - if (length_type==LengthType::Edge) { - for (std::vector<MeshLib::Element*>::const_iterator it(elements.cbegin()); - it != elements.cend(); ++it) { - std::size_t const n_edges((*it)->getNEdges()); - for (std::size_t k(0); k<n_edges; k++) { - auto edge = (*it)->getEdge(k); // allocation inside getEdge(). - double const len = edge->getContent(); - delete edge; - sum += len; - sum_of_sqr += len*len; - } - n_sampling += n_edges; - } - } else { - double min=0, max=0; - for (const MeshLib::Element* e : elements) { - e->computeSqrNodeDistanceRange(min, max, true); - sum += std::sqrt(min); - sum_of_sqr += min; - } - n_sampling = _mesh.getNElements(); - } + if (length_type==LengthType::Edge) { + for (std::vector<MeshLib::Element*>::const_iterator it(elements.cbegin()); + it != elements.cend(); ++it) { + std::size_t const n_edges((*it)->getNEdges()); + for (std::size_t k(0); k<n_edges; k++) { + auto edge = (*it)->getEdge(k); // allocation inside getEdge(). + double const len = edge->getContent(); + delete edge; + sum += len; + sum_of_sqr += len*len; + } + n_sampling += n_edges; + } + } else { + double min=0, max=0; + for (const MeshLib::Element* e : elements) { + e->computeSqrNodeDistanceRange(min, max, true); + sum += std::sqrt(min); + sum_of_sqr += min; + } + n_sampling = _mesh.getNElements(); + } - const double mean (sum/n_sampling); - const double variance ((sum_of_sqr - (sum*sum)/n_sampling)/(n_sampling-1)); + const double mean (sum/n_sampling); + const double variance ((sum_of_sqr - (sum*sum)/n_sampling)/(n_sampling-1)); - // Set the search length for the case of non-positive variance (which can - // happen due to numerics). - _search_length = mean/2; + // Set the search length for the case of non-positive variance (which can + // happen due to numerics). + _search_length = mean/2; - if (variance > 0) { - if (variance < mean*mean/4) - _search_length -= std::sqrt(variance); - else - _search_length = std::numeric_limits<double>::epsilon(); - } + if (variance > 0) { + if (variance < mean*mean/4) + _search_length -= std::sqrt(variance); + else + _search_length = std::numeric_limits<double>::epsilon(); + } - DBUG("[MeshNodeSearcher::MeshNodeSearcher] Calculated search length for mesh \"%s\" is %f.", - _mesh.getName().c_str(), _search_length); + DBUG("[MeshNodeSearcher::MeshNodeSearcher] Calculated search length for mesh \"%s\" is %f.", + _mesh.getName().c_str(), _search_length); } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/HeuristicSearchLength.h b/MeshGeoToolsLib/HeuristicSearchLength.h index a5d8538a96c3d23cb1e1cf030ee6ed80f4e87b00..0714e1c26cadd5143058eb8b4d1f994c3fff6e76 100644 --- a/MeshGeoToolsLib/HeuristicSearchLength.h +++ b/MeshGeoToolsLib/HeuristicSearchLength.h @@ -30,21 +30,21 @@ namespace MeshGeoToolsLib class HeuristicSearchLength : public SearchLength { public: - /// Type of length to be sampled - enum class LengthType - { - Edge, /// edge length of elements, which is recommended for meshes without nonlinear nodes - Node /// distance between nodes - }; - - /** - * Constructor - * @param mesh mesh object - * @param sampled_len length type to be sampled - */ - HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthType length_type = LengthType::Edge); + /// Type of length to be sampled + enum class LengthType + { + Edge, /// edge length of elements, which is recommended for meshes without nonlinear nodes + Node /// distance between nodes + }; + + /** + * Constructor + * @param mesh mesh object + * @param sampled_len length type to be sampled + */ + HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthType length_type = LengthType::Edge); private: - MeshLib::Mesh const& _mesh; + MeshLib::Mesh const& _mesh; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp index 5561a2b88cafb204a9a69ba783a49cdf8b230c0a..b4c1095d8e212d671439c9a7ca64b9fd7a3fda00 100644 --- a/MeshGeoToolsLib/MeshNodeSearcher.cpp +++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp @@ -35,129 +35,129 @@ std::vector<std::unique_ptr<MeshNodeSearcher>> MeshNodeSearcher::_mesh_node_sear MeshNodeSearcher::MeshNodeSearcher(MeshLib::Mesh const& mesh, - MeshGeoToolsLib::SearchLength const& search_length_algorithm, bool search_all_nodes) : - _mesh(mesh), _mesh_grid(_mesh.getNodes().cbegin(), _mesh.getNodes().cend()), - _search_length(0.0), _search_all_nodes(search_all_nodes) + MeshGeoToolsLib::SearchLength const& search_length_algorithm, bool search_all_nodes) : + _mesh(mesh), _mesh_grid(_mesh.getNodes().cbegin(), _mesh.getNodes().cend()), + _search_length(0.0), _search_all_nodes(search_all_nodes) { - DBUG("Constructing MeshNodeSearcher obj."); - _search_length = search_length_algorithm.getSearchLength(); + DBUG("Constructing MeshNodeSearcher obj."); + _search_length = search_length_algorithm.getSearchLength(); - DBUG("Calculated search length for mesh \"%s\" is %e.", - _mesh.getName().c_str(), _search_length); + DBUG("Calculated search length for mesh \"%s\" is %e.", + _mesh.getName().c_str(), _search_length); } MeshNodeSearcher::~MeshNodeSearcher() { - for (auto pointer : _mesh_nodes_on_points) - delete pointer; - for (auto pointer : _mesh_nodes_along_polylines) - delete pointer; - for (auto pointer : _mesh_nodes_along_surfaces) - delete pointer; + for (auto pointer : _mesh_nodes_on_points) + delete pointer; + for (auto pointer : _mesh_nodes_along_polylines) + delete pointer; + for (auto pointer : _mesh_nodes_along_surfaces) + delete pointer; } std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(GeoLib::GeoObject const& geoObj) { - std::vector<std::size_t> vec_nodes; - switch (geoObj.getGeoType()) { - case GeoLib::GEOTYPE::POINT: - { - vec_nodes = this->getMeshNodeIDsForPoint(*static_cast<const GeoLib::Point*>(&geoObj)); - break; - } - case GeoLib::GEOTYPE::POLYLINE: - vec_nodes = this->getMeshNodeIDsAlongPolyline(*static_cast<const GeoLib::Polyline*>(&geoObj)); - break; - case GeoLib::GEOTYPE::SURFACE: - vec_nodes = this->getMeshNodeIDsAlongSurface(*static_cast<const GeoLib::Surface*>(&geoObj)); - break; - default: - break; - } - return vec_nodes; + std::vector<std::size_t> vec_nodes; + switch (geoObj.getGeoType()) { + case GeoLib::GEOTYPE::POINT: + { + vec_nodes = this->getMeshNodeIDsForPoint(*static_cast<const GeoLib::Point*>(&geoObj)); + break; + } + case GeoLib::GEOTYPE::POLYLINE: + vec_nodes = this->getMeshNodeIDsAlongPolyline(*static_cast<const GeoLib::Polyline*>(&geoObj)); + break; + case GeoLib::GEOTYPE::SURFACE: + vec_nodes = this->getMeshNodeIDsAlongSurface(*static_cast<const GeoLib::Surface*>(&geoObj)); + break; + default: + break; + } + return vec_nodes; } std::vector<std::size_t> const& MeshNodeSearcher::getMeshNodeIDsForPoint(GeoLib::Point const& pnt) { - return getMeshNodesOnPoint(pnt).getNodeIDs(); + return getMeshNodesOnPoint(pnt).getNodeIDs(); } std::vector<std::size_t> const& MeshNodeSearcher::getMeshNodeIDsAlongPolyline( - GeoLib::Polyline const& ply) + GeoLib::Polyline const& ply) { - return getMeshNodesAlongPolyline(ply).getNodeIDs(); + return getMeshNodesAlongPolyline(ply).getNodeIDs(); } std::vector<std::size_t> const& MeshNodeSearcher::getMeshNodeIDsAlongSurface(GeoLib::Surface const& sfc) { - return getMeshNodesAlongSurface(sfc).getNodeIDs(); + return getMeshNodesAlongSurface(sfc).getNodeIDs(); } MeshNodesOnPoint& MeshNodeSearcher::getMeshNodesOnPoint(GeoLib::Point const& pnt) { - std::vector<MeshNodesOnPoint*>::const_iterator it(_mesh_nodes_on_points.begin()); - for (; it != _mesh_nodes_on_points.end(); ++it) { - if (&(*it)->getPoint() == &pnt) { - return *(*it); - } - } - - _mesh_nodes_on_points.push_back( - new MeshNodesOnPoint(_mesh, _mesh_grid, pnt, _search_length, _search_all_nodes)); - return *_mesh_nodes_on_points.back(); + std::vector<MeshNodesOnPoint*>::const_iterator it(_mesh_nodes_on_points.begin()); + for (; it != _mesh_nodes_on_points.end(); ++it) { + if (&(*it)->getPoint() == &pnt) { + return *(*it); + } + } + + _mesh_nodes_on_points.push_back( + new MeshNodesOnPoint(_mesh, _mesh_grid, pnt, _search_length, _search_all_nodes)); + return *_mesh_nodes_on_points.back(); } MeshNodesAlongPolyline& MeshNodeSearcher::getMeshNodesAlongPolyline(GeoLib::Polyline const& ply) { - std::vector<MeshNodesAlongPolyline*>::const_iterator it(_mesh_nodes_along_polylines.begin()); - for (; it != _mesh_nodes_along_polylines.end(); ++it) { - if (&(*it)->getPolyline() == &ply) { - // we calculated mesh nodes for this polyline already - return *(*it); - } - } - - // compute nodes (and supporting points) along polyline - _mesh_nodes_along_polylines.push_back( - new MeshNodesAlongPolyline(_mesh, ply, _search_length, _search_all_nodes)); - return *_mesh_nodes_along_polylines.back(); + std::vector<MeshNodesAlongPolyline*>::const_iterator it(_mesh_nodes_along_polylines.begin()); + for (; it != _mesh_nodes_along_polylines.end(); ++it) { + if (&(*it)->getPolyline() == &ply) { + // we calculated mesh nodes for this polyline already + return *(*it); + } + } + + // compute nodes (and supporting points) along polyline + _mesh_nodes_along_polylines.push_back( + new MeshNodesAlongPolyline(_mesh, ply, _search_length, _search_all_nodes)); + return *_mesh_nodes_along_polylines.back(); } MeshNodesAlongSurface& MeshNodeSearcher::getMeshNodesAlongSurface(GeoLib::Surface const& sfc) { - std::vector<MeshNodesAlongSurface*>::const_iterator it(_mesh_nodes_along_surfaces.begin()); - for (; it != _mesh_nodes_along_surfaces.end(); ++it) { - if (&(*it)->getSurface() == &sfc) { - // we calculated mesh nodes for this polyline already - return *(*it); - } - } - - // compute nodes (and supporting points) along polyline - _mesh_nodes_along_surfaces.push_back( - new MeshNodesAlongSurface(_mesh, sfc, _search_all_nodes)); - return *_mesh_nodes_along_surfaces.back(); + std::vector<MeshNodesAlongSurface*>::const_iterator it(_mesh_nodes_along_surfaces.begin()); + for (; it != _mesh_nodes_along_surfaces.end(); ++it) { + if (&(*it)->getSurface() == &sfc) { + // we calculated mesh nodes for this polyline already + return *(*it); + } + } + + // compute nodes (and supporting points) along polyline + _mesh_nodes_along_surfaces.push_back( + new MeshNodesAlongSurface(_mesh, sfc, _search_all_nodes)); + return *_mesh_nodes_along_surfaces.back(); } MeshNodeSearcher& MeshNodeSearcher::getMeshNodeSearcher(MeshLib::Mesh const& mesh) { - std::size_t const mesh_id = mesh.getID(); - if (_mesh_node_searchers.size() < mesh_id+1) - _mesh_node_searchers.resize(mesh_id+1); + std::size_t const mesh_id = mesh.getID(); + if (_mesh_node_searchers.size() < mesh_id+1) + _mesh_node_searchers.resize(mesh_id+1); - if (!_mesh_node_searchers[mesh_id]) - _mesh_node_searchers[mesh_id].reset( - new MeshGeoToolsLib::MeshNodeSearcher(mesh)); + if (!_mesh_node_searchers[mesh_id]) + _mesh_node_searchers[mesh_id].reset( + new MeshGeoToolsLib::MeshNodeSearcher(mesh)); - return *_mesh_node_searchers[mesh_id]; + return *_mesh_node_searchers[mesh_id]; } std::size_t MeshNodeSearcher::getMeshId() const { - return _mesh.getID(); + return _mesh.getID(); } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodeSearcher.h b/MeshGeoToolsLib/MeshNodeSearcher.h index f98de563a31a1520635b92a0dd90d6cfbc8db8f7..8957d80d690afc6676bfae9b9dc960457cd645a7 100644 --- a/MeshGeoToolsLib/MeshNodeSearcher.h +++ b/MeshGeoToolsLib/MeshNodeSearcher.h @@ -53,100 +53,100 @@ namespace MeshGeoToolsLib class MeshNodeSearcher { public: - /** - * Constructor for objects of class MeshNodeSearcher. It calculates - * internally used search length from the given MeshLib::Mesh. - * @param mesh Run search algorithm on this mesh. It is assumed - * that the mesh does not change its geometry. - */ - explicit MeshNodeSearcher(MeshLib::Mesh const& mesh, - MeshGeoToolsLib::SearchLength const& search_length_algorithm - = MeshGeoToolsLib::SearchLength(), - bool search_all_nodes = true); - - virtual ~MeshNodeSearcher(); - - /** - * Searches for the nearest mesh nodes on the given geometric object (point, polyline, surface). - * @param geoObj a GeoLib::GeoObject where the nearest mesh node is searched for - * @return a vector of mesh node ids - */ - std::vector<std::size_t> getMeshNodeIDs(GeoLib::GeoObject const& geoObj); - - /** - * Searches for the node nearest by the given point. If there are two nodes - * with the same distance the id of the one that was first found will be - * returned. The algorithm for the search is using GeoLib::Grid data - * structure. - * @param pnt a GeoLib::Point the nearest mesh node is searched for - * @return a vector of mesh node ids - */ - std::vector<std::size_t> const& getMeshNodeIDsForPoint(GeoLib::Point const& pnt); - - /** - * Searches for the nearest mesh nodes along a GeoLib::Polyline. - * The search for mesh nodes along a specific polyline will be performed - * only once. The result ids will be stored inside an object - * (@see class MeshGeoToolsLib::MeshNodesAlongPolyline). - * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for - * @return a vector of mesh node ids - */ - std::vector<std::size_t> const& getMeshNodeIDsAlongPolyline(GeoLib::Polyline const& ply); - - /** - * Searches for the nearest mesh nodes along a GeoLib::Surface. - * The search for mesh nodes along a specific surface will be performed - * only once. The result ids will be stored inside an object - * (@see class MeshGeoToolsLib::MeshNodesAlongSurface). - * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for - * @return a vector of mesh node ids - */ - std::vector<std::size_t> const& getMeshNodeIDsAlongSurface(GeoLib::Surface const& sfc); - - /** - * Return a MeshNodesOnPoint object for the given GeoLib::Point object. - * @param pnt the GeoLib::Point the nearest mesh nodes are searched for - * @return a reference to a MeshNodesOnPoint object - */ - MeshNodesOnPoint& getMeshNodesOnPoint(GeoLib::Point const& pnt); - - /** - * Return a MeshNodesAlongPolyline object for the given GeoLib::Polyline object. - * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for - * @return a reference to a MeshNodesAlongPolyline object - */ - MeshNodesAlongPolyline& getMeshNodesAlongPolyline(GeoLib::Polyline const& ply); - - /** - * Return a MeshNodesAlongSurface object for the given GeoLib::Surface object. - * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for - * @return a reference to a MeshNodesAlongSurface object - */ - MeshNodesAlongSurface& getMeshNodesAlongSurface(GeoLib::Surface const& sfc); - - /** - * Get the mesh this searcher operates on. - */ - std::size_t getMeshId() const; - - /** - * Returns a (possibly new) mesh node searcher for the mesh. - * A new one will be created, if it does not already exists. - */ - static MeshNodeSearcher& getMeshNodeSearcher(MeshLib::Mesh const& mesh); + /** + * Constructor for objects of class MeshNodeSearcher. It calculates + * internally used search length from the given MeshLib::Mesh. + * @param mesh Run search algorithm on this mesh. It is assumed + * that the mesh does not change its geometry. + */ + explicit MeshNodeSearcher(MeshLib::Mesh const& mesh, + MeshGeoToolsLib::SearchLength const& search_length_algorithm + = MeshGeoToolsLib::SearchLength(), + bool search_all_nodes = true); + + virtual ~MeshNodeSearcher(); + + /** + * Searches for the nearest mesh nodes on the given geometric object (point, polyline, surface). + * @param geoObj a GeoLib::GeoObject where the nearest mesh node is searched for + * @return a vector of mesh node ids + */ + std::vector<std::size_t> getMeshNodeIDs(GeoLib::GeoObject const& geoObj); + + /** + * Searches for the node nearest by the given point. If there are two nodes + * with the same distance the id of the one that was first found will be + * returned. The algorithm for the search is using GeoLib::Grid data + * structure. + * @param pnt a GeoLib::Point the nearest mesh node is searched for + * @return a vector of mesh node ids + */ + std::vector<std::size_t> const& getMeshNodeIDsForPoint(GeoLib::Point const& pnt); + + /** + * Searches for the nearest mesh nodes along a GeoLib::Polyline. + * The search for mesh nodes along a specific polyline will be performed + * only once. The result ids will be stored inside an object + * (@see class MeshGeoToolsLib::MeshNodesAlongPolyline). + * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for + * @return a vector of mesh node ids + */ + std::vector<std::size_t> const& getMeshNodeIDsAlongPolyline(GeoLib::Polyline const& ply); + + /** + * Searches for the nearest mesh nodes along a GeoLib::Surface. + * The search for mesh nodes along a specific surface will be performed + * only once. The result ids will be stored inside an object + * (@see class MeshGeoToolsLib::MeshNodesAlongSurface). + * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for + * @return a vector of mesh node ids + */ + std::vector<std::size_t> const& getMeshNodeIDsAlongSurface(GeoLib::Surface const& sfc); + + /** + * Return a MeshNodesOnPoint object for the given GeoLib::Point object. + * @param pnt the GeoLib::Point the nearest mesh nodes are searched for + * @return a reference to a MeshNodesOnPoint object + */ + MeshNodesOnPoint& getMeshNodesOnPoint(GeoLib::Point const& pnt); + + /** + * Return a MeshNodesAlongPolyline object for the given GeoLib::Polyline object. + * @param ply the GeoLib::Polyline the nearest mesh nodes are searched for + * @return a reference to a MeshNodesAlongPolyline object + */ + MeshNodesAlongPolyline& getMeshNodesAlongPolyline(GeoLib::Polyline const& ply); + + /** + * Return a MeshNodesAlongSurface object for the given GeoLib::Surface object. + * @param sfc the GeoLib::Surface the nearest mesh nodes are searched for + * @return a reference to a MeshNodesAlongSurface object + */ + MeshNodesAlongSurface& getMeshNodesAlongSurface(GeoLib::Surface const& sfc); + + /** + * Get the mesh this searcher operates on. + */ + std::size_t getMeshId() const; + + /** + * Returns a (possibly new) mesh node searcher for the mesh. + * A new one will be created, if it does not already exists. + */ + static MeshNodeSearcher& getMeshNodeSearcher(MeshLib::Mesh const& mesh); private: - MeshLib::Mesh const& _mesh; - GeoLib::Grid<MeshLib::Node> _mesh_grid; - double _search_length; - bool _search_all_nodes; - // with newer compiler we can omit to use a pointer here - std::vector<MeshNodesOnPoint*> _mesh_nodes_on_points; - std::vector<MeshNodesAlongPolyline*> _mesh_nodes_along_polylines; - std::vector<MeshNodesAlongSurface*> _mesh_nodes_along_surfaces; - - /// Mesh node searcher for the meshes indexed by the meshs' ids. - static std::vector<std::unique_ptr<MeshNodeSearcher>> _mesh_node_searchers; + MeshLib::Mesh const& _mesh; + GeoLib::Grid<MeshLib::Node> _mesh_grid; + double _search_length; + bool _search_all_nodes; + // with newer compiler we can omit to use a pointer here + std::vector<MeshNodesOnPoint*> _mesh_nodes_on_points; + std::vector<MeshNodesAlongPolyline*> _mesh_nodes_along_polylines; + std::vector<MeshNodesAlongSurface*> _mesh_nodes_along_surfaces; + + /// Mesh node searcher for the meshes indexed by the meshs' ids. + static std::vector<std::unique_ptr<MeshNodeSearcher>> _mesh_node_searchers; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp b/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp index e428ebfd6f602d60b2df1c25c81483253396a86e..439e2eae9216d2de71b970f8a08750c05b6c19fc 100644 --- a/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp +++ b/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp @@ -23,46 +23,46 @@ namespace MeshGeoToolsLib { MeshNodesAlongPolyline::MeshNodesAlongPolyline( - MeshLib::Mesh const& mesh, - GeoLib::Polyline const& ply, - double epsilon_radius, - bool search_all_nodes) : - _mesh(mesh), _ply(ply) + MeshLib::Mesh const& mesh, + GeoLib::Polyline const& ply, + double epsilon_radius, + bool search_all_nodes) : + _mesh(mesh), _ply(ply) { - assert(epsilon_radius > 0); - const std::size_t n_nodes (search_all_nodes ? _mesh.getNNodes() : _mesh.getNBaseNodes()); - auto &mesh_nodes = _mesh.getNodes(); - // loop over all nodes - for (std::size_t i = 0; i < n_nodes; i++) { - double dist = _ply.getDistanceAlongPolyline(*mesh_nodes[i], epsilon_radius); - if (dist >= 0.0) { - _msh_node_ids.push_back(mesh_nodes[i]->getID()); - _dist_of_proj_node_from_ply_start.push_back(dist); - } - } + assert(epsilon_radius > 0); + const std::size_t n_nodes (search_all_nodes ? _mesh.getNNodes() : _mesh.getNBaseNodes()); + auto &mesh_nodes = _mesh.getNodes(); + // loop over all nodes + for (std::size_t i = 0; i < n_nodes; i++) { + double dist = _ply.getDistanceAlongPolyline(*mesh_nodes[i], epsilon_radius); + if (dist >= 0.0) { + _msh_node_ids.push_back(mesh_nodes[i]->getID()); + _dist_of_proj_node_from_ply_start.push_back(dist); + } + } - // sort the nodes along the polyline according to their distances - BaseLib::quicksort<double> (_dist_of_proj_node_from_ply_start, 0, - _dist_of_proj_node_from_ply_start.size(), _msh_node_ids); + // sort the nodes along the polyline according to their distances + BaseLib::quicksort<double> (_dist_of_proj_node_from_ply_start, 0, + _dist_of_proj_node_from_ply_start.size(), _msh_node_ids); } MeshLib::Mesh const& MeshNodesAlongPolyline::getMesh () const { - return _mesh; + return _mesh; } std::vector<std::size_t> const& MeshNodesAlongPolyline::getNodeIDs () const { - return _msh_node_ids; + return _msh_node_ids; } GeoLib::Polyline const& MeshNodesAlongPolyline::getPolyline () const { - return _ply; + return _ply; } std::vector<double> const& MeshNodesAlongPolyline::getDistOfProjNodeFromPlyStart() const { - return _dist_of_proj_node_from_ply_start; + return _dist_of_proj_node_from_ply_start; } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesAlongPolyline.h b/MeshGeoToolsLib/MeshNodesAlongPolyline.h index d86b1c21aef555c77f3679e859037393bd381793..919f93dd5f039da576f678a175a86d9568b8ec1e 100644 --- a/MeshGeoToolsLib/MeshNodesAlongPolyline.h +++ b/MeshGeoToolsLib/MeshNodesAlongPolyline.h @@ -36,42 +36,42 @@ namespace MeshGeoToolsLib class MeshNodesAlongPolyline { public: - /** - * Constructor of object, that search mesh nodes along a - * GeoLib::Polyline polyline within a given search radius. So the polyline - * is something like a tube. - * @param mesh_nodes Nodes the search will be performed on. - * @param ply Along the GeoLib::Polyline ply the mesh nodes are searched. - * @param epsilon_radius Search / tube radius - */ - MeshNodesAlongPolyline(MeshLib::Mesh const& mesh, - GeoLib::Polyline const& ply, double epsilon_radius, bool search_all_nodes = true); + /** + * Constructor of object, that search mesh nodes along a + * GeoLib::Polyline polyline within a given search radius. So the polyline + * is something like a tube. + * @param mesh_nodes Nodes the search will be performed on. + * @param ply Along the GeoLib::Polyline ply the mesh nodes are searched. + * @param epsilon_radius Search / tube radius + */ + MeshNodesAlongPolyline(MeshLib::Mesh const& mesh, + GeoLib::Polyline const& ply, double epsilon_radius, bool search_all_nodes = true); - /// return the mesh object - MeshLib::Mesh const& getMesh() const; + /// return the mesh object + MeshLib::Mesh const& getMesh() const; - /** - * Access the vector of mesh node ids. - * @return The vector of mesh node ids calculated in the constructor - */ - std::vector<std::size_t> const& getNodeIDs () const; + /** + * Access the vector of mesh node ids. + * @return The vector of mesh node ids calculated in the constructor + */ + std::vector<std::size_t> const& getNodeIDs () const; - /** - * Deploying this method the user can get access to the underlying - * GeoLib::Polyline. - * @return the underlying GeoLib::Polyline - */ - GeoLib::Polyline const& getPolyline () const; + /** + * Deploying this method the user can get access to the underlying + * GeoLib::Polyline. + * @return the underlying GeoLib::Polyline + */ + GeoLib::Polyline const& getPolyline () const; - /// return a vector of node distance from the polyline start. The vector - /// corresponds to a vector returned in getNodeIDs() - std::vector<double> const & getDistOfProjNodeFromPlyStart() const; + /// return a vector of node distance from the polyline start. The vector + /// corresponds to a vector returned in getNodeIDs() + std::vector<double> const & getDistOfProjNodeFromPlyStart() const; private: - MeshLib::Mesh const& _mesh; - GeoLib::Polyline const& _ply; - std::vector<std::size_t> _msh_node_ids; - std::vector<double> _dist_of_proj_node_from_ply_start; + MeshLib::Mesh const& _mesh; + GeoLib::Polyline const& _ply; + std::vector<std::size_t> _msh_node_ids; + std::vector<double> _dist_of_proj_node_from_ply_start; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp index 3b9f0e890927ddc1bc9ed45535ab584631ed3100..173ff3c9cefd7a6b4bea0185b0fa0fc1dad40498 100644 --- a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp +++ b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp @@ -24,37 +24,37 @@ namespace MeshGeoToolsLib { MeshNodesAlongSurface::MeshNodesAlongSurface( - MeshLib::Mesh const& mesh, - GeoLib::Surface const& sfc, - bool search_all_nodes) : - _mesh(mesh), _sfc(sfc) + MeshLib::Mesh const& mesh, + GeoLib::Surface const& sfc, + bool search_all_nodes) : + _mesh(mesh), _sfc(sfc) { - auto& mesh_nodes = _mesh.getNodes(); - const std::size_t n_nodes (search_all_nodes ? _mesh.getNNodes() : _mesh.getNBaseNodes()); - // loop over all nodes - for (std::size_t i = 0; i < n_nodes; i++) { - auto* node = mesh_nodes[i]; - if (!sfc.isPntInBoundingVolume(*node)) - continue; - if (sfc.isPntInSfc(*node)) { - _msh_node_ids.push_back(node->getID()); - } - } + auto& mesh_nodes = _mesh.getNodes(); + const std::size_t n_nodes (search_all_nodes ? _mesh.getNNodes() : _mesh.getNBaseNodes()); + // loop over all nodes + for (std::size_t i = 0; i < n_nodes; i++) { + auto* node = mesh_nodes[i]; + if (!sfc.isPntInBoundingVolume(*node)) + continue; + if (sfc.isPntInSfc(*node)) { + _msh_node_ids.push_back(node->getID()); + } + } } MeshLib::Mesh const& MeshNodesAlongSurface::getMesh () const { - return _mesh; + return _mesh; } std::vector<std::size_t> const& MeshNodesAlongSurface::getNodeIDs () const { - return _msh_node_ids; + return _msh_node_ids; } GeoLib::Surface const& MeshNodesAlongSurface::getSurface () const { - return _sfc; + return _sfc; } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesAlongSurface.h b/MeshGeoToolsLib/MeshNodesAlongSurface.h index aeda4717223a4a137b924e4717698b3b2c51ee46..5707f221317d8df1e9b0861f434f76e2c6ab6367 100644 --- a/MeshGeoToolsLib/MeshNodesAlongSurface.h +++ b/MeshGeoToolsLib/MeshNodesAlongSurface.h @@ -33,34 +33,34 @@ namespace MeshGeoToolsLib class MeshNodesAlongSurface { public: - /** - * Constructor of object, that search mesh nodes along a - * GeoLib::Surface object within a given search radius. - * @param mesh_nodes Nodes the search will be performed on. - * @param sfc Along the GeoLib::Surface sfc the mesh nodes are searched. - */ - MeshNodesAlongSurface(MeshLib::Mesh const& mesh, GeoLib::Surface const& sfc, bool search_all_nodes = true); + /** + * Constructor of object, that search mesh nodes along a + * GeoLib::Surface object within a given search radius. + * @param mesh_nodes Nodes the search will be performed on. + * @param sfc Along the GeoLib::Surface sfc the mesh nodes are searched. + */ + MeshNodesAlongSurface(MeshLib::Mesh const& mesh, GeoLib::Surface const& sfc, bool search_all_nodes = true); - /// return the mesh object - MeshLib::Mesh const& getMesh() const; + /// return the mesh object + MeshLib::Mesh const& getMesh() const; - /** - * Access the vector of mesh node ids. - * @return The vector of mesh node ids calculated in the constructor - */ - std::vector<std::size_t> const& getNodeIDs () const; + /** + * Access the vector of mesh node ids. + * @return The vector of mesh node ids calculated in the constructor + */ + std::vector<std::size_t> const& getNodeIDs () const; - /** - * Deploying this method the user can get access to the underlying - * GeoLib::Surface. - * @return the underlying GeoLib::Surface - */ - GeoLib::Surface const& getSurface () const; + /** + * Deploying this method the user can get access to the underlying + * GeoLib::Surface. + * @return the underlying GeoLib::Surface + */ + GeoLib::Surface const& getSurface () const; private: - MeshLib::Mesh const& _mesh; - GeoLib::Surface const& _sfc; - std::vector<std::size_t> _msh_node_ids; + MeshLib::Mesh const& _mesh; + GeoLib::Surface const& _sfc; + std::vector<std::size_t> _msh_node_ids; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesOnPoint.cpp b/MeshGeoToolsLib/MeshNodesOnPoint.cpp index 6b9c00daf84b519c0abe22e8b82b1f576c0acf0a..db8d03f3a5de78f5af9b1e8f4758aae49c4a989c 100644 --- a/MeshGeoToolsLib/MeshNodesOnPoint.cpp +++ b/MeshGeoToolsLib/MeshNodesOnPoint.cpp @@ -17,14 +17,14 @@ MeshNodesOnPoint::MeshNodesOnPoint(MeshLib::Mesh const& mesh, GeoLib::Grid<MeshL GeoLib::Point const& pnt, double epsilon_radius, bool search_all_nodes) : _mesh(mesh), _pnt(pnt) { - std::vector<std::size_t> vec_ids(mesh_grid.getPointsInEpsilonEnvironment(pnt, epsilon_radius)); - if (search_all_nodes) - _msh_node_ids = vec_ids; - else { - for (auto id : vec_ids) - if (mesh.isBaseNode(id)) - _msh_node_ids.push_back(id); - } + std::vector<std::size_t> vec_ids(mesh_grid.getPointsInEpsilonEnvironment(pnt, epsilon_radius)); + if (search_all_nodes) + _msh_node_ids = vec_ids; + else { + for (auto id : vec_ids) + if (mesh.isBaseNode(id)) + _msh_node_ids.push_back(id); + } } } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesOnPoint.h b/MeshGeoToolsLib/MeshNodesOnPoint.h index 1b4d36c43c98f3901a603221fed384dc3ac10af4..3c8c6d9015712a27b1dbd24e9bbae108ee4674ae 100644 --- a/MeshGeoToolsLib/MeshNodesOnPoint.h +++ b/MeshGeoToolsLib/MeshNodesOnPoint.h @@ -29,38 +29,38 @@ namespace MeshGeoToolsLib class MeshNodesOnPoint { public: - /** - * Constructor of object, that search mesh nodes at a GeoLib::Point point - * within a given search radius. - * @param mesh Mesh object whose nodes are searched - * @param mesh_grid Grid object constructed with mesh nodes - * @param pnt a point - * @param epsilon_radius Search radius - * @param search_all_nodes whether this shearches all nodes or only base nodes - */ - MeshNodesOnPoint(MeshLib::Mesh const& mesh, GeoLib::Grid<MeshLib::Node> const &mesh_grid, - GeoLib::Point const& pnt, double epsilon_radius, bool search_all_nodes = true); + /** + * Constructor of object, that search mesh nodes at a GeoLib::Point point + * within a given search radius. + * @param mesh Mesh object whose nodes are searched + * @param mesh_grid Grid object constructed with mesh nodes + * @param pnt a point + * @param epsilon_radius Search radius + * @param search_all_nodes whether this shearches all nodes or only base nodes + */ + MeshNodesOnPoint(MeshLib::Mesh const& mesh, GeoLib::Grid<MeshLib::Node> const &mesh_grid, + GeoLib::Point const& pnt, double epsilon_radius, bool search_all_nodes = true); - /// return the mesh object - MeshLib::Mesh const& getMesh() const { return _mesh; } + /// return the mesh object + MeshLib::Mesh const& getMesh() const { return _mesh; } - /** - * Access the vector of mesh node ids. - * @return The vector of mesh node ids calculated in the constructor - */ - std::vector<std::size_t> const& getNodeIDs () const { return _msh_node_ids; } + /** + * Access the vector of mesh node ids. + * @return The vector of mesh node ids calculated in the constructor + */ + std::vector<std::size_t> const& getNodeIDs () const { return _msh_node_ids; } - /** - * Deploying this method the user can get access to the underlying - * GeoLib::Point. - * @return the underlying GeoLib::Point - */ - GeoLib::Point const& getPoint () const { return _pnt; } + /** + * Deploying this method the user can get access to the underlying + * GeoLib::Point. + * @return the underlying GeoLib::Point + */ + GeoLib::Point const& getPoint () const { return _pnt; } private: - MeshLib::Mesh const& _mesh; - GeoLib::Point const& _pnt; - std::vector<std::size_t> _msh_node_ids; + MeshLib::Mesh const& _mesh; + GeoLib::Point const& _pnt; + std::vector<std::size_t> _msh_node_ids; }; } // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/SearchLength.h b/MeshGeoToolsLib/SearchLength.h index 6caf37d6035a0e42151eb51814b6d8215cf39804..019bddd5dc20dc3b74dafdb4e3e78434a4453c01 100644 --- a/MeshGeoToolsLib/SearchLength.h +++ b/MeshGeoToolsLib/SearchLength.h @@ -22,23 +22,23 @@ namespace MeshGeoToolsLib class SearchLength { public: - /// Constructor for SearchLength object with a default search length - /// of 10 angstrom (\f$10^{-9}\f$ m) - explicit SearchLength(double search_length = 1e-9) - : _search_length(search_length) {} + /// Constructor for SearchLength object with a default search length + /// of 10 angstrom (\f$10^{-9}\f$ m) + explicit SearchLength(double search_length = 1e-9) + : _search_length(search_length) {} - SearchLength(SearchLength const&) = default; - SearchLength& operator=(SearchLength const&) = default; + SearchLength(SearchLength const&) = default; + SearchLength& operator=(SearchLength const&) = default; - virtual ~SearchLength() = default; + virtual ~SearchLength() = default; - virtual double getSearchLength() const - { - return _search_length; - } + virtual double getSearchLength() const + { + return _search_length; + } protected: - double _search_length; + double _search_length; }; } // end namespace MeshGeoToolsLib