OGS
Vtu2Grid.cpp
Go to the documentation of this file.
1
10#include <tclap/CmdLine.h>
11#include <vtkXMLUnstructuredGridReader.h>
12
13#include "InfoLib/GitInfo.h"
15#include "MeshLib/Mesh.h"
20
21#ifdef USE_PETSC
22#include <mpi.h>
23#endif
24int main(int argc, char* argv[])
25{
26 TCLAP::CmdLine cmd(
27 "Reads a 3D unstructured mesh and samples it onto a structured grid of "
28 "the same extent. Cell properties are mapped onto the grid (sampled at "
29 "the centre-points of each cube), node properties are ignored. Note, "
30 "that a large cube size may result in an undersampling of the original "
31 "mesh structure.\nCube sizes are defines by x/y/z-parameters. For "
32 "equilateral cubes, only the x-parameter needs to be set.\n\n"
33 "OpenGeoSys-6 software, version " +
35 ".\n"
36 "Copyright (c) 2012-2024, OpenGeoSys Community "
37 "(http://www.opengeosys.org)",
39
40 TCLAP::ValueArg<double> z_arg("z", "cellsize-z",
41 "edge length of cubes in z-direction (depth)",
42 false, 1000, "floating point number");
43 cmd.add(z_arg);
44
45 TCLAP::ValueArg<double> y_arg(
46 "y", "cellsize-y", "edge length of cubes in y-direction (latitude)",
47 false, 1000, "floating point number");
48 cmd.add(y_arg);
49
50 TCLAP::ValueArg<double> x_arg(
51 "x", "cellsize-x",
52 "edge length of cubes in x-direction (longitude) or all directions, if "
53 "y and z are not set",
54 true, 1000, "floating point number");
55 cmd.add(x_arg);
56
57 TCLAP::ValueArg<std::string> output_arg(
58 "o", "output", "the output grid (*.vtu)", true, "", "output.vtu");
59 cmd.add(output_arg);
60
61 TCLAP::ValueArg<std::string> input_arg("i", "input",
62 "the 3D input mesh (*.vtu, *.msh)",
63 true, "", "input.vtu");
64 cmd.add(input_arg);
65 cmd.parse(argc, argv);
66
67#ifdef USE_PETSC
68 MPI_Init(&argc, &argv);
69#endif
70
71 if ((y_arg.isSet() && !z_arg.isSet()) ||
72 ((!y_arg.isSet() && z_arg.isSet())))
73 {
74 ERR("For equilateral cubes, only x needs to be set. For unequal "
75 "cuboids, all three edge lengths (x/y/z) need to be specified.");
76#ifdef USE_PETSC
77 MPI_Finalize();
78#endif
79 return -1;
80 }
81 using namespace MeshToolsLib::MeshGenerator;
82
83 double const x_size = x_arg.getValue();
84 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
85 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
86
87 if (x_size <= 0 || y_size <= 0 || z_size <= 0)
88 {
89 ERR("A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
90 z_size);
91#ifdef USE_PETSC
92 MPI_Finalize();
93#endif
94 return -1;
95 }
96
97 std::array<double, 3> const cellsize = {x_size, y_size, z_size};
98
99 vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
100 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
101 reader->SetFileName(input_arg.getValue().c_str());
102 reader->Update();
103 vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
104
105 double* const bounds = mesh->GetBounds();
106 MathLib::Point3d const min(
107 std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
108 MathLib::Point3d const max(
109 std::array<double, 3>{bounds[1], bounds[3], bounds[5]});
110 std::array<double, 3> ranges = {max[0] - min[0], max[1] - min[1],
111 max[2] - min[2]};
112 if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
113 {
114 ERR("The range ({},{},{}) is not allowed to be < 0", ranges[0],
115 ranges[1], ranges[2]);
116#ifdef USE_PETSC
117 MPI_Finalize();
118#endif
119 return -1;
120 }
121 std::array<std::size_t, 3> const dims =
122 VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
123 std::unique_ptr<MeshLib::Mesh> grid(
125 dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
126 min, "grid"));
127
128 std::vector<int> const tmp_ids =
129 VoxelGridFromMesh::assignCellIds(mesh, min, dims, cellsize);
130 std::vector<int>& cell_ids =
131 *grid->getProperties().createNewPropertyVector<int>(
132 VoxelGridFromMesh::cell_id_name, MeshLib::MeshItemType::Cell, 1);
133 std::copy(tmp_ids.cbegin(), tmp_ids.cend(), std::back_inserter(cell_ids));
134
135 if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
136 {
137#ifdef USE_PETSC
138 MPI_Finalize();
139#endif
140 return EXIT_FAILURE;
141 }
142
143 VoxelGridFromMesh::mapMeshArraysOntoGrid(mesh, grid);
144
145 if (MeshLib::IO::writeMeshToFile(*grid, output_arg.getValue()) != 0)
146 {
147#ifdef USE_PETSC
148 MPI_Finalize();
149#endif
150 return EXIT_FAILURE;
151 }
152
153#ifdef USE_PETSC
154 MPI_Finalize();
155#endif
156 return EXIT_SUCCESS;
157}
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
Definition of the Mesh class.
int main(int argc, char *argv[])
Definition Vtu2Grid.cpp:24
GITINFOLIB_EXPORT const std::string ogs_version
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, std::set< std::string > variable_output_names)
MeshLib::Mesh * generateRegularHexMesh(const BaseLib::ISubdivision &div_x, const BaseLib::ISubdivision &div_y, const BaseLib::ISubdivision &div_z, MathLib::Point3d const &origin=MathLib::ORIGIN, std::string const &mesh_name="mesh")