OGS
Layers2Grid.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11
12#include <algorithm>
13#include <memory>
14#include <string>
15#include <vector>
16
18#include "BaseLib/Logging.h"
19#include "BaseLib/MPI.h"
21#include "GeoLib/AABB.h"
22#include "InfoLib/GitInfo.h"
23#include "MathLib/Point3d.h"
28
29int main(int argc, char* argv[])
30{
31 TCLAP::CmdLine cmd(
32 "Reads a list of 2D unstructured mesh layers and samples them onto a "
33 "structured grid of the same extent. Note, that a large cube size may "
34 "result in an undersampling of the original structure.\nCube sizes are "
35 "defines by x/y/z-parameters. For equilateral cubes, only the "
36 "x-parameter needs to be set.\n\n"
37 "OpenGeoSys-6 software, version " +
39 ".\n"
40 "Copyright (c) 2012-2025, OpenGeoSys Community "
41 "(http://www.opengeosys.org)",
43 TCLAP::SwitchArg dilate_arg(
44 "d", "dilate",
45 "assign mat IDs based on single nodes instead of a majority of nodes, "
46 "which can result in a slightly increased voxel grid extent",
47 false);
48 cmd.add(dilate_arg);
49
50 TCLAP::ValueArg<double> z_arg(
51 "z", "cellsize-z",
52 "edge length of cubes in z-direction (depth), "
53 "(min = 0)",
54 false, 1000, "CELLSIZE-Z");
55 cmd.add(z_arg);
56
57 TCLAP::ValueArg<double> y_arg(
58 "y", "cellsize-y",
59 "edge length of cubes in y-direction (latitude), "
60 "(min = 0)",
61 false, 1000, "CELLSIZE-Y");
62 cmd.add(y_arg);
63
64 TCLAP::ValueArg<double> x_arg(
65 "x", "cellsize-x",
66 "edge length of cubes in x-direction (longitude) or all directions, if "
67 "y and z are not set, (min = 0)",
68 true, 1000, "CELLSIZE-X");
69 cmd.add(x_arg);
70
71 TCLAP::ValueArg<std::string> output_arg(
72 "o", "output", "Output (.vtu). Name of output mesh file", true, "",
73 "OUTPUT_FILE");
74 cmd.add(output_arg);
75
76 TCLAP::ValueArg<std::string> input_arg(
77 "i", "input",
78 "Input (.vtu). Name of the input file list containing the paths the "
79 "all input layers "
80 "in correct order from top to bottom",
81 true, "", "INPUT_FILE_LIST");
82 cmd.add(input_arg);
83
84 auto log_level_arg = BaseLib::makeLogLevelArg();
85 cmd.add(log_level_arg);
86
87 cmd.parse(argc, argv);
88
89 BaseLib::MPI::Setup mpi_setup(argc, argv);
90 BaseLib::initOGSLogger(log_level_arg.getValue());
91
92 if ((y_arg.isSet() && !z_arg.isSet()) ||
93 ((!y_arg.isSet() && z_arg.isSet())))
94 {
95 ERR("For equilateral cubes, only x needs to be set. For unequal "
96 "cuboids, all three edge lengths (x/y/z) need to be specified.");
97 return EXIT_FAILURE;
98 }
99
100 double const x_size = x_arg.getValue();
101 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
102 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
103 std::array<double, 3> const cellsize = {x_size, y_size, z_size};
104
105 std::string const input_name = input_arg.getValue();
106 std::string const output_name = output_arg.getValue();
107 auto const layer_names = BaseLib::IO::readStringListFromFile(input_name);
108 if (layer_names.size() < 2)
109 {
110 ERR("At least two layers are required to create a 3D Mesh");
111 return EXIT_FAILURE;
112 }
113
114 std::vector<std::unique_ptr<MeshLib::Mesh>> layers;
115 layers.reserve(layer_names.size());
116 constexpr double minval = std::numeric_limits<double>::max();
117 constexpr double maxval = std::numeric_limits<double>::lowest();
118 std::pair<MathLib::Point3d, MathLib::Point3d> extent(
119 MathLib::Point3d{{minval, minval, minval}},
120 MathLib::Point3d{{maxval, maxval, maxval}});
121
122 for (auto const& layer : layer_names)
123 {
124 auto mesh(MeshLib::IO::readMeshFromFile(layer));
125 if (mesh == nullptr)
126 {
127 ERR("Input layer '{:s}' not found. Aborting...", layer);
128 return EXIT_FAILURE;
129 }
130 layers.emplace_back(mesh);
131 }
132 std::vector<MeshLib::Mesh const*> layers_ptr;
133 std::transform(std::begin(layers), std::end(layers),
134 std::back_inserter(layers_ptr),
135 [](auto const& layer) { return layer.get(); });
136 bool const dilate = dilate_arg.getValue();
138 createVoxelFromLayeredMesh(extent, layers_ptr, cellsize, dilate);
139 if (mesh == nullptr)
140 {
141 ERR("The VoxelGrid could not be created.");
142 return EXIT_FAILURE;
143 }
144 MeshLib::IO::VtuInterface vtu(mesh.get());
145 vtu.writeToFile(output_name);
146 return EXIT_SUCCESS;
147}
Definition of the AABB class.
Definition of the Element class.
Git information.
int main(int argc, char *argv[])
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
Definition of the Point3d class.
Implementation of the VtuInterface class.
Reads and writes VtkXMLUnstructuredGrid-files (vtu) to and from OGS data structures....
bool writeToFile(std::filesystem::path const &file_path)
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
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
std::unique_ptr< MeshLib::Mesh > createVoxelFromLayeredMesh(std::pair< MathLib::Point3d, MathLib::Point3d > &extent, std::vector< MeshLib::Mesh const * > const &layers, std::array< double, 3 > const cellsize, bool const dilate)
Definition of readMeshFromFile function.