New assign_co_dimensional_manifold_indicators function.

diff --git a/doc/news/changes/minor/20190319LucaHeltai_b b/doc/news/changes/minor/20190319LucaHeltai_b
new file mode 100644 (file)
index 0000000..c41dc60
--- /dev/null
+++ b/doc/news/changes/minor/20190319LucaHeltai_b
@@ -0,0 +1,5 @@
+New: The new method GridTools::assign_co_dimensional_manifold_indicators() propagate manifold ids from 
+cells to faces and edges, according to a disambiguation function that takes the set of manifold ids of
+the cells that share the given face or edge. 
+<br>
+(Luca Heltai, 2019/03/19)

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index df3711d54f44d8e3ce0329e6753d0ce6a2167517..5a341b9e9f2a517433eb6d6544783cc6aee83b49 100644 (file)
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -2632,7 +2632,36 @@ namespace GridTools
   copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                const bool compute_face_ids = false);
 
-
+  /**
+   * Propagate manifold indicators associated to the cells of the Triangulation
+   * @p tria to their co-dimension one and two objects.
+   *
+   * This function sets the @p manifold_id of shared faces and edges to the value
+   * returnd by the @p disambiguation_function method, called with the set of
+   * manifold indicators of the cells that share the same face or edge.
+   *
+   * By default, the @p disambiguation_function returns
+   * numbers::flat_manifold_id when the set has size greater than one (i.e.,
+   * when it is not possible to decide what manifold indicator a face or edge
+   * should have according to the manifold indicators of the adjacent cells)
+   * and it returns the manifold indicator contained in the set when it has
+   * dimension one (i.e., when all adjacent cells and faces have the same
+   * manifold indicator).
+   *
+   * @author Luca Heltai, 2019.
+   */
+  template <int dim, int spacedim>
+  void
+  assign_co_dimensional_manifold_indicators(
+    Triangulation<dim, spacedim> &            tria,
+    const std::function<types::manifold_id(
+      const std::set<types::manifold_id> &)> &disambiguation_function =
+      [](const std::set<types::manifold_id> &manifold_ids) {
+        if (manifold_ids.size() == 1)
+          return *manifold_ids.begin();
+        else
+          return numbers::invalid_manifold_id;
+      });
   /*@}*/
 
   /**

diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h
index 9698bb6038e02d22dceb1979f4bdc1b79feb5a36..47bd2399170ff2f4e429b57507a70a41d8c5fa7f 100644 (file)
--- a/include/deal.II/grid/tria_accessor.h
+++ b/include/deal.II/grid/tria_accessor.h
@@ -632,6 +632,24 @@ public:
   types::manifold_id
   manifold_id() const;
 
+  /**
+   * Dummy function that always returns numbers::invalid_unsigned_int.
+   */
+  unsigned int
+  user_index() const;
+
+  /**
+   * Dummy function that always throws.
+   */
+  void
+  set_user_index(const unsigned int p) const;
+
+  /**
+   * Dummy function that always throws.
+   */
+  void
+  set_manifold_id(const types::manifold_id) const;
+
   /**
    * Dummy function to extract vertices. Returns the origin.
    */

diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h
index bd9bc7b43cd2a786b32cfa53147ac243a3641fd6..3f74a7fdbc1130ae67cccb6e03ccb49152c7fe36 100644 (file)
--- a/include/deal.II/grid/tria_accessor.templates.h
+++ b/include/deal.II/grid/tria_accessor.templates.h
@@ -520,6 +520,40 @@ InvalidAccessor<structdim, dim, spacedim>::manifold_id() const
 }
 
 
+
+template <int structdim, int dim, int spacedim>
+unsigned int
+InvalidAccessor<structdim, dim, spacedim>::user_index() const
+{
+  return numbers::invalid_unsigned_int;
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+void
+InvalidAccessor<structdim, dim, spacedim>::set_user_index(
+  const unsigned int) const
+{
+  Assert(false,
+         ExcMessage("You are trying to set the user index of an "
+                    "invalid object."));
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+void
+InvalidAccessor<structdim, dim, spacedim>::set_manifold_id(
+  const types::manifold_id) const
+{
+  Assert(false,
+         ExcMessage("You are trying to set the manifold id of an "
+                    "invalid object."));
+}
+
+
+
 template <int structdim, int dim, int spacedim>
 inline Point<spacedim> &
 InvalidAccessor<structdim, dim, spacedim>::vertex(const unsigned int) const

diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index 44d49753d1661040521dd2e23e78e528c7fb3e14..219159c1facd5f3a64e6255f913b2f2e6597e8dc 100644 (file)
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -3756,6 +3756,79 @@ namespace GridTools
       }
   }
 
+
+  template <int dim, int spacedim>
+  void
+  assign_co_dimensional_manifold_indicators(
+    Triangulation<dim, spacedim> &            tria,
+    const std::function<types::manifold_id(
+      const std::set<types::manifold_id> &)> &disambiguation_function)
+  {
+    // Easy case first:
+    if (dim == 1)
+      return;
+    const unsigned int size =
+      dim == 2 ? tria.n_lines() : tria.n_lines() + tria.n_quads();
+
+    // If user index is zero, then it has not been set.
+    std::vector<std::set<types::manifold_id>> manifold_ids(size + 1);
+    std::vector<unsigned int>                 backup;
+    tria.save_user_indices(backup);
+
+    unsigned next_index = 1;
+    for (auto cell : tria.active_cell_iterators())
+      {
+        if (dim > 1)
+          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+            {
+              if (cell->line(l)->user_index() == 0)
+                {
+                  AssertIndexRange(next_index, size + 1);
+                  manifold_ids[next_index].insert(cell->manifold_id());
+                  cell->line(l)->set_user_index(next_index++);
+                }
+              else
+                manifold_ids[cell->line(l)->user_index()].insert(
+                  cell->manifold_id());
+            }
+        if (dim > 2)
+          for (unsigned int l = 0; l < GeometryInfo<dim>::quads_per_cell; ++l)
+            {
+              if (cell->quad(l)->user_index() == 0)
+                {
+                  AssertIndexRange(next_index, size + 1);
+                  manifold_ids[next_index].insert(cell->manifold_id());
+                  cell->quad(l)->set_user_index(next_index++);
+                }
+              else
+                manifold_ids[cell->quad(l)->user_index()].insert(
+                  cell->manifold_id());
+            }
+      }
+    for (auto cell : tria.active_cell_iterators())
+      {
+        if (dim > 1)
+          for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+            {
+              const auto id = cell->line(l)->user_index();
+              Assert(id != 0, ExcInternalError());
+              cell->line(l)->set_manifold_id(
+                disambiguation_function(manifold_ids[id]));
+            }
+        if (dim > 2)
+          for (unsigned int l = 0; l < GeometryInfo<dim>::quads_per_cell; ++l)
+            {
+              const auto id = cell->quad(l)->user_index();
+              Assert(id != 0, ExcInternalError());
+              cell->quad(l)->set_manifold_id(
+                disambiguation_function(manifold_ids[id]));
+            }
+      }
+    tria.load_user_indices(backup);
+  }
+
+
+
   template <int dim, int spacedim>
   std::pair<unsigned int, double>
   get_longest_direction(

diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index 23dd04b907b6debaae230217d13a17a10d5c1b03..4d04f4141da7a93861cdc3c85970800d902f9426 100644 (file)
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -375,6 +375,12 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         Triangulation<deal_II_dimension, deal_II_space_dimension> &,
         const std::vector<types::boundary_id> &);
 
+      template void
+      assign_co_dimensional_manifold_indicators(
+        Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+        const std::function<
+          types::manifold_id(const std::set<types::manifold_id> &)> &);
+
       template void
       regularize_corner_cells(
         Triangulation<deal_II_dimension, deal_II_space_dimension> &,

diff --git a/tests/grid/grid_tools_propagate_manifold_ids_01.cc b/tests/grid/grid_tools_propagate_manifold_ids_01.cc
new file mode 100644 (file)
index 0000000..d6ab8c2
--- /dev/null
+++ b/tests/grid/grid_tools_propagate_manifold_ids_01.cc
@@ -0,0 +1,74 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+// Validate GridTools::assign_co_dimensional_manifold_indicators
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+template <int dim, int spacedim>
+void
+test()
+{
+  deallog << "dim = " << dim << ", spacedim = " << spacedim << std::endl;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+  tria.begin_active()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  tria.begin_active()->set_manifold_id(1);
+
+  auto validate = [](const std::set<types::manifold_id> &set) {
+    if (set.size() > 1)
+      return 1;
+    else
+      return 0;
+  };
+
+  GridTools::assign_co_dimensional_manifold_indicators(tria, validate);
+  for (auto cell : tria.active_cell_iterators())
+    {
+      deallog << "Cell: " << (int)cell->manifold_id() << std::endl;
+      if (dim > 1)
+        for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
+          deallog << "Line " << l << ", " << (int)cell->line(l)->manifold_id()
+                  << std::endl;
+      if (dim > 2)
+        for (unsigned int l = 0; l < GeometryInfo<dim>::quads_per_cell; ++l)
+          deallog << "Quad " << l << ", " << (int)cell->quad(l)->manifold_id()
+                  << std::endl;
+    }
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<1, 1>();
+  test<1, 2>();
+  test<2, 2>();
+  test<2, 3>();
+  test<3, 3>();
+
+  return 0;
+}

diff --git a/tests/grid/grid_tools_propagate_manifold_ids_01.output b/tests/grid/grid_tools_propagate_manifold_ids_01.output
new file mode 100644 (file)
index 0000000..8b385d5
--- /dev/null
+++ b/tests/grid/grid_tools_propagate_manifold_ids_01.output
@@ -0,0 +1,367 @@
+
+DEAL::dim = 1, spacedim = 1
+DEAL::Cell: 1
+DEAL::Cell: -1
+DEAL::Cell: -1
+DEAL::dim = 1, spacedim = 2
+DEAL::Cell: 1
+DEAL::Cell: -1
+DEAL::Cell: -1
+DEAL::dim = 2, spacedim = 2
+DEAL::Cell: 1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 1
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 1
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::dim = 2, spacedim = 3
+DEAL::Cell: 1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 1
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 1
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::dim = 3, spacedim = 3
+DEAL::Cell: 1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 1
+DEAL::Line 4, 1
+DEAL::Line 5, 1
+DEAL::Line 6, 1
+DEAL::Line 7, 1
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 1
+DEAL::Line 11, 1
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 1
+DEAL::Quad 4, 0
+DEAL::Quad 5, 1
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 1
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 1
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 1
+DEAL::Line 7, 0
+DEAL::Line 8, 1
+DEAL::Line 9, 1
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 1
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 1
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 1
+DEAL::Line 1, 1
+DEAL::Line 2, 1
+DEAL::Line 3, 1
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 1
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 1
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0
+DEAL::Cell: -1
+DEAL::Line 0, 0
+DEAL::Line 1, 0
+DEAL::Line 2, 0
+DEAL::Line 3, 0
+DEAL::Line 4, 0
+DEAL::Line 5, 0
+DEAL::Line 6, 0
+DEAL::Line 7, 0
+DEAL::Line 8, 0
+DEAL::Line 9, 0
+DEAL::Line 10, 0
+DEAL::Line 11, 0
+DEAL::Quad 0, 0
+DEAL::Quad 1, 0
+DEAL::Quad 2, 0
+DEAL::Quad 3, 0
+DEAL::Quad 4, 0
+DEAL::Quad 5, 0