OGS
createLayeredMeshFromRasters.cpp File Reference
#include <tclap/CmdLine.h>
#include <algorithm>
#include <iterator>
#include <memory>
#include <string>
#include <vector>
#include <mpi.h>
#include "BaseLib/FileTools.h"
#include "BaseLib/IO/readStringListFromFile.h"
#include "GeoLib/IO/AsciiRasterInterface.h"
#include "InfoLib/GitInfo.h"
#include "MeshLib/IO/VtkIO/VtuInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshToolsLib/MeshGenerators/MeshLayerMapper.h"
Include dependency graph for createLayeredMeshFromRasters.cpp:

Go to the source code of this file.

Functions

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

Function Documentation

◆ main()

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

Definition at line 34 of file createLayeredMeshFromRasters.cpp.

35{
36 TCLAP::CmdLine cmd(
37 "Creates a layered 3D OGS mesh from an existing 2D OGS mesh and a list "
38 "of raster files representing subsurface layers. Supported raster "
39 "formats are ArcGIS ascii rasters (*.asc), Surfer Grids (*.grd), or "
40 "gridded XYZ rasters (*.xyz)."
41 "Only input meshes consisting of line and triangle elements are "
42 "currently supported as mapping of quads might result in invalid mesh "
43 "elements.\n\n"
44 "OpenGeoSys-6 software, version " +
46 ".\n"
47 "Copyright (c) 2012-2024, OpenGeoSys Community "
48 "(http://www.opengeosys.org)",
50
51 TCLAP::SwitchArg use_ascii_arg("", "ascii_output",
52 "Write VTU output in ASCII format.");
53 cmd.add(use_ascii_arg);
54
55 double min_thickness(std::numeric_limits<double>::epsilon());
56 TCLAP::ValueArg<double> min_thickness_arg(
57 "t", "thickness",
58 "The minimum thickness of a layer to be integrated at any given "
59 "location.",
60 false, min_thickness, "floating point number");
61 cmd.add(min_thickness_arg);
62
63 TCLAP::ValueArg<std::string> raster_path_arg(
64 "r", "raster-list",
65 "An ascii-file containing a list of raster files, starting from top "
66 "(DEM) to bottom.",
67 true, "", "file name");
68 cmd.add(raster_path_arg);
69
70 TCLAP::ValueArg<std::string> mesh_out_arg(
71 "o", "output-mesh-file", "The file name of the resulting 3D mesh.",
72 true, "", "file name");
73 cmd.add(mesh_out_arg);
74
75 TCLAP::ValueArg<std::string> mesh_arg("i", "input-mesh-file",
76 "The file name of the 2D input mesh.",
77 true, "", "file name");
78 cmd.add(mesh_arg);
79
80 cmd.parse(argc, argv);
81
82#ifdef USE_PETSC
83 MPI_Init(&argc, &argv);
84#endif
85
86 if (min_thickness_arg.isSet())
87 {
88 min_thickness = min_thickness_arg.getValue();
89 if (min_thickness < 0)
90 {
91 ERR("Minimum layer thickness must be non-negative value.");
92#ifdef USE_PETSC
93 MPI_Finalize();
94#endif
95 return EXIT_FAILURE;
96 }
97 }
98
99 INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
100 std::unique_ptr<MeshLib::Mesh> const sfc_mesh(
101 MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
102 if (!sfc_mesh)
103 {
104 ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
105#ifdef USE_PETSC
106 MPI_Finalize();
107#endif
108 return EXIT_FAILURE;
109 }
110 if (sfc_mesh->getDimension() != 2)
111 {
112 ERR("Input mesh must be a 2D mesh.");
113#ifdef USE_PETSC
114 MPI_Finalize();
115#endif
116 return EXIT_FAILURE;
117 }
118 INFO("done.");
119
120 std::vector<std::string> raster_paths =
121 BaseLib::IO::readStringListFromFile(raster_path_arg.getValue());
122 if (raster_paths.size() < 2)
123 {
124 ERR("At least two raster files needed to create 3D mesh.");
125#ifdef USE_PETSC
126 MPI_Finalize();
127#endif
128 return EXIT_FAILURE;
129 }
130 std::reverse(raster_paths.begin(), raster_paths.end());
131
133 if (auto const rasters = FileIO::readRasters(raster_paths))
134 {
135 if (!mapper.createLayers(*sfc_mesh, *rasters, min_thickness))
136 {
137#ifdef USE_PETSC
138 MPI_Finalize();
139#endif
140 return EXIT_FAILURE;
141 }
142 }
143 else
144 {
145 ERR("Reading raster files.");
146#ifdef USE_PETSC
147 MPI_Finalize();
148#endif
149 return EXIT_FAILURE;
150 }
151
152 std::string output_name(mesh_out_arg.getValue());
153 if (!BaseLib::hasFileExtension(".vtu", output_name))
154 {
155 output_name.append(".vtu");
156 }
157
158 INFO("Writing mesh '{:s}' ... ", output_name);
159 auto const result_mesh = mapper.getMesh("SubsurfaceMesh");
160 if (result_mesh == nullptr)
161 {
162 ERR("Mapper returned empty result for 'SubsurfaceMesh'.");
163#ifdef USE_PETSC
164 MPI_Finalize();
165#endif
166 return EXIT_FAILURE;
167 }
168
169 auto const data_mode =
170 use_ascii_arg.getValue() ? vtkXMLWriter::Ascii : vtkXMLWriter::Binary;
171
172 MeshLib::IO::writeVtu(*result_mesh, output_name, data_mode);
173 INFO("done.");
174
175#ifdef USE_PETSC
176 MPI_Finalize();
177#endif
178 return EXIT_SUCCESS;
179}
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
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.
Manipulating and adding prism element layers to an existing 2D mesh.
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
int writeVtu(MeshLib::Mesh const &mesh, std::string const &file_name, int const data_mode)
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)

References LayeredMeshGenerator::createLayers(), ERR(), LayeredMeshGenerator::getMesh(), BaseLib::hasFileExtension(), INFO(), GitInfoLib::GitInfo::ogs_version, MeshLib::IO::readMeshFromFile(), FileIO::readRasters(), BaseLib::IO::readStringListFromFile(), and MeshLib::IO::writeVtu().