]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Set boundary indicators in GridIn on vertices in 1d.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Aug 2015 20:16:11 +0000 (15:16 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Aug 2015 21:14:23 +0000 (16:14 -0500)
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
source/grid/grid_in.cc
tests/deal.II/grid_in_msh_02.cc
tests/deal.II/grid_in_msh_02.output
tests/grid/grid_in_msh_02.cc [new file with mode: 0644]
tests/grid/grid_in_msh_02.output [new file with mode: 0644]
tests/grid/grids/grid_in_msh_02.msh [new file with mode: 0644]

index 5965319d368d949b0f1f12ed5752c112d6ebafa8..c3f7f5cff3407694f14de8e96893e24ea9c7cbf6 100644 (file)
@@ -108,6 +108,12 @@ inconvenience this causes.
 
 
 <ol>
+  <li> Fixed: In 1d, GridIn::read_msh() ignored boundary indicators
+  associated with vertices. This is now fixed.
+  <br>
+  (Jan Strebel, Wolfgang Bangerth, 2015/08/25)
+  </li>
+
   <li> Improved: The interface and documentation for periodic boundary
   conditions have been restructured. A
   @ref GlossPeriodicConstraints "glossary entry" has been written.
index a516dd1e32d7f14f541838cc656fa9cc5f4e9fff..29c131f0fa4921cbc022c1006783bb79c5d15b76 100644 (file)
 
 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 <int spacedim>
+  void
+  assign_1d_boundary_indicators (const std::map<unsigned int, types::boundary_id> &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<GeometryInfo<1>::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 <int dim, int spacedim>
+  void
+  assign_1d_boundary_indicators (const std::map<unsigned int, types::boundary_id> &,
+                                 Triangulation<dim,spacedim> &)
+  {
+    // we shouldn't get here since boundary ids are not assigned to
+    // vertices except in 1d
+    Assert (dim != 1, ExcInternalError());
+  }
+}
+
 template <int dim, int spacedim>
 GridIn<dim, spacedim>::GridIn () :
   tria(0, typeid(*this).name()), default_format(ucd)
@@ -1244,9 +1287,12 @@ void GridIn<dim, spacedim>::read_msh (std::istream &in)
 
   in >> n_cells;
 
-  // set up array of cells
-  std::vector<CellData<dim> > 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<CellData<dim> >                cells;
+  SubCellData                                subcelldata;
+  std::map<unsigned int, types::boundary_id> boundary_ids_1d;
 
   for (unsigned int cell=0; cell<n_cells; ++cell)
     {
@@ -1289,8 +1335,9 @@ void GridIn<dim, spacedim>::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<dim, spacedim>::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<nod_num; ++i)
-                in >> 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<dim, spacedim>::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<dim, spacedim>::read_msh (std::istream &in)
     GridReordering<dim,spacedim>::invert_all_cells_of_negative_grid (vertices, cells);
   GridReordering<dim,spacedim>::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);
 }
 
 
index d3144ce76e28f3bc3283b2b917981dfc641e2bec..cfa014f4e34753a085fcafbdeefce12f9667aa1d 100644 (file)
 // ---------------------------------------------------------------------
 
 
-// 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 <deal.II/dofs/dof_handler.h>
 std::ofstream logfile("output");
 
 
-template<int dim>
-void check_file (const std::string name,
-                 typename GridIn<dim>::Format format)
+void check_file ()
 {
-  Triangulation<dim> tria;
-  GridIn<dim> 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 ();
 }
 
index 81a9950edde4dd9af4f01d32afb09d06caaafc26..2ebba38d85c8be09b42d3625ba12bf44f2c3f188 100644 (file)
@@ -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 (file)
index 0000000..b615e5c
--- /dev/null
@@ -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 <deal.II/base/logstream.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_in.h>
+
+#include <fstream>
+#include <cmath>
+
+
+template <int dim>
+void gmsh_grid (const char *name)
+{
+  Triangulation<dim> tria;
+  GridIn<dim> 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<dim>::active_cell_iterator c=tria.begin_active();
+       c!=tria.end(); ++c, ++index)
+    for (unsigned int i=0; i<GeometryInfo<dim>::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 (file)
index 0000000..e69de29
diff --git a/tests/grid/grids/grid_in_msh_02.msh b/tests/grid/grids/grid_in_msh_02.msh
new file mode 100644 (file)
index 0000000..3182d01
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.