VerticalSliceFromLayers.cpp File Reference

Detailed Description

Copyright

Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) Distributed under a Modified BSD License. See accompanying file LICENSE.txt or http://www.opengeosys.org/project/license

Definition in file VerticalSliceFromLayers.cpp.

#include <algorithm>
#include <cmath>
#include <memory>
#include <string>
#include <vector>
#include <tclap/CmdLine.h>
#include <QCoreApplication>
#include "Applications/FileIO/Gmsh/GMSHInterface.h"
#include "Applications/FileIO/Gmsh/GmshReader.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/IO/readStringListFromFile.h"
#include "GeoLib/AABB.h"
#include "GeoLib/AnalyticalGeometry.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
#include "GeoLib/Point.h"
#include "GeoLib/Polygon.h"
#include "GeoLib/Polyline.h"
#include "InfoLib/GitInfo.h"
#include "MathLib/MathTools.h"
#include "MathLib/Point3d.h"
#include "MeshGeoToolsLib/GeoMapper.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshEditing/MeshRevision.h"
#include "MeshLib/MeshEditing/RemoveMeshComponents.h"
#include "MeshLib/MeshSearch/ElementSearch.h"
#include "MeshLib/MeshSurfaceExtraction.h"
#include "MeshLib/Node.h"

Include dependency graph for VerticalSliceFromLayers.cpp:

Go to the source code of this file.

Functions
std::unique_ptr< std::vector< GeoLib::Point * > >	createPoints (MathLib::Point3d const &start, MathLib::Point3d const &end, std::size_t const n_intervals)
	creates a vector of sampling points based on the specified resolution More...
GeoLib::Polyline *	createPolyline (std::vector< GeoLib::Point * > const &points)
	creates a polyline to be mapped on a mesh layer More...
std::vector< std::string >	createGeometries (GeoLib::GEOObjects &geo, std::vector< std::string > const &layer_names, MathLib::Point3d const &pnt_start, MathLib::Point3d const &pnt_end, double const resolution)
	creates a mapped line for each of the mesh layers More...
void	mergeGeometries (GeoLib::GEOObjects &geo, std::vector< std::string > const &geo_names, std::string &merged_geo_name)
std::pair< Eigen::Matrix3d, double >	rotateGeometryToXY (std::vector< GeoLib::Point * > &points)
	rotates the merged geometry into the XY-plane More...
void	consolidateGeometry (GeoLib::GEOObjects &geo, std::string const &output_name, std::string &merged_geo_name, bool const keep_gml_file)
MeshLib::Mesh *	generateMesh (GeoLib::GEOObjects &geo, std::string const &geo_name, std::string const &output_name, double res)
	converts geometry into GMSH format and creates mesh More...
void	rotateMesh (MeshLib::Mesh &mesh, Eigen::Matrix3d const &rot_mat, double const z_shift)
	inverse rotation of the mesh, back into original position More...
MeshLib::Mesh *	removeLineElements (MeshLib::Mesh const &mesh)
	removes line elements from mesh such that only triangles remain More...
void	writeBoundary (MeshLib::Mesh const &mesh, std::vector< std::size_t > const &idx_array, std::string const &file_name)
void	extractBoundaries (MeshLib::Mesh const &mesh, std::string const &output_name, MathLib::Point3d const &pnt_start)
int	main (int argc, char *argv[])

Function Documentation

consolidateGeometry()

void consolidateGeometry	(	GeoLib::GEOObjects &	geo,
		std::string const &	output_name,
		std::string &	merged_geo_name,
		bool const	keep_gml_file
	)

This encapsulates a workaround: For unknown reasons, the GML->GEO converter will not work correctly when inputting the merged geometry directly. However, if the geometry is saved to a file and immedeately loaded again, everything works fine.

Definition at line 197 of file VerticalSliceFromLayers.cpp.

{
    std::string const filename(output_name + ".gml");
    GeoLib::IO::XmlGmlInterface xml(geo);
    xml.export_name = merged_geo_name;
    BaseLib::IO::writeStringToFile(xml.writeToString(), filename);
 
    geo.removePolylineVec(merged_geo_name);
    geo.removePointVec(merged_geo_name);
 
    xml.readFile(filename);
 
    if (!keep_gml_file)
    {
        BaseLib::removeFile(filename);
        BaseLib::removeFile(filename + ".md5");
    }
}

References BaseLib::IO::XMLInterface::export_name, GeoLib::IO::XmlGmlInterface::readFile(), BaseLib::removeFile(), GeoLib::GEOObjects::removePointVec(), GeoLib::GEOObjects::removePolylineVec(), BaseLib::IO::writeStringToFile(), and BaseLib::IO::Writer::writeToString().

Referenced by main().

createGeometries()

std::vector<std::string> createGeometries	(	GeoLib::GEOObjects &	geo,
		std::vector< std::string > const &	layer_names,
		MathLib::Point3d const &	pnt_start,
		MathLib::Point3d const &	pnt_end,
		double const	resolution
	)

creates a mapped line for each of the mesh layers

Definition at line 82 of file VerticalSliceFromLayers.cpp.

{
    std::vector<std::string> geo_name_list;
    std::size_t const n_layers = layer_names.size();
    for (std::size_t i = 0; i < n_layers; ++i)
    {
        std::unique_ptr<MeshLib::Mesh> const layer(
            MeshLib::IO::readMeshFromFile(layer_names[i]));
        if (layer == nullptr)
        {
            ERR("Could not read file {:s}. Skipping layer...", layer_names[i]);
            continue;
        }
        if (layer->getDimension() != 2)
        {
            ERR("Layer {:d} is not a 2D mesh. Skipping layer...", i);
            continue;
        }
 
        std::string geo_name(std::to_string(i));
        std::unique_ptr<std::vector<GeoLib::Point*>> points(
            createPoints(pnt_start, pnt_end, resolution));
        geo.addPointVec(std::move(points), geo_name);
 
        auto lines = std::make_unique<std::vector<GeoLib::Polyline*>>();
        lines->push_back(createPolyline(*geo.getPointVec(geo_name)));
        geo.addPolylineVec(std::move(lines), geo_name);
 
        MeshGeoToolsLib::GeoMapper mapper(geo, geo_name);
        mapper.mapOnMesh(layer.get());
        geo_name_list.push_back(geo_name);
    }
    return geo_name_list;
}

References GeoLib::GEOObjects::addPointVec(), GeoLib::GEOObjects::addPolylineVec(), createPoints(), createPolyline(), ERR(), GeoLib::GEOObjects::getPointVec(), MeshGeoToolsLib::GeoMapper::mapOnMesh(), and MeshLib::IO::readMeshFromFile().

Referenced by main().

createPoints()

std::unique_ptr<std::vector<GeoLib::Point*> > createPoints	(	MathLib::Point3d const &	start,
		MathLib::Point3d const &	end,
		std::size_t const	n_intervals
	)

creates a vector of sampling points based on the specified resolution

Definition at line 47 of file VerticalSliceFromLayers.cpp.

{
    auto points = std::make_unique<std::vector<GeoLib::Point*>>();
    points->push_back(new GeoLib::Point(start, 0));
    double const length = std::sqrt(MathLib::sqrDist(start, end));
    double const interval = length / n_intervals;
 
    GeoLib::Point const vec((end[0] - start[0]), (end[1] - start[1]),
                            (end[2] - start[2]));
    for (std::size_t i = 1; i < n_intervals; ++i)
    {
        points->push_back(new GeoLib::Point(
            start[0] + ((i * interval) / length * vec[0]),
            start[1] + ((i * interval) / length * vec[1]), 0, i));
    }
    points->push_back(new GeoLib::Point(end, n_intervals));
    return points;
}

References MathLib::sqrDist().

Referenced by createGeometries().

createPolyline()

GeoLib::Polyline* createPolyline

(

std::vector< GeoLib::Point * > const &

points

)

creates a polyline to be mapped on a mesh layer

Definition at line 70 of file VerticalSliceFromLayers.cpp.

{
    GeoLib::Polyline* line = new GeoLib::Polyline(points);
    std::size_t const length = points.size();
    for (std::size_t i = 0; i < length; ++i)
    {
        line->addPoint(i);
    }
    return line;
}

References GeoLib::Polyline::addPoint().

Referenced by createGeometries().

extractBoundaries()

void extractBoundaries	(	MeshLib::Mesh const &	mesh,
		std::string const &	output_name,
		MathLib::Point3d const &	pnt_start
	)

Definition at line 293 of file VerticalSliceFromLayers.cpp.

{
    double const eps = mesh.getMinEdgeLength() / 100.0;
    std::unique_ptr<MeshLib::Mesh> boundary_mesh(
        MeshLib::BoundaryExtraction::getBoundaryElementsAsMesh(
            mesh, "bulk_node_ids", "bulk_element_ids", "bulk_face_ids"));
 
    auto const& elems = boundary_mesh->getElements();
    std::vector<std::size_t> left_bound_idx, right_bound_idx, top_bound_idx,
        bottom_bound_idx;
    Eigen::Vector2d const anchor(pnt_start[0], pnt_start[1]);
    for (auto e : elems)
    {
        Eigen::Vector2d const n1((*e->getNode(0))[0], (*e->getNode(0))[1]);
        Eigen::Vector2d const n2((*e->getNode(1))[0], (*e->getNode(1))[1]);
        Eigen::Vector2d const dist1(n1 - anchor);
        Eigen::Vector2d const dist2(n2 - anchor);
        std::size_t const id = e->getID();
        // elements located at left or right side
        if ((dist1 - dist2).squaredNorm() < eps)
        {
            top_bound_idx.push_back(id);
            bottom_bound_idx.push_back(id);
 
            if (dist1.squaredNorm() < eps)
            {
                right_bound_idx.push_back(id);
            }
            else
            {
                left_bound_idx.push_back(id);
            }
            continue;
        }
        // elements located at top or bottom
        if (dist2.squaredNorm() < dist1.squaredNorm())
        {
            top_bound_idx.push_back(id);
            left_bound_idx.push_back(id);
            right_bound_idx.push_back(id);
        }
        else
        {
            bottom_bound_idx.push_back(id);
            left_bound_idx.push_back(id);
            right_bound_idx.push_back(id);
        }
    }
 
    writeBoundary(*boundary_mesh, left_bound_idx, output_name + "_left");
    writeBoundary(*boundary_mesh, right_bound_idx, output_name + "_right");
    writeBoundary(*boundary_mesh, top_bound_idx, output_name + "_top");
    writeBoundary(*boundary_mesh, bottom_bound_idx, output_name + "_bottom");
}

References MeshLib::BoundaryExtraction::getBoundaryElementsAsMesh(), MeshLib::Mesh::getMinEdgeLength(), and writeBoundary().

Referenced by main().

generateMesh()

MeshLib::Mesh* generateMesh	(	GeoLib::GEOObjects &	geo,
		std::string const &	geo_name,
		std::string const &	output_name,
		double	res
	)

converts geometry into GMSH format and creates mesh

Definition at line 220 of file VerticalSliceFromLayers.cpp.

{
    std::string const gmsh_geo_name(output_name + ".geo");
    std::vector<std::string> gmsh_geo;
    gmsh_geo.push_back(geo_name);
    FileIO::GMSH::GMSHInterface gmsh_io(
        geo, true, FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity, res, 0,
        0, gmsh_geo, false, false);
    gmsh_io.writePhysicalGroups(true);
    if (!BaseLib::IO::writeStringToFile(gmsh_io.writeToString(), gmsh_geo_name))
    {
        ERR("Writing gmsh geo file '{:s}' failed.", gmsh_geo_name);
    }
 
    std::string const gmsh_mesh_name = output_name + ".msh";
    std::string gmsh_command = "gmsh -2 -algo meshadapt " + gmsh_geo_name;
    gmsh_command += " -o " + gmsh_mesh_name + " -format msh22";
    int const return_value = std::system(gmsh_command.c_str());
    if (return_value != 0)
    {
        ERR("Execution of gmsh command returned non-zero status, %d",
            return_value);
    }
    return FileIO::GMSH::readGMSHMesh(gmsh_mesh_name);
}

References ERR(), FileIO::GMSH::FixedMeshDensity, FileIO::GMSH::readGMSHMesh(), FileIO::GMSH::GMSHInterface::writePhysicalGroups(), BaseLib::IO::writeStringToFile(), and BaseLib::IO::Writer::writeToString().

Referenced by main().

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 350 of file VerticalSliceFromLayers.cpp.

{
    QCoreApplication a(argc, argv);
    TCLAP::CmdLine cmd(
        "Creates a triangle-mesh of a vertical slice out of a list of input "
        "layer meshes. The slice is defined by a start- and end-point. In "
        "addition, the resolution for meshing the extracted slice (i.e. the "
        "maximum edge length of the domain discretisation) needs to be "
        "specified. The utility requires access to the meshing utility GMSH to "
        "work correctly.\n\n"
        "OpenGeoSys-6 software, version " +
            GitInfoLib::GitInfo::ogs_version +
            ".\n"
            "Copyright (c) 2012-2021, OpenGeoSys Community "
            "(http://www.opengeosys.org)",
        ' ', GitInfoLib::GitInfo::ogs_version);
    TCLAP::SwitchArg test_arg("t", "testdata", "keep test data", false);
    cmd.add(test_arg);
    TCLAP::SwitchArg bound_arg("b", "bounds", "save mesh boundaries", false);
    cmd.add(bound_arg);
    TCLAP::ValueArg<double> res_arg(
        "r", "resolution",
        "desired edge length of triangles in the resulting slice.", true, 0,
        "floating point number");
    cmd.add(res_arg);
    TCLAP::ValueArg<double> end_y_arg(
        "", "end-y", "y-coordinates of the end point defining the slice", true,
        0, "floating point number");
    cmd.add(end_y_arg);
    TCLAP::ValueArg<double> end_x_arg(
        "", "end-x", "x-coordinates of the end point defining the slice", true,
        0, "floating point number");
    cmd.add(end_x_arg);
    TCLAP::ValueArg<double> start_y_arg(
        "", "start-y", "y-coordinates of the start point defining the slice",
        true, 0, "floating point number");
    cmd.add(start_y_arg);
    TCLAP::ValueArg<double> start_x_arg(
        "", "start-x", "x-coordinates of the start point defining the slice",
        true, 0, "floating point number");
    cmd.add(start_x_arg);
    TCLAP::ValueArg<std::string> output_arg(
        "o", "output", "name of output mesh (*.vtu)", true, "", "string");
    cmd.add(output_arg);
    TCLAP::ValueArg<std::string> input_arg(
        "i", "input",
        "name of the input file list containing the paths the all input layers "
        "in correct order from top to bottom",
        true, "", "string");
    cmd.add(input_arg);
    cmd.parse(argc, argv);
 
    std::string const input_name = input_arg.getValue();
    std::string const output_name = output_arg.getValue();
 
    MathLib::Point3d const pnt_start{
        {start_x_arg.getValue(), start_y_arg.getValue(), 0.0}};
    MathLib::Point3d const pnt_end{
        {end_x_arg.getValue(), end_y_arg.getValue(), 0.0}};
    double const length = std::sqrt(MathLib::sqrDist(pnt_start, pnt_end));
 
    std::size_t const res = std::ceil(length / res_arg.getValue());
    double const interval_length = length / res;
 
    std::vector<std::string> const layer_names =
        BaseLib::IO::readStringListFromFile(input_name);
    if (layer_names.size() < 2)
    {
        ERR("At least two layers are required to extract a slice.");
        return EXIT_FAILURE;
    }
 
    GeoLib::GEOObjects geo;
    std::vector<std::string> const geo_name_list =
        createGeometries(geo, layer_names, pnt_start, pnt_end, res);
 
    if (geo_name_list.size() < 2)
    {
        ERR("Less than two geometries could be created from layers. Aborting "
            "extraction...");
        return EXIT_FAILURE;
    }
 
    std::string merged_geo_name = "merged_geometries";
    mergeGeometries(geo, geo_name_list, merged_geo_name);
    std::vector<GeoLib::Point*> points = *geo.getPointVec(merged_geo_name);
    auto const [rot_mat, z_shift] = rotateGeometryToXY(points);
    consolidateGeometry(geo, output_name, merged_geo_name, test_arg.getValue());
 
    std::unique_ptr<MeshLib::Mesh> mesh(
        generateMesh(geo, merged_geo_name, output_name, interval_length));
    if (!test_arg.getValue())
    {
        BaseLib::removeFile(output_name + ".geo");
        BaseLib::removeFile(output_name + ".msh");
    }
    if (mesh == nullptr)
    {
        ERR("Error generating mesh... (GMSH was unable to output mesh)");
        return EXIT_FAILURE;
    }
    rotateMesh(*mesh, rot_mat, z_shift);
    std::unique_ptr<MeshLib::Mesh> new_mesh(removeLineElements(*mesh));
    if (new_mesh == nullptr)
    {
        ERR("Error generating mesh... (GMSH created line mesh)");
        return EXIT_FAILURE;
    }
 
    // collapse all nodes that might have been created due to gmsh physically
    // separating layers
    MeshLib::MeshRevision rev(*new_mesh);
    std::unique_ptr<MeshLib::Mesh> revised_mesh(
        rev.simplifyMesh("RevisedMesh", new_mesh->getMinEdgeLength() / 100.0));
 
    if (bound_arg.getValue())
    {
        extractBoundaries(*revised_mesh, output_name, pnt_start);
    }
    MeshLib::IO::VtuInterface vtu(revised_mesh.get());
    vtu.writeToFile(output_name + ".vtu");
    return EXIT_SUCCESS;
}

References consolidateGeometry(), createGeometries(), ERR(), extractBoundaries(), generateMesh(), GeoLib::GEOObjects::getPointVec(), mergeGeometries(), GitInfoLib::GitInfo::ogs_version, BaseLib::IO::readStringListFromFile(), BaseLib::removeFile(), removeLineElements(), rotateGeometryToXY(), rotateMesh(), MeshLib::MeshRevision::simplifyMesh(), MathLib::sqrDist(), and MeshLib::IO::VtuInterface::writeToFile().

mergeGeometries()

void mergeGeometries	(	GeoLib::GEOObjects &	geo,
		std::vector< std::string > const &	geo_names,
		std::string &	merged_geo_name
	)

Merges all layer geometries into one. Each layer is specified by one polygon. This will ensure that GMSH deals with the subsequent assignment of material groups automatically.

Definition at line 125 of file VerticalSliceFromLayers.cpp.

{
    auto points = std::make_unique<std::vector<GeoLib::Point*>>();
    auto lines = std::make_unique<std::vector<GeoLib::Polyline*>>();
 
    auto layer_pnts = *geo.getPointVec(geo_names[0]);
    std::size_t const pnts_per_line = layer_pnts.size();
    std::size_t const n_layers = geo_names.size();
    std::vector<std::size_t> last_line_idx(pnts_per_line, 0);
 
    for (std::size_t i = 0; i < pnts_per_line; ++i)
    {
        std::size_t const idx = pnts_per_line - i - 1;
        points->push_back(new GeoLib::Point(*layer_pnts[i], idx));
        last_line_idx[i] = idx;
    }
    for (std::size_t j = 1; j < n_layers; ++j)
    {
        GeoLib::Polyline* line = new GeoLib::Polyline(*points);
        for (std::size_t i = 0; i < pnts_per_line; ++i)
        {
            line->addPoint(last_line_idx[i]);
        }
        layer_pnts = *geo.getPointVec(geo_names[j]);
        for (std::size_t i = 0; i < pnts_per_line; ++i)
        {
            // check if for current point the lower layer boundary is actually
            // located below the upper boundary
            std::size_t idx = last_line_idx[pnts_per_line - i - 1];
            if ((*(*points)[idx])[2] > (*layer_pnts[i])[2])
            {
                idx = points->size();
                points->push_back(new GeoLib::Point(*layer_pnts[i], idx));
                last_line_idx[pnts_per_line - i - 1] = idx;
            }
            line->addPoint(idx);
        }
        // close polygon
        line->addPoint(line->getPointID(0));
        lines->push_back(line);
    }
 
    geo.addPointVec(std::move(points), merged_geo_name);
    geo.addPolylineVec(std::move(lines), merged_geo_name);
}

References GeoLib::Polyline::addPoint(), GeoLib::GEOObjects::addPointVec(), GeoLib::GEOObjects::addPolylineVec(), GeoLib::Polyline::getPointID(), and GeoLib::GEOObjects::getPointVec().

Referenced by main().

removeLineElements()

MeshLib::Mesh* removeLineElements

(

MeshLib::Mesh const &

mesh

)

removes line elements from mesh such that only triangles remain

Definition at line 259 of file VerticalSliceFromLayers.cpp.

{
    std::vector<std::size_t> line_idx;
    std::vector<MeshLib::Element*> const& elems = mesh.getElements();
    std::for_each(elems.begin(), elems.end(),
                  [&](auto e)
                  {
                      if (e->getGeomType() == MeshLib::MeshElemType::LINE)
                      {
                          line_idx.push_back(e->getID());
                      }
                  });
    if (line_idx.size() == mesh.getNumberOfElements())
    {
        return nullptr;
    }
    return MeshLib::removeElements(mesh, line_idx, "mesh");
}

References MeshLib::Mesh::getElements().

Referenced by main().

rotateGeometryToXY()

std::pair<Eigen::Matrix3d, double> rotateGeometryToXY

(

std::vector< GeoLib::Point * > &

points

)

rotates the merged geometry into the XY-plane

Definition at line 174 of file VerticalSliceFromLayers.cpp.

{
    // compute the plane normal
    auto const [plane_normal, d] =
        GeoLib::getNewellPlane(points.begin(), points.end());
    // rotate points into x-y-plane
    Eigen::Matrix3d const rotation_matrix =
        GeoLib::computeRotationMatrixToXY(plane_normal);
    GeoLib::rotatePoints(rotation_matrix, points.begin(), points.end());
 
    GeoLib::AABB aabb(points.begin(), points.end());
    double const z_shift =
        (aabb.getMinPoint()[2] + aabb.getMaxPoint()[2]) / 2.0;
    std::for_each(points.begin(), points.end(),
                  [z_shift](GeoLib::Point* p) { (*p)[2] -= z_shift; });
    return {rotation_matrix, z_shift};
}

References GeoLib::computeRotationMatrixToXY(), GeoLib::AABB::getMaxPoint(), GeoLib::AABB::getMinPoint(), GeoLib::getNewellPlane(), MathLib::p, and GeoLib::rotatePoints().

Referenced by main().

rotateMesh()

void rotateMesh	(	MeshLib::Mesh &	mesh,
		Eigen::Matrix3d const &	rot_mat,
		double const	z_shift
	)

inverse rotation of the mesh, back into original position

Definition at line 249 of file VerticalSliceFromLayers.cpp.

{
    std::vector<MeshLib::Node*> const& nodes = mesh.getNodes();
    std::for_each(nodes.begin(), nodes.end(),
                  [z_shift](MeshLib::Node* n) { (*n)[2] += z_shift; });
    GeoLib::rotatePoints(rot_mat.transpose(), nodes.begin(), nodes.end());
}

References MeshLib::Mesh::getNodes(), and GeoLib::rotatePoints().

Referenced by main().

writeBoundary()

void writeBoundary	(	MeshLib::Mesh const &	mesh,
		std::vector< std::size_t > const &	idx_array,
		std::string const &	file_name
	)

Definition at line 278 of file VerticalSliceFromLayers.cpp.

{
    std::unique_ptr<MeshLib::Mesh> const boundary(
        MeshLib::removeElements(mesh, idx_array, "mesh"));
    if (boundary == nullptr)
    {
        ERR("Error extracting boundary '{:s}'", file_name);
        return;
    }
    MeshLib::IO::VtuInterface vtu(boundary.get());
    vtu.writeToFile(file_name + ".vtu");
}

References ERR(), MeshLib::removeElements(), and MeshLib::IO::VtuInterface::writeToFile().

Referenced by extractBoundaries().