From 1e9da15d130bf57677a7fab73c505a06bbcf2b57 Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Sat, 25 Jul 2015 18:55:21 +0200
Subject: [PATCH] Fixed documentation and a few const.

---
 doc/news/changes.h                            |  8 ++++++++
 include/deal.II/numerics/vector_tools.h       | 19 +++++++++---------
 .../deal.II/numerics/vector_tools.templates.h | 20 +++++++++++++++----
 3 files changed, 34 insertions(+), 13 deletions(-)

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 283865bdcd..5378458e53 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -80,6 +80,14 @@ inconvenience this causes.
 
 
 <ol>
+  <li> New: VectorTools::get_position_vector now works with arbitrary
+  FESystems, provided that the geometrical components are primitive, 
+  and that you provide a component mask to select what components of 
+  the finite element to use for the geometrical interpolation.
+  <br>
+  (Luca Heltai, 2015/07/25)
+  </li>
+
   <li> New: FESystem now does some work in parallel if your system
   has multiple processors.
   <br>
diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 142950baae..7d0608694a 100644
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -2345,17 +2345,15 @@ namespace VectorTools
    */
   //@{
   /**
-   * Given a DoFHandler containing at least a spacedim vector field, this
-   * function interpolates the Triangulation at the support points of a FE_Q()
-   * finite element of the same degree as dh->get_fe().degree.
+   * Given a DoFHandler containing at least a spacedim vector field,
+   * this function interpolates the Triangulation at the support
+   * points of a FE_Q() finite element of the same degree as the
+   * degree of the required components.
    *
    * Curved manifold are respected, and the resulting VECTOR will be
-   * geometrically consistent.
-   *
-   * The resulting map is guaranteed to be interpolatory at the support points
-   * of a FE_Q() finite element of the same degree as dh->get_fe().degree.
-   * Notice that this may or may not be meaningful, depending on the
-   * FiniteElement you have distributed in dh.
+   * geometrically consistent. The resulting map is guaranteed to be
+   * interpolatory at the support points of a FE_Q() finite element of
+   * the same degree as the degree of the required components.
    *
    * If the underlying finite element is an FE_Q(1)^spacedim, then the
    * resulting VECTOR is a finite element field representation of the vertices
@@ -2367,6 +2365,9 @@ namespace VectorTools
    * the first spacedim components of the FiniteElement are assumed to
    * represent the geometry of the problem.
    *
+   * This function is only implemented for FiniteElements where the
+   * specified components are primitive.
+   *
    * @author Luca Heltai, 2015
    */
   template<class DH, class VECTOR>
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 9a207bf594..746d8525df 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -7185,13 +7185,19 @@ namespace VectorTools
         for (unsigned int i=0; i<fe_mask.size(); ++i)
           if (fe_mask[i])
             {
-              unsigned int base_i = fe_system->component_to_base_index(i).first;
+              const unsigned int base_i = fe_system->component_to_base_index(i).first;
               Assert(degree == numbers::invalid_unsigned_int ||
                      degree == fe_system->base_element(base_i).degree,
                      ExcNotImplemented());
+              Assert(fe_system->base_element(base_i).is_primitive(),
+                     ExcNotImplemented());
               degree = fe_system->base_element(base_i).degree;
             }
 
+        // We create an intermediate FE_Q vector space, and then
+        // interpolate from that vector space to this one, by
+        // carefully selecting the right components.
+
         FESystem<dim,spacedim> feq(FE_Q<dim,spacedim>(degree), spacedim);
         DH dhq(dh.get_tria());
         dhq.distribute_dofs(feq);
@@ -7213,11 +7219,17 @@ namespace VectorTools
         // v_j phi_j(q_i) = w_k psi_k(q_i) = w_k delta_ki = w_i
         //
         // where psi_k are the basis functions for fe_q, and phi_i are
-        // the basis functions of the target space.
+        // the basis functions of the target space while q_i are the
+        // support points for the fe_q space. With this choice of
+        // interpolation points, on the matrices is the identity
+        // matrix, and we have to invert only one matrix. The
+        // resulting vector will be interpolatory at the support
+        // points of fe_q, even if the finite element does not have
+        // support points.
         //
         // Morally, we should invert the matrix T_ij = phi_i(q_j),
         // however in general this matrix is not invertible, since
-        // there may by components which do not contribute to the
+        // there may be components which do not contribute to the
         // displacement vector. Since we are not interested in those
         // components, we construct a square matrix with the same
         // number of components of the FE_Q system. The FE_Q system
@@ -7241,7 +7253,7 @@ namespace VectorTools
 
         for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
           {
-            unsigned int comp_j = fe.system_to_component_index(j).first;
+            const unsigned int comp_j = fe.system_to_component_index(j).first;
             if (fe_mask[comp_j])
               for (unsigned int i=0; i<points.size(); ++i)
                 {
-- 
2.39.5