diff --git a/Applications/Utils/MeshEdit/Vtu2Grid.cpp b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
index e25412b662aaa3756cc2ebc9c0a0826ea12778bc..8ee5777bfa76e358e2b992354ee87b124c69da78 100644
--- a/Applications/Utils/MeshEdit/Vtu2Grid.cpp
+++ b/Applications/Utils/MeshEdit/Vtu2Grid.cpp
@@ -8,14 +8,6 @@
*/
#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"
@@ -24,116 +16,11 @@
#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,10 +78,22 @@ int main(int argc, char* argv[])
#endif
return -1;
}
+ using namespace MeshToolsLib::MeshGenerator;
double const x_size = x_arg.getValue();
double const y_size = (y_arg.isSet()) ? y_arg.getValue() : x_arg.getValue();
double const z_size = (z_arg.isSet()) ? z_arg.getValue() : x_arg.getValue();
+
+ if (x_size <= 0 || y_size <= 0 || z_size <= 0)
+ {
+ ERR("A cellsize ({},{},{}) is not allowed to be <= 0", x_size, y_size,
+ z_size);
+#ifdef USE_PETSC
+ MPI_Finalize();
+#endif
+ return -1;
+ }
+
std::array<double, 3> const cellsize = {x_size, y_size, z_size};
vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
@@ -208,19 +107,32 @@ 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<double, 3> ranges = {max[0] - min[0], max[1] - min[1],
+ max[2] - min[2]};
+ if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
+ {
+ ERR("The range ({},{},{}) is not allowed to be < 0",
+ ranges[0], ranges[1], ranges[2]);
+#ifdef USE_PETSC
+ MPI_Finalize();
+#endif
+ return -1;
+ }
+ std::array<std::size_t, 3> const dims =
+ VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, 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);
+ VoxelGridFromMesh::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 +140,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 0000000000000000000000000000000000000000..a8d0284f27bf789bc2b3066cdaeb454a05c85624
--- /dev/null
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
@@ -0,0 +1,154 @@
+/**
+ * \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"
+
+#include <vtkAbstractArray.h>
+#include <vtkCellData.h>
+#include <vtkCellLocator.h>
+#include <vtkDataArray.h>
+#include <vtkDoubleArray.h>
+#include <vtkIntArray.h>
+#include <vtkXMLUnstructuredGridReader.h>
+
+#include <cassert>
+#include <cmath>
+
+#include "InfoLib/GitInfo.h"
+#include "MathLib/Point3d.h"
+#include "MeshLib/MeshSearch/ElementSearch.h"
+#include "MeshLib/MeshSearch/MeshElementGrid.h"
+#include "MeshToolsLib/MeshEditing/RemoveMeshComponents.h"
+#include "MeshToolsLib/MeshGenerators/MeshGenerator.h"
+
+namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
+{
+std::array<std::size_t, 3> getNumberOfVoxelPerDimension(
+ std::array<double, 3> const& ranges, std::array<double, 3> const& cellsize)
+{
+ if (cellsize[0] <= 0 || cellsize[1] <= 0 || cellsize[2] <= 0)
+ {
+ OGS_FATAL("A cellsize ({},{},{}) is not allowed to be <= 0",
+ cellsize[0], cellsize[1], cellsize[2]);
+ }
+ std::array<double, 3> numberOfVoxel = {ranges[0] / cellsize[0],
+ ranges[1] / cellsize[1],
+ ranges[2] / cellsize[2]};
+
+ if (ranges[0] < 0 || ranges[1] < 0 || ranges[2] < 0)
+ {
+ OGS_FATAL(
+ "The difference of max-min ({},{},{}) is not allowed to be < 0",
+ ranges[0], ranges[1], ranges[2]);
+ }
+ std::replace(numberOfVoxel.begin(), numberOfVoxel.end(), 0, 1);
+
+ return {static_cast<std::size_t>(std::lround(numberOfVoxel[0])),
+ static_cast<std::size_t>(std::lround(numberOfVoxel[1])),
+ static_cast<std::size_t>(std::lround(numberOfVoxel[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]};
+
+ double const start[3] = {min[0] + cellsize[0] / 2.,
+ min[1] + cellsize[1] / 2.,
+ min[2] + cellsize[2] / 2.};
+ double pnt[3];
+ for (pnt[2] = start[2]; pnt[2] < grid_max[2]; pnt[2] += cellsize[2])
+ {
+ for (pnt[1] = start[1]; pnt[1] < grid_max[1]; pnt[1] += cellsize[1])
+ {
+ for (pnt[0] = start[0]; pnt[0] < grid_max[0]; pnt[0] += cellsize[0])
+ {
+ 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;
+}
+
+// PropertyVector is required to contain parameter cell_id_name
+template <typename T, typename VTK_TYPE>
+void mapArray(MeshLib::Mesh& grid, VTK_TYPE vtk_arr, std::string const& name)
+{
+ auto const* cell_ids = grid.getProperties().getPropertyVector<int>(
+ cell_id_name, MeshLib::MeshItemType::Cell, 1);
+ if (cell_ids == nullptr)
+ {
+ // Error message
+ return;
+ }
+ auto& 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 MeshToolsLib::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 0000000000000000000000000000000000000000..aef731e64cd1886a5b41dfe92fb367be59c072f2
--- /dev/null
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
@@ -0,0 +1,52 @@
+/**
+ * \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 <vtkSmartPointer.h>
+#include <vtkUnstructuredGrid.h>
+
+#include <array>
+#include <memory>
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace MathLib
+{
+class Point3d;
+}
+namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
+{
+static std::string const cell_id_name = "CellIds";
+
+/// getNumberOfVoxelPerDimension is used to calculate how many voxel fit into
+/// a bounding box. For this calculation the difference of min and max point of
+/// the bounding box is divided by the cell size, for every dimension. The
+/// calculation is restricted to work only with positive values for the cell
+/// size. If the difference between min and max is zero, we assign one voxel for
+/// the respective dimension.
+std::array<std::size_t, 3> getNumberOfVoxelPerDimension(
+ std::array<double, 3> const& ranges, 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);
+
+// grid has to contain a PropertyVector with the name 'CellIds'
+bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid);
+
+void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
+ std::unique_ptr<MeshLib::Mesh>& grid);
+}; // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
\ No newline at end of file
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 3dbcdfaa98e98071aa4da060aa09fea772b4de10..374245558425d119dd84f584ccd31dbbcdeeda41 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -24,6 +24,7 @@ append_source_files(TEST_SOURCES GeoLib/IO)
append_source_files(TEST_SOURCES MaterialLib)
append_source_files(TEST_SOURCES MathLib)
append_source_files(TEST_SOURCES MeshLib)
+append_source_files(TEST_SOURCES MeshToolsLib)
append_source_files(TEST_SOURCES MeshGeoToolsLib)
append_source_files(TEST_SOURCES_NUMLIB NumLib)
# Disable Unity build for some files
@@ -93,6 +94,7 @@ target_link_libraries(
MaterialLib
MathLib
MeshLib
+ MeshToolsLib
NumLib
ParameterLib
$<$<NOT:$<BOOL:${OGS_BUILD_WHEEL}>>:ProcessLib>
diff --git a/Tests/MeshToolsLib/TestGetDimensions.cpp b/Tests/MeshToolsLib/TestGetDimensions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3d49aeace8cc600195caad18afa9f35e24d3b6dd
--- /dev/null
+++ b/Tests/MeshToolsLib/TestGetDimensions.cpp
@@ -0,0 +1,50 @@
+/**
+ * \file
+ * \date 2023-05-31
+ * \brief Tests for getNumberOfVoxelPerDimension()
+ *
+ * \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 <gtest/gtest.h>
+
+#include <array>
+
+#include "MathLib/Point3d.h"
+#include "MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h"
+
+TEST(MeshToolsLib, getNumberOfVoxelPerDimension)
+{
+ std::array<double, 3> const cellsize = {1, 2, 5};
+ std::array<double, 3> const ranges = {10, 10, 10};
+ std::array<std::size_t, 3> const dims = MeshToolsLib::MeshGenerator::
+ VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
+
+ ASSERT_EQ(10, dims[0]);
+ ASSERT_EQ(5, dims[1]);
+ ASSERT_EQ(2, dims[2]);
+}
+
+TEST(MeshToolsLib, getNumberOfVoxelPerDimensionMinMax)
+{
+ std::array<double, 3> const cellsize = {1, 2, 5};
+ std::array<double, 3> const ranges = {0, 0, 10};
+ std::array<std::size_t, 3> const dims = MeshToolsLib::MeshGenerator::
+ VoxelGridFromMesh::getNumberOfVoxelPerDimension(ranges, cellsize);
+
+ ASSERT_EQ(1, dims[0]);
+ ASSERT_EQ(1, dims[1]);
+ ASSERT_EQ(2, dims[2]);
+}
+
+TEST(MeshToolsLib, getNumberOfVoxelPerDimensionCellSize)
+{
+ std::array<double, 3> const cellsize = {0, 2, 5};
+ std::array<double, 3> const ranges = {1, 2, 5};
+ ASSERT_ANY_THROW(MeshToolsLib::MeshGenerator::VoxelGridFromMesh::
+ getNumberOfVoxelPerDimension(ranges, cellsize));
+}
\ No newline at end of file