10 #include <tclap/CmdLine.h>
11 #include <vtkAbstractArray.h>
12 #include <vtkCellData.h>
13 #include <vtkCellLocator.h>
14 #include <vtkDataArray.h>
15 #include <vtkDoubleArray.h>
16 #include <vtkIntArray.h>
17 #include <vtkSmartPointer.h>
18 #include <vtkUnstructuredGrid.h>
19 #include <vtkXMLUnstructuredGridReader.h>
32 std::array<double, 3>
const& cellsize)
35 static_cast<std::size_t
>(std::ceil((max[0] - min[0]) / cellsize[0])),
36 static_cast<std::size_t
>(std::ceil((max[1] - min[1]) / cellsize[1])),
37 static_cast<std::size_t
>(std::ceil((max[2] - min[2]) / cellsize[2]))};
40 std::vector<int>
assignCellIds(vtkSmartPointer<vtkUnstructuredGrid>
const& mesh,
42 std::array<std::size_t, 3>
const& dims,
43 std::array<double, 3>
const& cellsize)
45 vtkSmartPointer<vtkCellLocator> locator =
46 vtkSmartPointer<vtkCellLocator>::New();
47 locator->SetDataSet(mesh);
50 std::vector<int> cell_ids;
51 cell_ids.reserve(dims[0] * dims[1] * dims[2]);
52 std::array<double, 3>
const grid_max = {min[0] + dims[0] * cellsize[0],
53 min[1] + dims[1] * cellsize[1],
54 min[2] + dims[2] * cellsize[2]};
55 for (
double k = min[2] + (cellsize[2] / 2.0); k < grid_max[2];
58 for (
double j = min[1] + (cellsize[1] / 2.0); j < grid_max[1];
61 for (
double i = min[0] + (cellsize[0] / 2.0); i < grid_max[0];
64 double pnt[3] = {i, j, k};
65 cell_ids.push_back(
static_cast<int>(locator->FindCell(pnt)));
73 std::unique_ptr<MeshLib::Mesh>& grid)
77 cell_id_name, 0,
static_cast<int>(mesh->GetNumberOfCells()),
true);
79 if (n_elems_marked == grid->getNumberOfElements())
81 ERR(
"No valid elements found. Aborting...");
93 template <
typename T,
typename VTK_TYPE>
101 std::size_t
const n_elems = cell_ids.
size();
103 for (std::size_t j = 0; j < n_elems; ++j)
104 arr[j] = vtk_arr->GetValue(cell_ids[j]);
108 std::unique_ptr<MeshLib::Mesh>& grid)
110 vtkSmartPointer<vtkCellData>
const cell_data = mesh->GetCellData();
111 for (
int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
113 std::string
const&
name = cell_data->GetArrayName(i);
114 vtkSmartPointer<vtkDoubleArray>
const dbl_arr =
115 dynamic_cast<vtkDoubleArray*
>(cell_data->GetArray(
name.c_str()));
118 mapArray<double, vtkSmartPointer<vtkDoubleArray>>(*grid, dbl_arr,
122 vtkSmartPointer<vtkIntArray>
const int_arr =
123 dynamic_cast<vtkIntArray*
>(cell_data->GetArray(
name.c_str()));
126 mapArray<int, vtkSmartPointer<vtkIntArray>>(*grid, int_arr,
name);
129 WARN(
"Ignoring array '{:s}', array type not implemented...",
name);
133 int main(
int argc,
char* argv[])
136 "Reads a 3D unstructured mesh and samples it onto a structured grid of "
137 "the same extent. Cell properties are mapped onto the grid (sampled at "
138 "the centre-points of each cube), node properties are ignored. Note, "
139 "that a large cube size may result in an undersampling of the original "
140 "mesh structure.\nCube sizes are defines by x/y/z-parameters. For "
141 "equilateral cubes, only the x-parameter needs to be set.\n\n"
142 "OpenGeoSys-6 software, version " +
145 "Copyright (c) 2012-2021, OpenGeoSys Community "
146 "(http://www.opengeosys.org)",
149 TCLAP::ValueArg<double> z_arg(
"z",
"cellsize-z",
150 "edge length of cubes in z-direction (depth)",
151 false, 1000,
"floating point number");
154 TCLAP::ValueArg<double> y_arg(
155 "y",
"cellsize-y",
"edge length of cubes in y-direction (latitude)",
156 false, 1000,
"floating point number");
159 TCLAP::ValueArg<double> x_arg(
161 "edge length of cubes in x-direction (longitude) or all directions, if "
162 "y and z are not set",
163 true, 1000,
"floating point number");
166 TCLAP::ValueArg<std::string> output_arg(
167 "o",
"output",
"the output grid (*.vtu)",
true,
"",
"output.vtu");
170 TCLAP::ValueArg<std::string> input_arg(
"i",
"input",
171 "the 3D input mesh (*.vtu, *.msh)",
172 true,
"",
"input.vtu");
174 cmd.parse(argc, argv);
176 if ((y_arg.isSet() && !z_arg.isSet()) ||
177 ((!y_arg.isSet() && z_arg.isSet())))
179 ERR(
"For equilateral cubes, only x needs to be set. For unequal "
180 "cuboids, all three edge lengths (x/y/z) need to be specified.");
184 double const x_size = x_arg.getValue();
185 double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
186 double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
187 std::array<double, 3>
const cellsize = {x_size, y_size, z_size};
189 vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
190 vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
191 reader->SetFileName(input_arg.getValue().c_str());
193 vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
195 double*
const bounds = mesh->GetBounds();
197 std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
199 std::array<double, 3>{bounds[1], bounds[3], bounds[5]});
200 std::array<std::size_t, 3>
const dims =
getDimensions(min, max, cellsize);
201 std::unique_ptr<MeshLib::Mesh> grid(
203 dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
206 std::vector<int>
const tmp_ids =
assignCellIds(mesh, min, dims, cellsize);
207 std::vector<int>& cell_ids =
208 *grid->getProperties().createNewPropertyVector<
int>(
210 std::copy(tmp_ids.cbegin(), tmp_ids.cend(), std::back_inserter(cell_ids));
void ERR(char const *fmt, Args const &... args)
void WARN(char const *fmt, Args const &... args)
Definition of the Mesh class.
int main(int argc, char *argv[])
std::vector< int > assignCellIds(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, MathLib::Point3d const &min, std::array< std::size_t, 3 > const &dims, std::array< double, 3 > const &cellsize)
std::array< std::size_t, 3 > getDimensions(MathLib::Point3d const &min, MathLib::Point3d const &max, std::array< double, 3 > const &cellsize)
std::string const cell_id_name
void mapArray(MeshLib::Mesh &grid, VTK_TYPE vtk_arr, std::string const &name)
bool removeUnusedGridCells(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > &grid)
void mapMeshArraysOntoGrid(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > &grid)
std::size_t searchByPropertyValueRange(std::string const &property_name, PROPERTY_TYPE const min_property_value, PROPERTY_TYPE const max_property_value, bool outside_of)
const std::vector< std::size_t > & getSearchedElementIDs() const
return marked elements
Properties & getProperties()
PropertyVector< T > const * getPropertyVector(std::string const &name) const
PropertyVector< T > * createNewPropertyVector(std::string const &name, MeshItemType mesh_item_type, std::size_t n_components=1)
GITINFOLIB_EXPORT const std::string ogs_version
void copy(PETScVector const &x, PETScVector &y)
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, [[maybe_unused]] std::set< std::string > variable_output_names)
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")
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)