23{
24 TCLAP::CmdLine cmd(
25 "Reads a 3D unstructured mesh and samples it onto a structured grid of "
26 "the same extent. Cell properties are mapped onto the grid (sampled at "
27 "the centre-points of each cube), node properties are ignored. Note, "
28 "that a large cube size may result in an undersampling of the original "
29 "mesh structure.\nCube sizes are defines by x/y/z-parameters. For "
30 "equilateral cubes, only the x-parameter needs to be set.\n\n"
31 "OpenGeoSys-6 software, version " +
33 ".\n"
34 "Copyright (c) 2012-2024, OpenGeoSys Community "
35 "(http://www.opengeosys.org)",
37
38 TCLAP::ValueArg<double> z_arg("z", "cellsize-z",
39 "edge length of cubes in z-direction (depth)",
40 false, 1000, "floating point number");
41 cmd.add(z_arg);
42
43 TCLAP::ValueArg<double> y_arg(
44 "y", "cellsize-y", "edge length of cubes in y-direction (latitude)",
45 false, 1000, "floating point number");
46 cmd.add(y_arg);
47
48 TCLAP::ValueArg<double> x_arg(
49 "x", "cellsize-x",
50 "edge length of cubes in x-direction (longitude) or all directions, if "
51 "y and z are not set",
52 true, 1000, "floating point number");
53 cmd.add(x_arg);
54
55 TCLAP::ValueArg<std::string> output_arg(
56 "o", "output", "the output grid (*.vtu)", true, "", "output.vtu");
57 cmd.add(output_arg);
58
59 TCLAP::ValueArg<std::string> input_arg("i", "input",
60 "the 3D input mesh (*.vtu, *.msh)",
61 true, "", "input.vtu");
62 cmd.add(input_arg);
63 cmd.parse(argc, argv);
64
66
67 if ((y_arg.isSet() && !z_arg.isSet()) ||
68 ((!y_arg.isSet() && z_arg.isSet())))
69 {
70 ERR(
"For equilateral cubes, only x needs to be set. For unequal "
71 "cuboids, all three edge lengths (x/y/z) need to be specified.");
72 return -1;
73 }
75
76 double const x_size = x_arg.getValue();
77 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
78 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
79
80 if (x_size <= 0 || y_size <= 0 || z_size <= 0)
81 {
82 ERR(
"A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
83 z_size);
84 return -1;
85 }
86
87 std::array<double, 3> const cellsize = {x_size, y_size, z_size};
88
89 vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
90 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
91 reader->SetFileName(input_arg.getValue().c_str());
92 reader->Update();
93 vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
94
95 double* const bounds = mesh->GetBounds();
97 std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
99 std::array<double, 3>{bounds[1], bounds[3], bounds[5]});
100 std::array<double, 3> ranges = {max[0] - min[0], max[1] - min[1],
101 max[2] - min[2]};
102 if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
103 {
104 ERR(
"The range ({},{},{}) is not allowed to be < 0", ranges[0],
105 ranges[1], ranges[2]);
106 return -1;
107 }
108 std::array<std::size_t, 3> const dims =
109 VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
110 std::unique_ptr<MeshLib::Mesh> grid(
112 dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
113 min, "grid"));
114
115 std::vector<int> const tmp_ids =
116 VoxelGridFromMesh::assignCellIds(mesh, min, dims, cellsize);
117 std::vector<int>& cell_ids =
118 *grid->getProperties().createNewPropertyVector<int>(
120 std::copy(tmp_ids.cbegin(), tmp_ids.cend(), std::back_inserter(cell_ids));
121
122 if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
123 {
124 return EXIT_FAILURE;
125 }
126
127 VoxelGridFromMesh::mapMeshArraysOntoGrid(mesh, grid);
128
130 {
131 return EXIT_FAILURE;
132 }
133
134 return EXIT_SUCCESS;
135}
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
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)
1.11.0