From: Wolfgang Bangerth Date: Tue, 25 Aug 2015 20:16:11 +0000 (-0500) Subject: Set boundary indicators in GridIn on vertices in 1d. X-Git-Tag: v8.4.0-rc2~548^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e0a3d7b5a784122bc458429a915490d90626bcf;p=dealii.git 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. --- 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