]> https://gitweb.dealii.org/ - dealii.git/commitdiff
NonLocalDoFHandler.
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 23 Jul 2020 13:15:05 +0000 (15:15 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 30 Jul 2020 05:26:57 +0000 (07:26 +0200)
include/deal.II/dofs/dof_accessor.templates.h
include/deal.II/dofs/dof_handler.h
include/deal.II/dofs/non_local_dof_handler.h [new file with mode: 0644]
include/deal.II/fe/fe.h
source/dofs/dof_handler.cc
tests/non_local/fe_q1_nonlocal.h
tests/non_local/fe_q1_nonlocal_03.output
tests/non_local/fe_q1_nonlocal_04.cc [new file with mode: 0644]
tests/non_local/fe_q1_nonlocal_04.output [new file with mode: 0644]

index 9f2d6dd98cc6e1c9a9c46a5c7eceacf35cac9efe..08a9006f025059810488b9a655e70afd2b856fcb 100644 (file)
@@ -2252,88 +2252,6 @@ namespace internal
 
 
 
-      /**
-       * Implement setting dof indices on a cell. TO CHECK ZHAOWEI + LH
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static void
-      set_dof_indices(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
-        const std::vector<types::global_dof_index> &local_dof_indices)
-      {
-        Assert(accessor.has_children() == false, ExcInternalError());
-
-        const unsigned int dofs_per_vertex =
-                             accessor.get_fe().n_dofs_per_vertex(),
-                           dofs_per_line = accessor.get_fe().n_dofs_per_line(),
-                           dofs_per_quad = accessor.get_fe().n_dofs_per_quad(),
-                           dofs_per_hex  = accessor.get_fe().n_dofs_per_hex();
-
-        Assert(local_dof_indices.size() == accessor.get_fe().dofs_per_cell,
-               ExcInternalError());
-
-        unsigned int index = 0;
-
-        for (unsigned int vertex = 0;
-             vertex < GeometryInfo<dim>::vertices_per_cell;
-             ++vertex)
-          for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index)
-            accessor.set_vertex_dof_index(vertex,
-                                          d,
-                                          local_dof_indices[index],
-                                          accessor.active_fe_index());
-        // now copy dof numbers into the line. for lines in 3d with the
-        // wrong orientation, we have already made sure that we're ok
-        // by picking the correct vertices (this happens automatically
-        // in the vertex() function). however, if the line is in wrong
-        // orientation, we look at it in flipped orientation and we
-        // will have to adjust the shape function indices that we see
-        // to correspond to the correct (cell-local) ordering.
-        //
-        // of course, if dim<3, then there is nothing to adjust
-        for (unsigned int line = 0; line < GeometryInfo<dim>::lines_per_cell;
-             ++line)
-          for (unsigned int d = 0; d < dofs_per_line; ++d, ++index)
-            accessor.line(line)->set_dof_index(
-              dim < 3 ?
-                d :
-                accessor.get_fe().adjust_line_dof_index_for_line_orientation(
-                  d, accessor.line_orientation(line)),
-              local_dof_indices[index],
-              accessor.active_fe_index());
-        // now copy dof numbers into the face. for faces in 3d with the
-        // wrong orientation, we have already made sure that we're ok
-        // by picking the correct lines and vertices (this happens
-        // automatically in the line() and vertex()
-        // functions). however, if the face is in wrong orientation,
-        // we look at it in flipped orientation and we will have to
-        // adjust the shape function indices that we see to correspond
-        // to the correct (cell-local) ordering. The same applies, if
-        // the face_rotation or face_orientation is non-standard
-        //
-        // again, if dim<3, then there is nothing to adjust
-        for (unsigned int quad = 0; quad < GeometryInfo<dim>::quads_per_cell;
-             ++quad)
-          for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index)
-            accessor.quad(quad)->set_dof_index(
-              dim < 3 ?
-                d :
-                accessor.get_fe().adjust_quad_dof_index_for_face_orientation(
-                  d,
-                  accessor.face_orientation(quad),
-                  accessor.face_flip(quad),
-                  accessor.face_rotation(quad)),
-              local_dof_indices[index],
-              accessor.active_fe_index());
-        for (unsigned int d = 0; d < dofs_per_hex; ++d, ++index)
-          accessor.set_dof_index(d,
-                                 local_dof_indices[index],
-                                 accessor.active_fe_index());
-        Assert(index == accessor.get_fe().dofs_per_cell, ExcInternalError());
-      }
-
-
-
       /**
        * Implement setting non-local dof indices on a cell.
        */
@@ -2345,7 +2263,16 @@ namespace internal
       {
         Assert(accessor.has_children() == false, ExcInternalError());
 
-        const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell;
+        const unsigned int dofs_per_cell = accessor.get_fe().n_dofs_per_cell();
+        const unsigned int n_non_local_dofs =
+          accessor.get_fe().n_non_local_dofs_per_cell();
+        const unsigned int n_local_dofs = dofs_per_cell - n_non_local_dofs;
+
+        AssertDimension(local_non_local_dof_indices.size(), n_non_local_dofs);
+
+        // First the easy case.
+        if (n_non_local_dofs == 0)
+          return;
 
         types::global_dof_index *next_dof_index =
           const_cast<types::global_dof_index *>(
@@ -2353,22 +2280,11 @@ namespace internal
               get_cache_ptr(accessor.dof_handler,
                             accessor.present_level,
                             accessor.present_index,
-                            dofs_per_cell));
-
-        const unsigned int non_local_dofs =
-                             accessor.get_fe().n_non_local_dofs_per_cell(),
-                           n_dofs = accessor.get_fe().n_dofs_per_cell();
-
-        unsigned int index = 0;
+                            dofs_per_cell)) +
+          n_local_dofs;
 
-        for (unsigned int d = 0; d < n_dofs; ++d, ++next_dof_index)
-          if (d >= n_dofs - non_local_dofs)
-            {
-              *next_dof_index = local_non_local_dof_indices[index];
-              ++index;
-            }
-        Assert(index == accessor.get_fe().n_non_local_dofs_per_cell(),
-               ExcInternalError());
+        for (unsigned int d = 0; d < n_non_local_dofs; ++d, ++next_dof_index)
+          *next_dof_index = local_non_local_dof_indices[d];
       }
 
 
index f727ef1e27e91bbf06708779f957ca642ad42155..a0d8a2f79cdbe82bc7d44edfa0d4fb657a5e69ee 100644 (file)
@@ -592,6 +592,13 @@ public:
   void
   distribute_mg_dofs();
 
+  /**
+   * Distribute non local degrees of freedom. The local DoFs need to be
+   * distributed using distribute_dofs() before calling this function.
+   */
+  void
+  distribute_non_local_dofs();
+
   /**
    * This function returns whether this DoFHandler has DoFs distributed on
    * each multigrid level or in other words if distribute_mg_dofs() has been
diff --git a/include/deal.II/dofs/non_local_dof_handler.h b/include/deal.II/dofs/non_local_dof_handler.h
new file mode 100644 (file)
index 0000000..7cfe068
--- /dev/null
@@ -0,0 +1,149 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_non_local_dof_handler_h
+#define dealii_non_local_dof_handler_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/subscriptor.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declarations
+#ifndef DOXYGEN
+template <int dim, int spacedim, bool>
+class DoFCellAccessor;
+#endif
+
+/**
+ * This class is used to enumerate non local degrees of freedom. Its default
+ * implementation does nothing, since in general FiniteElement spaces only
+ * define degrees of freedom on vertices, edges, faces, or cells. There are
+ * cases, however, in which a FiniteElement space may define some non-local
+ * basis functions which are non-zero on a given cell, even if the basis
+ * function cannot be attributed to local subobjects of the given cell.
+ *
+ * In these cases, the DoFHandler class needs to know how to distribute these
+ * extra degrees of freedom, and since it cannot do so on its own, it asks
+ * the FiniteElement to return a NonLocalDoFHandler class that knows how to
+ * enumerate and distribute these non local degrees of freedom.
+ *
+ * This class is intended as a base class for all FiniteElement spaces that
+ * need to enumerate non local degrees of freedom.
+ *
+ * @ingroup dofs
+ */
+template <int dim, int spacedim = dim>
+class NonLocalDoFHandler : public Subscriptor
+{
+public:
+  /**
+   * Make the dimension available in function templates.
+   */
+  static const unsigned int dimension = dim;
+
+  /**
+   * Make the space dimension available in function templates.
+   */
+  static const unsigned int space_dimension = spacedim;
+
+  /**
+   * Standard constructor, not initializing any data.
+   */
+  NonLocalDoFHandler() = default;
+
+  /**
+   * Destructor.
+   */
+  virtual ~NonLocalDoFHandler() = default;
+
+  /**
+   * Copy operator. NonLocalDoFHandler objects may be large and expensive.
+   * They should not be copied, in particular not by accident, but
+   * rather deliberately constructed. As a consequence, this operator
+   * is explicitly removed from the interface of this class.
+   */
+  NonLocalDoFHandler &
+  operator=(const NonLocalDoFHandler &) = delete;
+
+  /**
+   * Return the non local dof indices associated to the current cell, for
+   * active cell iterators.
+   */
+  virtual std::vector<types::global_dof_index>
+  get_non_local_dof_indices(
+    const DoFCellAccessor<dim, spacedim, false> &) const;
+
+  /**
+   * Return the non local dof indices associated to the current cell, for
+   * level cell iterators.
+   */
+  virtual std::vector<types::global_dof_index>
+  get_non_local_dof_indices(const DoFCellAccessor<dim, spacedim, true> &) const;
+
+  /**
+   * Return the number of non local dof indices that are required in addition to
+   * the local ones.
+   */
+  virtual types::global_dof_index
+  n_additional_non_local_dofs() const;
+
+  /**
+   * At any given moment, only one NonLocalDoFHandler can be used. Throw an
+   * exception if two cells define an active FiniteElement for which
+   * get_non_local_dof_handler() return a different object.
+   *
+   * @ingroup Exceptions
+   */
+  DeclException0(ExcDifferentNonLocalDoFHandlers);
+};
+
+
+
+// ----------------------------------------------------------------------------
+template <int dim, int spacedim>
+inline std::vector<types::global_dof_index>
+NonLocalDoFHandler<dim, spacedim>::get_non_local_dof_indices(
+  const DoFCellAccessor<dim, spacedim, true> &) const
+{
+  // By default, we return an empty vector.
+  return std::vector<types::global_dof_index>();
+}
+
+
+
+template <int dim, int spacedim>
+inline std::vector<types::global_dof_index>
+NonLocalDoFHandler<dim, spacedim>::get_non_local_dof_indices(
+  const DoFCellAccessor<dim, spacedim, false> &) const
+{
+  // By default, we return an empty vector.
+  return std::vector<types::global_dof_index>();
+}
+
+
+
+template <int dim, int spacedim>
+inline types::global_dof_index
+NonLocalDoFHandler<dim, spacedim>::n_additional_non_local_dofs() const
+{
+  return 0;
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 586a3784ffa9efb875a41595b54ff48197109479..5ccb1589f9396320056d555f28c1136c1ef3c6df 100644 (file)
@@ -18,6 +18,8 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/dofs/non_local_dof_handler.h>
+
 #include <deal.II/fe/block_mask.h>
 #include <deal.II/fe/component_mask.h>
 #include <deal.II/fe/fe_base.h>
@@ -2260,6 +2262,17 @@ public:
 
   //@}
 
+  /**
+   * @name Non local dofs support
+   * @{
+   */
+  /**
+   * Return an object that knows how to handle non local dof indices.
+   */
+  virtual std::shared_ptr<const NonLocalDoFHandler<dim, spacedim>>
+  get_non_local_dof_handler() const;
+  //@}
+
   /**
    * Determine an estimate for the memory consumption (in bytes) of this
    * object.
@@ -3301,6 +3314,14 @@ FiniteElement<dim, spacedim>::get_associated_geometry_primitive(
 
 
 
+template <int dim, int spacedim>
+inline std::shared_ptr<const NonLocalDoFHandler<dim, spacedim>>
+FiniteElement<dim, spacedim>::get_non_local_dof_handler() const
+{
+  return std::make_shared<const NonLocalDoFHandler<dim, spacedim>>();
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 70a4c40cb17af2017360683218b926933206d2d3..1485977780ed7d3950e3341510e0356f9ab72e4e 100644 (file)
@@ -2565,6 +2565,8 @@ DoFHandler<dim, spacedim>::distribute_dofs(
             &*this->tria) == nullptr)
         this->block_info_object.initialize(*this, false, true);
     }
+  // now distribute non_local_dofs.
+  distribute_non_local_dofs();
 }
 
 
@@ -2601,6 +2603,31 @@ DoFHandler<dim, spacedim>::distribute_mg_dofs()
 
 
 
+template <int dim, int spacedim>
+void
+DoFHandler<dim, spacedim>::distribute_non_local_dofs()
+{
+  Assert(
+    this->object_dof_indices.size() > 0,
+    ExcMessage(
+      "Distribute active DoFs using distribute_dofs() before calling distribute_non_local_dofs()."));
+
+  auto non_local_dh = begin()->get_fe().get_non_local_dof_handler();
+  for (auto cell : active_cell_iterators())
+    {
+      // Make sure all cells use the same NonLocalDoFHandler.
+      Assert(non_local_dh == cell->get_fe().get_non_local_dof_handler(),
+             ExcMessage(
+               "You are trying to use different NonLocalDoFHandler objects"));
+      cell->set_non_local_dof_indices(
+        non_local_dh->get_non_local_dof_indices(*cell));
+    }
+  this->number_cache.n_global_dofs +=
+    non_local_dh->n_additional_non_local_dofs();
+}
+
+
+
 template <int dim, int spacedim>
 void
 DoFHandler<dim, spacedim>::initialize_local_block_info()
index 0e1a518af13413550209177ce281682805874f90..1d713624ac3841b51a98ce8534f981aa48035ead 100644 (file)
@@ -36,11 +36,45 @@ get_dpo_vector()
   return dpo;
 }
 
+
+/**
+ * Associates non local dofs according to the vertex index.
+ */
+template <int dim>
+class NonLocalQ1DoFHandler : public NonLocalDoFHandler<dim>
+{
+public:
+  virtual std::vector<types::global_dof_index>
+  get_non_local_dof_indices(
+    const DoFCellAccessor<dim, dim, false> &accessor) const override
+  {
+    if (!tria)
+      tria = &(accessor.get_triangulation());
+    std::vector<types::global_dof_index> dofs(accessor.n_vertices());
+    for (unsigned int i = 0; i < dofs.size(); ++i)
+      dofs[i] = accessor.vertex_index(i);
+    return dofs;
+  }
+
+  virtual types::global_dof_index
+  n_additional_non_local_dofs() const override
+  {
+    Assert(tria, ExcInternalError());
+    return tria->n_vertices();
+  }
+
+private:
+  mutable SmartPointer<const Triangulation<dim>> tria = nullptr;
+};
+
+
+
 template <int dim>
 class FE_Q1_Nonlocal : public FE_Q_Base<TensorProductPolynomials<dim>, dim, dim>
 {
 public:
-  FE_Q1_Nonlocal()
+  FE_Q1_Nonlocal(const std::shared_ptr<NonLocalQ1DoFHandler<dim>> &ptr =
+                   std::shared_ptr<NonLocalQ1DoFHandler<dim>>())
     : FE_Q_Base<TensorProductPolynomials<dim>, dim, dim>(
         TensorProductPolynomials<dim>(
           Polynomials::generate_complete_Lagrange_basis(
@@ -50,6 +84,7 @@ public:
                                1,
                                FiniteElementData<dim>::H1),
         std::vector<bool>(1, false))
+    , non_local_dh(ptr ? ptr : std::make_shared<NonLocalQ1DoFHandler<dim>>())
   {
     this->unit_support_points = QTrapez<dim>().get_points();
   }
@@ -57,7 +92,7 @@ public:
   virtual std::unique_ptr<FiniteElement<dim>>
   clone() const override
   {
-    return std::make_unique<FE_Q1_Nonlocal<dim>>();
+    return std::make_unique<FE_Q1_Nonlocal<dim>>(non_local_dh);
   }
 
   virtual std::string
@@ -65,4 +100,30 @@ public:
   {
     return "FE_Q_Nonlocal<dim>";
   }
+
+  virtual std::shared_ptr<const NonLocalDoFHandler<dim>>
+  get_non_local_dof_handler() const override
+  {
+    return non_local_dh;
+  }
+
+  virtual void
+  convert_generalized_support_point_values_to_dof_values(
+    const std::vector<Vector<double>> &support_point_values,
+    std::vector<double> &              nodal_values) const override
+  {
+    AssertDimension(support_point_values.size(),
+                    this->get_unit_support_points().size());
+    AssertDimension(support_point_values.size(), nodal_values.size());
+    AssertDimension(this->n_dofs_per_cell(), nodal_values.size());
+
+    for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i)
+      {
+        AssertDimension(support_point_values[i].size(), 1);
+        nodal_values[i] = support_point_values[i](0);
+      }
+  }
+
+private:
+  std::shared_ptr<NonLocalQ1DoFHandler<dim>> non_local_dh;
 };
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..fb71de2867cdd56077fa628bd1498d2d19382019 100644 (file)
@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/non_local/fe_q1_nonlocal_04.cc b/tests/non_local/fe_q1_nonlocal_04.cc
new file mode 100644 (file)
index 0000000..cebe016
--- /dev/null
@@ -0,0 +1,64 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 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.
+//
+// ---------------------------------------------------------------------
+
+// build an FE_Q1_Nonlocal finite element, distribute dofs on a simple
+// Triangulation, and interpolate a scalar function.
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../non_local/fe_q1_nonlocal.h"
+
+template <int dim>
+void
+test()
+{
+  FE_Q1_Nonlocal<dim> fe;
+  Triangulation<dim>  tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+  DoFHandler<dim> dh(tria);
+  dh.distribute_dofs(fe);
+
+  Tensor<1, dim> ones;
+  for (unsigned int d = 0; d < dim; ++d)
+    ones[d] = 1;
+
+  Vector<double> solution(dh.n_dofs());
+
+  VectorTools::interpolate(dh, Functions::Monomial<dim>(ones), solution);
+
+  deallog << solution << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  test<1>();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/non_local/fe_q1_nonlocal_04.output b/tests/non_local/fe_q1_nonlocal_04.output
new file mode 100644 (file)
index 0000000..e25d0a7
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::0.00000 1.00000 0.500000
+DEAL::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.500000 0.500000 0.250000
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.500000 0.500000 0.00000 0.00000 0.00000 0.250000 0.250000 0.250000 0.125000

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.