From 7b7abdd1a25afa52cde0399d315eb3583295a06e Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 5 Jun 2023 14:27:23 -0400 Subject: [PATCH] Fix dim = 1 mesh loading with gmsh. --- source/grid/grid_in.cc | 126 +++++++++++++++++++----------- tests/grid/grid_in_gmsh_04.cc | 75 ++++++++++++++++++ tests/grid/grid_in_gmsh_04.msh | 113 +++++++++++++++++++++++++++ tests/grid/grid_in_gmsh_04.output | 7 ++ 4 files changed, 277 insertions(+), 44 deletions(-) create mode 100644 tests/grid/grid_in_gmsh_04.cc create mode 100644 tests/grid/grid_in_gmsh_04.msh create mode 100644 tests/grid/grid_in_gmsh_04.output diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 03a3ae8208..24e456b551 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -64,35 +64,41 @@ namespace * boundary indicators on vertices after the triangulation has already been * created. * - * TODO: Fix this properly via SubcellData + * This function ignores all non-boundary vertices. Boundary ids can only be + * assigned to internal vertices - however, the grid fixup routines (e.g., + * combining duplicated vertices) may switch vertices from boundary to + * non-boundary. Get around this in two ways: + * + * 1. Ignore any attempt to set a boundary id on a non-boundary vertex. + * 2. Work around possible vertex renumberings by using Point + * instead of the corresponding vertex number (which may be out of date) + * + * TODO: fix this properly via SubCellData */ template void assign_1d_boundary_ids( - const std::map &boundary_ids, - Triangulation<1, spacedim> & triangulation) + const std::vector, types::boundary_id>> + & boundary_ids, + Triangulation<1, spacedim> &triangulation) { - if (boundary_ids.size() > 0) - for (const auto &cell : triangulation.active_cell_iterators()) - for (unsigned int f : GeometryInfo<1>::face_indices()) - if (boundary_ids.find(cell->vertex_index(f)) != boundary_ids.end()) - { - AssertThrow( - cell->at_boundary(f), - ExcMessage( - "You are trying to prescribe boundary ids on the face " - "of a 1d cell (i.e., on a vertex), but this face is not actually at " - "the boundary of the mesh. This is not allowed.")); - cell->face(f)->set_boundary_id( - boundary_ids.find(cell->vertex_index(f))->second); - } + for (auto &cell : triangulation.active_cell_iterators()) + for (unsigned int face_no : ReferenceCells::Line.face_indices()) + if (cell->face(face_no)->at_boundary()) + for (const auto &pair : boundary_ids) + if (cell->face(face_no)->vertex(0) == pair.first) + { + cell->face(face_no)->set_boundary_id(pair.second); + break; + } } template void - assign_1d_boundary_ids(const std::map &, - Triangulation &) + assign_1d_boundary_ids( + const std::vector, types::boundary_id>> &, + Triangulation &) { // we shouldn't get here since boundary ids are not assigned to // vertices except in 1d @@ -2319,6 +2325,12 @@ GridIn::read_msh(std::istream &in) SubCellData subcelldata; std::map boundary_ids_1d; + // Track the number of times each vertex is used in 1D. This determines + // whether or not we can assign a boundary id to a vertex. This is necessary + // because sometimes gmsh saves internal vertices in the $ELEM list in codim + // 1 or codim 2. + std::map vertex_counts; + { static constexpr std::array local_vertex_numbering = { {0, 1, 5, 4, 2, 3, 7, 6}}; @@ -2485,22 +2497,18 @@ GridIn::read_msh(std::istream &in) // allocate and read indices cells.emplace_back(); - cells.back().vertices.resize(vertices_per_cell); + CellData &cell = cells.back(); + cell.vertices.resize(vertices_per_cell); for (unsigned int i = 0; i < vertices_per_cell; ++i) { // hypercube cells need to be reordered if (vertices_per_cell == GeometryInfo::vertices_per_cell) - { - in >> cells.back() - .vertices[dim == 3 ? + in >> cell.vertices[dim == 3 ? local_vertex_numbering[i] : GeometryInfo::ucd_to_deal[i]]; - } else - { - in >> cells.back().vertices[i]; - } + in >> cell.vertices[i]; } // to make sure that the cast won't fail @@ -2514,21 +2522,21 @@ GridIn::read_msh(std::istream &in) // numbers::invalid_material_id-1 AssertIndexRange(material_id, numbers::invalid_material_id); - cells.back().material_id = material_id; + cell.material_id = material_id; // transform from gmsh to consecutive numbering for (unsigned int i = 0; i < vertices_per_cell; ++i) { - AssertThrow( - vertex_indices.find(cells.back().vertices[i]) != - vertex_indices.end(), - ExcInvalidVertexIndexGmsh(cell_per_entity, - elm_number, - cells.back().vertices[i])); - - // vertex with this index exists - cells.back().vertices[i] = - vertex_indices[cells.back().vertices[i]]; + AssertThrow(vertex_indices.find(cell.vertices[i]) != + vertex_indices.end(), + ExcInvalidVertexIndexGmsh(cell_per_entity, + elm_number, + cell.vertices[i])); + + const auto vertex = vertex_indices[cell.vertices[i]]; + if (dim == 1) + vertex_counts[vertex] += 1u; + cell.vertices[i] = vertex; } } else if ((cell_type == 1) && @@ -2662,13 +2670,23 @@ GridIn::read_msh(std::istream &in) ExcGmshNoCellInformation(subcelldata.boundary_lines.size(), subcelldata.boundary_quads.size())); + // apply_grid_fixup_functions() may invalidate the vertex indices (since it + // will delete duplicated or unused vertices). Get around this by storing + // Points directly in that case so that the comparisons are valid. + std::vector, types::boundary_id>> boundary_id_pairs; + if (dim == 1) + for (const auto &pair : vertex_counts) + if (pair.second == 1u) + boundary_id_pairs.emplace_back(vertices[pair.first], + boundary_ids_1d[pair.first]); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); // in 1d, we also have to attach boundary ids to vertices, which does not - // currently work through the call above + // currently work through the call above. if (dim == 1) - assign_1d_boundary_ids(boundary_ids_1d, *tria); + assign_1d_boundary_ids(boundary_id_pairs, *tria); } @@ -2699,6 +2717,12 @@ GridIn::read_msh(const std::string &fname) SubCellData subcelldata; std::map boundary_ids_1d; + // Track the number of times each vertex is used in 1D. This determines + // whether or not we can assign a boundary id to a vertex. This is necessary + // because sometimes gmsh saves internal vertices in the $ELEM list in codim + // 1 or codim 2. + std::map vertex_counts; + gmsh::initialize(); gmsh::option::setNumber("General.Verbosity", 0); gmsh::open(fname); @@ -2832,8 +2856,12 @@ GridIn::read_msh(const std::string &fname) cells.emplace_back(n_vertices); auto &cell = cells.back(); for (unsigned int v = 0; v < n_vertices; ++v) - cell.vertices[v] = - nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1; + { + cell.vertices[v] = + nodes[n_vertices * j + gmsh_to_dealii[type][v]] - 1; + if (dim == 1) + vertex_counts[cell.vertices[v]] += 1u; + } cell.manifold_id = manifold_id; cell.material_id = boundary_id; } @@ -2870,13 +2898,23 @@ GridIn::read_msh(const std::string &fname) } } + // apply_grid_fixup_functions() may invalidate the vertex indices (since it + // will delete duplicated or unused vertices). Get around this by storing + // Points directly in that case so that the comparisons are valid. + std::vector, types::boundary_id>> boundary_id_pairs; + if (dim == 1) + for (const auto &pair : vertex_counts) + if (pair.second == 1u) + boundary_id_pairs.emplace_back(vertices[pair.first], + boundary_ids_1d[pair.first]); + apply_grid_fixup_functions(vertices, cells, subcelldata); tria->create_triangulation(vertices, cells, subcelldata); // in 1d, we also have to attach boundary ids to vertices, which does not - // currently work through the call above + // currently work through the call above. if (dim == 1) - assign_1d_boundary_ids(boundary_ids_1d, *tria); + assign_1d_boundary_ids(boundary_id_pairs, *tria); gmsh::clear(); gmsh::finalize(); diff --git a/tests/grid/grid_in_gmsh_04.cc b/tests/grid/grid_in_gmsh_04.cc new file mode 100644 index 0000000000..f8eac02427 --- /dev/null +++ b/tests/grid/grid_in_gmsh_04.cc @@ -0,0 +1,75 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test that we can read dim = 1, codim = 2 meshes with boundary ids correctly. + +#include +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +template +void +print_mesh_info(const Triangulation &triangulation, + const std::string & filename) +{ + deallog << "Mesh info:" << std::endl + << " dimension: " << dim << std::endl + << " no. of cells: " << triangulation.n_active_cells() << std::endl; + + { + std::map boundary_count; + for (const auto &cell : triangulation.active_cell_iterators()) + for (unsigned int face_no : cell->face_indices()) + if (cell->face(face_no)->at_boundary()) + boundary_count[cell->face(face_no)->boundary_id()]++; + + deallog << " boundary indicators: "; + for (const std::pair &pair : + boundary_count) + { + deallog << pair.first << "(" << pair.second << " times) "; + } + deallog << std::endl; + } + + // Finally, produce a graphical representation of the mesh to an output + // file: + std::ofstream out(filename); + GridOut grid_out; + grid_out.write_vtu(triangulation, out); + deallog << " written to " << filename << std::endl << std::endl; +} + + +int +main() +{ + initlog(); + std::ifstream in(SOURCE_DIR "/grid_in_gmsh_04.msh"); + + Triangulation<1, 2> tria; + GridIn<1, 2> gi; + gi.attach_triangulation(tria); + gi.read_msh(in); + + print_mesh_info(tria, "grid_in_gmsh_04.vtu"); +} diff --git a/tests/grid/grid_in_gmsh_04.msh b/tests/grid/grid_in_gmsh_04.msh new file mode 100644 index 0000000000..8582326fa0 --- /dev/null +++ b/tests/grid/grid_in_gmsh_04.msh @@ -0,0 +1,113 @@ +$NOD +49 +1 -6 -4 0 +2 -6 -4 0 +3 0 0 0 +4 6 -4 0 +5 -4 -4 0 +6 -6 0 0 +7 6 0 0 +8 4 -4 0 +9 4 -2 0 +10 -4 -2 0 +11 -4 -2.5 0 +12 -4 -3 0 +13 -4 -3.5 0 +14 -4.5 -4 0 +15 -5 -4 0 +16 -5.5 -4 0 +17 -6 -3 0 +18 -6 -2 0 +19 -6 -1 0 +20 -5.25 0 0 +21 -4.5 0 0 +22 -3.75 0 0 +23 -3 0 0 +24 -2.25 0 0 +25 -1.5 0 0 +26 -0.75 0 0 +27 0.75 0 0 +28 1.5 0 0 +29 2.25 0 0 +30 3 0 0 +31 3.75 0 0 +32 4.5 0 0 +33 5.25 0 0 +34 6 -1 0 +35 6 -2 0 +36 6 -3 0 +37 5.5 -4 0 +38 5 -4 0 +39 4.5 -4 0 +40 4 -3.5 0 +41 4 -3 0 +42 4 -2.5 0 +43 3 -2 0 +44 2 -2 0 +45 1 -2 0 +46 0 -2 0 +47 -1 -2 0 +48 -2 -2 0 +49 -3 -2 0 +$ENDNOD +$ELM +58 +1 15 0 1 1 1 +2 15 0 2 1 2 +3 15 0 3 1 3 +4 15 0 4 1 4 +5 15 0 6 1 5 +6 15 0 7 1 6 +7 15 0 8 1 7 +8 15 0 9 1 8 +9 15 0 10 1 9 +10 15 0 11 1 10 +11 1 0 1 2 10 11 +12 1 0 1 2 11 12 +13 1 0 1 2 12 13 +14 1 0 1 2 13 5 +15 1 0 2 2 5 14 +16 1 0 2 2 14 15 +17 1 0 2 2 15 16 +18 1 0 2 2 16 1 +19 1 0 3 2 1 17 +20 1 0 3 2 17 18 +21 1 0 3 2 18 19 +22 1 0 3 2 19 6 +23 1 0 4 2 6 20 +24 1 0 4 2 20 21 +25 1 0 4 2 21 22 +26 1 0 4 2 22 23 +27 1 0 4 2 23 24 +28 1 0 4 2 24 25 +29 1 0 4 2 25 26 +30 1 0 4 2 26 3 +31 1 0 5 2 3 27 +32 1 0 5 2 27 28 +33 1 0 5 2 28 29 +34 1 0 5 2 29 30 +35 1 0 5 2 30 31 +36 1 0 5 2 31 32 +37 1 0 5 2 32 33 +38 1 0 5 2 33 7 +39 1 0 6 2 7 34 +40 1 0 6 2 34 35 +41 1 0 6 2 35 36 +42 1 0 6 2 36 4 +43 1 0 7 2 4 37 +44 1 0 7 2 37 38 +45 1 0 7 2 38 39 +46 1 0 7 2 39 8 +47 1 0 8 2 8 40 +48 1 0 8 2 40 41 +49 1 0 8 2 41 42 +50 1 0 8 2 42 9 +51 1 0 9 2 9 43 +52 1 0 9 2 43 44 +53 1 0 9 2 44 45 +54 1 0 9 2 45 46 +55 1 0 9 2 46 47 +56 1 0 9 2 47 48 +57 1 0 9 2 48 49 +58 1 0 9 2 49 10 +$ENDELM diff --git a/tests/grid/grid_in_gmsh_04.output b/tests/grid/grid_in_gmsh_04.output new file mode 100644 index 0000000000..5b718101cc --- /dev/null +++ b/tests/grid/grid_in_gmsh_04.output @@ -0,0 +1,7 @@ + +DEAL::Mesh info: +DEAL:: dimension: 1 +DEAL:: no. of cells: 48 +DEAL:: boundary indicators: +DEAL:: written to grid_in_gmsh_04.vtu +DEAL:: -- 2.39.5