24int main(
int argc,
char* argv[])
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 " +
36 "Copyright (c) 2012-2025, OpenGeoSys Community "
37 "(http://www.opengeosys.org)",
40 TCLAP::ValueArg<double> z_arg(
"z",
"cellsize-z",
41 "edge length of cubes in z-direction "
43 false, 1000,
"CELLSIZE_Z");
46 TCLAP::ValueArg<double> y_arg(
"y",
"cellsize-y",
47 "edge length of cubes in y-direction "
48 "(latitude), (min = 0)",
49 false, 1000,
"CELLSIZE_Y");
52 TCLAP::ValueArg<double> x_arg(
54 "edge length of cubes in x-direction (longitude) or all directions, if "
55 "y and z are not set, (min = 0)",
56 true, 1000,
"CELLSIZE_X");
59 TCLAP::ValueArg<std::string> output_arg(
60 "o",
"output",
"Output (.vtu). The output grid file",
true,
"",
64 TCLAP::ValueArg<std::string> input_arg(
65 "i",
"input",
"Input (.vtu | .msh). The 3D input mesh file",
true,
"",
69 cmd.add(log_level_arg);
70 cmd.parse(argc, argv);
75 if ((y_arg.isSet() && !z_arg.isSet()) ||
76 ((!y_arg.isSet() && z_arg.isSet())))
78 ERR(
"For equilateral cubes, only x needs to be set. For unequal "
79 "cuboids, all three edge lengths (x/y/z) need to be specified.");
84 double const x_size = x_arg.getValue();
85 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
86 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
88 if (x_size <= 0 || y_size <= 0 || z_size <= 0)
90 ERR(
"A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
95 std::array<double, 3>
const cellsize = {x_size, y_size, z_size};
97 vtkSmartPointer<vtkUnstructuredGrid> mesh =
99 input_arg.getValue());
105 double*
const bounds = mesh->GetBounds();
107 std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
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],
112 if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
114 ERR(
"The range ({},{},{}) is not allowed to be < 0", ranges[0],
115 ranges[1], ranges[2]);
118 std::array<std::size_t, 3>
const dims =
119 VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
120 std::unique_ptr<MeshLib::Mesh> grid(
122 dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
125 std::vector<int>
const tmp_ids =
126 VoxelGridFromMesh::assignCellIds(mesh, min, dims, cellsize);
127 auto*
const cell_ids = grid->getProperties().createNewPropertyVector<
int>(
130 cell_ids->assign(tmp_ids);
132 if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
137 VoxelGridFromMesh::mapMeshArraysOntoGrid(mesh, grid);