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-2025, 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 "
40 "(depth), (min = 0)",
41 false, 1000, "CELLSIZE_Z");
42 cmd.add(z_arg);
43
44 TCLAP::ValueArg<double> y_arg("y", "cellsize-y",
45 "edge length of cubes in y-direction "
46 "(latitude), (min = 0)",
47 false, 1000, "CELLSIZE_Y");
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, (min = 0)",
54 true, 1000, "CELLSIZE_X");
55 cmd.add(x_arg);
56
57 TCLAP::ValueArg<std::string> output_arg(
58 "o", "output", "Output (.vtu). The output grid file", true, "",
59 "OUTPUT_FILE");
60 cmd.add(output_arg);
61
62 TCLAP::ValueArg<std::string> input_arg(
63 "i", "input", "Input (.vtu | .msh). The 3D input mesh file", true, "",
64 "INPUT_FILE");
65 cmd.add(input_arg);
66 cmd.parse(argc, argv);
67
69
70 if ((y_arg.isSet() && !z_arg.isSet()) ||
71 ((!y_arg.isSet() && z_arg.isSet())))
72 {
73 ERR(
"For equilateral cubes, only x needs to be set. For unequal "
74 "cuboids, all three edge lengths (x/y/z) need to be specified.");
75 return -1;
76 }
78
79 double const x_size = x_arg.getValue();
80 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
81 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
82
83 if (x_size <= 0 || y_size <= 0 || z_size <= 0)
84 {
85 ERR(
"A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
86 z_size);
87 return -1;
88 }
89
90 std::array<double, 3> const cellsize = {x_size, y_size, z_size};
91
92 vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
93 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
94 reader->SetFileName(input_arg.getValue().c_str());
95 reader->Update();
96 vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
97
98 double* const bounds = mesh->GetBounds();
100 std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
102 std::array<double, 3>{bounds[1], bounds[3], bounds[5]});
103 std::array<double, 3> ranges = {max[0] - min[0], max[1] - min[1],
104 max[2] - min[2]};
105 if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
106 {
107 ERR(
"The range ({},{},{}) is not allowed to be < 0", ranges[0],
108 ranges[1], ranges[2]);
109 return -1;
110 }
111 std::array<std::size_t, 3> const dims =
112 VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
113 std::unique_ptr<MeshLib::Mesh> grid(
115 dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
116 min, "grid"));
117
118 std::vector<int> const tmp_ids =
119 VoxelGridFromMesh::assignCellIds(mesh, min, dims, cellsize);
120 auto* const cell_ids = grid->getProperties().createNewPropertyVector<int>(
122 assert(cell_ids);
123 cell_ids->assign(tmp_ids);
124
125 if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
126 {
127 return EXIT_FAILURE;
128 }
129
130 VoxelGridFromMesh::mapMeshArraysOntoGrid(mesh, grid);
131
133 {
134 return EXIT_FAILURE;
135 }
136
137 return EXIT_SUCCESS;
138}
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.12.0