OGS
createLayeredMeshFromRasters.cpp
Go to the documentation of this file.
1/*
2 * \file
3 * \date 2016-02-11
4 * \brief Creates a layered mesh from a 2D mesh and a bunch of raster files.
5 *
6 * \copyright
7 * Copyright (c) 2012-2025, OpenGeoSys Community (http://www.opengeosys.org)
8 * Distributed under a Modified BSD License.
9 * See accompanying file LICENSE.txt or
10 * http://www.opengeosys.org/project/license
11 */
12
13#include <tclap/CmdLine.h>
14
15#include <algorithm>
16#include <iterator>
17#include <memory>
18#include <range/v3/algorithm/copy.hpp>
19#include <range/v3/view/join.hpp>
20#include <range/v3/view/repeat_n.hpp>
21#include <string>
22#include <vector>
23
24#include "BaseLib/FileTools.h"
26#include "BaseLib/Logging.h"
27#include "BaseLib/MPI.h"
30#include "InfoLib/GitInfo.h"
33#include "MeshLib/Mesh.h"
35
37 MeshLib::Mesh* subsurface_mesh)
38{
39 auto const* surface_material_ids = materialIDs(*sfc_mesh);
40 auto* subsurface_material_ids = materialIDs(*subsurface_mesh);
41 // reset the material ids
42 if (!surface_material_ids)
43 {
44 OGS_FATAL("Surface mesh does not contain material IDs");
45 }
46 if (surface_material_ids->empty())
47 {
48 OGS_FATAL("Surface mesh material IDs doesn't contain any values");
49 }
50 if (!subsurface_material_ids)
51 {
52 OGS_FATAL("Subsurface mesh does not contain material IDs");
53 }
54 if (subsurface_material_ids->size() % surface_material_ids->size() != 0)
55 {
56 OGS_FATAL("Could not determine the number of subsurface layers.");
57 }
58 int const number_of_layers =
59 subsurface_material_ids->size() / surface_material_ids->size();
60
61 ranges::copy(
62 ranges::views::repeat_n(*surface_material_ids, number_of_layers) |
63 ranges::views::join,
64 subsurface_material_ids->begin());
65}
66
67int main(int argc, char* argv[])
68{
69 TCLAP::CmdLine cmd(
70 "Creates a layered 3D OGS mesh from an existing 2D OGS mesh and a list "
71 "of raster files representing subsurface layers. Supported raster "
72 "formats are ArcGIS ascii rasters (*.asc), Surfer Grids (*.grd), or "
73 "gridded XYZ rasters (*.xyz)."
74 "Only input meshes consisting of line and triangle elements are "
75 "currently supported as mapping of quads might result in invalid mesh "
76 "elements.\n\n"
77 "OpenGeoSys-6 software, version " +
79 ".\n"
80 "Copyright (c) 2012-2025, OpenGeoSys Community "
81 "(http://www.opengeosys.org)",
83
84 TCLAP::SwitchArg use_ascii_arg("", "ascii_output",
85 "Write VTU output in ASCII format.");
86 cmd.add(use_ascii_arg);
87
88 double min_thickness(std::numeric_limits<double>::epsilon());
89 TCLAP::ValueArg<double> min_thickness_arg(
90 "t", "thickness",
91 "The minimum thickness of a layer to be integrated at any given "
92 "location, (min = 0)",
93 false, min_thickness, "MIN_THICKNESS");
94 cmd.add(min_thickness_arg);
95
96 TCLAP::ValueArg<std::string> raster_path_arg(
97 "r", "raster-list",
98 "Input (.vtu). An ascii-file containing a list of input "
99 "raster files, starting from"
100 "top (DEM) to bottom",
101 true, "", "INPUT_FILE_LIST");
102 cmd.add(raster_path_arg);
103
104 TCLAP::SwitchArg keep_materials_arg(
105 "", "keep-surface-material-ids",
106 "if the argument is present the materials defined in the surface mesh "
107 "are used to set the material information for the subsurface cells",
108 false);
109 cmd.add(keep_materials_arg);
110
111 TCLAP::ValueArg<std::string> mesh_out_arg(
112 "o", "output-mesh-file",
113 "Output (.vtu). The file name of the resulting 3D mesh", true, "",
114 "OUTPUT_FILE");
115 cmd.add(mesh_out_arg);
116
117 TCLAP::ValueArg<std::string> mesh_arg(
118 "i", "input-mesh-file",
119 "Input (.vtu). The file name of the 2D input mesh", true, "",
120 "INPUT_FILE");
121 cmd.add(mesh_arg);
122 auto log_level_arg = BaseLib::makeLogLevelArg();
123 cmd.add(log_level_arg);
124
125 cmd.parse(argc, argv);
126
127 BaseLib::MPI::Setup mpi_setup(argc, argv);
128 BaseLib::initOGSLogger(log_level_arg.getValue());
129
130 if (min_thickness_arg.isSet())
131 {
132 min_thickness = min_thickness_arg.getValue();
133 if (min_thickness < 0)
134 {
135 ERR("Minimum layer thickness must be non-negative value.");
136 return EXIT_FAILURE;
137 }
138 }
139
140 INFO("Reading mesh '{:s}' ... ", mesh_arg.getValue());
141 std::unique_ptr<MeshLib::Mesh> const sfc_mesh(
142 MeshLib::IO::readMeshFromFile(mesh_arg.getValue()));
143 if (!sfc_mesh)
144 {
145 ERR("Error reading mesh '{:s}'.", mesh_arg.getValue());
146 return EXIT_FAILURE;
147 }
148 if (sfc_mesh->getDimension() != 2)
149 {
150 ERR("Input mesh must be a 2D mesh.");
151 return EXIT_FAILURE;
152 }
153 INFO("done.");
154
155 std::vector<std::string> raster_paths =
156 BaseLib::IO::readStringListFromFile(raster_path_arg.getValue());
157 if (raster_paths.size() < 2)
158 {
159 ERR("At least two raster files needed to create 3D mesh.");
160 return EXIT_FAILURE;
161 }
162 std::reverse(raster_paths.begin(), raster_paths.end());
163
165 if (auto const rasters = FileIO::readRasters(raster_paths))
166 {
167 if (!mapper.createLayers(*sfc_mesh, *rasters, min_thickness))
168 {
169 return EXIT_FAILURE;
170 }
171 }
172 else
173 {
174 ERR("Reading raster files.");
175 return EXIT_FAILURE;
176 }
177
178 auto const result_mesh = mapper.getMesh("SubsurfaceMesh");
179 if (result_mesh == nullptr)
180 {
181 ERR("Mapper returned empty result for 'SubsurfaceMesh'.");
182 return EXIT_FAILURE;
183 }
184
185 if (keep_materials_arg.getValue())
186 {
188 result_mesh.get());
189 }
190
191 std::string output_name(mesh_out_arg.getValue());
192 if (!BaseLib::hasFileExtension(".vtu", output_name))
193 {
194 output_name.append(".vtu");
195 }
196
197 INFO("Writing mesh '{:s}' ... ", output_name);
198 auto const data_mode =
199 use_ascii_arg.getValue() ? vtkXMLWriter::Ascii : vtkXMLWriter::Binary;
200
201 MeshLib::IO::writeVtu(*result_mesh, output_name, data_mode);
202 INFO("done.");
203
204 return EXIT_SUCCESS;
205}
Definition of the AsciiRasterInterface class.
#define OGS_FATAL(...)
Definition Error.h:26
Filename manipulation routines.
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
Definition of the MeshLayerMapper class.
Definition of the Mesh class.
Implementation of the VtuInterface class.
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.
int main(int argc, char *argv[])
void assignSurfaceMaterialIDsToSubsurfaceLayers(MeshLib::Mesh *sfc_mesh, MeshLib::Mesh *subsurface_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.
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
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)
Definition of readMeshFromFile function.