Loading [MathJax]/extensions/tex2jax.js
OGS
VoxelGridFromMesh.cpp
Go to the documentation of this file.
1
9#include "VoxelGridFromMesh.h"
10
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 <vtkLongLongArray.h>
18#include <vtkXMLUnstructuredGridReader.h>
19
20#include <cassert>
21#include <cmath>
22
23#include "InfoLib/GitInfo.h"
24#include "MathLib/Point3d.h"
29
31{
32std::array<std::size_t, 3> getNumberOfVoxelPerDimension(
33 std::array<double, 3> const& ranges, std::array<double, 3> const& cellsize)
34{
35 if (cellsize[0] <= 0 || cellsize[1] <= 0 || cellsize[2] <= 0)
36 {
37 OGS_FATAL("A cellsize ({},{},{}) is not allowed to be <= 0",
38 cellsize[0], cellsize[1], cellsize[2]);
39 }
40 std::array<double, 3> numberOfVoxel = {ranges[0] / cellsize[0],
41 ranges[1] / cellsize[1],
42 ranges[2] / cellsize[2]};
43
44 if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
45 {
47 "The difference of max-min ({},{},{}) is not allowed to be < 0",
48 ranges[0], ranges[1], ranges[2]);
49 }
50 std::replace(numberOfVoxel.begin(), numberOfVoxel.end(), 0, 1);
51
52 return {static_cast<std::size_t>(std::lround(numberOfVoxel[0])),
53 static_cast<std::size_t>(std::lround(numberOfVoxel[1])),
54 static_cast<std::size_t>(std::lround(numberOfVoxel[2]))};
55}
56
57std::vector<int> assignCellIds(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
58 MathLib::Point3d const& min,
59 std::array<std::size_t, 3> const& dims,
60 std::array<double, 3> const& cellsize)
61{
62 vtkSmartPointer<vtkCellLocator> locator =
63 vtkSmartPointer<vtkCellLocator>::New();
64 locator->SetDataSet(mesh);
65 locator->Update();
66
67 std::vector<int> cell_ids;
68 cell_ids.reserve(dims[0] * dims[1] * dims[2]);
69 std::array<double, 3> const grid_max = {min[0] + dims[0] * cellsize[0],
70 min[1] + dims[1] * cellsize[1],
71 min[2] + dims[2] * cellsize[2]};
72
73 double const start[3] = {min[0] + cellsize[0] / 2.,
74 min[1] + cellsize[1] / 2.,
75 min[2] + cellsize[2] / 2.};
76 double pnt[3];
77 for (pnt[2] = start[2]; pnt[2] < grid_max[2]; pnt[2] += cellsize[2])
78 {
79 for (pnt[1] = start[1]; pnt[1] < grid_max[1]; pnt[1] += cellsize[1])
80 {
81 for (pnt[0] = start[0]; pnt[0] < grid_max[0]; pnt[0] += cellsize[0])
82 {
83 cell_ids.push_back(static_cast<int>(locator->FindCell(pnt)));
84 }
85 }
86 }
87 return cell_ids;
88}
89
90bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
91 std::unique_ptr<MeshLib::Mesh>& grid)
92{
93 MeshLib::ElementSearch search(*grid);
94 std::size_t const n_elems_marked = search.searchByPropertyValueRange<int>(
95 cell_id_name, 0, static_cast<int>(mesh->GetNumberOfCells()), true);
96
97 if (n_elems_marked == grid->getNumberOfElements())
98 {
99 ERR("No valid elements found. Aborting...");
100 return false;
101 }
102
103 if (n_elems_marked)
104 {
106 *grid, search.getSearchedElementIDs(), "trimmed_grid"));
107 }
108 return true;
109}
110
111// PropertyVector is required to contain parameter cell_id_name
112template <typename T, typename VTK_TYPE>
113void mapArray(MeshLib::Mesh& grid, VTK_TYPE vtk_arr, std::string const& name)
114{
115 auto const* cell_ids = grid.getProperties().getPropertyVector<int>(
117 if (cell_ids == nullptr)
118 {
119 // Error message
120 return;
121 }
122 std::size_t const n_elems = cell_ids->size();
123 auto& arr = *grid.getProperties().createNewPropertyVector<T>(
124 name, MeshLib::MeshItemType::Cell, n_elems,
125 vtk_arr->GetNumberOfComponents());
126 for (std::size_t j = 0; j < n_elems; ++j)
127 {
128 arr[j] = vtk_arr->GetValue((*cell_ids)[j]);
129 }
130}
131
132template <typename T>
133// check whether dynamic_cast of cell data is possible the type of cell data to
134// map them on voxelgrid
136 vtkSmartPointer<vtkCellData> const cell_data,
137 char const* const name)
138{
139 using DataArrayType = vtkAOSDataArrayTemplate<T>;
140 vtkSmartPointer<DataArrayType> const arr =
141 dynamic_cast<DataArrayType*>(cell_data->GetArray(name));
142 if (!arr)
143 {
144 return false;
145 }
147 return true;
148}
149
150void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
151 std::unique_ptr<MeshLib::Mesh> const& grid)
152{
153 assert(mesh != nullptr);
154 assert(grid != nullptr);
155 vtkSmartPointer<vtkCellData> const cell_data = mesh->GetCellData();
156 for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
157 {
158 auto const name = cell_data->GetArrayName(i);
159
160 if (!(checkDyncast<double>(*grid, cell_data, name) ||
161 checkDyncast<long>(*grid, cell_data, name) ||
162 checkDyncast<long long>(*grid, cell_data, name) ||
163 checkDyncast<int>(*grid, cell_data, name)))
164 {
165 WARN("Ignoring array '{:s}', array type {:s} not implemented...",
166 name,
167 cell_data->GetArray(name)->GetDataTypeAsString());
168 }
169 }
170}
171} // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
#define OGS_FATAL(...)
Definition Error.h:26
Git information.
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Definition of the Point3d class.
Element search class.
const std::vector< std::size_t > & getSearchedElementIDs() const
return marked elements
std::size_t searchByPropertyValueRange(std::string const &property_name, PROPERTY_TYPE const min_property_value, PROPERTY_TYPE const max_property_value, bool outside_of)
Properties & getProperties()
Definition Mesh.h:136
PropertyVector< T > * createNewPropertyVector(std::string_view name, MeshItemType mesh_item_type, std::size_t n_components=1)
Definition Properties.h:17
PropertyVector< T > const * getPropertyVector(std::string_view name) const
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)
void mapArray(MeshLib::Mesh &grid, VTK_TYPE vtk_arr, std::string const &name)
void mapMeshArraysOntoGrid(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > const &grid)
bool removeUnusedGridCells(vtkSmartPointer< vtkUnstructuredGrid > const &mesh, std::unique_ptr< MeshLib::Mesh > &grid)
bool checkDyncast(MeshLib::Mesh &mesh, vtkSmartPointer< vtkCellData > const cell_data, char const *const name)
std::array< std::size_t, 3 > getNumberOfVoxelPerDimension(std::array< double, 3 > const &ranges, std::array< double, 3 > const &cellsize)
MeshLib::Mesh * removeElements(const MeshLib::Mesh &mesh, const std::vector< std::size_t > &removed_element_ids, const std::string &new_mesh_name)