]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add hp to FEInterfaceValues 14547/head
authorMarco Feder <marco.feder@sissa.it>
Thu, 8 Dec 2022 11:03:27 +0000 (12:03 +0100)
committerMarco Feder <marco.feder@sissa.it>
Tue, 20 Dec 2022 10:25:07 +0000 (11:25 +0100)
include/deal.II/fe/fe_interface_values.h
tests/feinterface/fe_interface_values_11.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_11.output [new file with mode: 0644]
tests/feinterface/fe_interface_values_12.cc [new file with mode: 0644]
tests/feinterface/fe_interface_values_12.output [new file with mode: 0644]

index 52f430d44eac352b61bd8f3f63ecaf8d64d88035..70287b46610afd4f791d3abf1354ee0ccd395c2e 100644 (file)
@@ -23,6 +23,8 @@
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/mapping.h>
 
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/fe_values.h>
 #include <deal.II/hp/q_collection.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -1430,6 +1432,23 @@ public:
                     const Quadrature<dim - 1> &         quadrature,
                     const UpdateFlags                   update_flags);
 
+  /**
+   * Construct the FEInterfaceValues object with different FiniteElements
+   * assigned to different faces.
+   */
+  FEInterfaceValues(
+    const hp::MappingCollection<dim, spacedim> &mapping_collection,
+    const hp::FECollection<dim, spacedim> &     fe_collection,
+    const hp::QCollection<dim - 1> &            quadrature_collection,
+    const UpdateFlags                           update_flags);
+
+  /**
+   * Same as above, but using a Q1 mapping.
+   */
+  FEInterfaceValues(const hp::FECollection<dim, spacedim> &fe_collection,
+                    const hp::QCollection<dim - 1> &quadrature_collection,
+                    const UpdateFlags               update_flags);
+
   /**
    * Re-initialize this object to be used on a new interface given by two faces
    * of two neighboring cells. The `cell` and `cell_neighbor` cells will be
@@ -2117,22 +2136,23 @@ private:
   /**
    * The FEFaceValues object for the current cell.
    */
-  FEFaceValues<dim, spacedim> internal_fe_face_values;
+  std::unique_ptr<FEFaceValues<dim, spacedim>> internal_fe_face_values;
 
   /**
    * The FEFaceValues object for the current cell if the cell is refined.
    */
-  FESubfaceValues<dim, spacedim> internal_fe_subface_values;
+  std::unique_ptr<FESubfaceValues<dim, spacedim>> internal_fe_subface_values;
 
   /**
    * The FEFaceValues object for the neighboring cell.
    */
-  FEFaceValues<dim, spacedim> internal_fe_face_values_neighbor;
+  std::unique_ptr<FEFaceValues<dim, spacedim>> internal_fe_face_values_neighbor;
 
   /**
    * The FEFaceValues object for the neighboring cell if the cell is refined.
    */
-  FESubfaceValues<dim, spacedim> internal_fe_subface_values_neighbor;
+  std::unique_ptr<FESubfaceValues<dim, spacedim>>
+    internal_fe_subface_values_neighbor;
 
   /**
    * Pointer to internal_fe_face_values or internal_fe_subface_values,
@@ -2147,8 +2167,37 @@ private:
    */
   FEFaceValuesBase<dim, spacedim> *fe_face_values_neighbor;
 
-  /* Make the view classes friends of this class, since they access internal
-   * data.
+  /**
+   * An hp::FEValues object for the FEFaceValues on the
+   * present cell.
+   */
+  std::unique_ptr<hp::FEFaceValues<dim, spacedim>> internal_hp_fe_face_values;
+
+  /**
+   * An hp::FEValues object for the FESubfaceValues on the
+   * present cell.
+   */
+  std::unique_ptr<hp::FESubfaceValues<dim, spacedim>>
+    internal_hp_fe_subface_values;
+
+  /**
+   * An hp::FEValues object for the FEFaceValues on the
+   * neighbor of the present cell.
+   */
+  std::unique_ptr<hp::FEFaceValues<dim, spacedim>>
+    internal_hp_fe_face_values_neighbor;
+
+  /**
+   * An hp::FEValues object for the FESubfaceValues on the
+   * neighboring cell.
+   */
+  std::unique_ptr<hp::FESubfaceValues<dim, spacedim>>
+    internal_hp_fe_subface_values_neighbor;
+
+
+  /*
+   * Make the view classes friends of this class, since they
+   * access internal data.
    */
   template <int, int>
   friend class FEInterfaceViews::Scalar;
@@ -2169,10 +2218,26 @@ FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
   const Quadrature<dim - 1> &         quadrature,
   const UpdateFlags                   update_flags)
   : n_quadrature_points(quadrature.size())
-  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
-  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
-  , internal_fe_face_values_neighbor(mapping, fe, quadrature, update_flags)
-  , internal_fe_subface_values_neighbor(mapping, fe, quadrature, update_flags)
+  , internal_fe_face_values(
+      std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+                                                    fe,
+                                                    quadrature,
+                                                    update_flags))
+  , internal_fe_subface_values(
+      std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+                                                       fe,
+                                                       quadrature,
+                                                       update_flags))
+  , internal_fe_face_values_neighbor(
+      std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+                                                    fe,
+                                                    quadrature,
+                                                    update_flags))
+  , internal_fe_subface_values_neighbor(
+      std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+                                                       fe,
+                                                       quadrature,
+                                                       update_flags))
   , fe_face_values(nullptr)
   , fe_face_values_neighbor(nullptr)
 {}
@@ -2186,45 +2251,138 @@ FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
   const hp::QCollection<dim - 1> &    quadrature,
   const UpdateFlags                   update_flags)
   : n_quadrature_points(quadrature.max_n_quadrature_points())
-  , internal_fe_face_values(mapping, fe, quadrature, update_flags)
-  , internal_fe_subface_values(mapping, fe, quadrature, update_flags)
-  , internal_fe_face_values_neighbor(mapping, fe, quadrature[0], update_flags)
-  , internal_fe_subface_values_neighbor(mapping,
-                                        fe,
-                                        quadrature[0],
-                                        update_flags)
+
+  , internal_fe_face_values(
+      std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+                                                    fe,
+                                                    quadrature,
+                                                    update_flags))
+  , internal_fe_subface_values(
+      std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+                                                       fe,
+                                                       quadrature,
+                                                       update_flags))
+  , internal_fe_face_values_neighbor(
+      std::make_unique<FEFaceValues<dim, spacedim>>(mapping,
+                                                    fe,
+                                                    quadrature[0],
+                                                    update_flags))
+  , internal_fe_subface_values_neighbor(
+      std::make_unique<FESubfaceValues<dim, spacedim>>(mapping,
+                                                       fe,
+                                                       quadrature[0],
+                                                       update_flags))
   , fe_face_values(nullptr)
   , fe_face_values_neighbor(nullptr)
 {}
 
 
 
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+  const hp::MappingCollection<dim, spacedim> &mapping_collection,
+  const hp::FECollection<dim, spacedim> &     fe_collection,
+  const hp::QCollection<dim - 1> &            quadrature_collection,
+  const UpdateFlags                           update_flags)
+  : n_quadrature_points(quadrature_collection.max_n_quadrature_points())
+  , fe_face_values(nullptr)
+  , fe_face_values_neighbor(nullptr)
+  , internal_hp_fe_face_values(
+      std::make_unique<hp::FEFaceValues<dim, spacedim>>(mapping_collection,
+                                                        fe_collection,
+                                                        quadrature_collection,
+                                                        update_flags))
+  , internal_hp_fe_subface_values(
+      std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+        mapping_collection,
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+  , internal_hp_fe_face_values_neighbor(
+      std::make_unique<hp::FEFaceValues<dim, spacedim>>(mapping_collection,
+                                                        fe_collection,
+                                                        quadrature_collection,
+                                                        update_flags))
+  , internal_hp_fe_subface_values_neighbor(
+      std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+        mapping_collection,
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+{
+  AssertDimension(dim, spacedim);
+}
+
+
+
+template <int dim, int spacedim>
+FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
+  const hp::FECollection<dim, spacedim> &fe_collection,
+  const hp::QCollection<dim - 1> &       quadrature_collection,
+  const UpdateFlags                      update_flags)
+  : n_quadrature_points(quadrature_collection.max_n_quadrature_points())
+  , fe_face_values(nullptr)
+  , fe_face_values_neighbor(nullptr)
+  , internal_hp_fe_face_values(
+      std::make_unique<hp::FEFaceValues<dim, spacedim>>(
+        fe_collection.get_reference_cell_default_linear_mapping(),
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+  , internal_hp_fe_subface_values(
+      std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+        fe_collection.get_reference_cell_default_linear_mapping(),
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+  , internal_hp_fe_face_values_neighbor(
+      std::make_unique<hp::FEFaceValues<dim, spacedim>>(
+        fe_collection.get_reference_cell_default_linear_mapping(),
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+  , internal_hp_fe_subface_values_neighbor(
+      std::make_unique<hp::FESubfaceValues<dim, spacedim>>(
+        fe_collection.get_reference_cell_default_linear_mapping(),
+        fe_collection,
+        quadrature_collection,
+        update_flags))
+{
+  AssertDimension(dim, spacedim);
+}
+
+
+
 template <int dim, int spacedim>
 FEInterfaceValues<dim, spacedim>::FEInterfaceValues(
   const FiniteElement<dim, spacedim> &fe,
   const Quadrature<dim - 1> &         quadrature,
   const UpdateFlags                   update_flags)
   : n_quadrature_points(quadrature.size())
-  , internal_fe_face_values(
+  , internal_fe_face_values(std::make_unique<FEFaceValues<dim, spacedim>>(
       fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
       fe,
       quadrature,
-      update_flags)
-  , internal_fe_subface_values(
+      update_flags))
+  , internal_fe_subface_values(std::make_unique<FESubfaceValues<dim, spacedim>>(
       fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
       fe,
       quadrature,
-      update_flags)
+      update_flags))
   , internal_fe_face_values_neighbor(
-      fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
-      fe,
-      quadrature,
-      update_flags)
+      std::make_unique<FEFaceValues<dim, spacedim>>(
+        fe.reference_cell()
+          .template get_default_linear_mapping<dim, spacedim>(),
+        fe,
+        quadrature,
+        update_flags))
   , internal_fe_subface_values_neighbor(
-      fe.reference_cell().template get_default_linear_mapping<dim, spacedim>(),
-      fe,
-      quadrature,
-      update_flags)
+      std::make_unique<FESubfaceValues<dim, spacedim>>(
+        fe.reference_cell()
+          .template get_default_linear_mapping<dim, spacedim>(),
+        fe,
+        quadrature,
+        update_flags))
   , fe_face_values(nullptr)
   , fe_face_values_neighbor(nullptr)
 {}
@@ -2242,34 +2400,94 @@ FEInterfaceValues<dim, spacedim>::reinit(
   const unsigned int              face_no_neighbor,
   const unsigned int              sub_face_no_neighbor)
 {
-  if (sub_face_no == numbers::invalid_unsigned_int)
-    {
-      internal_fe_face_values.reinit(cell, face_no);
-      fe_face_values = &internal_fe_face_values;
-    }
-  else
-    {
-      internal_fe_subface_values.reinit(cell, face_no, sub_face_no);
-      fe_face_values = &internal_fe_subface_values;
-    }
-  if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+  Assert(internal_fe_face_values || internal_hp_fe_face_values,
+         ExcNotInitialized());
+
+  if (internal_fe_face_values)
     {
-      internal_fe_face_values_neighbor.reinit(cell_neighbor, face_no_neighbor);
-      fe_face_values_neighbor = &internal_fe_face_values_neighbor;
+      if (sub_face_no == numbers::invalid_unsigned_int)
+        {
+          internal_fe_face_values->reinit(cell, face_no);
+          fe_face_values = internal_fe_face_values.get();
+        }
+      else
+        {
+          internal_fe_subface_values->reinit(cell, face_no, sub_face_no);
+          fe_face_values = internal_fe_subface_values.get();
+        }
+      if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+        {
+          internal_fe_face_values_neighbor->reinit(cell_neighbor,
+                                                   face_no_neighbor);
+          fe_face_values_neighbor = internal_fe_face_values_neighbor.get();
+        }
+      else
+        {
+          internal_fe_subface_values_neighbor->reinit(cell_neighbor,
+                                                      face_no_neighbor,
+                                                      sub_face_no_neighbor);
+          fe_face_values_neighbor = internal_fe_subface_values_neighbor.get();
+        }
+
+      AssertDimension(fe_face_values->n_quadrature_points,
+                      fe_face_values_neighbor->n_quadrature_points);
+
+      const_cast<unsigned int &>(this->n_quadrature_points) =
+        fe_face_values->n_quadrature_points;
     }
-  else
+  else if (internal_hp_fe_face_values)
     {
-      internal_fe_subface_values_neighbor.reinit(cell_neighbor,
-                                                 face_no_neighbor,
-                                                 sub_face_no_neighbor);
-      fe_face_values_neighbor = &internal_fe_subface_values_neighbor;
-    }
+      const unsigned int dominated_fe_index =
+        internal_hp_fe_face_values->get_fe_collection().find_dominated_fe(
+          {cell->active_fe_index(), cell_neighbor->active_fe_index()});
+
+      // Same as if above, but when hp is enabled.
+      if (sub_face_no == numbers::invalid_unsigned_int)
+        {
+          internal_hp_fe_face_values->reinit(cell,
+                                             face_no,
+                                             dominated_fe_index,
+                                             dominated_fe_index);
+          fe_face_values = &const_cast<FEFaceValues<dim, spacedim> &>(
+            internal_hp_fe_face_values->get_present_fe_values());
+        }
+      else
+        {
+          internal_hp_fe_subface_values->reinit(
+            cell, face_no, sub_face_no, dominated_fe_index, dominated_fe_index);
+
+          fe_face_values = &const_cast<FESubfaceValues<dim, spacedim> &>(
+            internal_hp_fe_subface_values->get_present_fe_values());
+        }
+      if (sub_face_no_neighbor == numbers::invalid_unsigned_int)
+        {
+          internal_hp_fe_face_values_neighbor->reinit(cell_neighbor,
+                                                      face_no_neighbor,
+                                                      dominated_fe_index,
+                                                      dominated_fe_index);
+
+          fe_face_values_neighbor = &const_cast<FEFaceValues<dim, spacedim> &>(
+            internal_hp_fe_face_values_neighbor->get_present_fe_values());
+        }
+      else
+        {
+          internal_hp_fe_subface_values_neighbor->reinit(cell_neighbor,
+                                                         face_no_neighbor,
+                                                         sub_face_no_neighbor,
+                                                         dominated_fe_index,
+                                                         dominated_fe_index);
+
+          fe_face_values_neighbor =
+            &const_cast<FESubfaceValues<dim, spacedim> &>(
+              internal_hp_fe_subface_values_neighbor->get_present_fe_values());
+        }
 
-  AssertDimension(fe_face_values->n_quadrature_points,
-                  fe_face_values_neighbor->n_quadrature_points);
+      AssertDimension(fe_face_values->n_quadrature_points,
+                      fe_face_values_neighbor->n_quadrature_points);
 
-  const_cast<unsigned int &>(this->n_quadrature_points) =
-    fe_face_values->n_quadrature_points;
+      const_cast<unsigned int &>(this->n_quadrature_points) =
+        fe_face_values->n_quadrature_points;
+    }
 
   // Set up dof mapping and remove duplicates (for continuous elements).
   {
@@ -2323,15 +2541,30 @@ void
 FEInterfaceValues<dim, spacedim>::reinit(const CellIteratorType &cell,
                                          const unsigned int      face_no)
 {
-  internal_fe_face_values.reinit(cell, face_no);
-  fe_face_values          = &internal_fe_face_values;
-  fe_face_values_neighbor = nullptr;
+  Assert(internal_fe_face_values || internal_hp_fe_face_values,
+         ExcNotInitialized());
 
-  interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
-  cell->get_active_or_mg_dof_indices(interface_dof_indices);
+  if (internal_fe_face_values)
+    {
+      internal_fe_face_values->reinit(cell, face_no);
+      fe_face_values          = internal_fe_face_values.get();
+      fe_face_values_neighbor = nullptr;
 
-  dofmap.resize(interface_dof_indices.size());
+      interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+      cell->get_active_or_mg_dof_indices(interface_dof_indices);
+    }
+  else if (internal_hp_fe_face_values)
+    {
+      internal_hp_fe_face_values->reinit(cell, face_no);
+      fe_face_values = &const_cast<FEFaceValues<dim> &>(
+        internal_hp_fe_face_values->get_present_fe_values());
+      fe_face_values_neighbor = nullptr;
 
+      interface_dof_indices.resize(fe_face_values->get_fe().n_dofs_per_cell());
+      cell->get_active_or_mg_dof_indices(interface_dof_indices);
+    }
+
+  dofmap.resize(interface_dof_indices.size());
   for (unsigned int i = 0; i < interface_dof_indices.size(); ++i)
     {
       dofmap[i] = {{i, numbers::invalid_unsigned_int}};
@@ -2377,7 +2610,7 @@ template <int dim, int spacedim>
 const Mapping<dim, spacedim> &
 FEInterfaceValues<dim, spacedim>::get_mapping() const
 {
-  return internal_fe_face_values.get_mapping();
+  return internal_fe_face_values->get_mapping();
 }
 
 
@@ -2386,7 +2619,7 @@ template <int dim, int spacedim>
 const FiniteElement<dim, spacedim> &
 FEInterfaceValues<dim, spacedim>::get_fe() const
 {
-  return internal_fe_face_values.get_fe();
+  return internal_fe_face_values->get_fe();
 }
 
 
@@ -2395,7 +2628,7 @@ template <int dim, int spacedim>
 const Quadrature<dim - 1> &
 FEInterfaceValues<dim, spacedim>::get_quadrature() const
 {
-  return internal_fe_face_values.get_quadrature();
+  return internal_fe_face_values->get_quadrature();
 }
 
 
@@ -2424,7 +2657,7 @@ template <int dim, int spacedim>
 UpdateFlags
 FEInterfaceValues<dim, spacedim>::get_update_flags() const
 {
-  return internal_fe_face_values.get_update_flags();
+  return internal_fe_face_values->get_update_flags();
 }
 
 
diff --git a/tests/feinterface/fe_interface_values_11.cc b/tests/feinterface/fe_interface_values_11.cc
new file mode 100644 (file)
index 0000000..52b09de
--- /dev/null
@@ -0,0 +1,196 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2022 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 basic properties of FEInterfaceValues in the hp case.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+#include <deal.II/hp/q_collection.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+inspect_fiv(FEInterfaceValues<dim> &fiv)
+{
+  deallog << "at_boundary(): " << fiv.at_boundary() << "\n"
+          << "n_current_interface_dofs(): " << fiv.n_current_interface_dofs()
+          << "\n";
+
+  std::vector<types::global_dof_index> indices =
+    fiv.get_interface_dof_indices();
+  Assert(indices.size() == fiv.n_current_interface_dofs(), ExcInternalError());
+
+  deallog << "interface_dof_indices: ";
+  for (auto i : indices)
+    deallog << i << ' ';
+  deallog << "\n";
+
+
+  unsigned int idx = 0;
+  for (auto v : indices)
+    {
+      deallog << "  index " << idx << " global_dof_index:" << v << ":\n";
+
+      const auto pair = fiv.interface_dof_to_dof_indices(idx);
+      deallog << "    dof indices: " << static_cast<int>(pair[0]) << " | "
+              << static_cast<int>(pair[1]) << "\n";
+
+      ++idx;
+    }
+
+  deallog << std::endl;
+}
+
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void
+make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void
+make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(const unsigned int p)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  DoFHandler<dim>          dofh(tria);
+  hp::FECollection<dim>    fe_collection;
+  hp::QCollection<dim - 1> q_collection;
+
+  fe_collection.push_back(FE_DGQ<dim>(p));
+  fe_collection.push_back(FE_DGQ<dim>(p + 1));
+
+  q_collection.push_back(QGauss<dim - 1>(p));
+  q_collection.push_back(QGauss<dim - 1>(p + 1));
+
+  // Set different finite elements spaces on the two cells.
+  unsigned int fe_index = 0;
+  for (const auto &cell : dofh.active_cell_iterators())
+    {
+      cell->set_active_fe_index(fe_index);
+      ++fe_index;
+    }
+
+  dofh.distribute_dofs(fe_collection);
+
+  UpdateFlags update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+
+  FEInterfaceValues<dim> fiv(fe_collection, q_collection, update_flags);
+
+
+  auto cell = dofh.begin();
+
+  deallog << "** interface between cell 0 and 1 **\n";
+
+  for (const unsigned int f : GeometryInfo<dim>::face_indices())
+    if (!cell->at_boundary(f))
+      {
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   cell->neighbor_of_neighbor(f),
+                   numbers::invalid_unsigned_int);
+
+        Assert(fiv.get_fe_face_values(0).get_cell() == cell,
+               ExcInternalError());
+        Assert(fiv.get_fe_face_values(1).get_cell() == cell->neighbor(f),
+               ExcInternalError());
+        Assert(fiv.n_current_interface_dofs() ==
+                 fe_collection[0].n_dofs_per_cell() +
+                   fe_collection[1].n_dofs_per_cell(),
+               ExcInternalError());
+        Assert(!fiv.at_boundary(), ExcInternalError());
+
+        auto mycell = cell;
+        for (unsigned int c = 0; c < 2; ++c)
+          {
+            std::vector<types::global_dof_index> indices(
+              fe_collection[c].n_dofs_per_cell());
+            mycell->get_dof_indices(indices);
+            deallog << "cell " << c << ": ";
+            for (auto i : indices)
+              deallog << i << ' ';
+            deallog << "\n";
+            ++mycell;
+          }
+
+        inspect_fiv(fiv);
+      }
+
+  deallog << "** boundary interface on cell 0 **\n" << std::endl;
+
+  {
+    fiv.reinit(cell, 0);
+    Assert(fiv.get_fe_face_values(0).get_cell() == cell, ExcInternalError());
+    Assert(fiv.n_current_interface_dofs() == fe_collection[0].n_dofs_per_cell(),
+           ExcInternalError());
+    Assert(fiv.at_boundary(), ExcInternalError());
+    inspect_fiv(fiv);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  for (unsigned int p : {1, 2})
+    test<2>(p);
+  for (unsigned int p : {1})
+    test<3>(p);
+}
diff --git a/tests/feinterface/fe_interface_values_11.output b/tests/feinterface/fe_interface_values_11.output
new file mode 100644 (file)
index 0000000..0817989
--- /dev/null
@@ -0,0 +1,228 @@
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 
+cell 1: 4 5 6 7 8 9 10 11 12 
+at_boundary(): 0
+n_current_interface_dofs(): 13
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+  index 4 global_dof_index:4:
+    dof indices: -1 | 0
+  index 5 global_dof_index:5:
+    dof indices: -1 | 1
+  index 6 global_dof_index:6:
+    dof indices: -1 | 2
+  index 7 global_dof_index:7:
+    dof indices: -1 | 3
+  index 8 global_dof_index:8:
+    dof indices: -1 | 4
+  index 9 global_dof_index:9:
+    dof indices: -1 | 5
+  index 10 global_dof_index:10:
+    dof indices: -1 | 6
+  index 11 global_dof_index:11:
+    dof indices: -1 | 7
+  index 12 global_dof_index:12:
+    dof indices: -1 | 8
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 4
+interface_dof_indices: 0 1 2 3 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 4 5 6 7 8 
+cell 1: 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
+at_boundary(): 0
+n_current_interface_dofs(): 25
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+  index 4 global_dof_index:4:
+    dof indices: 4 | -1
+  index 5 global_dof_index:5:
+    dof indices: 5 | -1
+  index 6 global_dof_index:6:
+    dof indices: 6 | -1
+  index 7 global_dof_index:7:
+    dof indices: 7 | -1
+  index 8 global_dof_index:8:
+    dof indices: 8 | -1
+  index 9 global_dof_index:9:
+    dof indices: -1 | 0
+  index 10 global_dof_index:10:
+    dof indices: -1 | 1
+  index 11 global_dof_index:11:
+    dof indices: -1 | 2
+  index 12 global_dof_index:12:
+    dof indices: -1 | 3
+  index 13 global_dof_index:13:
+    dof indices: -1 | 4
+  index 14 global_dof_index:14:
+    dof indices: -1 | 5
+  index 15 global_dof_index:15:
+    dof indices: -1 | 6
+  index 16 global_dof_index:16:
+    dof indices: -1 | 7
+  index 17 global_dof_index:17:
+    dof indices: -1 | 8
+  index 18 global_dof_index:18:
+    dof indices: -1 | 9
+  index 19 global_dof_index:19:
+    dof indices: -1 | 10
+  index 20 global_dof_index:20:
+    dof indices: -1 | 11
+  index 21 global_dof_index:21:
+    dof indices: -1 | 12
+  index 22 global_dof_index:22:
+    dof indices: -1 | 13
+  index 23 global_dof_index:23:
+    dof indices: -1 | 14
+  index 24 global_dof_index:24:
+    dof indices: -1 | 15
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 9
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+  index 4 global_dof_index:4:
+    dof indices: 4 | -1
+  index 5 global_dof_index:5:
+    dof indices: 5 | -1
+  index 6 global_dof_index:6:
+    dof indices: 6 | -1
+  index 7 global_dof_index:7:
+    dof indices: 7 | -1
+  index 8 global_dof_index:8:
+    dof indices: 8 | -1
+
+DEAL::** interface between cell 0 and 1 **
+cell 0: 0 1 2 3 4 5 6 7 
+cell 1: 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 
+at_boundary(): 0
+n_current_interface_dofs(): 35
+interface_dof_indices: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+  index 4 global_dof_index:4:
+    dof indices: 4 | -1
+  index 5 global_dof_index:5:
+    dof indices: 5 | -1
+  index 6 global_dof_index:6:
+    dof indices: 6 | -1
+  index 7 global_dof_index:7:
+    dof indices: 7 | -1
+  index 8 global_dof_index:8:
+    dof indices: -1 | 0
+  index 9 global_dof_index:9:
+    dof indices: -1 | 1
+  index 10 global_dof_index:10:
+    dof indices: -1 | 2
+  index 11 global_dof_index:11:
+    dof indices: -1 | 3
+  index 12 global_dof_index:12:
+    dof indices: -1 | 4
+  index 13 global_dof_index:13:
+    dof indices: -1 | 5
+  index 14 global_dof_index:14:
+    dof indices: -1 | 6
+  index 15 global_dof_index:15:
+    dof indices: -1 | 7
+  index 16 global_dof_index:16:
+    dof indices: -1 | 8
+  index 17 global_dof_index:17:
+    dof indices: -1 | 9
+  index 18 global_dof_index:18:
+    dof indices: -1 | 10
+  index 19 global_dof_index:19:
+    dof indices: -1 | 11
+  index 20 global_dof_index:20:
+    dof indices: -1 | 12
+  index 21 global_dof_index:21:
+    dof indices: -1 | 13
+  index 22 global_dof_index:22:
+    dof indices: -1 | 14
+  index 23 global_dof_index:23:
+    dof indices: -1 | 15
+  index 24 global_dof_index:24:
+    dof indices: -1 | 16
+  index 25 global_dof_index:25:
+    dof indices: -1 | 17
+  index 26 global_dof_index:26:
+    dof indices: -1 | 18
+  index 27 global_dof_index:27:
+    dof indices: -1 | 19
+  index 28 global_dof_index:28:
+    dof indices: -1 | 20
+  index 29 global_dof_index:29:
+    dof indices: -1 | 21
+  index 30 global_dof_index:30:
+    dof indices: -1 | 22
+  index 31 global_dof_index:31:
+    dof indices: -1 | 23
+  index 32 global_dof_index:32:
+    dof indices: -1 | 24
+  index 33 global_dof_index:33:
+    dof indices: -1 | 25
+  index 34 global_dof_index:34:
+    dof indices: -1 | 26
+
+DEAL::** boundary interface on cell 0 **
+
+DEAL::at_boundary(): 1
+n_current_interface_dofs(): 8
+interface_dof_indices: 0 1 2 3 4 5 6 7 
+  index 0 global_dof_index:0:
+    dof indices: 0 | -1
+  index 1 global_dof_index:1:
+    dof indices: 1 | -1
+  index 2 global_dof_index:2:
+    dof indices: 2 | -1
+  index 3 global_dof_index:3:
+    dof indices: 3 | -1
+  index 4 global_dof_index:4:
+    dof indices: 4 | -1
+  index 5 global_dof_index:5:
+    dof indices: 5 | -1
+  index 6 global_dof_index:6:
+    dof indices: 6 | -1
+  index 7 global_dof_index:7:
+    dof indices: 7 | -1
+
diff --git a/tests/feinterface/fe_interface_values_12.cc b/tests/feinterface/fe_interface_values_12.cc
new file mode 100644 (file)
index 0000000..e32d1f0
--- /dev/null
@@ -0,0 +1,173 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// evaluate jump_in_shape_values(), average_of_shape_values(), shape_value() of
+// FEInterfaceValues on an adaptive mesh in the hp scenario.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim>
+void
+make_2_cells(Triangulation<dim> &tria);
+
+template <>
+void
+make_2_cells<2>(Triangulation<2> &tria)
+{
+  const unsigned int        dim         = 2;
+  std::vector<unsigned int> repetitions = {2, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+template <>
+void
+make_2_cells<3>(Triangulation<3> &tria)
+{
+  const unsigned int        dim         = 3;
+  std::vector<unsigned int> repetitions = {2, 1, 1};
+  Point<dim>                p1;
+  Point<dim>                p2(2.0, 1.0, 1.0);
+
+  GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2);
+}
+
+
+template <int dim>
+void
+test(const unsigned int fe_degree0, const unsigned int fe_degree1 = 0)
+{
+  Triangulation<dim> tria;
+  make_2_cells(tria);
+
+  tria.begin()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+
+
+  UpdateFlags update_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+  DoFHandler<dim>          dofh(tria);
+  hp::FECollection<dim>    fe_collection;
+  hp::QCollection<dim - 1> q_collection;
+  fe_collection.push_back(FE_DGQ<dim>(fe_degree0));
+  fe_collection.push_back(FE_DGQ<dim>(fe_degree1));
+
+  q_collection.push_back(QGauss<dim - 1>(fe_degree0 + 1));
+  // Set different finite elements spaces on the two cells.
+  unsigned int fe_index = 0;
+  for (const auto &cell : dofh.active_cell_iterators_on_level(0))
+    {
+      cell->set_active_fe_index(fe_index);
+      ++fe_index;
+    }
+
+  deallog << fe_collection[0].get_name() << "-" << fe_collection[1].get_name()
+          << std::endl;
+  dofh.distribute_dofs(fe_collection);
+
+  FEInterfaceValues<dim> fiv(fe_collection, q_collection, update_flags);
+
+  auto cell = dofh.begin(1);
+  ++cell;
+
+  for (const unsigned int f : GeometryInfo<dim>::face_indices())
+    if (!cell->at_boundary(f))
+      {
+        if (!cell->neighbor_is_coarser(f))
+          continue;
+
+        auto nn = cell->neighbor_of_coarser_neighbor(f);
+        fiv.reinit(cell,
+                   f,
+                   numbers::invalid_unsigned_int,
+                   cell->neighbor(f),
+                   nn.first,
+                   nn.second);
+
+        const unsigned int n_dofs = fiv.n_current_interface_dofs();
+        Vector<double>     cell_vector(n_dofs);
+
+        const auto &q_points = fiv.get_quadrature_points();
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          deallog << "qpoint " << qpoint << ": " << q_points[qpoint]
+                  << std::endl;
+
+        for (unsigned int idx = 0; idx < n_dofs; ++idx)
+          {
+            const auto pair = fiv.interface_dof_to_dof_indices(idx);
+            deallog << "  idx: " << idx
+                    << " global: " << fiv.get_interface_dof_indices()[idx]
+                    << " dof indices: " << static_cast<int>(pair[0]) << " | "
+                    << static_cast<int>(pair[1]) << std::endl;
+          }
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(true, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(true): " << cell_vector << std::endl;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) +=
+              fiv.shape_value(false, i, qpoint) * fiv.get_JxW_values()[qpoint];
+        deallog << "shape_value(false): " << cell_vector << std::endl;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) += fiv.jump_in_shape_values(i, qpoint) *
+                              fiv.get_JxW_values()[qpoint];
+        deallog << "jump_in_shape_values(): " << cell_vector << std::endl;
+
+        cell_vector = 0.0;
+        for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            cell_vector(i) += fiv.average_of_shape_values(i, qpoint) *
+                              fiv.get_JxW_values()[qpoint];
+        deallog << "average_of_shape_values(): " << cell_vector << std::endl;
+      }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>(0);
+  test<2>(1, 2);
+  test<3>(0);
+  test<3>(1, 2);
+}
diff --git a/tests/feinterface/fe_interface_values_12.output b/tests/feinterface/fe_interface_values_12.output
new file mode 100644 (file)
index 0000000..93d406c
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL::FE_DGQ<2>(0)-FE_DGQ<2>(0)
+DEAL::qpoint 0: 1.00000 0.250000
+DEAL::  idx: 0 global: 0 dof indices: -1 | 0
+DEAL::  idx: 1 global: 2 dof indices: 0 | -1
+DEAL::shape_value(true): 0.00000 0.500000
+DEAL::shape_value(false): 0.500000 0.00000
+DEAL::jump_in_shape_values(): -0.500000 0.500000
+DEAL::average_of_shape_values(): 0.250000 0.250000
+DEAL::FE_DGQ<2>(1)-FE_DGQ<2>(2)
+DEAL::qpoint 0: 1.00000 0.105662
+DEAL::qpoint 1: 1.00000 0.394338
+DEAL::  idx: 0 global: 0 dof indices: -1 | 0
+DEAL::  idx: 1 global: 1 dof indices: -1 | 1
+DEAL::  idx: 2 global: 2 dof indices: -1 | 2
+DEAL::  idx: 3 global: 3 dof indices: -1 | 3
+DEAL::  idx: 4 global: 8 dof indices: 0 | -1
+DEAL::  idx: 5 global: 9 dof indices: 1 | -1
+DEAL::  idx: 6 global: 10 dof indices: 2 | -1
+DEAL::  idx: 7 global: 11 dof indices: 3 | -1
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 0.250000
+DEAL::shape_value(false): 0.375000 0.00000 0.125000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::jump_in_shape_values(): -0.375000 0.00000 -0.125000 0.00000 0.00000 0.250000 0.00000 0.250000
+DEAL::average_of_shape_values(): 0.187500 0.00000 0.0625000 0.00000 0.00000 0.125000 0.00000 0.125000
+DEAL::FE_DGQ<3>(0)-FE_DGQ<3>(0)
+DEAL::qpoint 0: 1.00000 0.250000 0.250000
+DEAL::  idx: 0 global: 0 dof indices: -1 | 0
+DEAL::  idx: 1 global: 2 dof indices: 0 | -1
+DEAL::shape_value(true): 0.00000 0.250000
+DEAL::shape_value(false): 0.250000 0.00000
+DEAL::jump_in_shape_values(): -0.250000 0.250000
+DEAL::average_of_shape_values(): 0.125000 0.125000
+DEAL::FE_DGQ<3>(1)-FE_DGQ<3>(2)
+DEAL::qpoint 0: 1.00000 0.105662 0.105662
+DEAL::qpoint 1: 1.00000 0.394338 0.105662
+DEAL::qpoint 2: 1.00000 0.105662 0.394338
+DEAL::qpoint 3: 1.00000 0.394338 0.394338
+DEAL::  idx: 0 global: 0 dof indices: -1 | 0
+DEAL::  idx: 1 global: 1 dof indices: -1 | 1
+DEAL::  idx: 2 global: 2 dof indices: -1 | 2
+DEAL::  idx: 3 global: 3 dof indices: -1 | 3
+DEAL::  idx: 4 global: 4 dof indices: -1 | 4
+DEAL::  idx: 5 global: 5 dof indices: -1 | 5
+DEAL::  idx: 6 global: 6 dof indices: -1 | 6
+DEAL::  idx: 7 global: 7 dof indices: -1 | 7
+DEAL::  idx: 8 global: 16 dof indices: 0 | -1
+DEAL::  idx: 9 global: 17 dof indices: 1 | -1
+DEAL::  idx: 10 global: 18 dof indices: 2 | -1
+DEAL::  idx: 11 global: 19 dof indices: 3 | -1
+DEAL::  idx: 12 global: 20 dof indices: 4 | -1
+DEAL::  idx: 13 global: 21 dof indices: 5 | -1
+DEAL::  idx: 14 global: 22 dof indices: 6 | -1
+DEAL::  idx: 15 global: 23 dof indices: 7 | -1
+DEAL::shape_value(true): 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000
+DEAL::shape_value(false): 0.140625 0.00000 0.0468750 0.00000 0.0468750 0.00000 0.0156250 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::jump_in_shape_values(): -0.140625 0.00000 -0.0468750 0.00000 -0.0468750 0.00000 -0.0156250 0.00000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000 0.00000 0.0625000
+DEAL::average_of_shape_values(): 0.0703125 0.00000 0.0234375 0.00000 0.0234375 0.00000 0.00781250 0.00000 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500 0.00000 0.0312500

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.