MeshGeoToolsLib Namespace Reference

Classes
class	BoundaryElementsAlongPolyline
class	BoundaryElementsAtPoint
	This class collects point elements located at a given point elements. More...
class	BoundaryElementsOnSurface
class	BoundaryElementsSearcher
class	GeoMapper
	A set of tools for mapping the elevation of geometric objects. More...
class	HeuristicSearchLength
class	MeshNodesAlongPolyline
class	MeshNodesAlongSurface
class	MeshNodeSearcher
class	MeshNodesOnPoint
class	SearchLength

Enumerations
enum class	SearchAllNodes : bool { Yes = true , No = false }

Functions
std::unique_ptr< MeshLib::Mesh >	appendLinesAlongPolylines (const MeshLib::Mesh &mesh, const GeoLib::PolylineVec &ply_vec)
template<typename CacheType, typename GeometryType>
std::vector< MeshLib::Element * > const &	getBoundaryElements (std::vector< CacheType * > &cached_elements, std::function< GeometryType(CacheType const &)> getCachedItem, GeometryType const &item, MeshLib::Mesh const &mesh, MeshNodeSearcher const &mesh_node_searcher, bool const multiple_nodes_allowed)
std::string	meshNameFromGeometry (std::string const &geometrical_set_name, std::string const &geometry_name)
template<typename GeometryVec>
std::vector< std::unique_ptr< MeshLib::Mesh > >	constructAdditionalMeshesFromGeometries (std::vector< GeometryVec * > const &geometries, MeshGeoToolsLib::BoundaryElementsSearcher &boundary_element_searcher, bool const multiple_nodes_allowed)
std::vector< std::unique_ptr< MeshLib::Mesh > >	constructAdditionalMeshesFromGeoObjects (GeoLib::GEOObjects const &geo_objects, MeshLib::Mesh const &mesh, std::unique_ptr< SearchLength > search_length_algorithm, bool const multiple_nodes_allowed)
std::unique_ptr< MeshGeoToolsLib::SearchLength >	createSearchLengthAlgorithm (BaseLib::ConfigTree const &external_config, MeshLib::Mesh const &mesh)
static MeshLib::Element const *	findElementContainingPointXY (std::vector< MeshLib::Element const * > const &elements, MathLib::Point3d const &p)
static std::vector< MathLib::Point3d >	computeElementSegmentIntersections (MeshLib::Element const &elem, GeoLib::LineSegment const &segment)
static std::vector< GeoLib::LineSegment >	createSubSegmentsForElement (std::vector< MathLib::Point3d > const &intersections, MeshLib::Element const *const beg_elem, MeshLib::Element const *const end_elem, MathLib::Point3d const &beg_pnt, MathLib::Point3d const &end_pnt, MeshLib::Element const *const elem)
static std::vector< GeoLib::LineSegment >	mapLineSegment (GeoLib::LineSegment const &segment, std::vector< MeshLib::Element const * > const &surface_elements, MeshLib::Element const *const beg_elem, MeshLib::Element const *const end_elem)
static void	mapPointOnSurfaceElement (MeshLib::Element const &elem, MathLib::Point3d &q)
static std::vector< MeshLib::Element const * >	getCandidateElementsForLineSegmentIntersection (MeshLib::MeshElementGrid const &mesh_element_grid, GeoLib::LineSegment const &segment)
static bool	snapPointToElementNode (MathLib::Point3d &p, MeshLib::Element const &elem, double rel_eps)
static void	insertSubSegments (GeoLib::Polyline &ply, GeoLib::PointVec &points, GeoLib::Polyline::SegmentIterator &segment_it, std::vector< GeoLib::LineSegment > const &sub_segments)
static void	mapPolylineOnSurfaceMesh (GeoLib::Polyline &ply, GeoLib::PointVec &orig_points, MeshLib::MeshElementGrid const &mesh_element_grid)
void	identifySubdomainMesh (MeshLib::Mesh &subdomain_mesh, MeshLib::Mesh const &bulk_mesh, MeshNodeSearcher const &mesh_node_searcher, bool const force_overwrite=false)
std::vector< bool >	markNodesOutSideOfPolygon (std::vector< MeshLib::Node * > const &nodes, GeoLib::Polygon const &polygon)
template<typename PT>
void	resetMeshElementProperty (MeshLib::Mesh &mesh, GeoLib::Polygon const &polygon, std::string const &property_name, PT new_property_value, int restrict_to_material_id, bool const any_of)
template<typename CacheType, typename GeometryType>
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, SearchAllNodes const search_all_nodes)

Enumeration Type Documentation

SearchAllNodes

enum class MeshGeoToolsLib::SearchAllNodes : bool

strong

Enumerator
Yes
No

Definition at line 8 of file SearchAllNodes.h.

{
    Yes = true,
    No = false
};

Function Documentation

appendLinesAlongPolylines()

std::unique_ptr< MeshLib::Mesh > MeshGeoToolsLib::appendLinesAlongPolylines	(	const MeshLib::Mesh &	mesh,
		const GeoLib::PolylineVec &	ply_vec )

Add line elements to a copy of a given mesh

The function creates line elements from nodes located along user-provided polylines. New elements will have a distinct material ID for each polyline.

Remarks

The function allows creation of duplicated line elements.

Line elements may not be placed along edges of existing elements.

Parameters

mesh	original mesh
ply_vec	polyline vector whose nodes are used to create line elements

Returns

a new mesh which is copied from a given mesh and additionally includes line elements

Definition at line 23 of file AppendLinesAlongPolyline.cpp.

{
    // 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);
 
    auto const material_ids = materialIDs(mesh);
    assert(material_ids != nullptr);
    int const max_matID = material_ids ? ranges::max(*material_ids) : 0;
 
    auto const edgeLengths = minMaxEdgeLength(mesh.getElements());
    double const min_edge = edgeLengths.first;
 
    std::vector<int> new_mat_ids;
    const std::size_t n_ply(ply_vec.size());
    // for each polyline
    for (std::size_t k(0); k < n_ply; k++)
    {
        auto const* const ply = ply_vec.getVector()[k];
 
        // search nodes on the polyline
        MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly(
            mesh, *ply, min_edge * 0.5, MeshGeoToolsLib::SearchAllNodes::Yes);
        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);
            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);
        }
    }
 
    // generate a mesh
    const std::string name = mesh.getName() + "_with_lines";
    auto new_mesh =
        std::make_unique<MeshLib::Mesh>(name, vec_new_nodes, vec_new_eles,
                                        true /* compute_element_neighbors */);
    auto new_material_ids =
        new_mesh->getProperties().createNewPropertyVector<int>(
            "MaterialIDs", MeshLib::MeshItemType::Cell);
    if (!new_material_ids)
    {
        OGS_FATAL("Could not create MaterialIDs cell vector in new mesh.");
    }
 
    auto initial_values =
        material_ids ? ranges::any_view<int const>(*material_ids)
                     : ranges::views::repeat_n(0, mesh.getNumberOfElements());
 
    new_material_ids->assign(ranges::views::common(
        ranges::views::concat(initial_values, new_mat_ids)));
 
    return new_mesh;
}

References MeshLib::Cell, MeshLib::copyElementVector(), MeshLib::copyNodeVector(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getName(), GeoLib::TemplateVec< T >::getNameOfElementByID(), MeshGeoToolsLib::MeshNodesAlongPolyline::getNodeIDs(), MeshLib::Mesh::getNodes(), MeshLib::Mesh::getNumberOfElements(), GeoLib::TemplateVec< T >::getVector(), INFO(), OGS_FATAL, GeoLib::TemplateVec< T >::size(), and Yes.

Referenced by main().

computeElementSegmentIntersections()

std::vector< MathLib::Point3d > MeshGeoToolsLib::computeElementSegmentIntersections ( MeshLib::Element const & elem, GeoLib::LineSegment const & segment )

static

Definition at line 284 of file GeoMapper.cpp.

{
    std::vector<MathLib::Point3d> element_intersections;
    for (std::size_t k(0); k < elem.getNumberOfEdges(); ++k)
    {
        auto const edge =
            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),
            new GeoLib::Point(*dynamic_cast<MathLib::Point3d*>(
                                  const_cast<MeshLib::Node*>(edge->getNode(1))),
                              0),
            true};
        std::vector<MathLib::Point3d> const intersections(
            GeoLib::lineSegmentIntersect2d(segment, elem_segment));
        element_intersections.insert(end(element_intersections),
                                     begin(intersections), end(intersections));
    }
    return element_intersections;
}

References computeElementSegmentIntersections(), MeshLib::Element::getEdge(), MeshLib::Element::getNumberOfEdges(), and GeoLib::lineSegmentIntersect2d().

Referenced by computeElementSegmentIntersections(), and mapLineSegment().

constructAdditionalMeshesFromGeometries()

template<typename GeometryVec>

std::vector< std::unique_ptr< MeshLib::Mesh > > MeshGeoToolsLib::constructAdditionalMeshesFromGeometries	(	std::vector< GeometryVec * > const &	geometries,
		MeshGeoToolsLib::BoundaryElementsSearcher &	boundary_element_searcher,
		bool const	multiple_nodes_allowed )

Definition at line 32 of file ConstructMeshesFromGeometries.cpp.

{
    std::vector<std::unique_ptr<MeshLib::Mesh>> additional_meshes;
 
#ifdef USE_PETSC
    // The subdomain_mesh is not yet a NodePartitionedMesh.
    // The bulk_mesh, which is a NodePartitionedMesh, is needed to construct the
    // subdomain NodePartitionedMesh.
    auto const* bulk_mesh = dynamic_cast<MeshLib::NodePartitionedMesh const*>(
        &boundary_element_searcher.mesh);
#endif
    for (GeometryVec* const geometry_vec : geometries)
    {
        // Each geometry_vec has a name, this is the first part of the full
        // name.
        auto const& vec_name = geometry_vec->getName();
 
        auto const& vec_data = geometry_vec->getVector();
 
        auto const vec_size = geometry_vec->size();
        for (std::size_t i = 0; i < vec_size; ++i)
        {
            // Each geometry has a name, this is the second part of the full
            // name.
            std::string geometry_name;
            bool const is_geometry_named =
                geometry_vec->getNameOfElementByID(i, geometry_name);
            if (!is_geometry_named)
            {
                continue;
            }
 
            auto const& geometry = *vec_data[i];
 
            DBUG("Creating mesh from geometry {:s} {:s}.", vec_name,
                 geometry_name);
 
            auto subdomain_mesh = createMeshFromElementSelection(
                meshNameFromGeometry(vec_name, geometry_name),
                MeshLib::cloneElements(
                    boundary_element_searcher.getBoundaryElements(
                        geometry, multiple_nodes_allowed)));
 
#ifdef USE_PETSC
            additional_meshes.push_back(
                MeshLib::transformMeshToNodePartitionedMesh(
                    bulk_mesh, subdomain_mesh.get()));
#else
            // Nothing special to be done in the serial case.
            additional_meshes.emplace_back(std::move(subdomain_mesh));
#endif
        }
    }
    return additional_meshes;
}

References MeshLib::cloneElements(), DBUG(), MeshGeoToolsLib::BoundaryElementsSearcher::getBoundaryElements(), MeshGeoToolsLib::BoundaryElementsSearcher::mesh, meshNameFromGeometry(), and MeshLib::transformMeshToNodePartitionedMesh().

Referenced by constructAdditionalMeshesFromGeoObjects().

constructAdditionalMeshesFromGeoObjects()

std::vector< std::unique_ptr< MeshLib::Mesh > > MeshGeoToolsLib::constructAdditionalMeshesFromGeoObjects	(	GeoLib::GEOObjects const &	geo_objects,
		MeshLib::Mesh const &	mesh,
		std::unique_ptr< SearchLength >	search_length_algorithm,
		bool const	multiple_nodes_allowed )

For each named geometry in the give geo_objects (defined on the given mesh) constructs a mesh corresponding to the geometry with mappings to the bulk mesh elements and nodes.

Definition at line 92 of file ConstructMeshesFromGeometries.cpp.

{
    std::vector<std::unique_ptr<MeshLib::Mesh>> additional_meshes;
 
    auto const& mesh_node_searcher =
        MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
            mesh, std::move(search_length_algorithm));
 
    MeshGeoToolsLib::BoundaryElementsSearcher boundary_element_searcher(
        mesh, mesh_node_searcher);
 
    //
    // Points
    //
    {
        auto point_meshes = constructAdditionalMeshesFromGeometries(
            geo_objects.getPoints(), boundary_element_searcher,
            multiple_nodes_allowed);
        std::move(begin(point_meshes), end(point_meshes),
                  std::back_inserter(additional_meshes));
    }
 
    //
    // Polylines
    //
    {
        auto polyline_meshes = constructAdditionalMeshesFromGeometries(
            geo_objects.getPolylines(), boundary_element_searcher,
            multiple_nodes_allowed);
        std::move(begin(polyline_meshes), end(polyline_meshes),
                  std::back_inserter(additional_meshes));
    }
 
    // Surfaces
    {
        auto surface_meshes = constructAdditionalMeshesFromGeometries(
            geo_objects.getSurfaces(), boundary_element_searcher,
            multiple_nodes_allowed);
        std::move(begin(surface_meshes), end(surface_meshes),
                  std::back_inserter(additional_meshes));
    }
 
    // Set axial symmetry for the additional meshes to the same value as the
    // "bulk" mesh.
    std::for_each(begin(additional_meshes), end(additional_meshes),
                  [axial_symmetry = mesh.isAxiallySymmetric()](auto& m)
                  { m->setAxiallySymmetric(axial_symmetry); });
    return additional_meshes;
}

References constructAdditionalMeshesFromGeometries(), MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(), GeoLib::GEOObjects::getPoints(), GeoLib::GEOObjects::getPolylines(), GeoLib::GEOObjects::getSurfaces(), and MeshLib::Mesh::isAxiallySymmetric().

Referenced by anonymous_namespace{ProjectData.cpp}::readMeshes().

createSearchLengthAlgorithm()

std::unique_ptr< MeshGeoToolsLib::SearchLength > MeshGeoToolsLib::createSearchLengthAlgorithm	(	BaseLib::ConfigTree const &	external_config,
		MeshLib::Mesh const &	mesh )

Creates a search length algorithm from the given config.

In case, that there is no tag for the search length algorithm found in the config, the default SearchLength algorithm is returned.

Input File Parameter

prj__search_length_algorithm

Input File Parameter

prj__search_length_algorithm__type

Input File Parameter

prj__search_length_algorithm__fixed

Input File Parameter

prj__search_length_algorithm__fixed__value

Input File Parameter

prj__search_length_algorithm__heuristic

Definition at line 13 of file CreateSearchLength.cpp.

{
    std::optional<BaseLib::ConfigTree> config =
        external_config.getConfigSubtreeOptional("search_length_algorithm");
 
    if (!config)
    {
        return std::make_unique<MeshGeoToolsLib::SearchLength>();
    }

    std::string const type = config->getConfigParameter<std::string>("type");

    if (type == "fixed")
    {
        auto const length = config->getConfigParameter<double>("value");
        return std::make_unique<MeshGeoToolsLib::SearchLength>(length);
    }
    if (type == "heuristic")
    {
        return std::make_unique<HeuristicSearchLength>(mesh);
    }
    OGS_FATAL("Unknown search length algorithm type '{:s}'.", type);
}

References BaseLib::ConfigTree::getConfigSubtreeOptional(), and OGS_FATAL.

Referenced by anonymous_namespace{ProjectData.cpp}::readMeshes().

createSubSegmentsForElement()

std::vector< GeoLib::LineSegment > MeshGeoToolsLib::createSubSegmentsForElement ( std::vector< MathLib::Point3d > const & intersections, MeshLib::Element const *const beg_elem, MeshLib::Element const *const end_elem, MathLib::Point3d const & beg_pnt, MathLib::Point3d const & end_pnt, MeshLib::Element const *const elem )

static

Definition at line 308 of file GeoMapper.cpp.

{
    std::vector<GeoLib::LineSegment> sub_segments;
    if (intersections.size() > 2)
    {
        std::stringstream out;
        out << "element with id " << elem->getID() << " and seg "
            << " intersecting at more than two edges\n";
        for (std::size_t k(0); k < intersections.size(); ++k)
        {
            out << k << " " << intersections[k] << "\n";
        }
        out << "Could not map segment on element. Aborting.\n";
        OGS_FATAL("{:s}", out.str());
    }
 
    if (intersections.size() == 1 && elem == beg_elem)
    {
        // The line segment intersects the element that contains the begin
        // point of the line segment. Here the first sub line segment is
        // added.
        if (MathLib::sqrDist(beg_pnt, intersections[0]) >
            std::numeric_limits<double>::epsilon())
        {
            sub_segments.emplace_back(new GeoLib::Point{beg_pnt, 0},
                                      new GeoLib::Point{intersections[0], 0},
                                      true);
        }
    }
 
    if (intersections.size() == 1 && elem == end_elem)
    {
        // The line segment intersects the element that contains the end
        // point of the line segment. Here the last sub line segment is
        // added.
        if (MathLib::sqrDist(end_pnt, intersections[0]) >
            std::numeric_limits<double>::epsilon())
        {
            sub_segments.emplace_back(new GeoLib::Point{intersections[0], 0},
                                      new GeoLib::Point{end_pnt, 0}, true);
        }
    }
 
    if (intersections.size() == 1 && (elem != beg_elem && elem != end_elem))
    {
        // Since the line segment enters and leaves the element in the same
        // point there isn't any need to insert a new sub line segment.
        return sub_segments;
    }
 
    // create sub segment for the current element
    if (intersections.size() == 2)
    {
        sub_segments.emplace_back(new GeoLib::Point{intersections[0], 0},
                                  new GeoLib::Point{intersections[1], 0}, true);
    }
    return sub_segments;
}

References createSubSegmentsForElement(), MeshLib::Element::getID(), OGS_FATAL, and MathLib::sqrDist().

Referenced by createSubSegmentsForElement(), and mapLineSegment().

findElementContainingPointXY()

MeshLib::Element const * MeshGeoToolsLib::findElementContainingPointXY ( std::vector< MeshLib::Element const * > const & elements, MathLib::Point3d const & p )

static

Find the 2d-element within the elements that contains the given point p.

The algorithm projects every element of the elements vector and the point p orthogonal to the \(x\)- \(y\) plane. In the \(x\)- \(y\) plane it is checked if the projected point is in the projected element.

Definition at line 249 of file GeoMapper.cpp.

{
    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)
        {
            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)
        {
            (*const_cast<MeshLib::Node*>(elem_2d->getNode(k)))[2] = 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)
            {
                delete elem_2d->getNode(k);
            }
            return elem;
        }
        // clean up the copied nodes
        for (std::size_t k(0); k < elem_2d->getNumberOfNodes(); ++k)
        {
            delete elem_2d->getNode(k);
        }
    }
    return nullptr;
}

Referenced by mapPolylineOnSurfaceMesh().

getBoundaryElements()

template<typename CacheType, typename GeometryType>

std::vector< MeshLib::Element * > const & MeshGeoToolsLib::getBoundaryElements	(	std::vector< CacheType * > &	cached_elements,
		std::function< GeometryType(CacheType const &)>	getCachedItem,
		GeometryType const &	item,
		MeshLib::Mesh const &	mesh,
		MeshNodeSearcher const &	mesh_node_searcher,
		bool const	multiple_nodes_allowed )

Definition at line 43 of file BoundaryElementsSearcher.cpp.

{
    if (auto const it = find_if(cbegin(cached_elements), cend(cached_elements),
                                [&](auto const& element)
                                { return getCachedItem(*element) == item; });
        it != cend(cached_elements))
    {
        return (*it)->getBoundaryElements();
    }
    // create new boundary elements
    if constexpr (std::is_convertible<GeometryType, GeoLib::Point>::value)
    {
        cached_elements.push_back(new CacheType(mesh, mesh_node_searcher, item,
                                                multiple_nodes_allowed));
    }
    else
    {
        cached_elements.push_back(
            new CacheType(mesh, mesh_node_searcher, item));
    }
    return cached_elements.back()->getBoundaryElements();
}

Referenced by MeshGeoToolsLib::BoundaryElementsSearcher::getBoundaryElements().

getCandidateElementsForLineSegmentIntersection()

std::vector< MeshLib::Element const * > MeshGeoToolsLib::getCandidateElementsForLineSegmentIntersection ( MeshLib::MeshElementGrid const & mesh_element_grid, GeoLib::LineSegment const & segment )

static

Definition at line 454 of file GeoMapper.cpp.

{
    GeoLib::LineSegment seg_deep_copy(
        new GeoLib::Point(segment.getBeginPoint()),
        new GeoLib::Point(segment.getEndPoint()), true);
    // modify z coordinates such that all surface elements around the line
    // segment are found
    seg_deep_copy.getBeginPoint()[2] = mesh_element_grid.getMinPoint()[2];
    seg_deep_copy.getEndPoint()[2] = mesh_element_grid.getMaxPoint()[2];
    std::array<MathLib::Point3d, 2> const pnts{
        {seg_deep_copy.getBeginPoint(), seg_deep_copy.getEndPoint()}};
    GeoLib::AABB aabb(pnts.cbegin(), pnts.cend());
 
    // TODO TF: remove after getElementsInVolume interface change
    auto convert_to_Point3d = [](Eigen::Vector3d const& v) {
        return MathLib::Point3d{std::array{v[0], v[1], v[2]}};
    };
 
    auto const min = convert_to_Point3d(aabb.getMinPoint());
    auto const max = convert_to_Point3d(aabb.getMaxPoint());
    auto candidate_elements = mesh_element_grid.getElementsInVolume(min, max);
 
    // make candidate elements unique
    BaseLib::makeVectorUnique(candidate_elements);
 
    return candidate_elements;
}

References GeoLib::LineSegment::getBeginPoint(), getCandidateElementsForLineSegmentIntersection(), MeshLib::MeshElementGrid::getElementsInVolume(), GeoLib::LineSegment::getEndPoint(), GeoLib::AABB::getMaxPoint(), MeshLib::MeshElementGrid::getMaxPoint(), GeoLib::AABB::getMinPoint(), MeshLib::MeshElementGrid::getMinPoint(), and BaseLib::makeVectorUnique().

Referenced by getCandidateElementsForLineSegmentIntersection(), and mapPolylineOnSurfaceMesh().

getMeshNodeIDs()

template<typename CacheType, typename GeometryType>

std::vector< std::size_t > const & MeshGeoToolsLib::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,
		SearchAllNodes const	search_all_nodes )

Definition at line 56 of file MeshNodeSearcher.cpp.

{
    if (auto const it = find_if(cbegin(cached_elements), cend(cached_elements),
                                [&](auto const& element)
                                { return getCachedItem(*element) == item; });
        it != cend(cached_elements))
    {
        return (*it)->getNodeIDs();
    }
    // search IDs for geometry object
    if constexpr (std::is_convertible<GeometryType, GeoLib::Point>::value)
    {
        cached_elements.push_back(new CacheType(
            mesh, mesh_grid, item, search_length, search_all_nodes));
    }
    else
    {
        cached_elements.push_back(
            new CacheType(mesh, item, search_length, search_all_nodes));
    }
    return cached_elements.back()->getNodeIDs();
}

Referenced by MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeIDs().

identifySubdomainMesh()

void MeshGeoToolsLib::identifySubdomainMesh	(	MeshLib::Mesh &	subdomain_mesh,
		MeshLib::Mesh const &	bulk_mesh,
		MeshNodeSearcher const &	mesh_node_searcher,
		bool const	force_overwrite = false )

Geometrically finds nodes and elements of the subdomain mesh in the bulk mesh, and updates or verifies the corresponding bulk_node_ids and bulk_element_ids properties.

In case of unexpected results OGS_FATAL is called.

Definition at line 232 of file IdentifySubdomainMesh.cpp.

{
    BaseLib::RunTime time;
    time.start();
    auto const& bulk_node_ids =
        identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
    INFO("identifySubdomainMesh(): identifySubdomainMeshNodes took {:g} s",
         time.elapsed());
 
    updateOrCheckExistingSubdomainProperty(
        subdomain_mesh, MeshLib::getBulkIDString(MeshLib::MeshItemType::Node),
        bulk_node_ids, MeshLib::MeshItemType::Node, force_overwrite);
 
    time.start();
    auto const& bulk_element_ids =
        identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
    INFO("identifySubdomainMesh(): identifySubdomainMeshElements took {:g} s",
         time.elapsed());
 
    // The bulk_element_ids could be of two types: one element per entry---this
    // is the expected case for the boundary meshes; multiple elements per
    // entry---this happens if the subdomain mesh lies inside the bulk mesh and
    // has lower dimension.
    // First find out the type, then add/check the CellData or FieldData.
    bool const all_single_elements =
        all_of(begin(bulk_element_ids), end(bulk_element_ids),
               [](std::vector<std::size_t> const& v) { return v.size() == 1; });
    if (all_single_elements)
    {
        // All vectors are of size 1, so the data can be flattened and
        // stored in CellData or compared to existing CellData.
        std::vector<std::size_t> unique_bulk_element_ids;
        unique_bulk_element_ids.reserve(bulk_element_ids.size());
        transform(begin(bulk_element_ids), end(bulk_element_ids),
                  back_inserter(unique_bulk_element_ids),
                  [](std::vector<std::size_t> const& v) { return v[0]; });
 
        updateOrCheckExistingSubdomainProperty(
            subdomain_mesh,
            MeshLib::getBulkIDString(MeshLib::MeshItemType::Cell),
            unique_bulk_element_ids, MeshLib::MeshItemType::Cell,
            force_overwrite);
    }
    else
    {
        // Some of the boundary elements are connected to multiple bulk
        // elements; Store the array in FieldData with additional CellData array
        // for the number of elements, which also provides the offsets.
        std::vector<std::size_t> flat_bulk_element_ids;
        flat_bulk_element_ids.reserve(2 * bulk_element_ids.size());  // Guess.
        std::vector<std::size_t> number_of_bulk_element_ids;
        number_of_bulk_element_ids.reserve(bulk_element_ids.size());
 
        for (auto const& v : bulk_element_ids)
        {
            // The element ids are already sorted descending by
            // findElementsInMesh(), making the ordering deterministic across
            // compilers and parallel execution.
            // TODO: sort the bulk element ids by the positions of the
            // corresponding subdomain element to the subdomain element, e.g.
            // left-to-right for 1D subdomains in 2D bulk meshes or 2D
            // subdomains in 3D bulk meshes.
            number_of_bulk_element_ids.push_back(v.size());
            flat_bulk_element_ids.insert(end(flat_bulk_element_ids), begin(v),
                                         end(v));
        }
 
        updateOrCheckExistingSubdomainProperty(
            subdomain_mesh, "number_bulk_elements", number_of_bulk_element_ids,
            MeshLib::MeshItemType::Cell, force_overwrite);
        updateOrCheckExistingSubdomainProperty(
            subdomain_mesh,
            MeshLib::getBulkIDString(MeshLib::MeshItemType::Cell),
            flat_bulk_element_ids, MeshLib::MeshItemType::IntegrationPoint,
            force_overwrite);
    }
}

References MeshLib::Cell, BaseLib::RunTime::elapsed(), MeshLib::getBulkIDString(), INFO(), MeshLib::IntegrationPoint, MeshLib::Node, and BaseLib::RunTime::start().

insertSubSegments()

void MeshGeoToolsLib::insertSubSegments ( GeoLib::Polyline & ply, GeoLib::PointVec & points, GeoLib::Polyline::SegmentIterator & segment_it, std::vector< GeoLib::LineSegment > const & sub_segments )

static

Definition at line 513 of file GeoMapper.cpp.

{
    std::size_t const j(segment_it.getSegmentNumber());
    std::size_t new_pnts_cnt(0);
    for (auto const& segment : sub_segments)
    {
        auto const begin_id(points.push_back(
            new GeoLib::Point(segment.getBeginPoint(), points.size())));
        if (ply.insertPoint(j + new_pnts_cnt + 1, begin_id))
        {
            new_pnts_cnt++;
        }
        auto const end_id(points.push_back(
            new GeoLib::Point(segment.getEndPoint(), points.size())));
        if (ply.insertPoint(j + new_pnts_cnt + 1, end_id))
        {
            new_pnts_cnt++;
        }
    }
    std::advance(segment_it, new_pnts_cnt);
}

References GeoLib::Polyline::SegmentIterator::getSegmentNumber(), GeoLib::Polyline::insertPoint(), insertSubSegments(), GeoLib::PointVec::push_back(), and GeoLib::TemplateVec< T >::size().

Referenced by insertSubSegments(), and mapPolylineOnSurfaceMesh().

mapLineSegment()

std::vector< GeoLib::LineSegment > MeshGeoToolsLib::mapLineSegment ( GeoLib::LineSegment const & segment, std::vector< MeshLib::Element const * > const & surface_elements, MeshLib::Element const *const beg_elem, MeshLib::Element const *const end_elem )

static

Definition at line 371 of file GeoMapper.cpp.

{
    std::vector<GeoLib::LineSegment> sub_segments;
    MathLib::Point3d const& beg_pnt(segment.getBeginPoint());
    MathLib::Point3d const& end_pnt(segment.getEndPoint());
 
    for (auto const elem : surface_elements)
    {
        // compute element-segment-intersections (2d in x-y-plane)
        std::vector<MathLib::Point3d> element_intersections(
            computeElementSegmentIntersections(*elem, segment));
        if (element_intersections.empty())
        {
            continue;
        }
 
        BaseLib::makeVectorUnique(element_intersections);
 
        std::vector<GeoLib::LineSegment> sub_seg_elem(
            createSubSegmentsForElement(element_intersections, beg_elem,
                                        end_elem, beg_pnt, end_pnt, elem));
        sub_segments.insert(sub_segments.end(), sub_seg_elem.begin(),
                            sub_seg_elem.end());
    }
 
    // beg_elem == nullptr means there isn't any element corresponding to the
    // beg_pnt and as a consequence the above algorithm doesn't insert a sub
    // segment
    if (beg_elem == nullptr)
    {
        auto min_dist_segment = std::min_element(
            sub_segments.begin(), sub_segments.end(),
            [&beg_pnt](GeoLib::LineSegment const& seg0,
                       GeoLib::LineSegment const& seg1)
            {
                // min dist for segment 0
                const double d0(
                    std::min(MathLib::sqrDist(beg_pnt, seg0.getBeginPoint()),
                             MathLib::sqrDist(beg_pnt, seg0.getEndPoint())));
                // min dist for segment 1
                const double d1(
                    std::min(MathLib::sqrDist(beg_pnt, seg1.getBeginPoint()),
                             MathLib::sqrDist(beg_pnt, seg1.getEndPoint())));
                return d0 < d1;
            });
        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()}
                : new GeoLib::Point{min_dist_segment->getEndPoint()}};
        sub_segments.emplace_back(new GeoLib::Point{beg_pnt, 0}, pnt, true);
    }
    // sort all sub segments for the given segment (beg_pnt, end_pnt)
    GeoLib::sortSegments(beg_pnt, sub_segments);
 
    sub_segments.erase(std::unique(sub_segments.begin(), sub_segments.end()),
                       sub_segments.end());
 
    return sub_segments;
}

References computeElementSegmentIntersections(), createSubSegmentsForElement(), GeoLib::LineSegment::getBeginPoint(), GeoLib::LineSegment::getEndPoint(), BaseLib::makeVectorUnique(), mapLineSegment(), GeoLib::sortSegments(), and MathLib::sqrDist().

Referenced by mapLineSegment(), and mapPolylineOnSurfaceMesh().

mapPointOnSurfaceElement()

void MeshGeoToolsLib::mapPointOnSurfaceElement ( MeshLib::Element const & elem, MathLib::Point3d & q )

static

Definition at line 436 of file GeoMapper.cpp.

{
    // create plane equation: n*p = d
    auto const& p = elem.getNode(0)->asEigenVector3d();
    Eigen::Vector3d const n(MeshLib::FaceRule::getSurfaceNormal(elem));
    if (n[2] == 0.0)
    {  // vertical plane, z coordinate is arbitrary
        q[2] = p[2];
    }
    else
    {
        double const d(n.dot(p));
        q[2] = (d - n[0] * q[0] - n[1] * q[1]) / n[2];
    }
}

References MathLib::Point3d::asEigenVector3d(), MeshLib::Element::getNode(), MeshLib::FaceRule::getSurfaceNormal(), and mapPointOnSurfaceElement().

Referenced by mapPointOnSurfaceElement(), and mapPolylineOnSurfaceMesh().

mapPolylineOnSurfaceMesh()

void MeshGeoToolsLib::mapPolylineOnSurfaceMesh ( GeoLib::Polyline & ply, GeoLib::PointVec & orig_points, MeshLib::MeshElementGrid const & mesh_element_grid )

static

Definition at line 538 of file GeoMapper.cpp.

{
    // for each segment ...
    for (auto segment_it(ply.begin()); segment_it != ply.end(); ++segment_it)
    {
        auto candidate_elements(getCandidateElementsForLineSegmentIntersection(
            mesh_element_grid, *segment_it));
 
        auto mapPoint = [&candidate_elements](MathLib::Point3d& p)
        {
            auto const* elem(
                findElementContainingPointXY(candidate_elements, p));
            if (elem)
            {
                if (!snapPointToElementNode(p, *elem, 1e-3))
                {
                    mapPointOnSurfaceElement(*elem, p);
                }
            }
            return elem;
        };
 
        // map segment begin and end point
        auto const* beg_elem(mapPoint((*segment_it).getBeginPoint()));
        auto const* end_elem(mapPoint((*segment_it).getEndPoint()));
 
        // Since the mapping of the segment begin and end points the coordinates
        // changed. The internal data structures of PointVec are possibly
        // invalid and hence it is necessary to re-create them.
        orig_points.resetInternalDataStructures();
 
        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:
            // 1. beg_elem == end_elem and the segment intersects elements
            // 2. beg_elem == end_elem and the segment does not intersect any
            // element, i.e., the segment is located outside of the mesh area
            //
            // Case 1 needs additional work.
            continue;
        }
 
        // map the line segment (and if necessary for the mapping partition it)
        std::vector<GeoLib::LineSegment> sub_segments(mapLineSegment(
            *segment_it, candidate_elements, beg_elem, end_elem));
 
        if (sub_segments.empty())
        {
            continue;
        }
 
        // The case sub_segment.size() == 1 is already handled above.
 
        if (sub_segments.size() > 1)
        {
            insertSubSegments(ply, orig_points, segment_it, sub_segments);
        }
    }
}

References GeoLib::Polyline::begin(), GeoLib::Polyline::end(), findElementContainingPointXY(), getCandidateElementsForLineSegmentIntersection(), insertSubSegments(), mapLineSegment(), mapPointOnSurfaceElement(), mapPolylineOnSurfaceMesh(), GeoLib::PointVec::resetInternalDataStructures(), and snapPointToElementNode().

Referenced by MeshGeoToolsLib::GeoMapper::advancedMapOnMesh(), and mapPolylineOnSurfaceMesh().

markNodesOutSideOfPolygon()

std::vector< bool > MeshGeoToolsLib::markNodesOutSideOfPolygon	(	std::vector< MeshLib::Node * > const &	nodes,
		GeoLib::Polygon const &	polygon )

Definition at line 16 of file MarkNodesOutsideOfPolygon.h.

{
    // *** rotate polygon points to xy-plane
    auto [rotated_polygon_points, normal] =
        GeoLib::rotatePolygonPointsToXY(polygon);
 
    // *** rotate mesh nodes to xy-plane
    // 1 copy all mesh nodes to GeoLib::Points
    std::vector<GeoLib::Point*> rotated_nodes;
    for (auto node : nodes)
    {
        rotated_nodes.push_back(new GeoLib::Point(*node, node->getID()));
    }
    // 2 rotate the Points
    Eigen::Matrix3d const rot_mat = GeoLib::computeRotationMatrixToXY(normal);
    GeoLib::rotatePoints(rot_mat, rotated_nodes);
    // 3 set z coord to zero
    std::for_each(rotated_nodes.begin(), rotated_nodes.end(),
                  [](GeoLib::Point* p) { (*p)[2] = 0.0; });
 
    std::vector<bool> outside(rotated_nodes.size(), true);
    // *** mark rotated nodes inside rotated polygon
    {
        // create new polygon using the rotated points
        GeoLib::Polyline rotated_polyline(rotated_polygon_points);
        for (std::size_t k(0); k < polygon.getNumberOfPoints(); k++)
        {
            rotated_polyline.addPoint(k);
        }
        rotated_polyline.addPoint(0);
        GeoLib::Polygon const rotated_polygon(rotated_polyline);
 
        for (std::size_t k(0); k < rotated_nodes.size(); k++)
        {
            if (rotated_polygon.isPntInPolygon(*(rotated_nodes[k])))
            {
                outside[k] = false;
            }
        }
    }
 
    for (auto& rotated_node : rotated_nodes)
    {
        delete rotated_node;
    }
 
    for (auto& rot_polygon_pnt : rotated_polygon_points)
    {
        delete rot_polygon_pnt;
    }
 
    return outside;
}

References GeoLib::Polyline::addPoint(), GeoLib::computeRotationMatrixToXY(), GeoLib::Polyline::getNumberOfPoints(), GeoLib::Polygon::isPntInPolygon(), GeoLib::rotatePoints(), and GeoLib::rotatePolygonPointsToXY().

Referenced by resetMeshElementProperty().

meshNameFromGeometry()

std::string MeshGeoToolsLib::meshNameFromGeometry	(	std::string const &	geometrical_set_name,
		std::string const &	geometry_name )

Definition at line 24 of file ConstructMeshesFromGeometries.cpp.

{
    return geometrical_set_name + "_" + geometry_name;
}

Referenced by constructAdditionalMeshesFromGeometries(), and anonymous_namespace{ProcessVariable.cpp}::findMeshInConfig().

resetMeshElementProperty()

template<typename PT>

void MeshGeoToolsLib::resetMeshElementProperty	(	MeshLib::Mesh &	mesh,
		GeoLib::Polygon const &	polygon,
		std::string const &	property_name,
		PT	new_property_value,
		int	restrict_to_material_id,
		bool const	any_of )

Definition at line 26 of file ResetMeshElementProperty.h.

{
    auto* const pv = MeshLib::getOrCreateMeshProperty<PT>(
        mesh, property_name, MeshLib::MeshItemType::Cell, 1);
 
    if (pv->getMeshItemType() != MeshLib::MeshItemType::Cell)
    {
        ERR("Values of the PropertyVector are not assigned to cells.");
        return;
    }
    auto const* material_ids =
        mesh.getProperties().getPropertyVector<int>("MaterialIDs");
 
    if (restrict_to_material_id != -1 && !material_ids)
    {
        OGS_FATAL(
            "Restriction of resetting a property in a polygonal region "
            "requires that a MaterialIDs data array is available in the "
            "mesh.");
    }
 
    auto has_element_required_material_id = [&](int const element_id)
    {
        return restrict_to_material_id == -1 ||
               (*material_ids)[element_id] == restrict_to_material_id;
    };
 
    auto const outside = markNodesOutSideOfPolygon(mesh.getNodes(), polygon);
 
    auto is_node_outside = [&outside](std::size_t const node_id)
    { return outside[node_id]; };
 
    auto is_element_outside = [&](MeshLib::Element const& e)
    {
        auto ids = e.nodes() | MeshLib::views::ids;
        return any_of ? ranges::all_of(ids, is_node_outside)
                      : ranges::any_of(ids, is_node_outside);
    };
 
    auto is_valid_element = [&](MeshLib::Element const& e)
    {
        return !is_element_outside(e) &&
               has_element_required_material_id(e.getID());
    };
 
    auto compute_value = [&](MeshLib::Element const* const e) -> PT
    { return is_valid_element(*e) ? new_property_value : (*pv)[e->getID()]; };
 
    pv->assign(mesh.getElements() | ranges::views::transform(compute_value));
}

References MeshLib::Cell, ERR(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNodes(), MeshLib::getOrCreateMeshProperty(), MeshLib::Mesh::getProperties(), MeshLib::Properties::getPropertyVector(), MeshLib::views::ids, markNodesOutSideOfPolygon(), and OGS_FATAL.

Referenced by main().

snapPointToElementNode()

bool MeshGeoToolsLib::snapPointToElementNode ( MathLib::Point3d & p, MeshLib::Element const & elem, double rel_eps )

static

Definition at line 484 of file GeoMapper.cpp.

{
    // 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)
    {
        auto const& node(*elem.getNode(k));
        double const sqr_dist_2d(MathLib::sqrDist2d(p, node));
        if (sqr_dist_2d < sqr_eps)
        {
#ifdef DEBUG_GEOMAPPER
            std::stringstream out;
            out.precision(std::numeric_limits<double>::max_digits10);
            out << "Segment point snapped from " << p;
#endif
            p = node;
#ifdef DEBUG_GEOMAPPER
            out << "to " << p;
            DBUG("{:s}", out.str());
#endif
            return true;
        }
    }
    return false;
}

References MeshLib::computeSqrNodeDistanceRange(), DBUG(), MeshLib::Element::getNode(), MeshLib::Element::getNumberOfNodes(), snapPointToElementNode(), and MathLib::sqrDist2d().

Referenced by mapPolylineOnSurfaceMesh(), and snapPointToElementNode().