MeshToolsLib::MeshLayerMapper Class Reference

Detailed Description

Manipulating and adding prism element layers to an existing 2D mesh.

Definition at line 25 of file MeshLayerMapper.h.

#include <MeshLayerMapper.h>

Inheritance diagram for MeshToolsLib::MeshLayerMapper:

Collaboration diagram for MeshToolsLib::MeshLayerMapper:

Public Member Functions
	~MeshLayerMapper () override=default
bool	createRasterLayers (MeshLib::Mesh const &mesh, std::vector< GeoLib::Raster const * > const &rasters, double minimum_thickness, double noDataReplacementValue=0.0) override
Public Member Functions inherited from LayeredMeshGenerator
virtual bool	createLayers (MeshLib::Mesh const &mesh, std::vector< GeoLib::Raster const * > const &rasters, double minimum_thickness, double noDataReplacementValue=0.0) final
std::unique_ptr< MeshLib::Mesh >	getMesh (std::string const &mesh_name) const
	Returns a mesh of the subsurface representation.

Static Public Member Functions
static MeshLib::Mesh *	createStaticLayers (MeshLib::Mesh const &mesh, std::vector< float > const &layer_thickness_vector, std::string const &mesh_name="SubsurfaceMesh")
static bool	layerMapping (MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double nodata_replacement=0.0, bool const ignore_nodata=false)
static bool	mapToStaticValue (MeshLib::Mesh const &mesh, double value)
	Maps the elevation of all mesh nodes to the specified static value.

Private Member Functions
void	addLayerToMesh (const MeshLib::Mesh &dem_mesh, unsigned layer_id, GeoLib::Raster const &raster) override
	Adds another layer to a subsurface mesh.

Additional Inherited Members
Protected Member Functions inherited from LayeredMeshGenerator
	LayeredMeshGenerator ()=default
virtual	~LayeredMeshGenerator ()=default
MeshLib::Node *	getNewLayerNode (MeshLib::Node const &dem_node, MeshLib::Node const &last_layer_node, GeoLib::Raster const &raster, std::size_t new_node_id) const
Static Protected Member Functions inherited from LayeredMeshGenerator
static double	calcEpsilon (GeoLib::Raster const &low, GeoLib::Raster const &high)
	Calculates a data-dependent epsilon value.
Protected Attributes inherited from LayeredMeshGenerator
double	_elevation_epsilon {0.0001}
double	_minimum_thickness {std::numeric_limits<float>::epsilon()}
std::vector< int >	_materials
std::vector< MeshLib::Node * >	_nodes
std::vector< MeshLib::Element * >	_elements

Constructor & Destructor Documentation

~MeshLayerMapper()

MeshToolsLib::MeshLayerMapper::~MeshLayerMapper ( )

overridedefault

Member Function Documentation

addLayerToMesh()

void MeshToolsLib::MeshLayerMapper::addLayerToMesh ( const MeshLib::Mesh & dem_mesh, unsigned layer_id, GeoLib::Raster const & raster )

overrideprivatevirtual

Adds another layer to a subsurface mesh.

Implements LayeredMeshGenerator.

Definition at line 206 of file MeshLayerMapper.cpp.

{
    const unsigned pyramid_base[3][4] = {
        {1, 3, 4, 2},  // Point 4 missing
        {2, 4, 3, 0},  // Point 5 missing
        {0, 3, 4, 1},  // Point 6 missing
    };
 
    std::size_t const nNodes = dem_mesh.getNumberOfNodes();
    std::vector<MeshLib::Node*> const& top_nodes = dem_mesh.getNodes();
    int const last_layer_node_offset = layer_id * nNodes;
 
    // add nodes for new layer
    for (std::size_t i = 0; i < nNodes; ++i)
    {
        _nodes.push_back(getNewLayerNode(*top_nodes[i],
                                         *_nodes[last_layer_node_offset + i],
                                         raster, _nodes.size()));
    }
 
    std::vector<MeshLib::Element*> const& elems = dem_mesh.getElements();
    std::size_t const nElems(dem_mesh.getNumberOfElements());
 
    for (std::size_t i = 0; i < nElems; ++i)
    {
        MeshLib::Element* elem(elems[i]);
        if (elem->getGeomType() != MeshLib::MeshElemType::TRIANGLE)
        {
            continue;
        }
        unsigned node_counter(3);
        unsigned missing_idx(0);
        std::array<MeshLib::Node*, 6> new_elem_nodes{};
        for (unsigned j = 0; j < 3; ++j)
        {
            new_elem_nodes[j] =
                _nodes[_nodes[last_layer_node_offset + getNodeIndex(*elem, j)]
                           ->getID()];
            new_elem_nodes[node_counter] =
                (_nodes[last_layer_node_offset + getNodeIndex(*elem, j) +
                        nNodes]);
            if (new_elem_nodes[j]->getID() !=
                new_elem_nodes[node_counter]->getID())
            {
                node_counter++;
            }
            else
            {
                missing_idx = j;
            }
        }
 
        switch (node_counter)
        {
            case 6:
                _elements.push_back(new MeshLib::Prism(new_elem_nodes));
                _materials.push_back(layer_id);
                break;
            case 5:
            {
                std::array<MeshLib::Node*, 5> pyramid_nodes{};
                pyramid_nodes[0] = new_elem_nodes[pyramid_base[missing_idx][0]];
                pyramid_nodes[1] = new_elem_nodes[pyramid_base[missing_idx][1]];
                pyramid_nodes[2] = new_elem_nodes[pyramid_base[missing_idx][2]];
                pyramid_nodes[3] = new_elem_nodes[pyramid_base[missing_idx][3]];
                pyramid_nodes[4] = new_elem_nodes[missing_idx];
                _elements.push_back(new MeshLib::Pyramid(pyramid_nodes));
                _materials.push_back(layer_id);
                break;
            }
            case 4:
            {
                std::array<MeshLib::Node*, 4> tet_nodes{};
                std::copy(new_elem_nodes.begin(),
                          new_elem_nodes.begin() + node_counter,
                          tet_nodes.begin());
                _elements.push_back(new MeshLib::Tet(tet_nodes));
                _materials.push_back(layer_id);
                break;
            }
            default:
                continue;
        }
    }
}

References LayeredMeshGenerator::_elements, LayeredMeshGenerator::_materials, LayeredMeshGenerator::_nodes, MeshLib::Mesh::getElements(), MeshLib::Element::getGeomType(), LayeredMeshGenerator::getNewLayerNode(), MeshLib::Mesh::getNodes(), MeshLib::Mesh::getNumberOfElements(), MeshLib::Mesh::getNumberOfNodes(), and MeshLib::TRIANGLE.

Referenced by createRasterLayers().

createRasterLayers()

bool MeshToolsLib::MeshLayerMapper::createRasterLayers ( MeshLib::Mesh const & mesh, std::vector< GeoLib::Raster const * > const & rasters, double minimum_thickness, double noDataReplacementValue = 0.0 )

overridevirtual

Based on a 2D triangle mesh this method creates a 3D mesh with a given number of prism-layers. Note: While this method would technically also work with quad meshes, this is discouraged as quad elements will most likely not be coplanar after the mapping process which result in invaled mesh elements.

Parameters

mesh	The 2D triangle mesh that is the basis for the new 3D prism mesh
rasters	Containing all the raster-data for the subsurface layers from bottom to top (starting with the bottom of the oldest layer and ending with the DEM)
minimum_thickness	Minimum thickness of each of the newly created layers (i.e. nodes with a vertical distance smaller than this will be collapsed)
noDataReplacementValue	Default z-coordinate if there are mesh nodes not located on the DEM raster (i.e. raster_paths[0])

Returns

A mesh with the requested number of layers of prism elements (also including Tet- & Pyramid-elements in case of degenerated prisms)

Implements LayeredMeshGenerator.

Definition at line 154 of file MeshLayerMapper.cpp.

{
    const std::size_t nLayers(rasters.size());
    if (nLayers < 2 || mesh.getDimension() != 2)
    {
        ERR("MeshLayerMapper::createRasterLayers(): A 2D mesh and at least two "
            "rasters required as input.");
        return false;
    }
 
    auto top = std::make_unique<MeshLib::Mesh>(mesh);
    if (!layerMapping(*top, *rasters.back(), noDataReplacementValue))
    {
        return false;
    }
 
    auto bottom = std::make_unique<MeshLib::Mesh>(mesh);
    if (!layerMapping(*bottom, *rasters[0], 0))
    {
        return false;
    }
 
    this->_minimum_thickness = minimum_thickness;
    std::size_t const nNodes = mesh.getNumberOfNodes();
    _nodes.reserve(nLayers * nNodes);
 
    // number of triangles in the original mesh
    std::size_t const nElems(std::count_if(
        mesh.getElements().begin(), mesh.getElements().end(),
        [](MeshLib::Element const* elem)
        { return (elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE); }));
    _elements.reserve(nElems * (nLayers - 1));
    _materials.reserve(nElems * (nLayers - 1));
 
    // add bottom layer
    std::vector<MeshLib::Node*> const& nodes = bottom->getNodes();
    std::transform(nodes.begin(), nodes.end(), std::back_inserter(_nodes),
                   [](auto const* node) { return new MeshLib::Node(*node); });
 
    // add the other layers
    for (std::size_t i = 0; i < nLayers - 1; ++i)
    {
        addLayerToMesh(*top, i, *rasters[i + 1]);
    }
 
    return true;
}

References LayeredMeshGenerator::_elements, LayeredMeshGenerator::_materials, LayeredMeshGenerator::_minimum_thickness, LayeredMeshGenerator::_nodes, addLayerToMesh(), ERR(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElements(), MeshLib::Mesh::getNumberOfNodes(), and layerMapping().

createStaticLayers()

MeshLib::Mesh * MeshToolsLib::MeshLayerMapper::createStaticLayers ( MeshLib::Mesh const & mesh, std::vector< float > const & layer_thickness_vector, std::string const & mesh_name = "SubsurfaceMesh" )

static

Based on a 2D triangle-or quad mesh this method creates a 3D mesh with a given number of prism- or hex-layers

Parameters

mesh	The triangle/quad mesh that is the basis for the new prism/hex mesh
layer_thickness_vector	The size of the vector equals the number of layers of prism/hex elements that will be extruded from the triangle/quad elements of the original mesh. The thickness of the $i$-th layer corresponds to the $i$ entry of the vector.
mesh_name	The name of the newly created mesh

Returns

A mesh with the requested number of layers of prism/hex elements

Definition at line 35 of file MeshLayerMapper.cpp.

{
    std::vector<float> thickness;
    for (std::size_t i = 0; i < layer_thickness_vector.size(); ++i)
    {
        if (layer_thickness_vector[i] > std::numeric_limits<float>::epsilon())
        {
            thickness.push_back(layer_thickness_vector[i]);
        }
        else
        {
            WARN("Ignoring layer {:d} with thickness {:f}.", i,
                 layer_thickness_vector[i]);
        }
    }
 
    const std::size_t nLayers(thickness.size());
    if (nLayers < 1 || mesh.getDimension() != 2)
    {
        ERR("MeshLayerMapper::createStaticLayers(): A 2D mesh with nLayers > 0 "
            "is required as input.");
        return nullptr;
    }
 
    const std::size_t nNodes = mesh.getNumberOfNodes();
    // count number of 2d elements in the original mesh
    const std::size_t nElems(
        std::count_if(mesh.getElements().begin(), mesh.getElements().end(),
                      [](MeshLib::Element const* elem)
                      { return (elem->getDimension() == 2); }));
 
    const std::size_t nOrgElems(mesh.getNumberOfElements());
    const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
    const std::vector<MeshLib::Element*>& elems = mesh.getElements();
    std::vector<MeshLib::Node*> new_nodes(nNodes + (nLayers * nNodes));
    std::vector<MeshLib::Element*> new_elems;
    new_elems.reserve(nElems * nLayers);
    MeshLib::Properties properties;
    auto* const materials = properties.createNewPropertyVector<int>(
        "MaterialIDs", MeshLib::MeshItemType::Cell);
    if (!materials)
    {
        ERR("Could not create PropertyVector object 'MaterialIDs'.");
        return nullptr;
    }
 
    auto layer_to_material_id = [&](unsigned const layer_id)
    { return ranges::views::repeat_n(nLayers - layer_id, nOrgElems); };
    materials->assign(ranges::views::iota(0u, nLayers + 1) |
                      ranges::views::transform(layer_to_material_id) |
                      ranges::views::join);
 
    double z_offset(0.0);
 
    for (unsigned layer_id = 0; layer_id <= nLayers; ++layer_id)
    {
        // add nodes for new layer
        unsigned node_offset(nNodes * layer_id);
        if (layer_id > 0)
        {
            z_offset += thickness[layer_id - 1];
        }
 
        std::transform(nodes.cbegin(), nodes.cend(),
                       new_nodes.begin() + node_offset,
                       [&z_offset](MeshLib::Node* node) {
                           return new MeshLib::Node((*node)[0], (*node)[1],
                                                    (*node)[2] - z_offset);
                       });
 
        // starting with 2nd layer create prism or hex elements connecting the
        // last layer with the current one
        if (layer_id == 0)
        {
            continue;
        }
 
        node_offset -= nNodes;
 
        for (unsigned i = 0; i < nOrgElems; ++i)
        {
            const MeshLib::Element* sfc_elem(elems[i]);
            if (sfc_elem->getDimension() < 2)
            {  // ignore line-elements
                continue;
            }
 
            const unsigned nElemNodes(sfc_elem->getNumberOfBaseNodes());
            auto** e_nodes = new MeshLib::Node*[2 * nElemNodes];
 
            for (unsigned j = 0; j < nElemNodes; ++j)
            {
                const unsigned node_id =
                    sfc_elem->getNode(j)->getID() + node_offset;
                e_nodes[j] = new_nodes[node_id + nNodes];
                e_nodes[j + nElemNodes] = new_nodes[node_id];
            }
            if (sfc_elem->getGeomType() == MeshLib::MeshElemType::TRIANGLE)
            {
                // extrude triangles to prism
                new_elems.push_back(new MeshLib::Prism(e_nodes));
            }
            else if (sfc_elem->getGeomType() == MeshLib::MeshElemType::QUAD)
            {
                // extrude quads to hexes
                new_elems.push_back(new MeshLib::Hex(e_nodes));
            }
            else
            {
                OGS_FATAL("MeshLayerMapper: Unknown element type to extrude.");
            }
        }
    }
    return new MeshLib::Mesh(mesh_name, new_nodes, new_elems,
                             true /* compute_element_neighbors */, properties);
}

References MeshLib::Cell, MeshLib::Properties::createNewPropertyVector(), ERR(), MeshLib::Element::getDimension(), MeshLib::Mesh::getDimension(), MeshLib::Mesh::getElements(), MeshLib::Element::getGeomType(), MathLib::Point3dWithID::getID(), MeshLib::Element::getNode(), MeshLib::Mesh::getNodes(), MeshLib::Element::getNumberOfBaseNodes(), MeshLib::Mesh::getNumberOfElements(), MeshLib::Mesh::getNumberOfNodes(), OGS_FATAL, MeshLib::QUAD, MeshLib::TRIANGLE, and WARN().

Referenced by MeshLayerEditDialog::createPrismMesh(), and MeshLayerEditDialog::createTetMesh().

layerMapping()

bool MeshToolsLib::MeshLayerMapper::layerMapping ( MeshLib::Mesh const & mesh, const GeoLib::Raster & raster, double nodata_replacement = 0.0, bool const ignore_nodata = false )

static

Maps the elevation of nodes of a given 2D mesh according to the raster. At locations where no information is given, node elevation is set to noDataReplacementValue.

Definition at line 294 of file MeshLayerMapper.cpp.

{
    if (mesh.getDimension() != 2)
    {
        ERR("MeshLayerMapper::layerMapping() - requires 2D mesh");
        return false;
    }
 
    GeoLib::RasterHeader const& header(raster.getHeader());
    const double x0(header.origin[0]);
    const double y0(header.origin[1]);
    const double delta(header.cell_size);
 
    const std::pair<double, double> xDim(
        x0, x0 + header.n_cols * delta);  // extension in x-dimension
    const std::pair<double, double> yDim(
        y0, y0 + header.n_rows * delta);  // extension in y-dimension
 
    const std::size_t nNodes(mesh.getNumberOfNodes());
    const std::vector<MeshLib::Node*>& nodes = mesh.getNodes();
    for (unsigned i = 0; i < nNodes; ++i)
    {
        if (!ignore_nodata && !raster.isPntOnRaster(*nodes[i]))
        {
            // use either default value or elevation from layer above
            (*nodes[i])[2] = nodata_replacement;
            continue;
        }
 
        double elevation(raster.getValueAtPoint(*nodes[i]));
        constexpr double eps = std::numeric_limits<double>::epsilon();
        if (std::abs(elevation - header.no_data) < eps)
        {
            if (ignore_nodata)
            {
                continue;
            }
            elevation = nodata_replacement;
        }
        else
        {
            elevation = raster.interpolateValueAtPoint(*nodes[i]);
        }
        (*nodes[i])[2] = elevation;
    }
 
    return true;
}

References GeoLib::RasterHeader::cell_size, ERR(), MeshLib::Mesh::getDimension(), GeoLib::Raster::getHeader(), MeshLib::Mesh::getNodes(), MeshLib::Mesh::getNumberOfNodes(), GeoLib::Raster::getValueAtPoint(), GeoLib::Raster::interpolateValueAtPoint(), GeoLib::Raster::isPntOnRaster(), GeoLib::RasterHeader::n_cols, GeoLib::RasterHeader::n_rows, GeoLib::RasterHeader::no_data, and GeoLib::RasterHeader::origin.

Referenced by LayeredVolume::createRasterLayers(), createRasterLayers(), main(), and MeshView::openMap2dMeshDialog().

mapToStaticValue()

bool MeshToolsLib::MeshLayerMapper::mapToStaticValue ( MeshLib::Mesh const & mesh, double value )

static

Maps the elevation of all mesh nodes to the specified static value.

Definition at line 346 of file MeshLayerMapper.cpp.

{
    if (mesh.getDimension() != 2)
    {
        ERR("MshLayerMapper::mapToStaticValue() - requires 2D mesh");
        return false;
    }
 
    std::vector<MeshLib::Node*> const& nodes(mesh.getNodes());
    for (MeshLib::Node* node : nodes)
    {
        (*node)[2] = value;
    }
    return true;
}

References ERR(), MeshLib::Mesh::getDimension(), and MeshLib::Mesh::getNodes().

Referenced by main(), and MeshView::openMap2dMeshDialog().

---

The documentation for this class was generated from the following files:

• MeshToolsLib/MeshGenerators/MeshLayerMapper.h
• MeshToolsLib/MeshGenerators/MeshLayerMapper.cpp