diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index b278259..81f7b02 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -10,7 +10,7 @@ jobs: debug: runs-on: ubuntu-latest container: - image: dealii/dealii:v9.6.0-noble + image: dealii/dealii:v9.7.1-noble options: --user root steps: diff --git a/backends/dealii/include/register_types.h b/backends/dealii/include/register_types.h index 1573140..f5c98c4 100644 --- a/backends/dealii/include/register_types.h +++ b/backends/dealii/include/register_types.h @@ -20,6 +20,7 @@ #include "coral_network.h" #include "laplace.h" #include "poisson.h" +#include "utilities.h" /** \cond INTERNAL */ namespace nlohmann @@ -121,6 +122,12 @@ namespace coral "grid_generator_function_name", "grid_generator_function_arguments"}); + NodeObject::register_function(read_grid, + {"dealii::read_grid<" + + Utilities::dim_string(dim, spacedim) + ">", + "file_name", + "triangulation"}); + NodeObject:: register_method, void, unsigned int>( &Triangulation::refine_global, diff --git a/backends/dealii/include/utilities.h b/backends/dealii/include/utilities.h new file mode 100644 index 0000000..259fdb6 --- /dev/null +++ b/backends/dealii/include/utilities.h @@ -0,0 +1,364 @@ +#ifndef CORAL_BACKENDS_UTILITIES_H +#define CORAL_BACKENDS_UTILITIES_H + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#ifdef DEAL_II_WITH_VTK +# include +# include +# include +# include +# include +# include +# include +# include +#endif + +namespace dealii +{ + namespace internal + { + inline std::string + lowercase_extension(const std::string &file_name) + { + std::string extension = + std::filesystem::path(file_name).extension().string(); + std::transform(extension.begin(), + extension.end(), + extension.begin(), + [](const unsigned char c) { + return static_cast(std::tolower(c)); + }); + return extension; + } + +#ifdef DEAL_II_WITH_VTK + inline auto + get_optional_vtk_array(vtkCellData *cell_data, const char *name) + -> vtkDataArray * + { + return cell_data != nullptr ? cell_data->GetArray(name) : nullptr; + } + + + + inline long long + get_vtk_id(vtkDataArray *array, const vtkIdType cell_index) + { + return array != nullptr ? + static_cast(array->GetComponent(cell_index, 0)) : + 0LL; + } + + + + template + inline void + fill_vertices(CellData &cell_data, vtkCell *cell) + { + for (unsigned int j = 0; j < cell_data.vertices.size(); ++j) + cell_data.vertices[j] = static_cast(cell->GetPointId(j)); + } + + + + template + inline void + reorder_vtk_vertices(CellData &cell_data, const int vtk_cell_type) + { + if constexpr (dim == 2) + { + if (vtk_cell_type == VTK_QUAD) + std::swap(cell_data.vertices[2], cell_data.vertices[3]); + } + else if constexpr (dim == 3) + { + if (vtk_cell_type == VTK_HEXAHEDRON) + { + std::swap(cell_data.vertices[2], cell_data.vertices[3]); + std::swap(cell_data.vertices[6], cell_data.vertices[7]); + } + } + } + + + + template + inline void + set_subcell_metadata(CellData &subcell, + vtkDataArray *boundary_or_material_ids, + vtkDataArray *manifold_ids, + const vtkIdType cell_index) + { + subcell.boundary_id = + boundary_or_material_ids != nullptr ? + static_cast( + get_vtk_id(boundary_or_material_ids, cell_index)) : + types::boundary_id{0}; + subcell.manifold_id = manifold_ids != nullptr ? + static_cast( + get_vtk_id(manifold_ids, cell_index)) : + numbers::flat_manifold_id; + } + + + + template + inline void + append_vtk_cell(vtkCell *cell, + const vtkIdType cell_index, + vtkDataArray *boundary_or_material_ids, + vtkDataArray *manifold_ids, + std::vector> &cells, + SubCellData &subcell_data) + { + const int vtk_cell_type = cell->GetCellType(); + + if constexpr (dim == 1) + { + AssertThrow(vtk_cell_type == VTK_LINE, + ExcMessage("Unsupported cell type in 1D VTK file: only " + "VTK_LINE is supported.")); + AssertThrow(cell->GetNumberOfPoints() == 2, + ExcMessage( + "Only line cells with 2 points are supported.")); + + CellData<1> cell_data(2); + fill_vertices(cell_data, cell); + cell_data.material_id = static_cast( + get_vtk_id(boundary_or_material_ids, cell_index)); + cell_data.manifold_id = manifold_ids != nullptr ? + static_cast( + get_vtk_id(manifold_ids, cell_index)) : + numbers::flat_manifold_id; + cells.push_back(cell_data); + } + else if constexpr (dim == 2) + { + if (vtk_cell_type == VTK_QUAD || vtk_cell_type == VTK_TRIANGLE) + { + const unsigned int n_vertices = + vtk_cell_type == VTK_QUAD ? 4U : 3U; + AssertThrow(static_cast( + cell->GetNumberOfPoints()) == n_vertices, + ExcMessage("Unexpected number of vertices for a 2D " + "VTK cell.")); + + CellData<2> cell_data(n_vertices); + fill_vertices(cell_data, cell); + reorder_vtk_vertices(cell_data, vtk_cell_type); + cell_data.material_id = static_cast( + get_vtk_id(boundary_or_material_ids, cell_index)); + cell_data.manifold_id = + manifold_ids != nullptr ? + static_cast( + get_vtk_id(manifold_ids, cell_index)) : + numbers::flat_manifold_id; + cells.push_back(cell_data); + return; + } + + AssertThrow(vtk_cell_type == VTK_LINE, + ExcMessage("Unsupported cell type in 2D VTK file: only " + "VTK_QUAD, VTK_TRIANGLE, and VTK_LINE are " + "supported.")); + AssertThrow(cell->GetNumberOfPoints() == 2, + ExcMessage("Only line subcells with 2 points are " + "supported.")); + + CellData<1> line_data(2); + fill_vertices(line_data, cell); + set_subcell_metadata(line_data, + boundary_or_material_ids, + manifold_ids, + cell_index); + subcell_data.boundary_lines.push_back(line_data); + } + else if constexpr (dim == 3) + { + switch (vtk_cell_type) + { + case VTK_HEXAHEDRON: + case VTK_TETRA: + case VTK_WEDGE: + case VTK_PYRAMID: + { + const unsigned int n_vertices = + vtk_cell_type == VTK_HEXAHEDRON ? 8U : + vtk_cell_type == VTK_TETRA ? 4U : + vtk_cell_type == VTK_WEDGE ? 6U : + 5U; + + AssertThrow(static_cast( + cell->GetNumberOfPoints()) == n_vertices, + ExcMessage( + "Unexpected number of vertices for a 3D VTK " + "cell.")); + + CellData<3> cell_data(n_vertices); + fill_vertices(cell_data, cell); + reorder_vtk_vertices(cell_data, vtk_cell_type); + cell_data.material_id = static_cast( + get_vtk_id(boundary_or_material_ids, cell_index)); + cell_data.manifold_id = + manifold_ids != nullptr ? + static_cast( + get_vtk_id(manifold_ids, cell_index)) : + numbers::flat_manifold_id; + cells.push_back(cell_data); + return; + } + + case VTK_QUAD: + case VTK_TRIANGLE: + { + const unsigned int n_vertices = + vtk_cell_type == VTK_QUAD ? 4U : 3U; + AssertThrow( + static_cast(cell->GetNumberOfPoints()) == + n_vertices, + ExcMessage("Unexpected number of vertices for a 3D face.")); + + CellData<2> face_data(n_vertices); + fill_vertices(face_data, cell); + set_subcell_metadata(face_data, + boundary_or_material_ids, + manifold_ids, + cell_index); + subcell_data.boundary_quads.push_back(face_data); + return; + } + + case VTK_LINE: + { + AssertThrow(cell->GetNumberOfPoints() == 2, + ExcMessage("Only line subcells with 2 points are " + "supported.")); + + CellData<1> line_data(2); + fill_vertices(line_data, cell); + set_subcell_metadata(line_data, + boundary_or_material_ids, + manifold_ids, + cell_index); + subcell_data.boundary_lines.push_back(line_data); + return; + } + + default: + AssertThrow(false, + ExcMessage("Unsupported cell type in 3D VTK file: " + "only VTK_HEXAHEDRON, VTK_TETRA, " + "VTK_WEDGE, VTK_PYRAMID, VTK_QUAD, " + "VTK_TRIANGLE, and VTK_LINE are " + "supported.")); + } + } + else + { + AssertThrow(false, ExcMessage("Unsupported dimension.")); + } + } +#endif + } // namespace internal + + + + template + void + read_grid(const std::string &file_name, + Triangulation &triangulation) + { +#ifndef DEAL_II_WITH_VTK + GridIn grid_in(triangulation); + grid_in.read(file_name); +#else + std::ifstream file(file_name); + AssertThrow(file.good(), ExcMessage("VTK file not found: " + file_name)); + + vtkSmartPointer grid; + const auto extension = internal::lowercase_extension(file_name); + + if (extension == ".vtk") + { + auto reader = vtkSmartPointer::New(); + reader->SetFileName(file_name.c_str()); + reader->Update(); + grid = reader->GetOutput(); + } + else if (extension == ".vtu") + { + auto reader = vtkSmartPointer::New(); + reader->SetFileName(file_name.c_str()); + reader->Update(); + grid = reader->GetOutput(); + } + else + { + AssertThrow(false, + ExcMessage("Unsupported VTK grid extension: " + extension)); + } + + AssertThrow(grid != nullptr, ExcMessage("Failed to read VTK grid.")); + + vtkPoints *vtk_points = grid->GetPoints(); + AssertThrow(vtk_points != nullptr, + ExcMessage("VTK grid does not contain point data.")); + + const vtkIdType n_points = vtk_points->GetNumberOfPoints(); + std::vector> points(static_cast(n_points)); + for (vtkIdType i = 0; i < n_points; ++i) + { + std::array coords{{0.0, 0.0, 0.0}}; + vtk_points->GetPoint(i, coords.data()); + + for (unsigned int d = 0; d < spacedim; ++d) + points[static_cast(i)][d] = coords[d]; + + for (unsigned int d = spacedim; d < 3; ++d) + AssertThrow(coords[d] == 0.0, + ExcMessage("VTK grid has non-zero coordinate in an " + "unused dimension.")); + } + + vtkCellData *vtk_cell_data = grid->GetCellData(); + vtkDataArray *boundary_or_material_ids = + internal::get_optional_vtk_array(vtk_cell_data, "MaterialID"); + if (boundary_or_material_ids == nullptr) + boundary_or_material_ids = + internal::get_optional_vtk_array(vtk_cell_data, "MaterialID"); + vtkDataArray *manifold_ids = + internal::get_optional_vtk_array(vtk_cell_data, "ManifoldID"); + + std::vector> cells; + SubCellData subcell_data; + + const vtkIdType n_cells = grid->GetNumberOfCells(); + cells.reserve(static_cast(n_cells)); + + for (vtkIdType i = 0; i < n_cells; ++i) + internal::append_vtk_cell(grid->GetCell(i), + i, + boundary_or_material_ids, + manifold_ids, + cells, + subcell_data); + + triangulation.create_triangulation(points, cells, subcell_data); +#endif + } +} // namespace dealii + +#endif diff --git a/backends/dealii/tests/dealii_examples.cc b/backends/dealii/tests/dealii_examples.cc index f19e8d6..b49b298 100644 --- a/backends/dealii/tests/dealii_examples.cc +++ b/backends/dealii/tests/dealii_examples.cc @@ -407,6 +407,51 @@ TEST(dealiiExamples, PoissonSolver) file.close(); } +TEST(dealiiExamples, ReadGridVtu) +{ + Triangulation<2> triangulation; + + read_grid<2>((std::filesystem::path(SOURCE_DIR) / "test_files" / + "test_grid.vtu") + .string(), + triangulation); + + EXPECT_EQ(triangulation.n_active_cells(), 16); + EXPECT_EQ(triangulation.n_used_vertices(), 25); + + std::set cell_manifold_ids; + std::set face_manifold_ids; + std::set material_ids; + std::set boundary_ids; + for (const auto &cell : triangulation.active_cell_iterators()) + { + material_ids.insert(cell->material_id()); + cell_manifold_ids.insert(cell->manifold_id()); + for (const auto &f : cell->face_iterators()) + { + if (f->at_boundary()) + boundary_ids.insert(f->boundary_id()); + face_manifold_ids.insert(f->manifold_id()); + } + } + + + EXPECT_EQ(material_ids.size(), 1u); + EXPECT_TRUE(material_ids.find(types::material_id{0}) != material_ids.end()); + + EXPECT_EQ(boundary_ids.size(), 2u); + EXPECT_TRUE(boundary_ids.find(types::boundary_id{0}) != boundary_ids.end()); + EXPECT_TRUE(boundary_ids.find(types::boundary_id{1}) != boundary_ids.end()); + + EXPECT_EQ(cell_manifold_ids.size(), 1u); + EXPECT_TRUE(cell_manifold_ids.find(numbers::flat_manifold_id) != + cell_manifold_ids.end()); + + EXPECT_EQ(face_manifold_ids.size(), 1u); + EXPECT_TRUE(face_manifold_ids.find(numbers::flat_manifold_id) != + face_manifold_ids.end()); +} + TEST(dealiiExamples, NetworkFromJsonPoissonSolverSolution) { ScopedTestOutputDir output_dir( diff --git a/test_files/test_grid.vtu b/test_files/test_grid.vtu new file mode 100644 index 0000000..d4ff4bd --- /dev/null +++ b/test_files/test_grid.vtu @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + _AQAAAACAAACAAAAAEwAAAA==eJxjYKAOYESiGfEpRAMAARQABA==AQAAAACAAACAAAAADAAAAA==eJz7/39gAQAiIH+BAQAAAACAAABYAgAAUgAAAA==eJyNkLERADAIArNZ9u8cwZIyI6Syec8QOwERXetVZz/pxldPPI0P+Zonnsaf+hj2yuT5nasc1BdOfQz5Ze6ifxp//kHmP/TXgMewVw2/qdItKQ==AQAAAACAAAAAAwAAigAAAA==eJx1kUcOhEAMBMlxF9hA+P9LOeC6lOS5tNwqexyK4nlj6B46y29CPwlH/Att5ZM3iKPuFHqIxy9D3wlH/Nc/+OR14uiTfs7Ql3zqfROOeA2t5JPXi6Mu/Vzi8el7STjiTf/gk1eL43HnOfG5UyvO9/O9fDffgf1R1/urxHl/3ovn99ye/wax4QSrAQAAAACAAAAAAQAAUwAAAA==eJwtxaEORAAAAFBBEARBEARBEAThgiDY7WZmdruZ2c3M/P9XCN4rLwwekWMnTp05d+HSlWs3frl1595vfzx49OTZX/+8ePXmv3cfPn35BroeBzE=AQAAAACAAAAgAAAADgAAAA==eJzj5EQFzGgAAA+AAME= + +