From 9e0a3d7b5a784122bc458429a915490d90626bcf Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 25 Aug 2015 15:16:11 -0500 Subject: [PATCH] Set boundary indicators in GridIn on vertices in 1d. We did not keep track of boundary indicators on vertices at all, even though these are faces in 1d. Fix this. --- doc/news/changes.h | 6 ++ source/grid/grid_in.cc | 81 ++++++++++++++--- tests/deal.II/grid_in_msh_02.cc | 41 ++++----- tests/deal.II/grid_in_msh_02.output | 129 +--------------------------- tests/grid/grid_in_msh_02.cc | 93 ++++++++++++++++++++ tests/grid/grid_in_msh_02.output | 0 tests/grid/grids/grid_in_msh_02.msh | 32 +++++++ 7 files changed, 222 insertions(+), 160 deletions(-) create mode 100644 tests/grid/grid_in_msh_02.cc create mode 100644 tests/grid/grid_in_msh_02.output create mode 100644 tests/grid/grids/grid_in_msh_02.msh diff --git a/doc/news/changes.h b/doc/news/changes.h index 5965319d36..c3f7f5cff3 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -108,6 +108,12 @@ inconvenience this causes.
    +
  1. Fixed: In 1d, GridIn::read_msh() ignored boundary indicators + associated with vertices. This is now fixed. +
    + (Jan Strebel, Wolfgang Bangerth, 2015/08/25) +
    2. +
  2. Improved: The interface and documentation for periodic boundary conditions have been restructured. A @ref GlossPeriodicConstraints "glossary entry" has been written. diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index a516dd1e32..29c131f0fa 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -37,6 +37,49 @@ DEAL_II_NAMESPACE_OPEN + +namespace +{ + /** + * In 1d, boundary indicators are associated with vertices, but this is not + * currently passed through the SubcellData structure. This function sets + * boundary indicators on vertices after the triangulation has already been + * created. + * + * TODO: This this properly via SubcellData + */ + template + void + assign_1d_boundary_indicators (const std::map &boundary_ids, + Triangulation<1,spacedim> &triangulation) + { + if (boundary_ids.size() > 0) + for (typename Triangulation<1,spacedim>::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + 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); + } + } + + + template + void + assign_1d_boundary_indicators (const std::map &, + Triangulation &) + { + // we shouldn't get here since boundary ids are not assigned to + // vertices except in 1d + Assert (dim != 1, ExcInternalError()); + } +} + template GridIn::GridIn () : tria(0, typeid(*this).name()), default_format(ucd) @@ -1244,9 +1287,12 @@ void GridIn::read_msh (std::istream &in) in >> n_cells; - // set up array of cells - std::vector > cells; - SubCellData subcelldata; + // set up array of cells and subcells (faces). In 1d, there is currently no + // standard way in deal.II to pass boundary indicators attached to individual + // vertices, so do this by hand via the boundary_ids_1d array + std::vector > cells; + SubCellData subcelldata; + std::map boundary_ids_1d; for (unsigned int cell=0; cell::read_msh (std::istream &in) case 2: { - // read the tags; ignore - // all but the first one + // read the tags; ignore all but the first one which we will + // interpret as the material_id (for cells) or boundary_id + // (for faces) unsigned int n_tags; in >> n_tags; if (n_tags > 0) @@ -1431,29 +1478,32 @@ void GridIn::read_msh (std::istream &in) subcelldata.boundary_quads.back().vertices[i])); subcelldata.boundary_quads.back().vertices[i] = numbers::invalid_unsigned_int; - }; + } } else if (cell_type == 15) { - // Ignore vertices - // but read the - // number of nodes - // given + // read the indices of nodes given + unsigned int node_index; switch (gmsh_file_format) { case 1: { for (unsigned int i=0; i> dummy; + in >> node_index; break; } case 2: { - in >> dummy; + in >> node_index; break; } } + + // we only care about boundary indicators assigned to individual + // vertices in 1d (because otherwise the vertices are not faces) + if (dim == 1) + boundary_ids_1d[vertex_indices[node_index]] = material_id; } else // cannot read this, so throw @@ -1474,7 +1524,7 @@ void GridIn::read_msh (std::istream &in) AssertThrow (false, ExcGmshUnsupportedGeometry(cell_type)); } - }; + } // Assert we reached the end of the block in >> line; @@ -1499,6 +1549,11 @@ void GridIn::read_msh (std::istream &in) GridReordering::invert_all_cells_of_negative_grid (vertices, cells); GridReordering::reorder_cells (cells); tria->create_triangulation_compatibility (vertices, cells, subcelldata); + + // in 1d, we also have to attach boundary ids to vertices, which does not + // currently work through the call above + if (dim == 1) + assign_1d_boundary_indicators (boundary_ids_1d, *tria); } diff --git a/tests/deal.II/grid_in_msh_02.cc b/tests/deal.II/grid_in_msh_02.cc index d3144ce76e..cfa014f4e3 100644 --- a/tests/deal.II/grid_in_msh_02.cc +++ b/tests/deal.II/grid_in_msh_02.cc @@ -14,9 +14,11 @@ // --------------------------------------------------------------------- -// read a file in the MSH format used by the GMSH program. this -// particular mesh has type-15 cells (nodes) with more than one -// associated vertex. we failed to read the vertices +// in 1d, we have to read vertex information to set boundary +// indicators +// +// test case by Jan Strebel + #include "../tests.h" #include @@ -37,25 +39,24 @@ std::ofstream logfile("output"); -template -void check_file (const std::string name, - typename GridIn::Format format) +void check_file () { - Triangulation tria; - GridIn gi; + Triangulation<1> tria; + GridIn<1> gi; gi.attach_triangulation (tria); - gi.read(name, format); - deallog << '\t' << tria.n_vertices() - << '\t' << tria.n_cells() - << std::endl; + std::ifstream in (SOURCE_DIR "/../grid/grids/grid_in_msh_02.msh"); + gi.read_msh(in); - GridOut grid_out; - grid_out.write_gnuplot (tria, deallog.get_file_stream()); -} - -void filename_resolution() -{ - check_file<2> (std::string(SOURCE_DIR "/grid_in_msh_02/mesh"), GridIn<2>::msh); + for (Triangulation<1>::active_cell_iterator cell = tria.begin_active(); cell != tria.end(); ++cell) + { + for (unsigned int face = 0; face < 2; ++face) + { + if (cell->at_boundary(face)) + deallog << "vertex " << cell->face_index(face) + << " has boundary indicator " << (int)cell->face(face)->boundary_indicator() + << std::endl; + } + } } @@ -67,6 +68,6 @@ int main () deallog.depth_console(0); deallog.threshold_double(1.e-10); - filename_resolution(); + check_file (); } diff --git a/tests/deal.II/grid_in_msh_02.output b/tests/deal.II/grid_in_msh_02.output index 81a9950edd..2ebba38d85 100644 --- a/tests/deal.II/grid_in_msh_02.output +++ b/tests/deal.II/grid_in_msh_02.output @@ -1,128 +1,3 @@ -DEAL:: 28 18 -2.0e+03 0.0 0 0 -2.0e+03 6.7e+02 0 0 -1.4e+03 5.0e+02 0 0 -1.5e+03 0.0 0 0 -2.0e+03 0.0 0 0 - - -1.5e+03 0.0 0 0 -1.4e+03 5.0e+02 0 0 -9.6e+02 4.0e+02 0 0 -1.0e+03 0.0 0 0 -1.5e+03 0.0 0 0 - - -1.0e+03 0.0 0 0 -9.6e+02 4.0e+02 0 0 -5.0e+02 3.3e+02 0 0 -5.0e+02 0.0 0 0 -1.0e+03 0.0 0 0 - - -2.0e+03 6.7e+02 0 0 -2.0e+03 1.3e+03 0 0 -1.4e+03 9.6e+02 0 0 -1.4e+03 5.0e+02 0 0 -2.0e+03 6.7e+02 0 0 - - -1.4e+03 5.0e+02 0 0 -1.4e+03 9.6e+02 0 0 -9.0e+02 7.8e+02 0 0 -9.6e+02 4.0e+02 0 0 -1.4e+03 5.0e+02 0 0 - - -9.6e+02 4.0e+02 0 0 -9.0e+02 7.8e+02 0 0 -5.0e+02 6.7e+02 0 0 -5.0e+02 3.3e+02 0 0 -9.6e+02 4.0e+02 0 0 - - -2.0e+03 1.3e+03 0 0 -2.0e+03 2.0e+03 0 0 -1.2e+03 1.3e+03 0 0 -1.4e+03 9.6e+02 0 0 -2.0e+03 1.3e+03 0 0 - - -1.4e+03 9.6e+02 0 0 -1.2e+03 1.3e+03 0 0 -7.9e+02 1.1e+03 0 0 -9.0e+02 7.8e+02 0 0 -1.4e+03 9.6e+02 0 0 - - -9.0e+02 7.8e+02 0 0 -7.9e+02 1.1e+03 0 0 -5.0e+02 1.0e+03 0 0 -5.0e+02 6.7e+02 0 0 -9.0e+02 7.8e+02 0 0 - - -2.0e+03 2.0e+03 0 0 -1.3e+03 2.0e+03 0 0 -8.4e+02 1.6e+03 0 0 -1.2e+03 1.3e+03 0 0 -2.0e+03 2.0e+03 0 0 - - -1.2e+03 1.3e+03 0 0 -8.4e+02 1.6e+03 0 0 -6.0e+02 1.4e+03 0 0 -7.9e+02 1.1e+03 0 0 -1.2e+03 1.3e+03 0 0 - - -7.9e+02 1.1e+03 0 0 -6.0e+02 1.4e+03 0 0 -4.3e+02 1.2e+03 0 0 -5.0e+02 1.0e+03 0 0 -7.9e+02 1.1e+03 0 0 - - -1.3e+03 2.0e+03 0 0 -6.7e+02 2.0e+03 0 0 -4.4e+02 1.8e+03 0 0 -8.4e+02 1.6e+03 0 0 -1.3e+03 2.0e+03 0 0 - - -8.4e+02 1.6e+03 0 0 -4.4e+02 1.8e+03 0 0 -3.3e+02 1.6e+03 0 0 -6.0e+02 1.4e+03 0 0 -8.4e+02 1.6e+03 0 0 - - -6.0e+02 1.4e+03 0 0 -3.3e+02 1.6e+03 0 0 -2.5e+02 1.4e+03 0 0 -4.3e+02 1.2e+03 0 0 -6.0e+02 1.4e+03 0 0 - - -6.7e+02 2.0e+03 0 0 -0.0 2.0e+03 0 0 -0.0 1.8e+03 0 0 -4.4e+02 1.8e+03 0 0 -6.7e+02 2.0e+03 0 0 - - -4.4e+02 1.8e+03 0 0 -0.0 1.8e+03 0 0 -0.0 1.7e+03 0 0 -3.3e+02 1.6e+03 0 0 -4.4e+02 1.8e+03 0 0 - - -3.3e+02 1.6e+03 0 0 -0.0 1.7e+03 0 0 -0.0 1.5e+03 0 0 -2.5e+02 1.4e+03 0 0 -3.3e+02 1.6e+03 0 0 - - +DEAL::vertex 0 has boundary indicator 10 +DEAL::vertex 1 has boundary indicator 20 diff --git a/tests/grid/grid_in_msh_02.cc b/tests/grid/grid_in_msh_02.cc new file mode 100644 index 0000000000..b615e5cba3 --- /dev/null +++ b/tests/grid/grid_in_msh_02.cc @@ -0,0 +1,93 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2014 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// check whether we can read in with the gmsh format + +#include "../tests.h" +#include + +#include +#include +#include +#include + +#include +#include + + +template +void gmsh_grid (const char *name) +{ + Triangulation tria; + GridIn grid_in; + grid_in.attach_triangulation (tria); + std::ifstream input_file(name); + grid_in.read_msh(input_file); + + deallog << " " << tria.n_active_cells() << " active cells" << std::endl; + + int hash = 0; + int index = 0; + for (typename Triangulation::active_cell_iterator c=tria.begin_active(); + c!=tria.end(); ++c, ++index) + for (unsigned int i=0; i::vertices_per_cell; ++i) + hash += (index * i * c->vertex_index(i)) % (tria.n_active_cells()+1); + deallog << " hash=" << hash << std::endl; +} + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + gmsh_grid<2> (SOURCE_DIR "/grids/grid_in_msh_01.2d.msh"); + gmsh_grid<2> (SOURCE_DIR "/grids/grid_in_msh_01.2da.msh"); + gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3d.msh"); + gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3da.msh"); + gmsh_grid<3> (SOURCE_DIR "/grids/grid_in_msh_01.3d_neg.msh"); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + return 0; +} diff --git a/tests/grid/grid_in_msh_02.output b/tests/grid/grid_in_msh_02.output new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/grid/grids/grid_in_msh_02.msh b/tests/grid/grids/grid_in_msh_02.msh new file mode 100644 index 0000000000..3182d01506 --- /dev/null +++ b/tests/grid/grids/grid_in_msh_02.msh @@ -0,0 +1,32 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +11 +1 0 0 0 +2 1 0 0 +3 0.1 0 0 +4 0.2 0 0 +5 0.3 0 0 +6 0.4 0 0 +7 0.5 0 0 +8 0.6 0 0 +9 0.7 0 0 +10 0.8 0 0 +11 0.9 0 0 +$EndNodes +$Elements +12 +1 15 2 10 1 1 +2 15 2 20 2 2 +3 1 2 100 1 1 3 +4 1 2 100 1 3 4 +5 1 2 100 1 4 5 +6 1 2 100 1 5 6 +7 1 2 100 1 6 7 +8 1 2 100 1 7 8 +9 1 2 100 1 8 9 +10 1 2 100 1 9 10 +11 1 2 100 1 10 11 +12 1 2 100 1 11 2 +$EndElements -- 2.39.5