* 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<spacedim>
+ * instead of the corresponding vertex number (which may be out of date)
+ *
+ * TODO: fix this properly via SubCellData
*/
template <int spacedim>
void
assign_1d_boundary_ids(
- const std::map<unsigned int, types::boundary_id> &boundary_ids,
- Triangulation<1, spacedim> & triangulation)
+ const std::vector<std::pair<Point<spacedim>, 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 <int dim, int spacedim>
void
- assign_1d_boundary_ids(const std::map<unsigned int, types::boundary_id> &,
- Triangulation<dim, spacedim> &)
+ assign_1d_boundary_ids(
+ const std::vector<std::pair<Point<spacedim>, types::boundary_id>> &,
+ Triangulation<dim, spacedim> &)
{
// we shouldn't get here since boundary ids are not assigned to
// vertices except in 1d
SubCellData subcelldata;
std::map<unsigned int, types::boundary_id> 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<unsigned int, unsigned int> vertex_counts;
+
{
static constexpr std::array<unsigned int, 8> local_vertex_numbering = {
{0, 1, 5, 4, 2, 3, 7, 6}};
// allocate and read indices
cells.emplace_back();
- cells.back().vertices.resize(vertices_per_cell);
+ CellData<dim> &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<dim>::vertices_per_cell)
- {
- in >> cells.back()
- .vertices[dim == 3 ?
+ in >> cell.vertices[dim == 3 ?
local_vertex_numbering[i] :
GeometryInfo<dim>::ucd_to_deal[i]];
- }
else
- {
- in >> cells.back().vertices[i];
- }
+ in >> cell.vertices[i];
}
// to make sure that the cast won't fail
// 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) &&
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<std::pair<Point<spacedim>, 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);
}
SubCellData subcelldata;
std::map<unsigned int, types::boundary_id> 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<unsigned int, unsigned int> vertex_counts;
+
gmsh::initialize();
gmsh::option::setNumber("General.Verbosity", 0);
gmsh::open(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;
}
}
}
+ // 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<std::pair<Point<spacedim>, 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();
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+#include <map>
+
+#include "../tests.h"
+
+template <int dim, int spacedim = dim>
+void
+print_mesh_info(const Triangulation<dim, spacedim> &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<types::boundary_id, unsigned int> 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<const types::boundary_id, unsigned int> &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");
+}
projects / dealii.git / commitdiff