OGS
Vtu2Grid.cpp
Go to the documentation of this file.
1 
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>
20 
21 #include "InfoLib/GitInfo.h"
23 #include "MeshLib/Mesh.h"
27 
28 std::string const cell_id_name = "CellIds";
29 
30 std::array<std::size_t, 3> getDimensions(MathLib::Point3d const& min,
31  MathLib::Point3d const& max,
32  std::array<double, 3> const& cellsize)
33 {
34  return {
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]))};
38 }
39 
40 std::vector<int> assignCellIds(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
41  MathLib::Point3d const& min,
42  std::array<std::size_t, 3> const& dims,
43  std::array<double, 3> const& cellsize)
44 {
45  vtkSmartPointer<vtkCellLocator> locator =
46  vtkSmartPointer<vtkCellLocator>::New();
47  locator->SetDataSet(mesh);
48  locator->Update();
49 
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];
56  k += cellsize[2])
57  {
58  for (double j = min[1] + (cellsize[1] / 2.0); j < grid_max[1];
59  j += cellsize[1])
60  {
61  for (double i = min[0] + (cellsize[0] / 2.0); i < grid_max[0];
62  i += cellsize[0])
63  {
64  double pnt[3] = {i, j, k};
65  cell_ids.push_back(static_cast<int>(locator->FindCell(pnt)));
66  }
67  }
68  }
69  return cell_ids;
70 }
71 
72 bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
73  std::unique_ptr<MeshLib::Mesh>& grid)
74 {
75  MeshLib::ElementSearch search(*grid);
76  std::size_t const n_elems_marked = search.searchByPropertyValueRange<int>(
77  cell_id_name, 0, static_cast<int>(mesh->GetNumberOfCells()), true);
78 
79  if (n_elems_marked == grid->getNumberOfElements())
80  {
81  ERR("No valid elements found. Aborting...");
82  return false;
83  }
84 
85  if (n_elems_marked)
86  {
87  grid.reset(MeshLib::removeElements(
88  *grid, search.getSearchedElementIDs(), "trimmed_grid"));
89  }
90  return true;
91 }
92 
93 template <typename T, typename VTK_TYPE>
94 void mapArray(MeshLib::Mesh& grid, VTK_TYPE vtk_arr, std::string const& name)
95 {
96  MeshLib::PropertyVector<int> const& cell_ids =
97  *grid.getProperties().getPropertyVector<int>(
99  std::vector<T>& arr = *grid.getProperties().createNewPropertyVector<T>(
100  name, MeshLib::MeshItemType::Cell, vtk_arr->GetNumberOfComponents());
101  std::size_t const n_elems = cell_ids.size();
102  arr.resize(n_elems);
103  for (std::size_t j = 0; j < n_elems; ++j)
104  arr[j] = vtk_arr->GetValue(cell_ids[j]);
105 }
106 
107 void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
108  std::unique_ptr<MeshLib::Mesh>& grid)
109 {
110  vtkSmartPointer<vtkCellData> const cell_data = mesh->GetCellData();
111  for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
112  {
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()));
116  if (dbl_arr)
117  {
118  mapArray<double, vtkSmartPointer<vtkDoubleArray>>(*grid, dbl_arr,
119  name);
120  continue;
121  }
122  vtkSmartPointer<vtkIntArray> const int_arr =
123  dynamic_cast<vtkIntArray*>(cell_data->GetArray(name.c_str()));
124  if (int_arr)
125  {
126  mapArray<int, vtkSmartPointer<vtkIntArray>>(*grid, int_arr, name);
127  continue;
128  }
129  WARN("Ignoring array '{:s}', array type not implemented...", name);
130  }
131 }
132 
133 int main(int argc, char* argv[])
134 {
135  TCLAP::CmdLine cmd(
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 " +
144  ".\n"
145  "Copyright (c) 2012-2021, OpenGeoSys Community "
146  "(http://www.opengeosys.org)",
148 
149  TCLAP::ValueArg<double> z_arg("z", "cellsize-z",
150  "edge length of cubes in z-direction (depth)",
151  false, 1000, "floating point number");
152  cmd.add(z_arg);
153 
154  TCLAP::ValueArg<double> y_arg(
155  "y", "cellsize-y", "edge length of cubes in y-direction (latitude)",
156  false, 1000, "floating point number");
157  cmd.add(y_arg);
158 
159  TCLAP::ValueArg<double> x_arg(
160  "x", "cellsize-x",
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");
164  cmd.add(x_arg);
165 
166  TCLAP::ValueArg<std::string> output_arg(
167  "o", "output", "the output grid (*.vtu)", true, "", "output.vtu");
168  cmd.add(output_arg);
169 
170  TCLAP::ValueArg<std::string> input_arg("i", "input",
171  "the 3D input mesh (*.vtu, *.msh)",
172  true, "", "input.vtu");
173  cmd.add(input_arg);
174  cmd.parse(argc, argv);
175 
176  if ((y_arg.isSet() && !z_arg.isSet()) ||
177  ((!y_arg.isSet() && z_arg.isSet())))
178  {
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.");
181  return -1;
182  }
183 
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};
188 
189  vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
190  vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
191  reader->SetFileName(input_arg.getValue().c_str());
192  reader->Update();
193  vtkSmartPointer<vtkUnstructuredGrid> mesh = reader->GetOutput();
194 
195  double* const bounds = mesh->GetBounds();
196  MathLib::Point3d const min(
197  std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
198  MathLib::Point3d const max(
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],
204  min, "grid"));
205 
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));
211 
212  if (!removeUnusedGridCells(mesh, grid))
213  return EXIT_FAILURE;
214 
215  mapMeshArraysOntoGrid(mesh, grid);
216 
217  if (MeshLib::IO::writeMeshToFile(*grid, output_arg.getValue()) != 0)
218  return EXIT_FAILURE;
219 
220  return EXIT_SUCCESS;
221 }
Git information.
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of the Mesh class.
int main(int argc, char *argv[])
Definition: Vtu2Grid.cpp:133
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)
Definition: Vtu2Grid.cpp:40
std::array< std::size_t, 3 > getDimensions(MathLib::Point3d const &min, MathLib::Point3d const &max, std::array< double, 3 > const &cellsize)
Definition: Vtu2Grid.cpp:30
std::string const cell_id_name
Definition: Vtu2Grid.cpp:28
void mapArray(MeshLib::Mesh &grid, VTK_TYPE vtk_arr, std::string const &name)
Definition: Vtu2Grid.cpp:94
bool removeUnusedGridCells(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > &grid)
Definition: Vtu2Grid.cpp:72
void mapMeshArraysOntoGrid(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > &grid)
Definition: Vtu2Grid.cpp:107
Element search class.
Definition: ElementSearch.h:28
std::size_t searchByPropertyValueRange(std::string const &property_name, PROPERTY_TYPE const min_property_value, PROPERTY_TYPE const max_property_value, bool outside_of)
Definition: ElementSearch.h:70
const std::vector< std::size_t > & getSearchedElementIDs() const
return marked elements
Definition: ElementSearch.h:33
Properties & getProperties()
Definition: Mesh.h:123
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)
std::size_t size() const
GITINFOLIB_EXPORT const std::string ogs_version
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
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)