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