]> https://gitweb.dealii.org/ - dealii.git/commitdiff
make_flux_sparsity_pattern with flux couplings on a subset of the faces 9407/head
authorSimon Sticko <simon@sticko.se>
Tue, 14 Jan 2020 10:29:15 +0000 (11:29 +0100)
committerSimon Sticko <simon@sticko.se>
Wed, 29 Jan 2020 08:11:33 +0000 (09:11 +0100)
Add an optional std::function input parameter to
make_flux_sparsity_pattern which describe over which faces there should
be flux couplings in the created sparsity pattern.

internal::make_flux_sparsity_pattern has one implementation for a
generic DoFHandlerType and a specialization for hp::DoFHandler. To
circument a bug where MSVC can't deduce which of these internal
functions should be called, change DoFHandlerType to
dealii::DoFHandler in the generic function.

To avoid another bug where MSVC can't construct a templated class,
call one of the overloaded versions of make_flux_sparsity_pattern with
always_couple_on_faces (even if this is already the default value).

doc/news/changes/minor/20200128SimonSticko [new file with mode: 0644]
include/deal.II/dofs/dof_tools.h
source/dofs/dof_tools_sparsity.cc
source/dofs/dof_tools_sparsity.inst.in
tests/dofs/dof_tools_25.cc [new file with mode: 0644]
tests/dofs/dof_tools_25.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20200128SimonSticko b/doc/news/changes/minor/20200128SimonSticko
new file mode 100644 (file)
index 0000000..9081573
--- /dev/null
@@ -0,0 +1,5 @@
+Add an optional std::function parameter to make_flux_sparsity_pattern which can
+be used to specify over which faces there should be a flux coupling in the
+created sparsity pattern.
+<br>
+(Simon Sticko, 2020/01/28)
index b5a51ae3cf5edbbf7c24c0a8cacb55453ac9ed5c..f83c001f6a33f6906418e6a095dea805a74712f6 100644 (file)
@@ -68,6 +68,28 @@ namespace GridTools
   template <typename CellIterator>
   struct PeriodicFacePair;
 }
+
+namespace DoFTools
+{
+  namespace internal
+  {
+    /*
+     * Default value of the face_has_flux_coupling parameter of
+     * make_flux_sparsity_pattern. Defined here (instead of using a default
+     * lambda in the parameter list) to avoid a bug in gcc where the same lambda
+     * gets defined multiple times.
+     */
+    template <typename DoFHandlerType>
+    inline bool
+    always_couple_on_faces(
+      const typename DoFHandlerType::active_cell_iterator &,
+      const unsigned int)
+    {
+      return true;
+    }
+  } // namespace internal
+} // namespace DoFTools
+
 #endif
 
 /**
@@ -775,18 +797,40 @@ namespace DoFTools
    * components of a finite element are continuous and some discontinuous,
    * allowing constraints to be imposed on the continuous part while also
    * building the flux terms needed for the discontinuous part.
+   *
+   * The optional @param face_has_flux_coupling can be used to specify on which
+   * faces flux couplings occur. This allows for creating a sparser pattern when
+   * using a bilinear form where flux terms only appear on a subset of the faces
+   * in the triangulation. By default flux couplings are added over all internal
+   * faces. @param face_has_flux_coupling should be a function that takes an
+   * active_cell_iterator and a face index and should return true if there is a
+   * flux coupling over the face. When using the ::dealii::DoFHandler we could,
+   * for example, use
+   *
+   * @code
+   *  auto face_has_flux_coupling =
+   *    [](const typename DoFHandler<dim>::active_cell_iterator &cell,
+   *       const unsigned int                                    face_index) {
+   *      const Point<dim> &face_center = cell->face(face_index)->center();
+   *      return 0 < face_center[0];
+   *    };
    */
   template <typename DoFHandlerType,
             typename SparsityPatternType,
             typename number>
   void
-  make_flux_sparsity_pattern(const DoFHandlerType &           dof,
-                             SparsityPatternType &            sparsity,
-                             const AffineConstraints<number> &constraints,
-                             const bool                keep_constrained_dofs,
-                             const Table<2, Coupling> &couplings,
-                             const Table<2, Coupling> &face_couplings,
-                             const types::subdomain_id subdomain_id);
+  make_flux_sparsity_pattern(
+    const DoFHandlerType &           dof,
+    SparsityPatternType &            sparsity,
+    const AffineConstraints<number> &constraints,
+    const bool                       keep_constrained_dofs,
+    const Table<2, Coupling> &       couplings,
+    const Table<2, Coupling> &       face_couplings,
+    const types::subdomain_id        subdomain_id,
+    const std::function<
+      bool(const typename DoFHandlerType::active_cell_iterator &,
+           const unsigned int)> &face_has_flux_coupling =
+      &internal::always_couple_on_faces<DoFHandlerType>);
 
   /**
    * Create the sparsity pattern for boundary matrices. See the general
index b3c05e9e931e0a265235d77a5eebf05637b88734..a4c3ddc967d964f3facac123ed43ff929b1d2f6c 100644 (file)
@@ -759,20 +759,24 @@ namespace DoFTools
     {
       // implementation of the same function in namespace DoFTools for
       // non-hp DoFHandlers
-      template <typename DoFHandlerType,
+      template <int dim,
+                int spacedim,
                 typename SparsityPatternType,
                 typename number>
       void
-      make_flux_sparsity_pattern(const DoFHandlerType &           dof,
-                                 SparsityPatternType &            sparsity,
-                                 const AffineConstraints<number> &constraints,
-                                 const bool keep_constrained_dofs,
-                                 const Table<2, Coupling> &int_mask,
-                                 const Table<2, Coupling> &flux_mask,
-                                 const types::subdomain_id subdomain_id)
+      make_flux_sparsity_pattern(
+        const DoFHandler<dim, spacedim> &dof,
+        SparsityPatternType &            sparsity,
+        const AffineConstraints<number> &constraints,
+        const bool                       keep_constrained_dofs,
+        const Table<2, Coupling> &       int_mask,
+        const Table<2, Coupling> &       flux_mask,
+        const types::subdomain_id        subdomain_id,
+        const std::function<
+          bool(const typename DoFHandler<dim, spacedim>::active_cell_iterator &,
+               const unsigned int)> &face_has_flux_coupling)
       {
-        const FiniteElement<DoFHandlerType::dimension,
-                            DoFHandlerType::space_dimension> &fe = dof.get_fe();
+        const FiniteElement<dim, spacedim> &fe = dof.get_fe();
 
         std::vector<types::global_dof_index> dofs_on_this_cell(
           fe.dofs_per_cell);
@@ -783,13 +787,10 @@ namespace DoFTools
           int_dof_mask  = dof_couplings_from_component_couplings(fe, int_mask),
           flux_dof_mask = dof_couplings_from_component_couplings(fe, flux_mask);
 
-        Table<2, bool> support_on_face(
-          fe.dofs_per_cell,
-          GeometryInfo<DoFHandlerType::dimension>::faces_per_cell);
+        Table<2, bool> support_on_face(fe.dofs_per_cell,
+                                       GeometryInfo<dim>::faces_per_cell);
         for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
-          for (unsigned int f = 0;
-               f < GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
-               ++f)
+          for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
             support_on_face(i, f) = fe.has_support_on_face(i, f);
 
         // Convert the int_dof_mask to bool_int_dof_mask so we can pass it
@@ -801,8 +802,9 @@ namespace DoFTools
             if (int_dof_mask(i, j) != none)
               bool_int_dof_mask(i, j) = true;
 
-        typename DoFHandlerType::active_cell_iterator cell = dof.begin_active(),
-                                                      endc = dof.end();
+        typename DoFHandler<dim, spacedim>::active_cell_iterator
+          cell = dof.begin_active(),
+          endc = dof.end();
         for (; cell != endc; ++cell)
           if (((subdomain_id == numbers::invalid_subdomain_id) ||
                (subdomain_id == cell->subdomain_id())) &&
@@ -816,12 +818,11 @@ namespace DoFTools
                                                       bool_int_dof_mask);
               // Loop over all interior neighbors
               for (unsigned int face_n = 0;
-                   face_n <
-                   GeometryInfo<DoFHandlerType::dimension>::faces_per_cell;
+                   face_n < GeometryInfo<dim>::faces_per_cell;
                    ++face_n)
                 {
-                  const typename DoFHandlerType::face_iterator cell_face =
-                    cell->face(face_n);
+                  const typename DoFHandler<dim, spacedim>::face_iterator
+                    cell_face = cell->face(face_n);
 
                   const bool periodic_neighbor =
                     cell->has_periodic_neighbor(face_n);
@@ -846,8 +847,11 @@ namespace DoFTools
                     }
                   else
                     {
-                      typename DoFHandlerType::level_cell_iterator neighbor =
-                        cell->neighbor_or_periodic_neighbor(face_n);
+                      if (!face_has_flux_coupling(cell, face_n))
+                        continue;
+
+                      typename DoFHandler<dim, spacedim>::level_cell_iterator
+                        neighbor = cell->neighbor_or_periodic_neighbor(face_n);
                       // If the cells are on the same level (and both are
                       // active, locally-owned cells) then only add to the
                       // sparsity pattern if the current cell is 'greater' in
@@ -886,7 +890,7 @@ namespace DoFTools
                       // if (neighbor->has_children()) section below. We need to
                       // do this since we otherwise iterate over the children of
                       // the face, which are always 0 in 1D.
-                      if (DoFHandlerType::dimension == 1)
+                      if (dim == 1)
                         while (neighbor->has_children())
                           neighbor = neighbor->child(face_n == 0 ? 1 : 0);
 
@@ -896,8 +900,8 @@ namespace DoFTools
                                sub_nr != cell_face->n_children();
                                ++sub_nr)
                             {
-                              const typename DoFHandlerType::level_cell_iterator
-                                sub_neighbor =
+                              const typename DoFHandler<dim, spacedim>::
+                                level_cell_iterator sub_neighbor =
                                   periodic_neighbor ?
                                     cell->periodic_neighbor_child_on_subface(
                                       face_n, sub_nr) :
@@ -1069,7 +1073,11 @@ namespace DoFTools
         const bool                                   keep_constrained_dofs,
         const Table<2, Coupling> &                   int_mask,
         const Table<2, Coupling> &                   flux_mask,
-        const types::subdomain_id                    subdomain_id)
+        const types::subdomain_id                    subdomain_id,
+        const std::function<
+          bool(const typename dealii::hp::DoFHandler<dim, spacedim>::
+                 active_cell_iterator &,
+               const unsigned int)> &face_has_flux_coupling)
       {
         // while the implementation above is quite optimized and caches a
         // lot of data (see e.g. the int/flux_dof_mask tables), this is no
@@ -1156,6 +1164,9 @@ namespace DoFTools
                         level_cell_iterator neighbor =
                           cell->neighbor_or_periodic_neighbor(face);
 
+                      if (!face_has_flux_coupling(cell, face))
+                        continue;
+
                       // Like the non-hp case: If the cells are on the same
                       // level (and both are active, locally-owned cells) then
                       // only add to the sparsity pattern if the current cell
@@ -1336,26 +1347,34 @@ namespace DoFTools
 
     const bool keep_constrained_dofs = true;
 
-    make_flux_sparsity_pattern(dof,
-                               sparsity,
-                               dummy,
-                               keep_constrained_dofs,
-                               int_mask,
-                               flux_mask,
-                               subdomain_id);
+    make_flux_sparsity_pattern(
+      dof,
+      sparsity,
+      dummy,
+      keep_constrained_dofs,
+      int_mask,
+      flux_mask,
+      subdomain_id,
+      internal::always_couple_on_faces<DoFHandlerType>);
   }
 
+
+
   template <typename DoFHandlerType,
             typename SparsityPatternType,
             typename number>
   void
-  make_flux_sparsity_pattern(const DoFHandlerType &           dof,
-                             SparsityPatternType &            sparsity,
-                             const AffineConstraints<number> &constraints,
-                             const bool                keep_constrained_dofs,
-                             const Table<2, Coupling> &int_mask,
-                             const Table<2, Coupling> &flux_mask,
-                             const types::subdomain_id subdomain_id)
+  make_flux_sparsity_pattern(
+    const DoFHandlerType &           dof,
+    SparsityPatternType &            sparsity,
+    const AffineConstraints<number> &constraints,
+    const bool                       keep_constrained_dofs,
+    const Table<2, Coupling> &       int_mask,
+    const Table<2, Coupling> &       flux_mask,
+    const types::subdomain_id        subdomain_id,
+    const std::function<
+      bool(const typename DoFHandlerType::active_cell_iterator &,
+           const unsigned int)> &face_has_flux_coupling)
   {
     // do the error checking and frame code here, and then pass on to more
     // specialized functions in the internal namespace
@@ -1390,16 +1409,22 @@ namespace DoFTools
              "associated DoF handler objects, asking for any subdomain other "
              "than the locally owned one does not make sense."));
 
+    Assert(
+      face_has_flux_coupling,
+      ExcMessage(
+        "The function which specifies if a flux coupling occurs over a given "
+        "face is empty."));
+
     internal::make_flux_sparsity_pattern(dof,
                                          sparsity,
                                          constraints,
                                          keep_constrained_dofs,
                                          int_mask,
                                          flux_mask,
-                                         subdomain_id);
+                                         subdomain_id,
+                                         face_has_flux_coupling);
   }
 
-
 } // end of namespace DoFTools
 
 
index f73035578f86fe2d8688451c1942558b74d8fef7..cce112d35e1f2a276ff59d3f2d564ea8532c3da7 100644 (file)
@@ -194,7 +194,10 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS;
       const bool,
       const Table<2, Coupling> &,
       const Table<2, Coupling> &,
-      const types::subdomain_id);
+      const types::subdomain_id,
+      const std::function<bool(
+        const typename DoFHandler<deal_II_dimension>::active_cell_iterator &,
+        const unsigned int)> &);
 
     template void
     DoFTools::make_flux_sparsity_pattern<DoFHandler<deal_II_dimension>, SP>(
@@ -220,7 +223,10 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS;
       const bool,
       const Table<2, Coupling> &,
       const Table<2, Coupling> &,
-      const types::subdomain_id);
+      const types::subdomain_id,
+      const std::function<bool(const typename hp::DoFHandler<
+                                 deal_II_dimension>::active_cell_iterator &,
+                               const unsigned int)> &);
 
 #if deal_II_dimension < 3
 
diff --git a/tests/dofs/dof_tools_25.cc b/tests/dofs/dof_tools_25.cc
new file mode 100644 (file)
index 0000000..7f67b75
--- /dev/null
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+// Test the version of make_flux_sparsity_pattern that takes an additional
+// std::function describing over which faces flux couplings occur.
+//
+// Set up a 1D-triangulation with 3 cells (left, middle and right) where the
+// face between left and middle is located at x=0:
+//
+// |-L-|-M-|-R-|
+//    x=0
+//
+// Distribute FE_Q<1>(1) elements over this grid. Then call
+// make_flux_sparsity_pattern with an std::function specifying that there should
+// only be a flux coupling between the left and middle cell. Verify that this
+// is the case for the created sparsity pattern.
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/dof_handler.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+
+#include "../tests.h"
+
+
+// Call the make_flux_sparsity_pattern with a lambda specifying that only the
+// face at x = 0 should have a flux coupling. Print the constructed sparsity
+// pattern to deallog.
+template <int dim, class DoFHandlerType>
+void
+create_and_print_pattern(const DoFHandlerType &dof_handler)
+{
+  DynamicSparsityPattern dynamic_pattern(dof_handler.n_dofs());
+
+  auto face_has_flux_coupling =
+    [](const typename DoFHandlerType::active_cell_iterator &cell,
+       const unsigned int                                   face_index) {
+      // Only add a flux coupling if the face is at x = 0.
+      const Point<dim> &center = cell->face(face_index)->center();
+      return std::abs(center[0]) < 1e-3;
+      return true;
+    };
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  const bool keep_constrained_dofs = true;
+
+  // Set all couplings to be connected.
+  const unsigned int           n_components = 1;
+  Table<2, DoFTools::Coupling> couplings(n_components, n_components);
+  couplings[0][0] = DoFTools::Coupling::always;
+  Table<2, DoFTools::Coupling> face_couplings(n_components, n_components);
+  face_couplings[0][0] = DoFTools::Coupling::always;
+
+  const types::subdomain_id subdomain_id = 0;
+
+  DoFTools::make_flux_sparsity_pattern(dof_handler,
+                                       dynamic_pattern,
+                                       constraints,
+                                       keep_constrained_dofs,
+                                       couplings,
+                                       face_couplings,
+                                       subdomain_id,
+                                       face_has_flux_coupling);
+
+
+  dynamic_pattern.print(deallog.get_file_stream());
+}
+
+
+
+// The implementation is different depending on which DoFHandler is used. Test
+// both of them.
+template <int dim>
+void
+test_with_both_dof_handlers(const Triangulation<dim> &triangulation)
+{
+  const FE_Q<dim> element(1);
+
+  deallog << "dealii::DoFHandler" << std::endl;
+  {
+    DoFHandler<dim> dof_handler(triangulation);
+    dof_handler.distribute_dofs(element);
+    create_and_print_pattern<dim>(dof_handler);
+  }
+
+  deallog << std::endl;
+
+  deallog << "hp::DoFHandler" << std::endl;
+  {
+    hp::FECollection<dim> fe_collection(element);
+    hp::DoFHandler<dim>   dof_handler(triangulation);
+    dof_handler.distribute_dofs(fe_collection);
+    create_and_print_pattern<dim>(dof_handler);
+  }
+}
+
+
+
+//  Create the 3-cell-triangulation described at the top.
+void create_3_cell_triangulation(Triangulation<1> &triangulation)
+{
+  const unsigned int n_elements = 3;
+  const double       left       = -1;
+  const double       right      = 2;
+
+  GridGenerator::subdivided_hyper_cube(triangulation, n_elements, left, right);
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  const int          dim = 1;
+  Triangulation<dim> triangulation;
+  create_3_cell_triangulation(triangulation);
+
+  test_with_both_dof_handlers(triangulation);
+  return 0;
+}
diff --git a/tests/dofs/dof_tools_25.output b/tests/dofs/dof_tools_25.output
new file mode 100644 (file)
index 0000000..c2353b3
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::dealii::DoFHandler
+[0,0,1,2]
+[1,0,1,2]
+[2,0,1,2,3]
+[3,2,3]
+DEAL::
+DEAL::hp::DoFHandler
+[0,0,1,2]
+[1,0,1,2]
+[2,0,1,2,3]
+[3,2,3]

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.