From 0bd6ea3d03c6c1d8d73c88748f1263a9da3f061a Mon Sep 17 00:00:00 2001
From: Julian Heinze <julian.heinze@ufz.de>
Date: Thu, 11 May 2023 09:37:06 +0200
Subject: [PATCH] Extract functions to create voxelgrid from mesh
Allowing to make the functions from Vtu2Grid visible
for accessing them in the DataExplorer.
---
Applications/Utils/MeshEdit/Vtu2Grid.cpp | 128 +-----------------
.../MeshGenerators/VoxelGridFromMesh.cpp | 117 ++++++++++++++++
.../MeshGenerators/VoxelGridFromMesh.h | 51 +++++++
3 files changed, 175 insertions(+), 121 deletions(-)
create mode 100644 MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
create mode 100644 MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
diff --git a/Applications/Utils/MeshEdit/Vtu2Grid.cpp b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
index e25412b662a..15ef9f3ce36 100644
--- a/Applications/Utils/MeshEdit/Vtu2Grid.cpp
+++ b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
@@ -8,132 +8,17 @@
*/
#include <tclap/CmdLine.h>
-#include <vtkAbstractArray.h>
-#include <vtkCellData.h>
-#include <vtkCellLocator.h>
-#include <vtkDataArray.h>
-#include <vtkDoubleArray.h>
-#include <vtkIntArray.h>
-#include <vtkSmartPointer.h>
-#include <vtkUnstructuredGrid.h>
-#include <vtkXMLUnstructuredGridReader.h>
-
-#include "InfoLib/GitInfo.h"
+
#include "MeshLib/IO/writeMeshToFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/MeshSearch/ElementSearch.h"
#include "MeshToolsLib/MeshEditing/RemoveMeshComponents.h"
#include "MeshToolsLib/MeshGenerators/MeshGenerator.h"
+#include "MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h"
#ifdef USE_PETSC
#include <mpi.h>
#endif
-
-std::string const cell_id_name = "CellIds";
-
-std::array<std::size_t, 3> getDimensions(MathLib::Point3d const& min,
- MathLib::Point3d const& max,
- std::array<double, 3> const& cellsize)
-{
- return {
- static_cast<std::size_t>(std::ceil((max[0] - min[0]) / cellsize[0])),
- static_cast<std::size_t>(std::ceil((max[1] - min[1]) / cellsize[1])),
- static_cast<std::size_t>(std::ceil((max[2] - min[2]) / cellsize[2]))};
-}
-
-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)
-{
- vtkSmartPointer<vtkCellLocator> locator =
- vtkSmartPointer<vtkCellLocator>::New();
- locator->SetDataSet(mesh);
- locator->Update();
-
- std::vector<int> cell_ids;
- cell_ids.reserve(dims[0] * dims[1] * dims[2]);
- std::array<double, 3> const grid_max = {min[0] + dims[0] * cellsize[0],
- min[1] + dims[1] * cellsize[1],
- min[2] + dims[2] * cellsize[2]};
- for (double k = min[2] + (cellsize[2] / 2.0); k < grid_max[2];
- k += cellsize[2])
- {
- for (double j = min[1] + (cellsize[1] / 2.0); j < grid_max[1];
- j += cellsize[1])
- {
- for (double i = min[0] + (cellsize[0] / 2.0); i < grid_max[0];
- i += cellsize[0])
- {
- double pnt[3] = {i, j, k};
- cell_ids.push_back(static_cast<int>(locator->FindCell(pnt)));
- }
- }
- }
- return cell_ids;
-}
-
-bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
- std::unique_ptr<MeshLib::Mesh>& grid)
-{
- MeshLib::ElementSearch search(*grid);
- std::size_t const n_elems_marked = search.searchByPropertyValueRange<int>(
- cell_id_name, 0, static_cast<int>(mesh->GetNumberOfCells()), true);
-
- if (n_elems_marked == grid->getNumberOfElements())
- {
- ERR("No valid elements found. Aborting...");
- return false;
- }
-
- if (n_elems_marked)
- {
- grid.reset(MeshToolsLib::removeElements(
- *grid, search.getSearchedElementIDs(), "trimmed_grid"));
- }
- return true;
-}
-
-template <typename T, typename VTK_TYPE>
-void mapArray(MeshLib::Mesh& grid, VTK_TYPE vtk_arr, std::string const& name)
-{
- MeshLib::PropertyVector<int> const& cell_ids =
- *grid.getProperties().getPropertyVector<int>(
- cell_id_name, MeshLib::MeshItemType::Cell, 1);
- std::vector<T>& arr = *grid.getProperties().createNewPropertyVector<T>(
- name, MeshLib::MeshItemType::Cell, vtk_arr->GetNumberOfComponents());
- std::size_t const n_elems = cell_ids.size();
- arr.resize(n_elems);
- for (std::size_t j = 0; j < n_elems; ++j)
- arr[j] = vtk_arr->GetValue(cell_ids[j]);
-}
-
-void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
- std::unique_ptr<MeshLib::Mesh>& grid)
-{
- vtkSmartPointer<vtkCellData> const cell_data = mesh->GetCellData();
- for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
- {
- std::string const& name = cell_data->GetArrayName(i);
- vtkSmartPointer<vtkDoubleArray> const dbl_arr =
- dynamic_cast<vtkDoubleArray*>(cell_data->GetArray(name.c_str()));
- if (dbl_arr)
- {
- mapArray<double, vtkSmartPointer<vtkDoubleArray>>(*grid, dbl_arr,
- name);
- continue;
- }
- vtkSmartPointer<vtkIntArray> const int_arr =
- dynamic_cast<vtkIntArray*>(cell_data->GetArray(name.c_str()));
- if (int_arr)
- {
- mapArray<int, vtkSmartPointer<vtkIntArray>>(*grid, int_arr, name);
- continue;
- }
- WARN("Ignoring array '{:s}', array type not implemented...", name);
- }
-}
-
int main(int argc, char* argv[])
{
TCLAP::CmdLine cmd(
@@ -191,6 +76,7 @@ int main(int argc, char* argv[])
#endif
return -1;
}
+ using namespace MeshLib::MeshGenerator;
double const x_size = x_arg.getValue();
double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
@@ -208,19 +94,19 @@ int main(int argc, char* argv[])
std::array<double, 3>{bounds[0], bounds[2], bounds[4]});
MathLib::Point3d const max(
std::array<double, 3>{bounds[1], bounds[3], bounds[5]});
- std::array<std::size_t, 3> const dims = getDimensions(min, max, cellsize);
+ std::array<std::size_t, 3> const dims = VoxelGridFromMesh::getDimensions(min, max, cellsize);
std::unique_ptr<MeshLib::Mesh> grid(
MeshToolsLib::MeshGenerator::generateRegularHexMesh(
dims[0], dims[1], dims[2], cellsize[0], cellsize[1], cellsize[2],
min, "grid"));
- std::vector<int> const tmp_ids = assignCellIds(mesh, min, dims, cellsize);
+ std::vector<int> const tmp_ids = VoxelGridFromMesh::assignCellIds(mesh, min, dims, cellsize);
std::vector<int>& cell_ids =
*grid->getProperties().createNewPropertyVector<int>(
cell_id_name, MeshLib::MeshItemType::Cell, 1);
std::copy(tmp_ids.cbegin(), tmp_ids.cend(), std::back_inserter(cell_ids));
- if (!removeUnusedGridCells(mesh, grid))
+ if (!VoxelGridFromMesh::removeUnusedGridCells(mesh, grid))
{
#ifdef USE_PETSC
MPI_Finalize();
@@ -228,7 +114,7 @@ int main(int argc, char* argv[])
return EXIT_FAILURE;
}
- mapMeshArraysOntoGrid(mesh, grid);
+ VoxelGridFromMesh::mapMeshArraysOntoGrid(mesh, grid);
if (MeshLib::IO::writeMeshToFile(*grid, output_arg.getValue()) != 0)
{
diff --git a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
new file mode 100644
index 00000000000..c2c9946edc8
--- /dev/null
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
@@ -0,0 +1,117 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "VoxelGridFromMesh.h"
+
+namespace MeshLib::MeshGenerator::VoxelGridFromMesh
+{
+
+std::array<std::size_t, 3> getDimensions(MathLib::Point3d const& min,
+ MathLib::Point3d const& max,
+ std::array<double, 3> const& cellsize)
+{
+ return {
+ static_cast<std::size_t>(std::ceil((max[0] - min[0]) / cellsize[0])),
+ static_cast<std::size_t>(std::ceil((max[1] - min[1]) / cellsize[1])),
+ static_cast<std::size_t>(std::ceil((max[2] - min[2]) / cellsize[2]))};
+}
+
+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)
+{
+ vtkSmartPointer<vtkCellLocator> locator =
+ vtkSmartPointer<vtkCellLocator>::New();
+ locator->SetDataSet(mesh);
+ locator->Update();
+
+ std::vector<int> cell_ids;
+ cell_ids.reserve(dims[0] * dims[1] * dims[2]);
+ std::array<double, 3> const grid_max = {min[0] + dims[0] * cellsize[0],
+ min[1] + dims[1] * cellsize[1],
+ min[2] + dims[2] * cellsize[2]};
+ for (double k = min[2] + (cellsize[2] / 2.0); k < grid_max[2];
+ k += cellsize[2])
+ {
+ for (double j = min[1] + (cellsize[1] / 2.0); j < grid_max[1];
+ j += cellsize[1])
+ {
+ for (double i = min[0] + (cellsize[0] / 2.0); i < grid_max[0];
+ i += cellsize[0])
+ {
+ double pnt[3] = {i, j, k};
+ cell_ids.push_back(static_cast<int>(locator->FindCell(pnt)));
+ }
+ }
+ }
+ return cell_ids;
+}
+
+bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid)
+{
+ MeshLib::ElementSearch search(*grid);
+ std::size_t const n_elems_marked = search.searchByPropertyValueRange<int>(
+ cell_id_name, 0, static_cast<int>(mesh->GetNumberOfCells()), true);
+
+ if (n_elems_marked == grid->getNumberOfElements())
+ {
+ ERR("No valid elements found. Aborting...");
+ return false;
+ }
+
+ if (n_elems_marked)
+ {
+ grid.reset(MeshLib::removeElements(
+ *grid, search.getSearchedElementIDs(), "trimmed_grid"));
+ }
+ return true;
+}
+
+template <typename T, typename VTK_TYPE>
+void mapArray(MeshLib::Mesh& grid, VTK_TYPE vtk_arr, std::string const& name)
+{
+ MeshLib::PropertyVector<int> const& cell_ids =
+ *grid.getProperties().getPropertyVector<int>(
+ cell_id_name, MeshLib::MeshItemType::Cell, 1);
+ std::vector<T>& arr = *grid.getProperties().createNewPropertyVector<T>(
+ name, MeshLib::MeshItemType::Cell, vtk_arr->GetNumberOfComponents());
+ std::size_t const n_elems = cell_ids.size();
+ arr.resize(n_elems);
+ for (std::size_t j = 0; j < n_elems; ++j)
+ arr[j] = vtk_arr->GetValue(cell_ids[j]);
+}
+
+void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid)
+{
+ vtkSmartPointer<vtkCellData> const cell_data = mesh->GetCellData();
+ for (int i = 0; i < cell_data->GetNumberOfArrays(); ++i)
+ {
+ std::string const& name = cell_data->GetArrayName(i);
+ vtkSmartPointer<vtkDoubleArray> const dbl_arr =
+ dynamic_cast<vtkDoubleArray*>(cell_data->GetArray(name.c_str()));
+ if (dbl_arr)
+ {
+ mapArray<double, vtkSmartPointer<vtkDoubleArray>>(*grid, dbl_arr,
+ name);
+ continue;
+ }
+ vtkSmartPointer<vtkIntArray> const int_arr =
+ dynamic_cast<vtkIntArray*>(cell_data->GetArray(name.c_str()));
+ if (int_arr)
+ {
+ mapArray<int, vtkSmartPointer<vtkIntArray>>(*grid, int_arr, name);
+ continue;
+ }
+ WARN("Ignoring array '{:s}', array type not implemented...", name);
+ }
+}
+} // namespace MeshLib::MeshGenerator::VoxelGridFromMesh
\ No newline at end of file
diff --git a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
new file mode 100644
index 00000000000..f62acdc37f1
--- /dev/null
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
@@ -0,0 +1,51 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <vtkAbstractArray.h>
+#include <vtkCellData.h>
+#include <vtkCellLocator.h>
+#include <vtkDataArray.h>
+#include <vtkDoubleArray.h>
+#include <vtkIntArray.h>
+#include <vtkSmartPointer.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkXMLUnstructuredGridReader.h>
+
+#include "InfoLib/GitInfo.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshEditing/RemoveMeshComponents.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/MeshSearch/ElementSearch.h"
+#include "MeshLib/MeshSearch/MeshElementGrid.h"
+
+static constexpr std::string const cell_id_name = "CellIds";
+
+namespace MeshLib::MeshGenerator::VoxelGridFromMesh
+{
+
+ std::array<std::size_t, 3> getDimensions(
+ MathLib::Point3d const& min,
+ MathLib::Point3d const& max,
+ std::array<double, 3> const& cellsize);
+
+ 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);
+
+ bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid);
+
+ void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid);
+};
\ No newline at end of file
--
GitLab