]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement getting and setting boundary indicators of triangulations in 1d as well.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 31 Aug 2011 01:19:35 +0000 (01:19 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 31 Aug 2011 01:19:35 +0000 (01:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@24224 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/grid/tria.h
deal.II/include/deal.II/grid/tria_accessor.h
deal.II/include/deal.II/grid/tria_accessor.templates.h
deal.II/source/grid/tria.cc
tests/deal.II/vertex_as_face_10.cc [new file with mode: 0644]
tests/deal.II/vertex_as_face_10/cmp/generic [new file with mode: 0644]
tests/deal.II/vertex_as_face_11.cc [new file with mode: 0644]
tests/deal.II/vertex_as_face_11/cmp/generic [new file with mode: 0644]

index 5ca92752991c8cd85eaa713aa2c51825be092ac6..7e18a31931e0d6e25fe1f530221ca490372565cb 100644 (file)
@@ -283,6 +283,11 @@ and DoF handlers embedded in higher dimensional space have been added.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: Code like <code>cell-@>face(1)-@>set_boundary_indicator(42);</code>
+now also works in 1d.
+<br>
+(Wolfgang Bangerth, 2011/08/30)
+
 <li> Fixed: The TimerOutput::print_summary() function changed the
 precision of output on the stream it prints to, but didn't restore
 the previous value. This is now fixed.
index ff8aa0c49a5fecd9e935b0109780f2d959152d57..6912b7d3f1feda7a84da38464b8b2bf88608516b 100644 (file)
@@ -38,7 +38,7 @@ template <int dim, int spacedim> class Boundary;
 template <int dim, int spacedim> class StraightBoundary;
 
 template <int, int, int> class TriaAccessor;
-template <int, int, int> class TriaAccessor;
+template <int spacedim> class TriaAccessor<0,1,spacedim>;
 
 namespace internal
 {
@@ -667,7 +667,7 @@ namespace internal
  *   refinement smoothes the mesh a bit.
  *
  *   The function will make sure that vertices on restricted faces (hanging
- *   nodes) will result in the correct place, i.e. in the middle of the two
+ *   nodes) will end up in the correct place, i.e. in the middle of the two
  *   other vertices of the mother line, and the analogue in higher space
  *   dimensions (vertices on the boundary are not corrected, so don't distort
  *   boundary vertices in more than two space dimension, i.e. in dimensions
@@ -945,7 +945,7 @@ namespace internal
  *
  *   <h3>Material and boundary information</h3>
  *
- *   Each line, quad, etc stores one byte of information denoting the
+ *   Each cell, face or edge stores information denoting the
  *   material or the part of the boundary that an object
  *   belongs to. The material of a cell may be used
  *   during matrix generation in order to implement different
@@ -986,16 +986,13 @@ namespace internal
  *   are separated by a vertex (in 2D) or a line (in 3D) such that each boundary
  *   line or quad has a unique boundary indicator.
  *
- *   Since in one dimension, no substructures of lower dimension exist to
- *   cells (of course apart from vertices, but these are handled
- *   in another way than the structures and substructures with dimension one
- *   or greater), there is no way to denote boundary indicators to boundary
- *   vertices (the endpoints). This is not a big thing, however, since you
- *   will normally not want to do a loop over all vertices, but rather work
- *   on the cells, which in turn have a possibility to find out whether they
- *   are at one of the endpoints. Only the handling of boundary values gets
- *   a bit more complicated, but this seems to be the price to be paid for
- *   the different handling of vertices from lines and quads.
+ *   By default (unless otherwise specified during creation of a
+ *   triangulation), all parts of the boundary have boundary indicator
+ *   zero. As a historical wart, this isn't true for 1d meshes, however: For
+ *   these, leftmost vertices have boundary indicator zero while rightmost
+ *   vertices have boundary indicator one. In either case, the boundary
+ *   indicator of a face can be changed using a call of the kind
+ *   <code>cell-@>face(1)-@>set_boundary_indicator(42);</code>.
  *
  *
  *
@@ -3744,6 +3741,31 @@ class Triangulation : public Subscriptor
                                      */
     internal::Triangulation::NumberCache<dim> number_cache;
 
+                                    /**
+                                     * A map that relates the number of a
+                                     * boundary vertex to the boundary
+                                     * indicator. This field is only used in
+                                     * 1d. We have this field because we
+                                     * store boundary indicator information
+                                     * with faces in 2d and higher where we
+                                     * have space in the structures that
+                                     * store data for faces, but in 1d there
+                                     * is no such space for faces.
+                                     *
+                                     * The field is declared as a pointer for
+                                     * a rather mundane reason: all other
+                                     * fields of this class that can be
+                                     * modified by the TriaAccessor hierarchy
+                                     * are pointers, and so these accessor
+                                     * classes store a const pointer to the
+                                     * triangulation. We could no longer do
+                                     * so for TriaAccessor<0,1,spacedim> if
+                                     * this field (that can be modified by
+                                     * TriaAccessor::set_boundary_indicator)
+                                     * were not a pointer.
+                                     */
+    std::map<unsigned int, unsigned char> *vertex_to_boundary_id_map_1d;
+
     /**
      * A map that correlates each refinement listener that has been added
      * through the outdated RefinementListener interface via
@@ -3763,6 +3785,8 @@ class Triangulation : public Subscriptor
                                     // friends
     template <int,int,int> friend class TriaAccessorBase;
     template <int,int,int> friend class TriaAccessor;
+    friend class TriaAccessor<0, 1, spacedim>;
+
     friend class CellAccessor<dim, spacedim>;
 
     friend struct internal::TriaAccessor::Implementation;
@@ -3874,6 +3898,9 @@ Triangulation<dim,spacedim>::save (Archive & ar,
   ar & number_cache;
 
   ar & check_for_distorted_cells;
+
+  if (dim == 1)
+    ar & vertex_to_boundary_id_map_1d;
 }
 
 
@@ -3909,6 +3936,9 @@ Triangulation<dim,spacedim>::load (Archive & ar,
                      "same setting with regard to reporting distorted "
                      "cell as the one previously stored."));
 
+  if (dim == 1)
+    ar & vertex_to_boundary_id_map_1d;
+
                                   // trigger the create signal to indicate
                                   // that new content has been imported into
                                   // the triangulation
index d6564f713450d59847274b4e05b74b884e8184c9..e5ba1977bcba8eddfc3558f6bba7e776aec3f813 100644 (file)
@@ -2020,17 +2020,40 @@ class TriaAccessor<0, 1, spacedim>
     int isotropic_child_index (const unsigned int i);
 
                                     /**
-                                     * @brief Do nothing and throw and error
+                                     * Set the boundary indicator.
+                                     * The same applies as for the
+                                     * <tt>boundary_indicator()</tt>
+                                     * function.
+                                     *
+                                     * @warning You should never set the
+                                     * boundary indicator of an interior face
+                                     * (a face not at the boundary of the
+                                     * domain), or set set the boundary
+                                     * indicator of an exterior face to 255
+                                     * (this value is reserved for another
+                                     * purpose). Algorithms may not work or
+                                     * produce very confusing results if
+                                     * boundary cells have a boundary
+                                     * indicator of 255 or if interior cells
+                                     * have boundary indicators other than
+                                     * 255. Unfortunately, the current object
+                                     * has no means of finding out whether it
+                                     * really is at the boundary of the
+                                     * domain and so cannot determine whether
+                                     * the value you are trying to set makes
+                                     * sense under the current circumstances.
                                      */
-    static
     void
     set_boundary_indicator (const unsigned char);
 
                                     /**
-                                     * @brief Do nothing and throw and error
+                                     * Since this object only represents a
+                                     * single vertex, call
+                                     * set_boundary_indicator with the same
+                                     * argument.
                                      */
-    static
-    void set_all_boundary_indicators (const unsigned char);
+    void
+    set_all_boundary_indicators (const unsigned char);
 
                                     /**
                                      * Return whether the vertex
index 7ac0529a8852497b490c428138a7b316160a42e5..eb0879bbb8a62a2e4c7443a35a0586b95a7d79a3 100644 (file)
@@ -2209,9 +2209,15 @@ TriaAccessor<0, 1, spacedim>::boundary_indicator () const
   switch (vertex_kind)
     {
       case left_vertex:
-           return 0;
       case right_vertex:
-           return 1;
+      {
+       Assert (tria->vertex_to_boundary_id_map_1d->find (this->vertex_index())
+               != tria->vertex_to_boundary_id_map_1d->end(),
+               ExcInternalError());
+
+       return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()];
+      }
+
       default:
            return 255;
     }
@@ -2334,22 +2340,22 @@ int TriaAccessor<0, 1, spacedim>::isotropic_child_index (const unsigned int)
 template <int spacedim>
 inline
 void
-TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const unsigned char)
+TriaAccessor<0, 1, spacedim>::set_boundary_indicator (const unsigned char b)
 {
-  Assert(false,
-        ExcMessage ("The boundary indicator of vertices is determined by their "
-                    "location within a triangulation and can not be set explicitly"));
+  Assert (tria->vertex_to_boundary_id_map_1d->find (this->vertex_index())
+         != tria->vertex_to_boundary_id_map_1d->end(),
+         ExcInternalError());
+
+  (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b;
 }
 
 
 
 template <int spacedim>
 inline
-void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const unsigned char)
+void TriaAccessor<0, 1, spacedim>::set_all_boundary_indicators (const unsigned char b)
 {
-  Assert(false,
-        ExcMessage ("The boundary indicator of vertices is determined by their "
-                    "location within a triangulation and can not be set explicitly"));
+  set_boundary_indicator (b);
 }
 
 
index 15431b627e47f438cd608cdf35bd0598677d4aab..705431cefed8a5db7b2ad1b8e307a0a2700b5f67 100644 (file)
@@ -1570,6 +1570,19 @@ namespace internal
                                lines_at_vertex[line->vertex_index(vertex)][0]);
                    line->set_neighbor (vertex, neighbor);
                  }
+
+                                            // finally set the
+                                            // vertex_to_boundary_id_map_1d
+                                            // map
+           triangulation.vertex_to_boundary_id_map_1d->clear();
+           for (typename Triangulation<dim,spacedim>::active_cell_iterator
+                  cell = triangulation.begin_active();
+                cell != triangulation.end(); ++cell)
+             for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+               if (cell->at_boundary(f))
+                 (*triangulation
+                  .vertex_to_boundary_id_map_1d)[cell->face(f)->vertex_index()]
+                   = f;
          }
 
 
@@ -8944,7 +8957,7 @@ namespace internal
        void
        prevent_distorted_boundary_cells (const Triangulation<1,spacedim> &);
 
-       
+
        template <int dim, int spacedim>
        static
        void
@@ -8954,7 +8967,7 @@ namespace internal
                                             // one, we cannot perform
                                             // this check yet.
            if(spacedim>dim) return;
-           
+
            for (typename Triangulation<dim,spacedim>::cell_iterator
                   cell=triangulation.begin(); cell!=triangulation.end(); ++cell)
              if (cell->at_boundary() &&
@@ -9378,13 +9391,18 @@ Triangulation (const MeshSmoothing smooth_grid,
                smooth_grid(smooth_grid),
                faces(NULL),
                anisotropic_refinement(false),
-               check_for_distorted_cells(check_for_distorted_cells)
+               check_for_distorted_cells(check_for_distorted_cells),
+               vertex_to_boundary_id_map_1d (0)
 {
                                   // set default boundary for all
                                   // possible components
   for (unsigned int i=0;i<255;++i)
     boundary[i] = &straight_boundary;
-  
+
+  if (dim == 1)
+    vertex_to_boundary_id_map_1d
+      = new std::map<unsigned int, unsigned char>();
+
   // connect the any_change signal to the other signals
   signals.create.connect (signals.any_change);
   signals.post_refinement.connect (signals.any_change);
@@ -9400,7 +9418,8 @@ Triangulation (const Triangulation<dim, spacedim> & other)
                                   // is an error!
                :
                Subscriptor(),
-               check_for_distorted_cells(other.check_for_distorted_cells)
+               check_for_distorted_cells(other.check_for_distorted_cells),
+               vertex_to_boundary_id_map_1d (0)
 {
   Assert (false, ExcInternalError());
 }
@@ -9414,6 +9433,14 @@ Triangulation<dim, spacedim>::~Triangulation ()
     delete levels[i];
   levels.clear ();
   delete faces;
+
+                                  // the vertex_to_boundary_id_map_1d field
+                                  // should be unused except in 1d
+  Assert ((dim == 1)
+         ||
+         (vertex_to_boundary_id_map_1d == 0),
+         ExcInternalError());
+  delete vertex_to_boundary_id_map_1d;
 }
 
 
@@ -9472,24 +9499,40 @@ template <int dim, int spacedim>
 std::vector<unsigned char>
 Triangulation<dim, spacedim>::get_boundary_indicators () const
 {
+                                  // in 1d, we store a map of all used
+                                  // boundary indicators. use it for our
+                                  // purposes
+  if (dim == 1)
+    {
+      std::vector<unsigned char> boundary_indicators;
+      for (std::map<unsigned int, unsigned char>::const_iterator
+            p = vertex_to_boundary_id_map_1d->begin();
+          p !=  vertex_to_boundary_id_map_1d->end();
+          ++p)
+       boundary_indicators.push_back (p->second);
+
+      return boundary_indicators;
+    }
+  else
+    {
+      std::vector<bool> bi_exists(255, false);
+      active_cell_iterator cell=begin_active();
+      for (; cell!=end(); ++cell)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         if (cell->at_boundary(face))
+           bi_exists[cell->face(face)->boundary_indicator()]=true;
 
-  std::vector<bool> bi_exists(255, false);
-  active_cell_iterator cell=begin_active();
-  for (; cell!=end(); ++cell)
-    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-      if (cell->at_boundary(face))
-       bi_exists[cell->face(face)->boundary_indicator()]=true;
-
-  const unsigned int n_bi=
-    std::count(bi_exists.begin(), bi_exists.end(), true);
+      const unsigned int n_bi=
+       std::count(bi_exists.begin(), bi_exists.end(), true);
 
-  std::vector<unsigned char> boundary_indicators(n_bi);
-  unsigned int bi_counter=0;
-  for (unsigned int i=0; i<bi_exists.size(); ++i)
-    if (bi_exists[i]==true)
-      boundary_indicators[bi_counter++]=i;
+      std::vector<unsigned char> boundary_indicators(n_bi);
+      unsigned int bi_counter=0;
+      for (unsigned int i=0; i<bi_exists.size(); ++i)
+       if (bi_exists[i]==true)
+         boundary_indicators[bi_counter++]=i;
 
-  return boundary_indicators;
+      return boundary_indicators;
+    }
 }
 
 
@@ -9530,6 +9573,14 @@ copy_triangulation (const Triangulation<dim, spacedim> &old_tria)
 
   number_cache = old_tria.number_cache;
 
+  if (dim == 1)
+    {
+      delete vertex_to_boundary_id_map_1d;
+      vertex_to_boundary_id_map_1d
+       = (new std::map<unsigned int, unsigned char>
+          (*old_tria.vertex_to_boundary_id_map_1d));
+    }
+
                                   // inform RefinementListeners of old_tria of
                                   // the copy operation
   old_tria.signals.copy (*this);
@@ -9594,7 +9645,7 @@ create_triangulation (const std::vector<Point<spacedim> >    &v,
     }
 
   compute_number_cache (*this, levels.size(), number_cache);
-  
+
                                   // now verify that there are indeed
                                   // no distorted cells. as per the
                                   // documentation of this class, we
@@ -12437,7 +12488,7 @@ Triangulation<dim, spacedim>::execute_coarsening_and_refinement ()
                                   // Inform RefinementListeners
                                    // about beginning of refinement.
   signals.pre_refinement();
-  
+
   execute_coarsening();
 
   const DistortedCellList
@@ -14507,8 +14558,8 @@ template<int dim, int spacedim>
 void
 Triangulation<dim, spacedim>::add_refinement_listener (RefinementListener &listener) const
 {
-  // in this compatibility mode with the old-style refinement listeners, an 
-  // external class presents itself as one that may or may not have 
+  // in this compatibility mode with the old-style refinement listeners, an
+  // external class presents itself as one that may or may not have
   // overloaded all of the functions that the RefinementListener
   // class has. consequently, we need to connect each of its functions
   // to the relevant signals. for those functions that haven't been
@@ -14516,25 +14567,25 @@ Triangulation<dim, spacedim>::add_refinement_listener (RefinementListener &liste
   // to the function in the RefinementListener base class which simply
   // does nothing
   std::vector<boost::signals2::connection> connections;
-  
-  connections.push_back 
+
+  connections.push_back
   (signals.create.connect (std_cxx1x::bind (&RefinementListener::create_notification,
                                            std_cxx1x::ref(listener),
                                            std_cxx1x::cref(*this))));
-  connections.push_back 
+  connections.push_back
   (signals.copy.connect (std_cxx1x::bind (&RefinementListener::copy_notification,
                                          std_cxx1x::ref(listener),
                                          std_cxx1x::cref(*this),
                                          std_cxx1x::_1)));
-  connections.push_back 
+  connections.push_back
   (signals.pre_refinement.connect (std_cxx1x::bind (&RefinementListener::pre_refinement_notification,
                                                    std_cxx1x::ref(listener),
                                                    std_cxx1x::cref(*this))));
-  connections.push_back 
+  connections.push_back
   (signals.post_refinement.connect (std_cxx1x::bind (&RefinementListener::post_refinement_notification,
                                                     std_cxx1x::ref(listener),
                                                    std_cxx1x::cref(*this))));
-  
+
   // now push the set of connections into the map
   refinement_listener_map.insert (std::make_pair(&listener, connections));
 }
@@ -14548,14 +14599,14 @@ Triangulation<dim, spacedim>::remove_refinement_listener (RefinementListener &li
   Assert (refinement_listener_map.find (&listener) != refinement_listener_map.end(),
          ExcMessage("You try to remove a refinement listener that does "
            "not appear to have been added previously."));
-  
+
   // get the element of the map, and terminate these
   // connections. then erase the element from the list
   std::vector<boost::signals2::connection> connections
   = refinement_listener_map.find(&listener)->second;
   for (unsigned int i=0; i<connections.size(); ++i)
     connections[i].disconnect ();
-  
+
   refinement_listener_map.erase (refinement_listener_map.find (&listener));
 }
 
diff --git a/tests/deal.II/vertex_as_face_10.cc b/tests/deal.II/vertex_as_face_10.cc
new file mode 100644 (file)
index 0000000..0f2f572
--- /dev/null
@@ -0,0 +1,55 @@
+//----------------------------  vertex_as_face_10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  vertex_as_face_10.cc  ---------------------------
+
+// verify that we can do things like cell->face() in 1d as well. here:
+// setting boundary indicators
+
+
+#include "../tests.h"
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+
+
+template <int spacedim>
+void test ()
+{
+  Triangulation<1,spacedim> tria;
+  GridGenerator::hyper_cube (tria);
+
+  tria.begin_active()->face(0)->set_boundary_indicator(2);
+  tria.begin_active()->face(1)->set_boundary_indicator(4);
+
+  std::vector<unsigned char>
+    boundary_ids = tria.get_boundary_indicators ();
+
+  for (unsigned int i=0; i<boundary_ids.size(); ++i)
+    deallog << (int)boundary_ids[i] << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("vertex_as_face_10/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1> ();
+  test<2> ();
+
+  return 0;
+}
diff --git a/tests/deal.II/vertex_as_face_10/cmp/generic b/tests/deal.II/vertex_as_face_10/cmp/generic
new file mode 100644 (file)
index 0000000..a9927d4
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::2
+DEAL::4
+DEAL::2
+DEAL::4
diff --git a/tests/deal.II/vertex_as_face_11.cc b/tests/deal.II/vertex_as_face_11.cc
new file mode 100644 (file)
index 0000000..534b552
--- /dev/null
@@ -0,0 +1,52 @@
+//----------------------------  vertex_as_face_11.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2010, 2011 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  vertex_as_face_11.cc  ---------------------------
+
+// verify that we can do things like cell->face() in 1d as well. here:
+// getting boundary indicators
+
+
+#include "../tests.h"
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <fstream>
+
+
+template <int spacedim>
+void test ()
+{
+  Triangulation<1,spacedim> tria;
+  GridGenerator::hyper_cube (tria);
+
+  tria.begin_active()->face(0)->set_boundary_indicator(2);
+  tria.begin_active()->face(1)->set_boundary_indicator(4);
+
+  deallog << (int)tria.begin_active()->face(0)->boundary_indicator() << std::endl;
+  deallog << (int)tria.begin_active()->face(1)->boundary_indicator() << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("vertex_as_face_11/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1> ();
+  test<2> ();
+
+  return 0;
+}
diff --git a/tests/deal.II/vertex_as_face_11/cmp/generic b/tests/deal.II/vertex_as_face_11/cmp/generic
new file mode 100644 (file)
index 0000000..a9927d4
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::2
+DEAL::4
+DEAL::2
+DEAL::4

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.