OGS
createTetgenSmeshFromRasters.cpp
Go to the documentation of this file.
1
13#include <tclap/CmdLine.h>
14
15#include <algorithm>
16#include <iterator>
17#include <memory>
18#include <string>
19#include <vector>
20
21#ifdef USE_PETSC
22#include <mpi.h>
23#endif
24
26#include "BaseLib/FileTools.h"
29#include "InfoLib/GitInfo.h"
32#include "MeshLib/Mesh.h"
35
36int main(int argc, char* argv[])
37{
38 TCLAP::CmdLine cmd(
39 "Creates a boundary representation for a layered 3D mesh in "
40 "*.smesh-format. This boundary representation can be used with the "
41 "TetGen Tetrahedral Mesh Generator to create 3D meshes of the given "
42 "geometry at arbitrary resolutions and with varying properties. "
43 "Details on command line switches and possible parametrisation can be "
44 "found in the TetGen User's Manual. Supported raster formats are "
45 "ArcGIS ascii rasters (*.asc), Surfer Grids (*.grd) and XYZ raster "
46 "files (*.xyz)."
47 "Only input meshes consisting of line and triangle elements are "
48 "currently supported as mapping of quads might result in invalid mesh "
49 "elements.\n\n"
50 "OpenGeoSys-6 software, version " +
52 ".\n"
53 "Copyright (c) 2012-2024, OpenGeoSys Community "
54 "(http://www.opengeosys.org)",
56
57 TCLAP::SwitchArg use_ascii_arg("", "ascii_output",
58 "Write VTU output in ASCII format.");
59 cmd.add(use_ascii_arg);
60
61 double min_thickness(std::numeric_limits<double>::epsilon());
62 TCLAP::ValueArg<double> min_thickness_arg(
63 "t", "thickness",
64 "The minimum thickness of a layer to be integrated at any given "
65 "location.",
66 false, min_thickness, "floating point number");
67 cmd.add(min_thickness_arg);
68
69 TCLAP::ValueArg<std::string> raster_path_arg(
70 "r", "raster-list",
71 "An ascii-file containing a list of raster files, starting from top "
72 "(DEM) to bottom.",
73 true, "", "file name");
74 cmd.add(raster_path_arg);
75
76 TCLAP::ValueArg<std::string> mesh_out_arg(
77 "o", "output-mesh-file", "The file name of the resulting 3D mesh.",
78 true, "", "file name");
79 cmd.add(mesh_out_arg);
80
81 TCLAP::ValueArg<std::string> mesh_arg("i", "input-mesh-file",
82 "The file name of the 2D input mesh.",
83 true, "", "file name");
84 cmd.add(mesh_arg);
85
86 cmd.parse(argc, argv);
87
88#ifdef USE_PETSC
89 MPI_Init(&argc, &argv);
90#endif
91
92 if (min_thickness_arg.isSet())
93 {
94 min_thickness = min_thickness_arg.getValue();
95 if (min_thickness < 0)
96 {
97 ERR("Minimum layer thickness must be non-negative value.");
98#ifdef USE_PETSC
99 MPI_Finalize();
100#endif
101 return EXIT_FAILURE;
102 }
103 }
104
105 INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
106 std::unique_ptr<MeshLib::Mesh> const sfc_mesh(MeshLib::IO::readMeshFromFile(
107 mesh_arg.getValue(), true /* compute_element_neighbors */));
108 if (!sfc_mesh)
109 {
110 ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
111#ifdef USE_PETSC
112 MPI_Finalize();
113#endif
114 return EXIT_FAILURE;
115 }
116 if (sfc_mesh->getDimension() != 2)
117 {
118 ERR("Input mesh must be a 2D mesh.");
119#ifdef USE_PETSC
120 MPI_Finalize();
121#endif
122 return EXIT_FAILURE;
123 }
124 INFO("done.");
125
126 std::vector<std::string> raster_paths =
127 BaseLib::IO::readStringListFromFile(raster_path_arg.getValue());
128 if (raster_paths.size() < 2)
129 {
130 ERR("At least two raster files needed to create 3D mesh.");
131#ifdef USE_PETSC
132 MPI_Finalize();
133#endif
134 return EXIT_FAILURE;
135 }
136 std::reverse(raster_paths.begin(), raster_paths.end());
137
138 std::string output_name(mesh_out_arg.getValue());
139 if (!BaseLib::hasFileExtension(".smesh", output_name))
140 {
141 output_name.append(".smesh");
142 }
143
144 auto const rasters = FileIO::readRasters(raster_paths);
145 LayeredVolume lv;
146 if (rasters)
147 {
148 if (!lv.createLayers(*sfc_mesh, *rasters, min_thickness))
149 {
150 ERR("Creating the layers was erroneous.");
151#ifdef USE_PETSC
152 MPI_Finalize();
153#endif
154 return EXIT_FAILURE;
155 }
156 }
157 else
158 {
159 ERR("The raster files are not accessible.");
160#ifdef USE_PETSC
161 MPI_Finalize();
162#endif
163 return EXIT_FAILURE;
164 }
165 std::unique_ptr<MeshLib::Mesh> tg_mesh =
166 lv.getMesh("BoundaryRepresentation");
167 if (tg_mesh != nullptr)
168 {
169 std::vector<MeshLib::Node> tg_attr(lv.getAttributePoints());
170 FileIO::TetGenInterface tetgen_interface;
171 tetgen_interface.writeTetGenSmesh(output_name, *tg_mesh, tg_attr);
172 INFO("Smesh was successfully written.");
173 }
174 else
175 {
176 ERR("The tetgen-smesh could not be created.");
177#ifdef USE_PETSC
178 MPI_Finalize();
179#endif
180 return EXIT_FAILURE;
181 }
182
183#ifdef USE_PETSC
184 MPI_Finalize();
185#endif
186 return EXIT_SUCCESS;
187}
Definition of the AsciiRasterInterface class.
Filename manipulation routines.
Git information.
Definition of the LayeredVolume class.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the MeshLayerMapper class.
Definition of the Mesh class.
Definition of the TetGenInterface class.
Implementation of the VtuInterface class.
static bool writeTetGenSmesh(const std::string &file_name, const GeoLib::GEOObjects &geo_objects, const std::string &geo_name, const std::vector< GeoLib::Point > &attribute_points)
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.
Creates a volume geometry from 2D mesh layers based on raster data.
std::vector< MeshLib::Node > getAttributePoints() const
int main(int argc, char *argv[])
std::vector< std::string > readStringListFromFile(std::string const &filename)
Reads non-empty lines from a list of strings from a file into a vector.
bool hasFileExtension(std::string const &extension, std::string const &filename)
std::optional< std::vector< GeoLib::Raster const * > > readRasters(std::vector< std::string > const &raster_paths)
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
Definition of readMeshFromFile function.