From f791b78827fa0444093b31dcb2f48ba6e1a32d82 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Sat, 31 Aug 2019 17:21:23 -0400
Subject: [PATCH] Permit reinitialization with iterators instead of (sub)face
 numbers.

This is convenient when paired with the new cell->face_iterators() function.
---
 doc/news/changes/minor/20191027DavidWells |  4 ++
 include/deal.II/fe/fe_values.h            | 66 ++++++++++++++++++
 source/fe/fe_values.cc                    | 54 +++++++++++++++
 source/fe/fe_values.inst.in               | 26 ++++++-
 tests/fe/shapes.h                         | 82 ++++++++++++++++++-----
 5 files changed, 212 insertions(+), 20 deletions(-)
 create mode 100644 doc/news/changes/minor/20191027DavidWells

diff --git a/doc/news/changes/minor/20191027DavidWells b/doc/news/changes/minor/20191027DavidWells
new file mode 100644
index 0000000000..fb70fa6b38
--- /dev/null
+++ b/doc/news/changes/minor/20191027DavidWells
@@ -0,0 +1,4 @@
+New: Added two new overloads of FEValues::reinit() that take face and subface
+iterators instead of face and subface numbers.
+<br>
+(David Wells, 2019/10/27)
diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h
index 8229d7b760..da9b132a57 100644
--- a/include/deal.II/fe/fe_values.h
+++ b/include/deal.II/fe/fe_values.h
@@ -3757,6 +3757,18 @@ public:
                                             level_dof_access>> &cell,
          const unsigned int                                     face_no);
 
+  /**
+   * Reinitialize the gradients, Jacobi determinants, etc for face @p face
+   * and cell @p cell.
+   *
+   * @note @p face must be one of @p cell's face iterators.
+   */
+  template <template <int, int> class DoFHandlerType, bool level_dof_access>
+  void
+  reinit(const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>,
+                                            level_dof_access>> &     cell,
+         const typename Triangulation<dim, spacedim>::face_iterator &face);
+
   /**
    * Reinitialize the gradients, Jacobi determinants, etc for the given face
    * on a given cell of type "iterator into a Triangulation object", and the
@@ -3774,6 +3786,25 @@ public:
   reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
          const unsigned int                                          face_no);
 
+  /*
+   * Reinitialize the gradients, Jacobi determinants, etc for the given face
+   * on a given cell of type "iterator into a Triangulation object", and the
+   * given finite element. Since iterators into a triangulation alone only
+   * convey information about the geometry of a cell, but not about degrees of
+   * freedom possibly associated with this cell, you will not be able to call
+   * some functions of this class if they need information about degrees of
+   * freedom. These functions are, above all, the
+   * <tt>get_function_value/gradients/hessians/third_derivatives</tt>
+   * functions. If you want to call these functions, you have to call the @p
+   * reinit variants that take iterators into DoFHandler or other DoF handler
+   * type objects.
+   *
+   * @note @p face must be one of @p cell's face iterators.
+   */
+  void
+  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+         const typename Triangulation<dim, spacedim>::face_iterator &face);
+
   /**
    * Return a reference to this very object.
    *
@@ -3877,6 +3908,17 @@ public:
          const unsigned int                                     face_no,
          const unsigned int                                     subface_no);
 
+  /**
+   * Alternative reinitialization function that takes, as arguments, iterators
+   * to the face and subface instead of their numbers.
+   */
+  template <template <int, int> class DoFHandlerType, bool level_dof_access>
+  void
+  reinit(const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>,
+                                            level_dof_access>> &     cell,
+         const typename Triangulation<dim, spacedim>::face_iterator &face,
+         const typename Triangulation<dim, spacedim>::face_iterator &subface);
+
   /**
    * Reinitialize the gradients, Jacobi determinants, etc for the given
    * subface on a given cell of type "iterator into a Triangulation object", and
@@ -3895,6 +3937,30 @@ public:
          const unsigned int                                          face_no,
          const unsigned int subface_no);
 
+  /**
+   * Reinitialize the gradients, Jacobi determinants, etc for the given
+   * subface on a given cell of type "iterator into a Triangulation object", and
+   * the given finite element. Since iterators into a triangulation alone only
+   * convey information about the geometry of a cell, but not about degrees of
+   * freedom possibly associated with this cell, you will not be able to call
+   * some functions of this class if they need information about degrees of
+   * freedom. These functions are, above all, the
+   * <tt>get_function_value/gradients/hessians/third_derivatives</tt>
+   * functions. If you want to call these functions, you have to call the @p
+   * reinit variants that take iterators into DoFHandler or other DoF handler
+   * type objects.
+   *
+   * This does the same thing as the previous function but takes iterators
+   * instead of numbers as arguments.
+   *
+   * @note @p face and @p subface must correspond to a face (and a subface of
+   * that face) of @p cell.
+   */
+  void
+  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+         const typename Triangulation<dim, spacedim>::face_iterator &face,
+         const typename Triangulation<dim, spacedim>::face_iterator &subface);
+
   /**
    * Return a reference to this very object.
    *
diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc
index 4208da6368..914733969f 100644
--- a/source/fe/fe_values.cc
+++ b/source/fe/fe_values.cc
@@ -4784,6 +4784,19 @@ FEFaceValues<dim, spacedim>::reinit(
 
 
 
+template <int dim, int spacedim>
+template <template <int, int> class DoFHandlerType, bool lda>
+void
+FEFaceValues<dim, spacedim>::reinit(
+  const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, lda>> &cell,
+  const typename Triangulation<dim, spacedim>::face_iterator &             face)
+{
+  const auto face_n = cell->face_iterator_to_index(face);
+  reinit(cell, face_n);
+}
+
+
+
 template <int dim, int spacedim>
 void
 FEFaceValues<dim, spacedim>::reinit(
@@ -4806,6 +4819,18 @@ FEFaceValues<dim, spacedim>::reinit(
 
 
 
+template <int dim, int spacedim>
+void
+FEFaceValues<dim, spacedim>::reinit(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const typename Triangulation<dim, spacedim>::face_iterator &face)
+{
+  const auto face_n = cell->face_iterator_to_index(face);
+  reinit(cell, face_n);
+}
+
+
+
 template <int dim, int spacedim>
 void
 FEFaceValues<dim, spacedim>::do_reinit(const unsigned int face_no)
@@ -4980,6 +5005,21 @@ FESubfaceValues<dim, spacedim>::reinit(
 
 
 
+template <int dim, int spacedim>
+template <template <int, int> class DoFHandlerType, bool lda>
+void
+FESubfaceValues<dim, spacedim>::reinit(
+  const TriaIterator<DoFCellAccessor<DoFHandlerType<dim, spacedim>, lda>> &cell,
+  const typename Triangulation<dim, spacedim>::face_iterator &             face,
+  const typename Triangulation<dim, spacedim>::face_iterator &subface)
+{
+  reinit(cell,
+         cell->face_iterator_to_index(face),
+         face->child_iterator_to_index(subface));
+}
+
+
+
 template <int dim, int spacedim>
 void
 FESubfaceValues<dim, spacedim>::reinit(
@@ -5005,6 +5045,20 @@ FESubfaceValues<dim, spacedim>::reinit(
 
 
 
+template <int dim, int spacedim>
+void
+FESubfaceValues<dim, spacedim>::reinit(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const typename Triangulation<dim, spacedim>::face_iterator &face,
+  const typename Triangulation<dim, spacedim>::face_iterator &subface)
+{
+  reinit(cell,
+         cell->face_iterator_to_index(face),
+         face->child_iterator_to_index(subface));
+}
+
+
+
 template <int dim, int spacedim>
 void
 FESubfaceValues<dim, spacedim>::do_reinit(const unsigned int face_no,
diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in
index 5eca3410a2..b2cfec8f4f 100644
--- a/source/fe/fe_values.inst.in
+++ b/source/fe/fe_values.inst.in
@@ -67,13 +67,35 @@ for (dof_handler : DOFHANDLER_TEMPLATES; deal_II_dimension : DIMENSIONS;
         DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>,
                         lda>> &,
       unsigned int);
+
+    template void
+    FEFaceValues<deal_II_dimension, deal_II_space_dimension>::reinit(
+      const TriaIterator<
+        DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>,
+                        lda>> &,
+      const TriaIterator<TriaAccessor<deal_II_dimension - 1,
+                                      deal_II_dimension,
+                                      deal_II_space_dimension>> &);
+
     template void
     FESubfaceValues<deal_II_dimension, deal_II_space_dimension>::reinit(
       const TriaIterator<
         DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>,
                         lda>> &,
-      unsigned int,
-      unsigned int);
+      const unsigned int,
+      const unsigned int);
+
+    template void
+    FESubfaceValues<deal_II_dimension, deal_II_space_dimension>::reinit(
+      const TriaIterator<
+        DoFCellAccessor<dof_handler<deal_II_dimension, deal_II_space_dimension>,
+                        lda>> &,
+      const TriaIterator<TriaAccessor<deal_II_dimension - 1,
+                                      deal_II_dimension,
+                                      deal_II_space_dimension>> &,
+      const TriaIterator<TriaAccessor<deal_II_dimension - 1,
+                                      deal_II_dimension,
+                                      deal_II_space_dimension>> &);
 #endif
   }
 
diff --git a/tests/fe/shapes.h b/tests/fe/shapes.h
index dbe5cda337..d61631e06d 100644
--- a/tests/fe/shapes.h
+++ b/tests/fe/shapes.h
@@ -157,6 +157,17 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
                            q,
                            UpdateFlags(uflags | update_quadrature_points));
 
+  // Test the iterator-based initialization routines too:
+  FEFaceValues<dim>    fe2(mapping,
+                        finel,
+                        q,
+                        UpdateFlags(uflags | update_quadrature_points));
+  FESubfaceValues<dim> sub2(mapping,
+                            finel,
+                            q,
+                            UpdateFlags(uflags | update_quadrature_points));
+
+
   sprintf(fname, "Face%dd-%s", dim, name);
   deallog.push(fname);
 
@@ -165,6 +176,7 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
       if (!c->face(f)->has_children())
         {
           fe.reinit(c, f);
+          fe2.reinit(c, c->face(f));
 
           unsigned int k = 0;
           for (unsigned int my = 0; my <= ((dim > 2) ? div : 0); ++my)
@@ -187,23 +199,47 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
                               c)
                             {
                               if (uflags & update_values)
-                                AssertThrow((fe.shape_value(i, k) ==
-                                             fe.shape_value_component(i, k, c)),
-                                            ExcInternalError());
+                                {
+                                  AssertThrow(
+                                    (fe.shape_value(i, k) ==
+                                     fe.shape_value_component(i, k, c)),
+                                    ExcInternalError());
+                                  AssertThrow((fe.shape_value(i, k) ==
+                                               fe2.shape_value(i, k)),
+                                              ExcInternalError());
+                                }
                               if (uflags & update_gradients)
-                                AssertThrow((fe.shape_grad(i, k) ==
-                                             fe.shape_grad_component(i, k, c)),
-                                            ExcInternalError());
+                                {
+                                  AssertThrow(
+                                    (fe.shape_grad(i, k) ==
+                                     fe.shape_grad_component(i, k, c)),
+                                    ExcInternalError());
+                                  AssertThrow((fe.shape_grad(i, k) ==
+                                               fe2.shape_grad(i, k)),
+                                              ExcInternalError());
+                                }
                               if (uflags & update_hessians)
-                                AssertThrow(
-                                  (fe.shape_hessian(i, k) ==
-                                   fe.shape_hessian_component(i, k, c)),
-                                  ExcInternalError());
+                                {
+                                  AssertThrow(
+                                    (fe.shape_hessian(i, k) ==
+                                     fe.shape_hessian_component(i, k, c)),
+                                    ExcInternalError());
+                                  AssertThrow((fe.shape_hessian(i, k) ==
+                                               fe2.shape_hessian(i, k)),
+                                              ExcInternalError());
+                                }
                               if (uflags & update_3rd_derivatives)
-                                AssertThrow(
-                                  (fe.shape_3rd_derivative(i, k) ==
-                                   fe.shape_3rd_derivative_component(i, k, c)),
-                                  ExcInternalError());
+                                {
+                                  AssertThrow(
+                                    (fe.shape_3rd_derivative(i, k) ==
+                                     fe.shape_3rd_derivative_component(i,
+                                                                       k,
+                                                                       c)),
+                                    ExcInternalError());
+                                  AssertThrow((fe.shape_3rd_derivative(i, k) ==
+                                               fe2.shape_3rd_derivative(i, k)),
+                                              ExcInternalError());
+                                }
                             }
                           else
                             {
@@ -239,6 +275,7 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
                ++s)
             {
               sub.reinit(c, f, s);
+              sub2.reinit(c, c->face(f), c->face(f)->child(s));
 
               unsigned int k = 0;
               for (unsigned int my = 0; my <= ((dim > 2) ? div : 0); ++my)
@@ -267,23 +304,29 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
                                       const double v1 = sub.shape_value(i, k),
                                                    v2 =
                                                      sub.shape_value_component(
-                                                       i, k, c);
+                                                       i, k, c),
+                                                   v3 = sub2.shape_value(i, k);
                                       Assert(v1 == v2, ExcInternalError());
+                                      Assert(v1 == v3, ExcInternalError());
                                     }
                                   if (uflags & update_gradients)
                                     {
                                       const Tensor<1, dim>
                                         g1 = sub.shape_grad(i, k),
-                                        g2 = sub.shape_grad_component(i, k, c);
+                                        g2 = sub.shape_grad_component(i, k, c),
+                                        g3 = sub2.shape_grad(i, k);
                                       Assert(g1 == g2, ExcInternalError());
+                                      Assert(g1 == g3, ExcInternalError());
                                     }
                                   if (uflags & update_hessians)
                                     {
                                       const Tensor<2, dim>
                                         s1 = sub.shape_hessian(i, k),
                                         s2 =
-                                          sub.shape_hessian_component(i, k, c);
+                                          sub.shape_hessian_component(i, k, c),
+                                        s3 = sub2.shape_hessian(i, k);
                                       Assert(s1 == s2, ExcInternalError());
+                                      Assert(s1 == s3, ExcInternalError());
                                     }
                                   if (uflags & update_3rd_derivatives)
                                     {
@@ -292,8 +335,11 @@ plot_face_shape_functions(Mapping<dim> &      mapping,
                                         t2 =
                                           sub.shape_3rd_derivative_component(i,
                                                                              k,
-                                                                             c);
+                                                                             c),
+                                        t3 = sub2.shape_3rd_derivative(i, k);
+
                                       Assert(t1 == t2, ExcInternalError());
+                                      Assert(t1 == t3, ExcInternalError());
                                     }
                                 }
                               else
-- 
2.39.5