diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp index 23cad184be33104d56fcd79381246ccdf36fa329..6065ad11ef871b6d838ff0513680c9dd0553ecb4 100644 --- a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp +++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp @@ -10,18 +10,15 @@ #include "AppendLinesAlongPolyline.h" #include "BaseLib/Logging.h" - #include "GeoLib/Polyline.h" #include "GeoLib/PolylineVec.h" - -#include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" -#include "MeshLib/Elements/Line.h" +#include "MeshGeoToolsLib/MeshNodesAlongPolyline.h" #include "MeshLib/Elements/Element.h" -#include "MeshLib/MeshEnums.h" +#include "MeshLib/Elements/Line.h" +#include "MeshLib/Mesh.h" #include "MeshLib/MeshEditing/DuplicateMeshComponents.h" - -#include "MeshGeoToolsLib/MeshNodesAlongPolyline.h" +#include "MeshLib/MeshEnums.h" +#include "MeshLib/Node.h" namespace MeshGeoToolsLib { @@ -29,8 +26,10 @@ 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); + 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); auto const material_ids = materialIDs(mesh); int const max_matID = @@ -39,7 +38,7 @@ std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines( : 0; std::vector<int> new_mat_ids; - const std::size_t n_ply (ply_vec.size()); + const std::size_t n_ply(ply_vec.size()); // for each polyline for (std::size_t k(0); k < n_ply; k++) { @@ -49,8 +48,9 @@ std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines( MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly( mesh, *ply, mesh.getMinEdgeLength() * 0.5, MeshGeoToolsLib::SearchAllNodes::Yes); - auto &vec_nodes_on_ply = mshNodesAlongPoly.getNodeIDs(); - if (vec_nodes_on_ply.empty()) { + 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); @@ -58,13 +58,14 @@ std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines( } // add line elements - for (std::size_t i=0; i<vec_nodes_on_ply.size()-1; i++) { + 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]]; + 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); + new_mat_ids.push_back(max_matID + k + 1); } } diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp index 75bab1c659ea5018601cc3d747790cfe555b3aff..9af0f991e76f96f24f9b6719070d5b99aa66fbd5 100644 --- a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp +++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp @@ -11,13 +11,12 @@ #include "GeoLib/Point.h" #include "MathLib/Point3d.h" +#include "MeshGeoToolsLib/MeshNodeSearcher.h" #include "MeshLib/Elements/Point.h" #include "MeshLib/Mesh.h" #include "MeshLib/MeshSearch/ElementSearch.h" #include "MeshLib/Node.h" -#include "MeshGeoToolsLib/MeshNodeSearcher.h" - namespace MeshGeoToolsLib { BoundaryElementsAtPoint::BoundaryElementsAtPoint( @@ -77,8 +76,7 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint( std::array<MeshLib::Node*, 1> const nodes = { {const_cast<MeshLib::Node*>(_mesh.getNode(nearest_node_id))}}; - _boundary_elements.push_back( - new MeshLib::Point{nodes, nearest_node_id}); + _boundary_elements.push_back(new MeshLib::Point{nodes, nearest_node_id}); } BoundaryElementsAtPoint::~BoundaryElementsAtPoint() diff --git a/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp b/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp index 5c9ca623821122bb5c94cd529950d978f674c9c1..49fa443a52390a3defeb4f647d1d488b695e1e5a 100644 --- a/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp +++ b/MeshGeoToolsLib/BoundaryElementsOnSurface.cpp @@ -10,13 +10,11 @@ #include "BoundaryElementsOnSurface.h" #include "GeoLib/Surface.h" - -#include "MeshLib/Mesh.h" +#include "MeshGeoToolsLib/MeshNodeSearcher.h" #include "MeshLib/Elements/Element.h" +#include "MeshLib/Mesh.h" #include "MeshLib/MeshSearch/ElementSearch.h" -#include "MeshGeoToolsLib/MeshNodeSearcher.h" - namespace MeshGeoToolsLib { BoundaryElementsOnSurface::BoundaryElementsOnSurface( @@ -28,10 +26,11 @@ BoundaryElementsOnSurface::BoundaryElementsOnSurface( auto node_ids_on_sfc = mshNodeSearcher.getMeshNodeIDs(sfc); MeshLib::ElementSearch es(_mesh); es.searchByNodeIDs(node_ids_on_sfc); - auto &ele_ids_near_sfc = es.getSearchedElementIDs(); + auto& ele_ids_near_sfc = es.getSearchedElementIDs(); // get a list of faces made of the nodes - for (auto ele_id : ele_ids_near_sfc) { + for (auto ele_id : ele_ids_near_sfc) + { auto* e = _mesh.getElement(ele_id); // skip internal elements if (!e->isBoundaryElement()) @@ -39,11 +38,13 @@ BoundaryElementsOnSurface::BoundaryElementsOnSurface( continue; } // find faces on surface - for (unsigned i=0; i<e->getNumberOfFaces(); i++) { + for (unsigned i = 0; i < e->getNumberOfFaces(); i++) + { auto* face = e->getFace(i); // check std::size_t cnt_match = 0; - for (std::size_t j=0; j<face->getNumberOfBaseNodes(); j++) { + for (std::size_t j = 0; j < face->getNumberOfBaseNodes(); j++) + { if (std::find(node_ids_on_sfc.begin(), node_ids_on_sfc.end(), face->getNodeIndex(j)) != node_ids_on_sfc.end()) { @@ -55,9 +56,10 @@ BoundaryElementsOnSurface::BoundaryElementsOnSurface( } } // update the list - if (cnt_match==face->getNumberOfBaseNodes()) + if (cnt_match == face->getNumberOfBaseNodes()) { - _boundary_elements.push_back(const_cast<MeshLib::Element*>(face)); + _boundary_elements.push_back( + const_cast<MeshLib::Element*>(face)); } else { @@ -75,5 +77,4 @@ BoundaryElementsOnSurface::~BoundaryElementsOnSurface() } } -} // end namespace MeshGeoToolsLib - +} // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp index 3f966b57e785e936d02a8178e0b75ec50f6a12bd..d54c54057ac922d15721ac3485c443995f04d051 100644 --- a/MeshGeoToolsLib/BoundaryElementsSearcher.cpp +++ b/MeshGeoToolsLib/BoundaryElementsSearcher.cpp @@ -12,24 +12,22 @@ #include "GeoLib/GeoObject.h" #include "GeoLib/Polyline.h" #include "GeoLib/Surface.h" - -#include "MeshLib/Mesh.h" -#include "MeshLib/Node.h" -#include "MeshLib/Elements/Element.h" -#include "MeshLib/Elements/Point.h" - -#include "MeshGeoToolsLib/MeshNodeSearcher.h" -#include "MeshGeoToolsLib/BoundaryElementsAtPoint.h" #include "MeshGeoToolsLib/BoundaryElementsAlongPolyline.h" +#include "MeshGeoToolsLib/BoundaryElementsAtPoint.h" #include "MeshGeoToolsLib/BoundaryElementsOnSurface.h" - +#include "MeshGeoToolsLib/MeshNodeSearcher.h" +#include "MeshLib/Elements/Element.h" +#include "MeshLib/Elements/Point.h" +#include "MeshLib/Mesh.h" +#include "MeshLib/Node.h" namespace MeshGeoToolsLib { BoundaryElementsSearcher::BoundaryElementsSearcher( MeshLib::Mesh const& mesh, MeshNodeSearcher const& mshNodeSearcher) : _mesh(mesh), _mshNodeSearcher(mshNodeSearcher) -{} +{ +} BoundaryElementsSearcher::~BoundaryElementsSearcher() { @@ -120,5 +118,4 @@ BoundaryElementsSearcher::getBoundaryElements(GeoLib::GeoObject const& geoObj, } } -} // end namespace MeshGeoToolsLib - +} // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp b/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp index 3a9b74e47fe349fca96d6efeffb3f1ef9ea4ffcb..7989e3f3380fdcca49fe7ac28489e3bafe1f5aa8 100644 --- a/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp +++ b/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp @@ -10,12 +10,10 @@ #include "ConstructMeshesFromGeometries.h" #include "BaseLib/Logging.h" - +#include "BoundaryElementsSearcher.h" #include "MeshLib/Elements/Element.h" #include "MeshLib/MeshEditing/DuplicateMeshComponents.h" #include "MeshLib/Node.h" - -#include "BoundaryElementsSearcher.h" #include "MeshNodeSearcher.h" namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/CreateSearchLength.cpp b/MeshGeoToolsLib/CreateSearchLength.cpp index 6115a1bf82c9bf76cb0817b4a7de3fb93efbb47b..ebeffa8286ad2bf6977049aaf83e85e7d2cbdd99 100644 --- a/MeshGeoToolsLib/CreateSearchLength.cpp +++ b/MeshGeoToolsLib/CreateSearchLength.cpp @@ -11,7 +11,6 @@ #include "BaseLib/ConfigTree.h" #include "BaseLib/Error.h" - #include "MeshGeoToolsLib/HeuristicSearchLength.h" #include "MeshGeoToolsLib/SearchLength.h" diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp index 9ec3af6fc341f9f629c6260fde90f159ab100450..1cef01614cfa2b9f175b7a96e1bda88f4e65ce92 100644 --- a/MeshGeoToolsLib/GeoMapper.cpp +++ b/MeshGeoToolsLib/GeoMapper.cpp @@ -15,31 +15,32 @@ #include "GeoMapper.h" #include <algorithm> -#include <sstream> #include <numeric> - -#include "BaseLib/Logging.h" +#include <sstream> #include "BaseLib/Algorithm.h" - +#include "BaseLib/Logging.h" #include "GeoLib/AABB.h" #include "GeoLib/AnalyticalGeometry.h" #include "GeoLib/GEOObjects.h" #include "GeoLib/Raster.h" #include "GeoLib/StationBorehole.h" - -#include "MeshLib/Mesh.h" #include "MeshLib/Elements/Element.h" #include "MeshLib/Elements/FaceRule.h" -#include "MeshLib/Node.h" -#include "MeshLib/MeshSurfaceExtraction.h" +#include "MeshLib/Mesh.h" #include "MeshLib/MeshSearch/MeshElementGrid.h" +#include "MeshLib/MeshSurfaceExtraction.h" +#include "MeshLib/Node.h" -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) +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) { } @@ -50,24 +51,31 @@ GeoMapper::~GeoMapper() void GeoMapper::mapOnDEM(std::unique_ptr<GeoLib::Raster const> raster) { - std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); - if (! pnts) { + std::vector<GeoLib::Point*> const* pnts( + _geo_objects.getPointVec(_geo_name)); + if (!pnts) + { ERR("Geometry '{:s}' does not exist.", _geo_name); return; } _raster = std::move(raster); - if (GeoLib::isStation((*pnts)[0])) { + if (GeoLib::isStation((*pnts)[0])) + { mapStationData(*pnts); - } else { + } + else + { mapPointDataToDEM(*pnts); } } -void GeoMapper::mapOnMesh(MeshLib::Mesh const*const mesh) +void GeoMapper::mapOnMesh(MeshLib::Mesh const* const mesh) { - std::vector<GeoLib::Point*> const* pnts(_geo_objects.getPointVec(_geo_name)); - if (! pnts) { + std::vector<GeoLib::Point*> const* pnts( + _geo_objects.getPointVec(_geo_name)); + if (!pnts) + { ERR("Geometry '{:s}' does not exist.", _geo_name); return; } @@ -82,26 +90,31 @@ void GeoMapper::mapOnMesh(MeshLib::Mesh const*const mesh) } else { - Eigen::Vector3d const dir(0,0,-1); + Eigen::Vector3d const 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::Point3d origin(std::array<double, 3>{{0, 0, 0}}); std::vector<MeshLib::Node> flat_nodes; flat_nodes.reserve(_surface_mesh->getNumberOfNodes()); // 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()) { + 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()); + _grid = + new GeoLib::Grid<MeshLib::Node>(flat_nodes.cbegin(), flat_nodes.cend()); - if (GeoLib::isStation((*pnts)[0])) { + if (GeoLib::isStation((*pnts)[0])) + { mapStationData(*pnts); - } else { + } + else + { mapPointDataToMeshSurface(*pnts); } @@ -110,13 +123,15 @@ void GeoMapper::mapOnMesh(MeshLib::Mesh const*const mesh) void GeoMapper::mapToConstantValue(double value) { - std::vector<GeoLib::Point*> const* points (this->_geo_objects.getPointVec(this->_geo_name)); + std::vector<GeoLib::Point*> const* points( + this->_geo_objects.getPointVec(this->_geo_name)); if (points == nullptr) { ERR("Geometry '{:s}' not found.", this->_geo_name); return; } - std::for_each(points->begin(), points->end(), [value](GeoLib::Point* pnt){ (*pnt)[2] = value; }); + std::for_each(points->begin(), points->end(), + [value](GeoLib::Point* pnt) { (*pnt)[2] = value; }); } void GeoMapper::mapStationData(std::vector<GeoLib::Point*> const& points) @@ -125,13 +140,13 @@ void GeoMapper::mapStationData(std::vector<GeoLib::Point*> const& points) double max_val(0); if (_surface_mesh) { - GeoLib::AABB bounding_box( - _surface_mesh->getNodes().begin(), _surface_mesh->getNodes().end()); + 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) + for (auto* pnt : points) { double offset = (_grid) @@ -143,8 +158,10 @@ void GeoMapper::mapStationData(std::vector<GeoLib::Point*> const& points) { continue; } - auto const& layers = static_cast<GeoLib::StationBorehole*>(pnt)->getProfile(); - for (auto * layer_pnt : layers) { + auto const& layers = + static_cast<GeoLib::StationBorehole*>(pnt)->getProfile(); + for (auto* layer_pnt : layers) + { (*layer_pnt)[2] = (*layer_pnt)[2] + offset; } } @@ -153,24 +170,26 @@ void GeoMapper::mapStationData(std::vector<GeoLib::Point*> const& points) void GeoMapper::mapPointDataToDEM( std::vector<GeoLib::Point*> const& points) const { - for (auto * pnt : points) + for (auto* pnt : points) { - GeoLib::Point &p(*pnt); + GeoLib::Point& p(*pnt); p[2] = getDemElevation(p); } } -void GeoMapper::mapPointDataToMeshSurface(std::vector<GeoLib::Point*> const& pnts) +void GeoMapper::mapPointDataToMeshSurface( + std::vector<GeoLib::Point*> const& pnts) { - GeoLib::AABB const aabb( - _surface_mesh->getNodes().cbegin(), _surface_mesh->getNodes().cend()); + 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) { + 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); + GeoLib::Point& p(*pnt); if (p[0] < aabb.getMinPoint()[0] || aabb.getMaxPoint()[0] < p[0]) { continue; @@ -186,7 +205,7 @@ void GeoMapper::mapPointDataToMeshSurface(std::vector<GeoLib::Point*> const& pnt float GeoMapper::getDemElevation(GeoLib::Point const& pnt) const { - double const elevation (_raster->getValueAtPoint(pnt)); + double const elevation(_raster->getValueAtPoint(pnt)); if (std::abs(elevation - _raster->getHeader().no_data) < std::numeric_limits<double>::epsilon()) { @@ -195,8 +214,8 @@ float GeoMapper::getDemElevation(GeoLib::Point const& pnt) const return static_cast<float>(elevation); } -double GeoMapper::getMeshElevation( - double x, double y, double min_val, double max_val) const +double GeoMapper::getMeshElevation(double x, double y, double min_val, + double max_val) const { const MeshLib::Node* pnt = _grid->getNearestPoint(MathLib::Point3d{{{x, y, 0}}}); @@ -204,7 +223,7 @@ double GeoMapper::getMeshElevation( _surface_mesh->getNode(pnt->getID())->getElements()); std::unique_ptr<GeoLib::Point> intersection; - for (auto const & element : elements) + for (auto const& element : elements) { if (intersection == nullptr && element->getGeomType() != MeshLib::MeshElemType::LINE) @@ -228,11 +247,13 @@ double GeoMapper::getMeshElevation( { return (*intersection)[2]; } - // if something goes wrong, simply take the elevation of the nearest mesh node + // if something goes wrong, simply take the elevation of the nearest mesh + // node return (*(_surface_mesh->getNode(pnt->getID())))[2]; } -/// Find the 2d-element within the \c elements that contains the given point \c p. +/// Find the 2d-element within the \c elements that contains the given point \c +/// p. /// /// The algorithm projects every element of the elements vector and the point /// \c p orthogonal to the \f$x\f$-\f$y\f$ plane. In the \f$x\f$-\f$y\f$ plane @@ -241,19 +262,24 @@ static MeshLib::Element const* findElementContainingPointXY( std::vector<MeshLib::Element const*> const& elements, MathLib::Point3d const& p) { - for (auto const elem : elements) { + for (auto const elem : elements) + { std::unique_ptr<MeshLib::Element> elem_2d(elem->clone()); // reset/copy the nodes - for (std::size_t k(0); k<elem_2d->getNumberOfNodes(); ++k) { + for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k) + { elem_2d->setNode(k, new MeshLib::Node(*elem_2d->getNode(k))); } // project to xy - for (std::size_t k(0); k<elem_2d->getNumberOfNodes(); ++k) { + for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k) + { (*const_cast<MeshLib::Node*>(elem_2d->getNode(k)))[2] = 0.0; } - if (elem_2d->isPntInElement(MathLib::Point3d{ {{p[0], p[1], 0.0}} })) { + if (elem_2d->isPntInElement(MathLib::Point3d{{{p[0], p[1], 0.0}}})) + { // clean up the copied nodes - for (std::size_t k(0); k<elem_2d->getNumberOfNodes(); ++k) { + for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k) + { delete elem_2d->getNode(k); } return elem; @@ -277,9 +303,11 @@ static std::vector<MathLib::Point3d> computeElementSegmentIntersections( std::unique_ptr<MeshLib::Element const>(elem.getEdge(k)); GeoLib::LineSegment elem_segment{ new GeoLib::Point(*dynamic_cast<MathLib::Point3d*>( - const_cast<MeshLib::Node*>(edge->getNode(0))), 0), + const_cast<MeshLib::Node*>(edge->getNode(0))), + 0), new GeoLib::Point(*dynamic_cast<MathLib::Point3d*>( - const_cast<MeshLib::Node*>(edge->getNode(1))), 0), + const_cast<MeshLib::Node*>(edge->getNode(1))), + 0), true}; std::vector<MathLib::Point3d> const intersections( GeoLib::lineSegmentIntersect2d(segment, elem_segment)); @@ -296,10 +324,11 @@ static std::vector<GeoLib::LineSegment> createSubSegmentsForElement( MathLib::Point3d const& end_pnt, MeshLib::Element const* const elem) { std::vector<GeoLib::LineSegment> sub_segments; - if (intersections.size() > 2) { + if (intersections.size() > 2) + { std::stringstream out; out << "element with id " << elem->getID() << " and seg " - << " intersecting at more than two edges\n"; + << " intersecting at more than two edges\n"; for (std::size_t k(0); k < intersections.size(); ++k) { out << k << " " << intersections[k] << "\n"; @@ -361,7 +390,8 @@ static std::vector<GeoLib::LineSegment> mapLineSegment( MathLib::Point3d const& beg_pnt(segment.getBeginPoint()); MathLib::Point3d const& end_pnt(segment.getEndPoint()); - for (auto const elem : surface_elements) { + for (auto const elem : surface_elements) + { // compute element-segment-intersections (2d in x-y-plane) std::vector<MathLib::Point3d> element_intersections( computeElementSegmentIntersections(*elem, segment)); @@ -387,8 +417,7 @@ static std::vector<GeoLib::LineSegment> mapLineSegment( auto min_dist_segment = std::min_element( sub_segments.begin(), sub_segments.end(), [&beg_pnt](GeoLib::LineSegment const& seg0, - GeoLib::LineSegment const& seg1) - { + GeoLib::LineSegment const& seg1) { // min dist for segment 0 const double d0( std::min(MathLib::sqrDist(beg_pnt, seg0.getBeginPoint()), @@ -399,7 +428,7 @@ static std::vector<GeoLib::LineSegment> mapLineSegment( MathLib::sqrDist(beg_pnt, seg1.getEndPoint()))); return d0 < d1; }); - GeoLib::Point * pnt{ + GeoLib::Point* pnt{ MathLib::sqrDist(beg_pnt, min_dist_segment->getBeginPoint()) < MathLib::sqrDist(beg_pnt, min_dist_segment->getEndPoint()) ? new GeoLib::Point{min_dist_segment->getBeginPoint()} @@ -422,11 +451,14 @@ static void mapPointOnSurfaceElement(MeshLib::Element const& elem, auto const p = Eigen::Map<Eigen::Vector3d const>(elem.getNode(0)->getCoords()); Eigen::Vector3d const n(MeshLib::FaceRule::getSurfaceNormal(&elem)); - if (n[2] == 0.0) { // vertical plane, z coordinate is arbitrary + if (n[2] == 0.0) + { // vertical plane, z coordinate is arbitrary q[2] = p[2]; - } else { + } + else + { double const d(n.dot(p)); - q[2] = (d - n[0]*q[0] - n[1]*q[1])/n[2]; + q[2] = (d - n[0] * q[0] - n[1] * q[1]) / n[2]; } } @@ -461,11 +493,13 @@ static bool snapPointToElementNode(MathLib::Point3d& p, // values will be initialized within computeSqrNodeDistanceRange auto const [sqr_min, sqr_max] = MeshLib::computeSqrNodeDistanceRange(elem); - double const sqr_eps(rel_eps*rel_eps * sqr_min); - for (std::size_t k(0); k<elem.getNumberOfNodes(); ++k) { + double const sqr_eps(rel_eps * rel_eps * sqr_min); + for (std::size_t k(0); k < elem.getNumberOfNodes(); ++k) + { auto const& node(*elem.getNode(k)); double const sqr_dist_2d(MathLib::sqrDist2d(p, node)); - if (sqr_dist_2d < sqr_eps) { + if (sqr_dist_2d < sqr_eps) + { #ifdef DEBUG_GEOMAPPER std::stringstream out; out.precision(std::numeric_limits<double>::digits10); @@ -540,7 +574,8 @@ static void mapPolylineOnSurfaceMesh( // invalid and hence it is necessary to re-create them. orig_points.resetInternalDataStructures(); - if (beg_elem == end_elem) { + if (beg_elem == end_elem) + { // TODO: handle cases: beg_elem == end_elem == nullptr // There are further checks necessary to determine which case we are // in: @@ -575,12 +610,15 @@ void GeoMapper::advancedMapOnMesh(MeshLib::Mesh const& mesh) // 1. extract surface delete _surface_mesh; - if (mesh.getDimension()<3) { + if (mesh.getDimension() < 3) + { _surface_mesh = new MeshLib::Mesh(mesh); - } else { + } + else + { Eigen::Vector3d const dir({0, 0, -1}); - _surface_mesh = - MeshLib::MeshSurfaceExtraction::getMeshSurface(mesh, dir, 90+1e-6); + _surface_mesh = MeshLib::MeshSurfaceExtraction::getMeshSurface( + mesh, dir, 90 + 1e-6); } // 2. compute mesh grid for surface @@ -595,4 +633,4 @@ void GeoMapper::advancedMapOnMesh(MeshLib::Mesh const& mesh) } } -} // end namespace MeshGeoToolsLib +} // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/HeuristicSearchLength.cpp b/MeshGeoToolsLib/HeuristicSearchLength.cpp index 318ae2c34523367824df1762b635e8c190ff8a62..096d7aa3c273fae8345c1d4b38eb99d52e7133a7 100644 --- a/MeshGeoToolsLib/HeuristicSearchLength.cpp +++ b/MeshGeoToolsLib/HeuristicSearchLength.cpp @@ -13,36 +13,39 @@ #include "HeuristicSearchLength.h" #include "BaseLib/Logging.h" - #include "MeshLib/Elements/Element.h" namespace MeshGeoToolsLib { - -HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthType length_type) -: _mesh(mesh) +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::size_t n_sampling(0); // total length of edges squared std::vector<MeshLib::Element*> const& elements(_mesh.getElements()); - if (length_type==LengthType::Edge) { + if (length_type == LengthType::Edge) + { for (auto element : elements) { std::size_t const n_edges(element->getNumberOfEdges()); - for (std::size_t k(0); k<n_edges; k++) { + for (std::size_t k(0); k < n_edges; k++) + { auto edge = element->getEdge(k); // allocation inside getEdge(). double const len = edge->getContent(); delete edge; sum += len; - sum_of_sqr += len*len; + sum_of_sqr += len * len; } n_sampling += n_edges; } - } else { + } + else + { for (const MeshLib::Element* e : elements) { auto const [min, max] = @@ -53,14 +56,16 @@ HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthTy n_sampling = _mesh.getNumberOfElements(); } - 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; + _search_length = mean / 2; - if (variance > 0) { + if (variance > 0) + { if (variance < mean * mean / 4) { _search_length -= std::sqrt(variance); @@ -77,4 +82,4 @@ HeuristicSearchLength::HeuristicSearchLength(MeshLib::Mesh const& mesh, LengthTy _mesh.getName(), _search_length); } -} // end namespace MeshGeoToolsLib +} // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp index 5246b2382a4235058122fd2dfb89f5e4fd2ca156..bf9ca67a1ea34984278f02c099ad87fcfe1b5ed3 100644 --- a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp +++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp @@ -13,7 +13,6 @@ #include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" - #include "MeshNodeSearcher.h" namespace diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp index cd3c90f89e2e578bf442c3c4b6a066e18c13ea33..7af112eb6b31b40c5927849d7483002534d572ef 100644 --- a/MeshGeoToolsLib/MeshNodeSearcher.cpp +++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp @@ -11,22 +11,19 @@ #include "MeshNodeSearcher.h" -#include <typeinfo> #include <sstream> - -#include "HeuristicSearchLength.h" -#include "MeshNodesAlongPolyline.h" -#include "MeshNodesAlongSurface.h" -#include "MeshNodesOnPoint.h" +#include <typeinfo> #include "BaseLib/Logging.h" - #include "GeoLib/Point.h" #include "GeoLib/Polyline.h" - +#include "HeuristicSearchLength.h" #include "MeshLib/Elements/Element.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" +#include "MeshNodesAlongPolyline.h" +#include "MeshNodesAlongSurface.h" +#include "MeshNodesOnPoint.h" namespace MeshGeoToolsLib { @@ -67,8 +64,7 @@ std::vector<std::size_t> const& getMeshNodeIDs( std::vector<CacheType*>& cached_elements, std::function<GeometryType(CacheType const&)> getCachedItem, GeometryType const& item, MeshLib::Mesh const& mesh, - GeoLib::Grid<MeshLib::Node> const& mesh_grid, - double const search_length, + GeoLib::Grid<MeshLib::Node> const& mesh_grid, double const search_length, SearchAllNodes const search_all_nodes) { if (auto const it = find_if(cbegin(cached_elements), cend(cached_elements), diff --git a/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp b/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp index 2c2765d394b94120318a12e1d71808d790bc7631..bc5e16aca4e025c7851f4ecefa3a3b21fd59a8f4 100644 --- a/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp +++ b/MeshGeoToolsLib/MeshNodesAlongPolyline.cpp @@ -14,8 +14,8 @@ #include <algorithm> #include "BaseLib/quicksort.h" -#include "MathLib/MathTools.h" #include "GeoLib/Polyline.h" +#include "MathLib/MathTools.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" @@ -31,32 +31,36 @@ MeshNodesAlongPolyline::MeshNodesAlongPolyline(MeshLib::Mesh const& mesh, const std::size_t n_nodes(search_all_nodes == SearchAllNodes::Yes ? _mesh.getNumberOfNodes() : _mesh.getNumberOfBaseNodes()); - auto &mesh_nodes = _mesh.getNodes(); + 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) { + 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); + 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 +MeshLib::Mesh const& MeshNodesAlongPolyline::getMesh() const { return _mesh; } -std::vector<std::size_t> const& MeshNodesAlongPolyline::getNodeIDs () const +std::vector<std::size_t> const& MeshNodesAlongPolyline::getNodeIDs() const { return _msh_node_ids; } -GeoLib::Polyline const& MeshNodesAlongPolyline::getPolyline () const +GeoLib::Polyline const& MeshNodesAlongPolyline::getPolyline() const { return _ply; } diff --git a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp index bafc4c2548641727d6d589ee86ddeea90bcd947d..2c39e13650cd09ad7cd1071a02f194372275606e 100644 --- a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp +++ b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp @@ -15,8 +15,8 @@ #include <algorithm> #include "BaseLib/quicksort.h" -#include "MathLib/MathTools.h" #include "GeoLib/Surface.h" +#include "MathLib/MathTools.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" @@ -33,31 +33,33 @@ MeshNodesAlongSurface::MeshNodesAlongSurface(MeshLib::Mesh const& mesh, ? _mesh.getNumberOfNodes() : _mesh.getNumberOfBaseNodes()); // loop over all nodes - for (std::size_t i = 0; i < n_nodes; i++) { + for (std::size_t i = 0; i < n_nodes; i++) + { auto* node = mesh_nodes[i]; if (!sfc.isPntInBoundingVolume(*node, epsilon_radius)) { continue; } - if (sfc.isPntInSfc(*node, epsilon_radius)) { + if (sfc.isPntInSfc(*node, epsilon_radius)) + { _msh_node_ids.push_back(node->getID()); } } } -MeshLib::Mesh const& MeshNodesAlongSurface::getMesh () const +MeshLib::Mesh const& MeshNodesAlongSurface::getMesh() const { return _mesh; } -std::vector<std::size_t> const& MeshNodesAlongSurface::getNodeIDs () const +std::vector<std::size_t> const& MeshNodesAlongSurface::getNodeIDs() const { return _msh_node_ids; } -GeoLib::Surface const& MeshNodesAlongSurface::getSurface () const +GeoLib::Surface const& MeshNodesAlongSurface::getSurface() const { return _sfc; } -} // end namespace MeshGeoToolsLib +} // end namespace MeshGeoToolsLib diff --git a/MeshGeoToolsLib/MeshNodesOnPoint.cpp b/MeshGeoToolsLib/MeshNodesOnPoint.cpp index d6b8d9fb1cdbb7bd4be04966d955310c989c50c5..30d0e93113fd8cd6606672f7e796d9fde8045f2a 100644 --- a/MeshGeoToolsLib/MeshNodesOnPoint.cpp +++ b/MeshGeoToolsLib/MeshNodesOnPoint.cpp @@ -39,4 +39,3 @@ MeshNodesOnPoint::MeshNodesOnPoint(MeshLib::Mesh const& mesh, } } // end namespace MeshGeoToolsLib -