From da57e2ebbee31b5740962a69df6ef0da6398be0a Mon Sep 17 00:00:00 2001
From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 18 Jan 2014 19:30:26 +0000
Subject: [PATCH] Implement more functionality for get_interpolated_dof_values.

git-svn-id: https://svn.dealii.org/trunk@32249 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/doc/news/changes.h                   | 11 ++++-
 deal.II/include/deal.II/dofs/dof_accessor.h  | 32 +++++++++----
 deal.II/source/dofs/dof_accessor_get.cc      | 50 ++++++++++++++++----
 deal.II/source/dofs/dof_accessor_get.inst.in | 18 ++++---
 deal.II/source/dofs/dof_accessor_set.cc      | 23 +++++----
 deal.II/source/dofs/dof_accessor_set.inst.in | 18 ++++---
 6 files changed, 112 insertions(+), 40 deletions(-)

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 8fd8bc22e3..255a8174ba 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -91,13 +91,20 @@ inconvenience this causes.
   (Timo Heister, 2014/01/17)
   </li>
 
-  <li> Fixed: DoFCellAccessor::set_dof_values_by_interpolation and
+  <li> Fixed/new: DoFCellAccessor::set_dof_values_by_interpolation and
   DoFCellAccessor::get_interpolated_dof_values could previously be
   called for hp::DoFHandler objects on cells that are non-active. This
   makes no sense since these cells have no associated finite element
   space. Doing so now raises an exception.
   <br>
-  (Wolfgang Bangerth, 2014/01/17)
+  However, there are legitimate cases where one may want to interpolate
+  from children to a parent's finite element space. Since in the hp
+  case no finite element space is naturally associated with an inactive
+  cell, it is now possible to pass an explicit finite element index
+  argument to the function specifying which element of an hp::FECollection
+  object describes the space onto which you want to interpolate.
+  <br>
+  (Wolfgang Bangerth, 2014/01/18)
   </li>
 
   <li> Fixed: The methods IndexSet::do_compress() and
diff --git a/deal.II/include/deal.II/dofs/dof_accessor.h b/deal.II/include/deal.II/dofs/dof_accessor.h
index 461d8d912a..885eeff0ee 100644
--- a/deal.II/include/deal.II/dofs/dof_accessor.h
+++ b/deal.II/include/deal.II/dofs/dof_accessor.h
@@ -1407,27 +1407,42 @@ public:
   /**
    * Return the interpolation of the given finite element function to
    * the present cell. In the simplest case, the cell is a terminal
-   * one, i.e. has no children; then, the returned value is the vector
-   * of nodal values on that cell. You could then as well get the
+   * one, i.e., it has no children; then, the returned value is the vector
+   * of nodal values on that cell. You could as well get the
    * desired values through the @p get_dof_values function. In the
    * other case, when the cell has children, we use the restriction
    * matrices provided by the finite element class to compute the
    * interpolation from the children to the present cell.
    *
-   * It is assumed that both vectors already have the right size
+   * If the cell is part of a hp::DoFHandler object, cells only have an
+   * associated finite element space if they are active. However, this
+   * function is supposed to also provide information on inactive cells
+   * with children. Consequently, it carries a third argument that can be
+   * used in the hp context that denotes the finite element space we are
+   * supposed to interpolate onto. If the cell is active, this function
+   * then obtains the finite element function from the <code>values</code>
+   * vector on this cell and interpolates it onto the space described
+   * by the <code>fe_index</code>th element of the hp::FECollection associated
+   * with the hp::DoFHandler of which this cell is a part of. If the cell
+   * is not active, then we first perform this interpolation on all of its
+   * terminal children and then interpolate this function down to the
+   * cell requested keeping the function space the same.
+   *
+   * It is assumed that both input vectors already have the right size
    * beforehand.
    *
-   * Unlike the get_dof_values() function, this function works on
-   * cells rather than to lines, quads, and hexes, since interpolation
+   * @note Unlike the get_dof_values() function, this function is only available
+   * on cells, rather than on lines, quads, and hexes, since interpolation
    * is presently only provided for cells by the finite element
    * classes.
    */
   template <class InputVector, typename number>
   void get_interpolated_dof_values (const InputVector &values,
-                                    Vector<number>    &interpolated_values) const;
+                                    Vector<number>    &interpolated_values,
+                                    const unsigned int fe_index = DH::default_fe_index) const;
 
   /**
-   * This, again, is the counterpart to get_interpolated_dof_values():
+   * This function is the counterpart to get_interpolated_dof_values():
    * you specify the dof values on a cell and these are interpolated
    * to the children of the present cell and set on the terminal
    * cells.
@@ -1476,7 +1491,8 @@ public:
    */
   template <class OutputVector, typename number>
   void set_dof_values_by_interpolation (const Vector<number> &local_values,
-                                        OutputVector         &values) const;
+                                        OutputVector         &values,
+                                        const unsigned int fe_index = DH::default_fe_index) const;
 
   /**
    * Distribute a local (cell based) vector to a global one by mapping
diff --git a/deal.II/source/dofs/dof_accessor_get.cc b/deal.II/source/dofs/dof_accessor_get.cc
index 33639cf5bc..427e7e9f2c 100644
--- a/deal.II/source/dofs/dof_accessor_get.cc
+++ b/deal.II/source/dofs/dof_accessor_get.cc
@@ -43,12 +43,37 @@ template <class InputVector, typename number>
 void
 DoFCellAccessor<DH,lda>::
 get_interpolated_dof_values (const InputVector &values,
-                             Vector<number>    &interpolated_values) const
+                             Vector<number>    &interpolated_values,
+                             const unsigned int fe_index) const
 {
   if (!this->has_children())
     // if this cell has no children: simply return the exact values on this
-    // cell
-    this->get_dof_values (values, interpolated_values);
+    // cell unless the finite element we need to interpolate to is different than
+    // the one we have on the current cell
+    {
+      if ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
+           (this->dof_handler)
+           != 0)
+          ||
+          (fe_index == this->active_fe_index())
+          ||
+          (fe_index == DH::default_fe_index))
+        this->get_dof_values (values, interpolated_values);
+      else
+        {
+          // well, here we need to first get the values from the current
+          // cell and then interpolate it to the element requested. this
+          // can clearly only happen for hp::DoFHandler objects
+          Vector<number> tmp (this->get_fe().dofs_per_cell);
+          this->get_dof_values (values, tmp);
+
+          FullMatrix<double> interpolation (this->dof_handler->get_fe()[fe_index].dofs_per_cell,
+                                            this->get_fe().dofs_per_cell);
+          this->dof_handler->get_fe()[fe_index].get_interpolation_matrix (this->get_fe(),
+                                                                          interpolation);
+          interpolation.vmult (interpolated_values, tmp);
+        }
+    }
   else
     // otherwise obtain them from the children
     {
@@ -57,17 +82,21 @@ get_interpolated_dof_values (const InputVector &values,
       // context, we can simply assume that the FE space to which we
       // want to interpolate is the same as for all elements in the
       // mesh). consequently, we cannot interpolate from children's FE
-      // space to this cell's (unknown) FE space
+      // space to this cell's (unknown) FE space unless an explicit
+      // fe_index is given
       Assert ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
 	       (this->dof_handler)
-	       != 0),
+	       != 0)
+              ||
+	      (fe_index != DH::default_fe_index),
 	      ExcMessage ("You cannot call this function on non-active cells "
-			  "of hp::DoFHandler objects because they do not have "
+			  "of hp::DoFHandler objects unless you provide an explicit "
+			  "finite element index because they do not have naturally "
 			  "associated finite element spaces associated: degrees "
 			  "of freedom are only distributed on active cells for which "
 			  "the active_fe_index has been set."));
 
-      const FiniteElement<dim,spacedim> &fe            = this->get_fe();
+      const FiniteElement<dim,spacedim> &fe            = this->get_dof_handler().get_fe()[fe_index];
       const unsigned int                 dofs_per_cell = fe.dofs_per_cell;
 
       Assert (this->dof_handler != 0,
@@ -117,9 +146,12 @@ get_interpolated_dof_values (const InputVector &values,
       for (unsigned int child=0; child<this->n_children(); ++child)
         {
           // get the values from the present child, if necessary by
-          // interpolation itself
+          // interpolation itself either from its own children or
+          // by interpolating from the finite element on an active
+          // child to the finite element space requested here
           this->child(child)->get_interpolated_dof_values (values,
-                                                           tmp1);
+                                                           tmp1,
+                                                           fe_index);
           // interpolate these to the mother cell
           fe.get_restriction_matrix(child, this->refinement_case()).vmult (tmp2, tmp1);
 
diff --git a/deal.II/source/dofs/dof_accessor_get.inst.in b/deal.II/source/dofs/dof_accessor_get.inst.in
index 75d1009f06..c2e4b79151 100644
--- a/deal.II/source/dofs/dof_accessor_get.inst.in
+++ b/deal.II/source/dofs/dof_accessor_get.inst.in
@@ -21,14 +21,16 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<DoFHandler<deal_II_dimension>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #if deal_II_dimension != 3
 
     template
       void
       DoFCellAccessor<DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -37,7 +39,8 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<DoFHandler<1,3>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -49,14 +52,16 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<hp::DoFHandler<deal_II_dimension>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #if deal_II_dimension != 3
 
     template
       void
       DoFCellAccessor<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -65,7 +70,8 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<hp::DoFHandler<1,3>, lda>::get_interpolated_dof_values
-      (const VEC&, Vector<SCALAR>&) const;
+      (const VEC&, Vector<SCALAR>&,
+       const unsigned int fe_index) const;
 
 #endif
   }
diff --git a/deal.II/source/dofs/dof_accessor_set.cc b/deal.II/source/dofs/dof_accessor_set.cc
index 4438c6b5b5..bfbedda11c 100644
--- a/deal.II/source/dofs/dof_accessor_set.cc
+++ b/deal.II/source/dofs/dof_accessor_set.cc
@@ -44,7 +44,8 @@ template <class OutputVector, typename number>
 void
 DoFCellAccessor<DH,lda>::
 set_dof_values_by_interpolation (const Vector<number> &local_values,
-                                 OutputVector         &values) const
+                                 OutputVector         &values,
+                                 const unsigned int fe_index) const
 {
   if (!this->has_children())
     // if this cell has no children: simply set the values on this cell
@@ -57,15 +58,19 @@ set_dof_values_by_interpolation (const Vector<number> &local_values,
       // context, we can simply assume that the FE space to from which we
       // want to interpolate is the same as for all elements in the
       // mesh). consequently, we cannot interpolate to children's FE
-      // space from this cell's (unknown) FE space
+      // space from this cell's (unknown) FE space unless an explicit
+      // fe_index is given
       Assert ((dynamic_cast<DoFHandler<DH::dimension,DH::space_dimension>*>
-	       (this->dof_handler)
-	       != 0),
-	      ExcMessage ("You cannot call this function on non-active cells "
-			  "of hp::DoFHandler objects because they do not have "
-			  "associated finite element spaces associated: degrees "
-			  "of freedom are only distributed on active cells for which "
-			  "the active_fe_index has been set."));
+               (this->dof_handler)
+               != 0)
+              ||
+              (fe_index != DH::default_fe_index),
+              ExcMessage ("You cannot call this function on non-active cells "
+                          "of hp::DoFHandler objects unless you provide an explicit "
+                          "finite element index because they do not have naturally "
+                          "associated finite element spaces associated: degrees "
+                          "of freedom are only distributed on active cells for which "
+                          "the active_fe_index has been set."));
 
       const unsigned int dofs_per_cell = this->get_fe().dofs_per_cell;
 
diff --git a/deal.II/source/dofs/dof_accessor_set.inst.in b/deal.II/source/dofs/dof_accessor_set.inst.in
index 74948318e8..4b54558287 100644
--- a/deal.II/source/dofs/dof_accessor_set.inst.in
+++ b/deal.II/source/dofs/dof_accessor_set.inst.in
@@ -21,14 +21,16 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<DoFHandler<deal_II_dimension>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #if deal_II_dimension != 3
 
     template
       void
       DoFCellAccessor<DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -37,7 +39,8 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<DoFHandler<1,3>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -49,14 +52,16 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<hp::DoFHandler<deal_II_dimension>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #if deal_II_dimension != 3
 
     template
       void
       DoFCellAccessor<hp::DoFHandler<deal_II_dimension,deal_II_dimension+1>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #endif
 
@@ -65,7 +70,8 @@ for (VEC : SERIAL_VECTORS; SCALAR : REAL_SCALARS; deal_II_dimension : DIMENSIONS
     template
       void
       DoFCellAccessor<hp::DoFHandler<1,3>, lda>::set_dof_values_by_interpolation
-      (const Vector<SCALAR>&, VEC&) const;
+      (const Vector<SCALAR>&, VEC&,
+       const unsigned int fe_index) const;
 
 #endif
   }
-- 
2.39.5