38 const std::string& geo_name)
39 : _geo_objects(geo_objects),
40 _geo_name(const_cast<std::string&>(geo_name)),
41 _surface_mesh(nullptr),
54 std::vector<GeoLib::Point*>
const* pnts(
75 std::vector<GeoLib::Point*>
const* pnts(
93 Eigen::Vector3d
const dir(0, 0, -1);
100 std::vector<MeshLib::Node> flat_nodes;
106 flat_nodes.emplace_back(*n_ptr);
107 flat_nodes.back()[2] = 0.0;
126 std::vector<GeoLib::Point*>
const* points(
128 if (points ==
nullptr)
133 std::for_each(points->begin(), points->end(),
149 for (
auto* pnt : points)
164 for (
auto* layer_pnt : layers)
166 (*layer_pnt)[2] = (*layer_pnt)[2] + offset;
172 std::vector<GeoLib::Point*>
const& points)
const
174 for (
auto* pnt : points)
182 std::vector<GeoLib::Point*>
const& pnts)
189 for (
auto* pnt : pnts)
209 double const elevation(
_raster->getValueAtPoint(pnt));
210 if (std::abs(elevation -
_raster->getHeader().no_data) <
211 std::numeric_limits<double>::epsilon())
215 return static_cast<float>(elevation);
219 double max_val)
const
227 for (
auto const& element : elements)
233 *element->getNode(0), *element->getNode(1),
242 *element->getNode(0), *element->getNode(2),
263 std::vector<MeshLib::Element const*>
const& elements,
266 for (
auto const elem : elements)
268 std::unique_ptr<MeshLib::Element> elem_2d(elem->clone());
270 for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k)
272 elem_2d->setNode(k,
new MeshLib::Node(*elem_2d->getNode(k)));
275 for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k)
282 for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k)
284 delete elem_2d->getNode(k);
289 for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k)
291 delete elem_2d->getNode(k);
300 std::vector<MathLib::Point3d> element_intersections;
304 std::unique_ptr<MeshLib::Element const>(elem.
getEdge(k));
313 std::vector<MathLib::Point3d>
const intersections(
315 element_intersections.insert(end(element_intersections),
316 begin(intersections), end(intersections));
318 return element_intersections;
322 std::vector<MathLib::Point3d>
const& intersections,
327 std::vector<GeoLib::LineSegment> sub_segments;
328 if (intersections.size() > 2)
330 std::stringstream out;
331 out <<
"element with id " << elem->
getID() <<
" and seg "
332 <<
" intersecting at more than two edges\n";
333 for (std::size_t k(0); k < intersections.size(); ++k)
335 out << k <<
" " << intersections[k] <<
"\n";
337 out <<
"Could not map segment on element. Aborting.\n";
341 if (intersections.size() == 1 && elem == beg_elem)
347 std::numeric_limits<double>::epsilon())
355 if (intersections.size() == 1 && elem == end_elem)
361 std::numeric_limits<double>::epsilon())
363 sub_segments.emplace_back(
new GeoLib::Point{intersections[0], 0},
368 if (intersections.size() == 1 && (elem != beg_elem && elem != end_elem))
376 if (intersections.size() == 2)
378 sub_segments.emplace_back(
new GeoLib::Point{intersections[0], 0},
386 std::vector<MeshLib::Element const*>
const& surface_elements,
390 std::vector<GeoLib::LineSegment> sub_segments;
394 for (
auto const elem : surface_elements)
397 std::vector<MathLib::Point3d> element_intersections(
399 if (element_intersections.empty())
406 std::vector<GeoLib::LineSegment> sub_seg_elem(
408 end_elem, beg_pnt, end_pnt, elem));
409 sub_segments.insert(sub_segments.end(), sub_seg_elem.begin(),
416 if (beg_elem ==
nullptr)
418 auto min_dist_segment = std::min_element(
419 sub_segments.begin(), sub_segments.end(),
425 std::min(MathLib::sqrDist(beg_pnt, seg0.getBeginPoint()),
426 MathLib::sqrDist(beg_pnt, seg0.getEndPoint())));
429 std::min(MathLib::sqrDist(beg_pnt, seg1.getBeginPoint()),
430 MathLib::sqrDist(beg_pnt, seg1.getEndPoint())));
438 sub_segments.emplace_back(
new GeoLib::Point{beg_pnt, 0}, pnt,
true);
443 sub_segments.erase(std::unique(sub_segments.begin(), sub_segments.end()),
462 double const d(n.dot(
p));
463 q[2] = (d - n[0] *
q[0] - n[1] *
q[1]) / n[2];
467 static std::vector<MeshLib::Element const*>
479 std::array<MathLib::Point3d, 2>
const pnts{
484 auto convert_to_Point3d = [](Eigen::Vector3d
const& v) {
488 auto const min = convert_to_Point3d(aabb.
getMinPoint());
489 auto const max = convert_to_Point3d(aabb.
getMaxPoint());
495 return candidate_elements;
504 double const sqr_eps(rel_eps * rel_eps * sqr_min);
507 auto const& node(*elem.
getNode(k));
509 if (sqr_dist_2d < sqr_eps)
511 #ifdef DEBUG_GEOMAPPER
512 std::stringstream out;
513 out.precision(std::numeric_limits<double>::digits10);
514 out <<
"Segment point snapped from " <<
p;
517 #ifdef DEBUG_GEOMAPPER
519 DBUG(
"{:s}", out.str());
530 std::vector<GeoLib::LineSegment>
const& sub_segments)
533 std::size_t new_pnts_cnt(0);
534 for (
auto const& segment : sub_segments)
538 if (ply.
insertPoint(j + new_pnts_cnt + 1, begin_id))
549 segment_it += new_pnts_cnt;
558 for (
auto segment_it(ply.
begin()); segment_it != ply.
end(); ++segment_it)
561 mesh_element_grid, *segment_it));
578 auto const* beg_elem(mapPoint((*segment_it).getBeginPoint()));
579 auto const* end_elem(mapPoint((*segment_it).getEndPoint()));
586 if (beg_elem == end_elem)
601 *segment_it, candidate_elements, beg_elem, end_elem));
603 if (sub_segments.empty())
610 if (sub_segments.size() > 1)
620 delete _surface_mesh;
628 Eigen::Vector3d
const dir({0, 0, -1});
630 mesh, dir, 90 + 1e-6);
637 auto org_lines(_geo_objects.getPolylineVec(_geo_name));
638 auto org_points(_geo_objects.getPointVecObj(_geo_name));
639 for (
auto org_line : *org_lines)
Definition of the AABB class.
Definition of analytical geometry functions.
Definition of the Element class.
Definition of the GEOObjects class.
Definition of the GeoMapper class.
void ERR(char const *fmt, Args const &... args)
void DBUG(char const *fmt, Args const &... args)
Definition of the Mesh class.
Definition of the Node class.
Definition of the GeoLib::Raster class.
Definition of the StationBorehole class.
Class AABB is an axis aligned bounding box around a given set of geometric points of (template) type ...
Eigen::Vector3d const & getMinPoint() const
Eigen::Vector3d const & getMaxPoint() const
Container class for geometric objects.
const std::vector< Point * > * getPointVec(const std::string &name) const
POINT * getNearestPoint(P const &pnt) const
GeoLib::Point const & getBeginPoint() const
GeoLib::Point const & getEndPoint() const
This class manages pointers to Points in a std::vector along with a name. It also handles the deletin...
std::size_t push_back(Point *pnt)
void resetInternalDataStructures()
std::size_t getSegmentNumber() const
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
virtual bool insertPoint(std::size_t pos, std::size_t pnt_id)
SegmentIterator begin() const
SegmentIterator end() const
A borehole as a geometric object.
std::size_t getID() const
const T * getCoords() const
virtual const Node * getNode(unsigned idx) const =0
virtual unsigned getNumberOfNodes() const =0
virtual unsigned getNumberOfEdges() const =0
Get the number of edges for this element.
virtual const Element * getEdge(unsigned i) const =0
Returns the i-th edge of the element.
virtual std::size_t getID() const final
Returns the ID of the element.
static Eigen::Vector3d getSurfaceNormal(Element const &e)
Returns the surface normal of a 2D element.
std::vector< MeshLib::Element const * > getElementsInVolume(POINT const &min, POINT const &max) const
Eigen::Vector3d const & getMinPoint() const
Eigen::Vector3d const & getMaxPoint() const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::size_t getNumberOfNodes() const
Get the number of nodes.
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
const Node * getNode(std::size_t idx) const
Get the node with the given index.
void makeVectorUnique(std::vector< T > &v)
bool isStation(GeoLib::Point const *pnt)
void sortSegments(MathLib::Point3d const &seg_beg_pnt, std::vector< GeoLib::LineSegment > &sub_segments)
bool isBorehole(GeoLib::Point const *pnt)
std::vector< MathLib::Point3d > lineSegmentIntersect2d(GeoLib::LineSegment const &ab, GeoLib::LineSegment const &cd)
std::unique_ptr< GeoLib::Point > triangleLineIntersection(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, MathLib::Point3d const &p, MathLib::Point3d const &q)
double sqrDist2d(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
std::pair< double, double > computeSqrNodeDistanceRange(MeshLib::Element const &element, bool const check_allnodes)
Compute the minimum and maximum node distances for this element.
TemplateElement< PointRule1 > Point
static std::vector< std::size_t > intersection(std::vector< std::size_t > const &vec1, std::vector< std::size_t > const &vec2)