From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 13 Mar 2017 23:41:06 +0000 (-0600)
Subject: Implement the new FiniteElement::convert_generalized_support_point_values_to_nodal_va... 
X-Git-Tag: v8.5.0-rc1~16^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cf3ac8ac32d96d9ac686e9e6b2ec63844cbed6f0;p=dealii.git

Implement the new FiniteElement::convert_generalized_support_point_values_to_nodal_values interface.

This is easily done by copying the implementation of the
existing FiniteElement::interpolate() functions in the
various derived classes, removing the 'offset' argument,
and reordering+renaming the other two arguments.
---

diff --git a/include/deal.II/fe/fe_abf.h b/include/deal.II/fe/fe_abf.h
index 24e3f0b8a2..8a6f0a7c16 100644
--- a/include/deal.II/fe/fe_abf.h
+++ b/include/deal.II/fe/fe_abf.h
@@ -123,6 +123,12 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double>                &local_dofs,
                            const std::vector<double>          &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_bdm.h b/include/deal.II/fe/fe_bdm.h
index e64dd95172..b6a1d384dc 100644
--- a/include/deal.II/fe/fe_bdm.h
+++ b/include/deal.II/fe/fe_bdm.h
@@ -72,6 +72,12 @@ public:
 
   virtual FiniteElement<dim> *clone () const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double>                &local_dofs,
                            const std::vector<double>          &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_nedelec.h b/include/deal.II/fe/fe_nedelec.h
index 52a3b0af9d..9ea468c122 100644
--- a/include/deal.II/fe/fe_nedelec.h
+++ b/include/deal.II/fe/fe_nedelec.h
@@ -255,6 +255,12 @@ public:
   get_prolongation_matrix (const unsigned int child,
                            const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate (std::vector<double> &local_dofs,
                             const std::vector<double> &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h
index 046d6fd676..27ca79a0b1 100644
--- a/include/deal.II/fe/fe_q_bubbles.h
+++ b/include/deal.II/fe/fe_q_bubbles.h
@@ -104,6 +104,12 @@ public:
    */
   virtual std::string get_name () const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double>       &local_dofs,
                            const std::vector<double> &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_q_dg0.h b/include/deal.II/fe/fe_q_dg0.h
index fa3467301f..ddb27bdbc8 100644
--- a/include/deal.II/fe/fe_q_dg0.h
+++ b/include/deal.II/fe/fe_q_dg0.h
@@ -258,6 +258,12 @@ public:
    */
   virtual std::string get_name () const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double>       &local_dofs,
                            const std::vector<double> &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_rannacher_turek.h b/include/deal.II/fe/fe_rannacher_turek.h
index 5c138583be..80c5b7a2ae 100644
--- a/include/deal.II/fe/fe_rannacher_turek.h
+++ b/include/deal.II/fe/fe_rannacher_turek.h
@@ -67,6 +67,12 @@ public:
   virtual std::string get_name() const;
   virtual FiniteElement<dim> *clone() const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double> &local_dofs,
                            const std::vector<double> &values) const DEAL_II_DEPRECATED;
 
diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h
index 1e76c5b4f0..bb7151844f 100644
--- a/include/deal.II/fe/fe_raviart_thomas.h
+++ b/include/deal.II/fe/fe_raviart_thomas.h
@@ -128,6 +128,12 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
+
   virtual void interpolate(std::vector<double>                &local_dofs,
                            const std::vector<double>          &values) const DEAL_II_DEPRECATED;
 
@@ -260,6 +266,11 @@ public:
 
   virtual FiniteElement<dim> *clone () const;
 
+  // documentation inherited from the base class
+  virtual
+  void
+  convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                            std::vector<double>                &nodal_values) const;
   virtual void interpolate(std::vector<double>                &local_dofs,
                            const std::vector<double> &values) const DEAL_II_DEPRECATED;
 
diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc
index 3daec17566..86d0b0b8ac 100644
--- a/source/fe/fe_abf.cc
+++ b/source/fe/fe_abf.cc
@@ -508,6 +508,68 @@ FE_ABF<dim>::has_support_on_face (const unsigned int shape_index,
 
 
 
+template <int dim>
+void
+FE_ABF<dim>::
+convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                          std::vector<double>                &nodal_values) const
+{
+  Assert (support_point_values.size() == this->generalized_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size()));
+  Assert (support_point_values[0].size() == this->n_components(),
+          ExcDimensionMismatch(support_point_values[0].size(), this->n_components()));
+  Assert (nodal_values.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell));
+
+  std::fill(nodal_values.begin(), nodal_values.end(), 0.);
+
+  const unsigned int n_face_points = boundary_weights.size(0);
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int k=0; k<n_face_points; ++k)
+      for (unsigned int i=0; i<boundary_weights.size(1); ++i)
+        {
+          nodal_values[i+face*this->dofs_per_face] += boundary_weights(k,i)
+                                                      * support_point_values[face*n_face_points+k][GeometryInfo<dim>::unit_normal_direction[face]];
+        }
+
+  const unsigned int start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+  const unsigned int start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+
+  for (unsigned int k=0; k<interior_weights.size(0); ++k)
+    for (unsigned int i=0; i<interior_weights.size(1); ++i)
+      for (unsigned int d=0; d<dim; ++d)
+        nodal_values[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * support_point_values[k+start_cell_points][d];
+
+  const unsigned int start_abf_dofs = start_cell_dofs + interior_weights.size(1) * dim;
+
+  // Cell integral of ABF terms
+  for (unsigned int k=0; k<interior_weights_abf.size(0); ++k)
+    for (unsigned int i=0; i<interior_weights_abf.size(1); ++i)
+      for (unsigned int d=0; d<dim; ++d)
+        nodal_values[start_abf_dofs+i] += interior_weights_abf(k,i,d) * support_point_values[k+start_cell_points][d];
+
+  // Face integral of ABF terms
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    {
+      double n_orient = (double) GeometryInfo<dim>::unit_normal_orientation[face];
+      for (unsigned int fp=0; fp < n_face_points; ++fp)
+        {
+          // TODO: Check what the face_orientation, face_flip and face_rotation have to be in 3D
+          unsigned int k = QProjector<dim>::DataSetDescriptor::face (face, false, false, false, n_face_points);
+          for (unsigned int i=0; i<boundary_weights_abf.size(1); ++i)
+            nodal_values[start_abf_dofs+i] += n_orient * boundary_weights_abf(k + fp, i)
+                                              * support_point_values[k + fp][GeometryInfo<dim>::unit_normal_direction[face]];
+        }
+    }
+
+  // TODO: Check if this "correction" can be removed.
+  for (unsigned int i=0; i<boundary_weights_abf.size(1); ++i)
+    if (std::fabs (nodal_values[start_abf_dofs+i]) < 1.0e-16)
+      nodal_values[start_abf_dofs+i] = 0.0;
+}
+
+
+
 template <int dim>
 void
 FE_ABF<dim>::interpolate(
@@ -580,6 +642,8 @@ FE_ABF<dim>::interpolate(
       local_dofs[start_abf_dofs+i] = 0.0;
 }
 
+
+
 template <int dim>
 void
 FE_ABF<dim>::interpolate(
@@ -641,6 +705,7 @@ FE_ABF<dim>::interpolate(
 }
 
 
+
 template <int dim>
 std::size_t
 FE_ABF<dim>::memory_consumption () const
diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index 120d2c8808..557c8fc7c3 100644
--- a/source/fe/fe_bdm.cc
+++ b/source/fe/fe_bdm.cc
@@ -122,6 +122,77 @@ FE_BDM<dim>::clone() const
 
 
 
+template <int dim>
+void
+FE_BDM<dim>::
+convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                          std::vector<double> &nodal_values) const
+{
+  Assert (support_point_values.size() == this->generalized_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size()));
+  AssertDimension (support_point_values[0].size(), dim);
+  Assert (nodal_values.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell));
+
+  // First do interpolation on faces. There, the component evaluated
+  // depends on the face direction and orientation.
+
+  // The index of the first dof on this face or the cell
+  unsigned int dbase = 0;
+  // The index of the first generalized support point on this face or the cell
+  unsigned int pbase = 0;
+  for (unsigned int f = 0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+    {
+      // Old version with no moments in 2D. See comment below in
+      // initialize_support_points()
+      if (test_values_face.size() == 0)
+        {
+          for (unsigned int i=0; i<this->dofs_per_face; ++i)
+            nodal_values[dbase+i] = support_point_values[pbase+i][GeometryInfo<dim>::unit_normal_direction[f]];
+          pbase += this->dofs_per_face;
+        }
+      else
+        {
+          for (unsigned int i=0; i<this->dofs_per_face; ++i)
+            {
+              double s = 0.;
+              for (unsigned int k=0; k<test_values_face.size(); ++k)
+                s += support_point_values[pbase+k][GeometryInfo<dim>::unit_normal_direction[f]] * test_values_face[k][i];
+              nodal_values[dbase+i] = s;
+            }
+          pbase += test_values_face.size();
+        }
+      dbase += this->dofs_per_face;
+    }
+
+  AssertDimension (dbase, this->dofs_per_face * GeometryInfo<dim>::faces_per_cell);
+  AssertDimension (pbase, this->generalized_support_points.size() - test_values_cell.size());
+
+  // Done for BDM1
+  if (dbase == this->dofs_per_cell) return;
+
+  // What's missing are the interior
+  // degrees of freedom. In each
+  // point, we take all components of
+  // the solution.
+  Assert ((this->dofs_per_cell - dbase) % dim == 0, ExcInternalError());
+
+  for (unsigned int d=0; d<dim; ++d, dbase += test_values_cell[0].size())
+    {
+      for (unsigned int i=0; i<test_values_cell[0].size(); ++i)
+        {
+          double s = 0.;
+          for (unsigned int k=0; k<test_values_cell.size(); ++k)
+            s += support_point_values[pbase+k][d] * test_values_cell[k][i];
+          nodal_values[dbase+i] = s;
+        }
+    }
+
+  Assert (dbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
 template <int dim>
 void
 FE_BDM<dim>::interpolate(
diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index 761cd84c2c..39fef89974 100644
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -3086,6 +3086,934 @@ FE_Nedelec<dim>
   return this->restriction[refinement_case-1][child];
 }
 
+
+// Interpolate a function, which is given by
+// its values at the generalized support
+// points in the finite element space on the
+// reference cell.
+// This is done as usual by projection-based
+// interpolation.
+template <int dim>
+void
+FE_Nedelec<dim>::
+convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
+                                                          std::vector<double> &nodal_values) const
+{
+  const unsigned int deg = this->degree-1;
+  Assert (support_point_values.size () == this->generalized_support_points.size (),
+          ExcDimensionMismatch (support_point_values.size (),
+                                this->generalized_support_points.size ()));
+  Assert (support_point_values[0].size () == this->n_components (),
+          ExcDimensionMismatch (support_point_values[0].size (), this->n_components ()));
+  Assert (nodal_values.size () == this->dofs_per_cell,
+          ExcDimensionMismatch (nodal_values.size (), this->dofs_per_cell));
+  std::fill (nodal_values.begin (), nodal_values.end (), 0.0);
+
+  switch (dim)
+    {
+    case 2:
+    {
+      // Let us begin with the
+      // interpolation part.
+      const QGauss<dim - 1> reference_edge_quadrature (this->degree);
+      const unsigned int &
+      n_edge_points = reference_edge_quadrature.size ();
+
+      for (unsigned int i = 0; i < 2; ++i)
+        for (unsigned int j = 0; j < 2; ++j)
+          {
+            for (unsigned int q_point = 0; q_point < n_edge_points;
+                 ++q_point)
+              nodal_values[(i + 2 * j) * this->degree]
+              += reference_edge_quadrature.weight (q_point)
+                 * support_point_values[q_point + (i + 2 * j) * n_edge_points][1 - j];
+
+            // Add the computed support_point_values to the resulting vector
+            // only, if they are not
+            // too small.
+            if (std::abs (nodal_values[(i + 2 * j) * this->degree]) < 1e-14)
+              nodal_values[(i + 2 * j) * this->degree] = 0.0;
+          }
+
+      // If the degree is greater
+      // than 0, then we have still
+      // some higher order edge
+      // shape functions to
+      // consider.
+      // Here the projection part
+      // starts. The dof support_point_values
+      // are obtained by solving
+      // a linear system of
+      // equations.
+      if (this->degree-1 > 1)
+        {
+          // We start with projection
+          // on the higher order edge
+          // shape function.
+          const std::vector<Polynomials::Polynomial<double> > &
+          lobatto_polynomials
+            = Polynomials::Lobatto::generate_complete_basis
+              (this->degree);
+          FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
+          std::vector<Polynomials::Polynomial<double> >
+          lobatto_polynomials_grad (this->degree);
+
+          for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+               ++i)
+            lobatto_polynomials_grad[i]
+              = lobatto_polynomials[i + 1].derivative ();
+
+          // Set up the system matrix.
+          // This can be used for all
+          // edges.
+          for (unsigned int i = 0; i < system_matrix.m (); ++i)
+            for (unsigned int j = 0; j < system_matrix.n (); ++j)
+              for (unsigned int q_point = 0; q_point < n_edge_points;
+                   ++q_point)
+                system_matrix (i, j)
+                += boundary_weights (q_point, j)
+                   * lobatto_polynomials_grad[i + 1].value
+                   (this->generalized_face_support_points[q_point]
+                    (1));
+
+          FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
+
+          system_matrix_inv.invert (system_matrix);
+
+          const unsigned int
+          line_coordinate[GeometryInfo<2>::lines_per_cell]
+            = {1, 1, 0, 0};
+          Vector<double> system_rhs (system_matrix.m ());
+          Vector<double> solution (system_rhs.size ());
+
+          for (unsigned int line = 0;
+               line < GeometryInfo<dim>::lines_per_cell; ++line)
+            {
+              // Set up the right hand side.
+              system_rhs = 0;
+
+              for (unsigned int q_point = 0; q_point < n_edge_points;
+                   ++q_point)
+                {
+                  const double tmp
+                    = support_point_values[line * n_edge_points + q_point][line_coordinate[line]]
+                      - nodal_values[line * this->degree]
+                      * this->shape_value_component
+                      (line * this->degree,
+                       this->generalized_support_points[line
+                                                        * n_edge_points
+                                                        + q_point],
+                       line_coordinate[line]);
+
+                  for (unsigned int i = 0; i < system_rhs.size (); ++i)
+                    system_rhs (i) += boundary_weights (q_point, i) * tmp;
+                }
+
+              system_matrix_inv.vmult (solution, system_rhs);
+
+              // Add the computed support_point_values
+              // to the resulting vector
+              // only, if they are not
+              // too small.
+              for (unsigned int i = 0; i < solution.size (); ++i)
+                if (std::abs (solution (i)) > 1e-14)
+                  nodal_values[line * this->degree + i + 1] = solution (i);
+            }
+
+          // Then we go on to the
+          // interior shape
+          // functions. Again we
+          // set up the system
+          // matrix and use it
+          // for both, the
+          // horizontal and the
+          // vertical, interior
+          // shape functions.
+          const QGauss<dim> reference_quadrature (this->degree);
+          const unsigned int &
+          n_interior_points = reference_quadrature.size ();
+          const std::vector<Polynomials::Polynomial<double> > &
+          legendre_polynomials
+            = Polynomials::Legendre::generate_complete_basis (this->degree-1);
+
+          system_matrix.reinit ((this->degree-1) * this->degree,
+                                (this->degree-1) * this->degree);
+          system_matrix = 0;
+
+          for (unsigned int i = 0; i < this->degree; ++i)
+            for (unsigned int j = 0; j < this->degree-1; ++j)
+              for (unsigned int k = 0; k < this->degree; ++k)
+                for (unsigned int l = 0; l < this->degree-1; ++l)
+                  for (unsigned int q_point = 0;
+                       q_point < n_interior_points; ++q_point)
+                    system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+                    += reference_quadrature.weight (q_point)
+                       * legendre_polynomials[i].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points]
+                        (0))
+                       * lobatto_polynomials[j + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points]
+                        (1))
+                       * lobatto_polynomials_grad[k].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points]
+                        (0))
+                       * lobatto_polynomials[l + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points]
+                        (1));
+
+          system_matrix_inv.reinit (system_matrix.m (),
+                                    system_matrix.m ());
+          system_matrix_inv.invert (system_matrix);
+          // Set up the right hand side
+          // for the horizontal shape
+          // functions.
+          system_rhs.reinit (system_matrix_inv.m ());
+          system_rhs = 0;
+
+          for (unsigned int q_point = 0; q_point < n_interior_points;
+               ++q_point)
+            {
+              double tmp
+                = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
+                                       * n_edge_points][0];
+
+              for (unsigned int i = 0; i < 2; ++i)
+                for (unsigned int j = 0; j <= deg; ++j)
+                  tmp -= nodal_values[(i + 2) * this->degree + j]
+                         * this->shape_value_component
+                         ((i + 2) * this->degree + j,
+                          this->generalized_support_points[q_point
+                                                           + GeometryInfo<dim>::lines_per_cell
+                                                           * n_edge_points],
+                          0);
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  system_rhs (i * deg + j)
+                  += reference_quadrature.weight (q_point) * tmp
+                     * lobatto_polynomials_grad[i].value
+                     (this->generalized_support_points[q_point
+                                                       + GeometryInfo<dim>::lines_per_cell
+                                                       * n_edge_points]
+                      (0))
+                     * lobatto_polynomials[j + 2].value
+                     (this->generalized_support_points[q_point
+                                                       + GeometryInfo<dim>::lines_per_cell
+                                                       * n_edge_points]
+                      (1));
+            }
+
+          solution.reinit (system_matrix.m ());
+          system_matrix_inv.vmult (solution, system_rhs);
+
+          // Add the computed support_point_values
+          // to the resulting vector
+          // only, if they are not
+          // too small.
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              if (std::abs (solution (i * deg + j)) > 1e-14)
+                nodal_values[(i + GeometryInfo<dim>::lines_per_cell) * deg
+                             + j + GeometryInfo<dim>::lines_per_cell]
+                  = solution (i * deg + j);
+
+          system_rhs = 0;
+          // Set up the right hand side
+          // for the vertical shape
+          // functions.
+
+          for (unsigned int q_point = 0; q_point < n_interior_points;
+               ++q_point)
+            {
+              double tmp
+                = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
+                                       * n_edge_points][1];
+
+              for (unsigned int i = 0; i < 2; ++i)
+                for (unsigned int j = 0; j <= deg; ++j)
+                  tmp -= nodal_values[i * this->degree + j]
+                         * this->shape_value_component
+                         (i * this->degree + j,
+                          this->generalized_support_points[q_point
+                                                           + GeometryInfo<dim>::lines_per_cell
+                                                           * n_edge_points],
+                          1);
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  system_rhs (i * deg + j)
+                  += reference_quadrature.weight (q_point) * tmp
+                     * lobatto_polynomials_grad[i].value
+                     (this->generalized_support_points[q_point
+                                                       + GeometryInfo<dim>::lines_per_cell
+                                                       * n_edge_points]
+                      (1))
+                     * lobatto_polynomials[j + 2].value
+                     (this->generalized_support_points[q_point
+                                                       + GeometryInfo<dim>::lines_per_cell
+                                                       * n_edge_points]
+                      (0));
+            }
+
+          system_matrix_inv.vmult (solution, system_rhs);
+
+          // Add the computed support_point_values
+          // to the resulting vector
+          // only, if they are not
+          // too small.
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              if (std::abs (solution (i * deg + j)) > 1e-14)
+                nodal_values[i + (j + GeometryInfo<dim>::lines_per_cell
+                                  + deg) * this->degree]
+                  = solution (i * deg + j);
+        }
+
+      break;
+    }
+
+    case 3:
+    {
+      // Let us begin with the
+      // interpolation part.
+      const QGauss<1> reference_edge_quadrature (this->degree);
+      const unsigned int &
+      n_edge_points = reference_edge_quadrature.size ();
+
+      for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
+        {
+          for (unsigned int i = 0; i < 4; ++i)
+            nodal_values[(i + 8) * this->degree]
+            += reference_edge_quadrature.weight (q_point)
+               * support_point_values[q_point + (i + 8) * n_edge_points][2];
+
+          for (unsigned int i = 0; i < 2; ++i)
+            for (unsigned int j = 0; j < 2; ++j)
+              for (unsigned int k = 0; k < 2; ++k)
+                nodal_values[(i + 2 * (2 * j + k)) * this->degree]
+                += reference_edge_quadrature.weight (q_point)
+                   * support_point_values[q_point + (i + 2 * (2 * j + k))
+                                          * n_edge_points][1 - k];
+        }
+
+      // Add the computed support_point_values
+      // to the resulting vector
+      // only, if they are not
+      // too small.
+      for (unsigned int i = 0; i < 4; ++i)
+        if (std::abs (nodal_values[(i + 8) * this->degree]) < 1e-14)
+          nodal_values[(i + 8) * this->degree] = 0.0;
+
+      for (unsigned int i = 0; i < 2; ++i)
+        for (unsigned int j = 0; j < 2; ++j)
+          for (unsigned int k = 0; k < 2; ++k)
+            if (std::abs (nodal_values[(i + 2 * (2 * j + k)) * this->degree])
+                < 1e-14)
+              nodal_values[(i + 2 * (2 * j + k)) * this->degree] = 0.0;
+
+      // If the degree is greater
+      // than 0, then we have still
+      // some higher order shape
+      // functions to consider.
+      // Here the projection part
+      // starts. The dof support_point_values
+      // are obtained by solving
+      // a linear system of
+      // equations.
+      if (this->degree > 1)
+        {
+          // We start with projection
+          // on the higher order edge
+          // shape function.
+          const std::vector<Polynomials::Polynomial<double> > &
+          lobatto_polynomials
+            = Polynomials::Lobatto::generate_complete_basis
+              (this->degree);
+          FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
+          std::vector<Polynomials::Polynomial<double> >
+          lobatto_polynomials_grad (this->degree);
+
+          for (unsigned int i = 0; i < lobatto_polynomials_grad.size ();
+               ++i)
+            lobatto_polynomials_grad[i]
+              = lobatto_polynomials[i + 1].derivative ();
+
+          // Set up the system matrix.
+          // This can be used for all
+          // edges.
+          for (unsigned int i = 0; i < system_matrix.m (); ++i)
+            for (unsigned int j = 0; j < system_matrix.n (); ++j)
+              for (unsigned int q_point = 0; q_point < n_edge_points;
+                   ++q_point)
+                system_matrix (i, j)
+                += boundary_weights (q_point, j)
+                   * lobatto_polynomials_grad[i + 1].value
+                   (this->generalized_face_support_points[q_point]
+                    (1));
+
+          FullMatrix<double> system_matrix_inv (this->degree-1, this->degree-1);
+
+          system_matrix_inv.invert (system_matrix);
+
+          const unsigned int
+          line_coordinate[GeometryInfo<3>::lines_per_cell]
+            = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
+          Vector<double> system_rhs (system_matrix.m ());
+          Vector<double> solution (system_rhs.size ());
+
+          for (unsigned int line = 0;
+               line < GeometryInfo<dim>::lines_per_cell; ++line)
+            {
+              // Set up the right hand side.
+              system_rhs = 0;
+
+              for (unsigned int q_point = 0; q_point < this->degree; ++q_point)
+                {
+                  const double tmp
+                    = support_point_values[line * this->degree
+                                           + q_point][line_coordinate[line]]
+                      - nodal_values[line * this->degree]
+                      * this->shape_value_component
+                      (line * this->degree,
+                       this->generalized_support_points[line
+                                                        * this->degree
+                                                        + q_point],
+                       line_coordinate[line]);
+
+                  for (unsigned int i = 0; i < system_rhs.size (); ++i)
+                    system_rhs (i) += boundary_weights (q_point, i)
+                                      * tmp;
+                }
+
+              system_matrix_inv.vmult (solution, system_rhs);
+
+              // Add the computed values
+              // to the resulting vector
+              // only, if they are not
+              // too small.
+              for (unsigned int i = 0; i < solution.size (); ++i)
+                if (std::abs (solution (i)) > 1e-14)
+                  nodal_values[line * this->degree + i + 1] = solution (i);
+            }
+
+          // Then we go on to the
+          // face shape functions.
+          // Again we set up the
+          // system matrix and
+          // use it for both, the
+          // horizontal and the
+          // vertical, shape
+          // functions.
+          const std::vector<Polynomials::Polynomial<double> > &
+          legendre_polynomials
+            = Polynomials::Legendre::generate_complete_basis (this->degree-1);
+          const unsigned int n_face_points = n_edge_points * n_edge_points;
+
+          system_matrix.reinit ((this->degree-1) * this->degree,
+                                (this->degree-1) * this->degree);
+          system_matrix = 0;
+
+          for (unsigned int i = 0; i < this->degree; ++i)
+            for (unsigned int j = 0; j < this->degree-1; ++j)
+              for (unsigned int k = 0; k < this->degree; ++k)
+                for (unsigned int l = 0; l < this->degree-1; ++l)
+                  for (unsigned int q_point = 0; q_point < n_face_points;
+                       ++q_point)
+                    system_matrix (i * (this->degree-1) + j, k * (this->degree-1) + l)
+                    += boundary_weights (q_point + n_edge_points,
+                                         2 * (k * (this->degree-1) + l))
+                       * legendre_polynomials[i].value
+                       (this->generalized_face_support_points[q_point
+                                                              + 4
+                                                              * n_edge_points]
+                        (0))
+                       * lobatto_polynomials[j + 2].value
+                       (this->generalized_face_support_points[q_point
+                                                              + 4
+                                                              * n_edge_points]
+                        (1));
+
+          system_matrix_inv.reinit (system_matrix.m (),
+                                    system_matrix.m ());
+          system_matrix_inv.invert (system_matrix);
+          solution.reinit (system_matrix.m ());
+          system_rhs.reinit (system_matrix.m ());
+
+          const unsigned int
+          face_coordinates[GeometryInfo<3>::faces_per_cell][2]
+          = {{1, 2}, {1, 2}, {2, 0}, {2, 0}, {0, 1}, {0, 1}};
+          const unsigned int
+          edge_indices[GeometryInfo<3>::faces_per_cell][GeometryInfo<3>::lines_per_face]
+          = {{0, 4, 8, 10}, {1, 5, 9, 11}, {8, 9, 2, 6},
+            {10, 11, 3, 7}, {2, 3, 0, 1}, {6, 7, 4, 5}
+          };
+
+          for (unsigned int face = 0;
+               face < GeometryInfo<dim>::faces_per_cell; ++face)
+            {
+              // Set up the right hand side
+              // for the horizontal shape
+              // functions.
+              system_rhs = 0;
+
+              for (unsigned int q_point = 0; q_point < n_face_points;
+                   ++q_point)
+                {
+                  double tmp
+                    = support_point_values[q_point
+                                           + GeometryInfo<dim>::lines_per_cell
+                                           * n_edge_points][face_coordinates[face][0]];
+
+                  for (unsigned int i = 0; i < 2; ++i)
+                    for (unsigned int j = 0; j <= deg; ++j)
+                      tmp -= nodal_values[edge_indices[face][i]
+                                          * this->degree + j]
+                             * this->shape_value_component
+                             (edge_indices[face][i] * this->degree + j,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points],
+                              face_coordinates[face][0]);
+
+                  for (unsigned int i = 0; i <= deg; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      system_rhs (i * deg + j)
+                      += boundary_weights (q_point + n_edge_points,
+                                           2 * (i * deg + j)) * tmp;
+                }
+
+              system_matrix_inv.vmult (solution, system_rhs);
+
+              // Add the computed support_point_values
+              // to the resulting vector
+              // only, if they are not
+              // too small.
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                    nodal_values[(2 * face * this->degree + i
+                                  + GeometryInfo<dim>::lines_per_cell) * deg
+                                 + j + GeometryInfo<dim>::lines_per_cell]
+                      = solution (i * deg + j);
+
+              // Set up the right hand side
+              // for the vertical shape
+              // functions.
+              system_rhs = 0;
+
+              for (unsigned int q_point = 0; q_point < n_face_points;
+                   ++q_point)
+                {
+                  double tmp
+                    = support_point_values[q_point
+                                           + GeometryInfo<dim>::lines_per_cell
+                                           * n_edge_points][face_coordinates[face][1]];
+
+                  for (int i = 2; i < (int) GeometryInfo<dim>::lines_per_face; ++i)
+                    for (unsigned int j = 0; j <= deg; ++j)
+                      tmp -= nodal_values[edge_indices[face][i]
+                                          * this->degree + j]
+                             * this->shape_value_component
+                             (edge_indices[face][i] * this->degree + j,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points],
+                              face_coordinates[face][1]);
+
+                  for (unsigned int i = 0; i <= deg; ++i)
+                    for (unsigned int j = 0; j < deg; ++j)
+                      system_rhs (i * deg + j)
+                      += boundary_weights (q_point + n_edge_points,
+                                           2 * (i * deg + j) + 1)
+                         * tmp;
+                }
+
+              system_matrix_inv.vmult (solution, system_rhs);
+
+              // Add the computed support_point_values
+              // to the resulting vector
+              // only, if they are not
+              // too small.
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  if (std::abs (solution (i * deg + j)) > 1e-14)
+                    nodal_values[((2 * face + 1) * deg + j + GeometryInfo<dim>::lines_per_cell)
+                                 * this->degree + i]
+                      = solution (i * deg + j);
+            }
+
+          // Finally we project
+          // the remaining parts
+          // of the function on
+          // the interior shape
+          // functions.
+          const QGauss<dim> reference_quadrature (this->degree);
+          const unsigned int
+          n_interior_points = reference_quadrature.size ();
+
+          // We create the
+          // system matrix.
+          system_matrix.reinit (this->degree * deg * deg,
+                                this->degree * deg * deg);
+          system_matrix = 0;
+
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              for (unsigned int k = 0; k < deg; ++k)
+                for (unsigned int l = 0; l <= deg; ++l)
+                  for (unsigned int m = 0; m < deg; ++m)
+                    for (unsigned int n = 0; n < deg; ++n)
+                      for (unsigned int q_point = 0;
+                           q_point < n_interior_points; ++q_point)
+                        system_matrix ((i * deg + j) * deg + k,
+                                       (l * deg + m) * deg + n)
+                        += reference_quadrature.weight (q_point)
+                           * legendre_polynomials[i].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (0))
+                           * lobatto_polynomials[j + 2].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (1))
+                           * lobatto_polynomials[k + 2].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (2))
+                           * lobatto_polynomials_grad[l].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (0))
+                           * lobatto_polynomials[m + 2].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (1))
+                           * lobatto_polynomials[n + 2].value
+                           (this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points]
+                            (2));
+
+          system_matrix_inv.reinit (system_matrix.m (),
+                                    system_matrix.m ());
+          system_matrix_inv.invert (system_matrix);
+          // Set up the right hand side.
+          system_rhs.reinit (system_matrix.m ());
+          system_rhs = 0;
+
+          for (unsigned int q_point = 0; q_point < n_interior_points;
+               ++q_point)
+            {
+              double tmp
+                = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
+                                       * n_edge_points
+                                       + GeometryInfo<dim>::faces_per_cell
+                                       * n_face_points][0];
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                {
+                  for (unsigned int j = 0; j < 2; ++j)
+                    for (unsigned int k = 0; k < 2; ++k)
+                      tmp -= nodal_values[i + (j + 4 * k + 2) * this->degree]
+                             * this->shape_value_component
+                             (i + (j + 4 * k + 2) * this->degree,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              0);
+
+                  for (unsigned int j = 0; j < deg; ++j)
+                    for (unsigned int k = 0; k < 4; ++k)
+                      tmp -= nodal_values[(i + 2 * (k + 2) * this->degree
+                                           + GeometryInfo<dim>::lines_per_cell)
+                                          * deg + j
+                                          + GeometryInfo<dim>::lines_per_cell]
+                             * this->shape_value_component
+                             ((i + 2 * (k + 2) * this->degree
+                               + GeometryInfo<dim>::lines_per_cell)
+                              * deg + j
+                              + GeometryInfo<dim>::lines_per_cell,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              0);
+                }
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  for (unsigned int k = 0; k < deg; ++k)
+                    system_rhs ((i * deg + j) * deg + k)
+                    += reference_quadrature.weight (q_point) * tmp
+                       * lobatto_polynomials_grad[i].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (0))
+                       * lobatto_polynomials[j + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (1))
+                       * lobatto_polynomials[k + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (2));
+            }
+
+          solution.reinit (system_rhs.size ());
+          system_matrix_inv.vmult (solution, system_rhs);
+
+          // Add the computed values
+          // to the resulting vector
+          // only, if they are not
+          // too small.
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              for (unsigned int k = 0; k < deg; ++k)
+                if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
+                  nodal_values[((i + 2 * GeometryInfo<dim>::faces_per_cell)
+                                * deg + j + GeometryInfo<dim>::lines_per_cell
+                                + 2 * GeometryInfo<dim>::faces_per_cell)
+                               * deg + k + GeometryInfo<dim>::lines_per_cell]
+                    = solution ((i * deg + j) * deg + k);
+
+          // Set up the right hand side.
+          system_rhs = 0;
+
+          for (unsigned int q_point = 0; q_point < n_interior_points;
+               ++q_point)
+            {
+              double tmp
+                = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
+                                       * n_edge_points
+                                       + GeometryInfo<dim>::faces_per_cell
+                                       * n_face_points][1];
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < 2; ++j)
+                  {
+                    for (unsigned int k = 0; k < 2; ++k)
+                      tmp -= nodal_values[i + (4 * j + k) * this->degree]
+                             * this->shape_value_component
+                             (i + (4 * j + k) * this->degree,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              1);
+
+                    for (unsigned int k = 0; k < deg; ++k)
+                      tmp -= nodal_values[(i + 2 * j * this->degree
+                                           + GeometryInfo<dim>::lines_per_cell)
+                                          * deg + k
+                                          + GeometryInfo<dim>::lines_per_cell]
+                             * this->shape_value_component
+                             ((i + 2 * j * this->degree
+                               + GeometryInfo<dim>::lines_per_cell)
+                              * deg + k
+                              + GeometryInfo<dim>::lines_per_cell,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              1)
+                             + nodal_values[i + ((2 * j + 9) * deg + k
+                                                 + GeometryInfo<dim>::lines_per_cell)
+                                            * this->degree]
+                             * this->shape_value_component
+                             (i + ((2 * j + 9) * deg + k
+                                   + GeometryInfo<dim>::lines_per_cell)
+                              * this->degree,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              1);
+                  }
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  for (unsigned int k = 0; k < deg; ++k)
+                    system_rhs ((i * deg + j) * deg + k)
+                    += reference_quadrature.weight (q_point) * tmp
+                       * lobatto_polynomials_grad[i].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (1))
+                       * lobatto_polynomials[j + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (0))
+                       * lobatto_polynomials[k + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (2));
+            }
+
+          system_matrix_inv.vmult (solution, system_rhs);
+
+          // Add the computed support_point_values
+          // to the resulting vector
+          // only, if they are not
+          // too small.
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              for (unsigned int k = 0; k < deg; ++k)
+                if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
+                  nodal_values[((i + this->degree + 2
+                                 * GeometryInfo<dim>::faces_per_cell) * deg
+                                + j + GeometryInfo<dim>::lines_per_cell + 2
+                                * GeometryInfo<dim>::faces_per_cell) * deg
+                               + k + GeometryInfo<dim>::lines_per_cell]
+                    = solution ((i * deg + j) * deg + k);
+
+          // Set up the right hand side.
+          system_rhs = 0;
+
+          for (unsigned int q_point = 0; q_point < n_interior_points;
+               ++q_point)
+            {
+              double tmp
+                = support_point_values[q_point + GeometryInfo<dim>::lines_per_cell
+                                       * n_edge_points
+                                       + GeometryInfo<dim>::faces_per_cell
+                                       * n_face_points][2];
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < 4; ++j)
+                  {
+                    tmp -= nodal_values[i + (j + 8) * this->degree]
+                           * this->shape_value_component
+                           (i + (j + 8) * this->degree,
+                            this->generalized_support_points[q_point
+                                                             + GeometryInfo<dim>::lines_per_cell
+                                                             * n_edge_points
+                                                             + GeometryInfo<dim>::faces_per_cell
+                                                             * n_face_points],
+                            2);
+
+                    for (unsigned int k = 0; k < deg; ++k)
+                      tmp -= nodal_values[i + ((2 * j + 1) * deg + k
+                                               + GeometryInfo<dim>::lines_per_cell)
+                                          * this->degree]
+                             * this->shape_value_component
+                             (i + ((2 * j + 1) * deg + k
+                                   + GeometryInfo<dim>::lines_per_cell)
+                              * this->degree,
+                              this->generalized_support_points[q_point
+                                                               + GeometryInfo<dim>::lines_per_cell
+                                                               * n_edge_points
+                                                               + GeometryInfo<dim>::faces_per_cell
+                                                               * n_face_points],
+                              2);
+                  }
+
+              for (unsigned int i = 0; i <= deg; ++i)
+                for (unsigned int j = 0; j < deg; ++j)
+                  for (unsigned int k = 0; k < deg; ++k)
+                    system_rhs ((i * deg + j) * deg + k)
+                    += reference_quadrature.weight (q_point) * tmp
+                       * lobatto_polynomials_grad[i].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (2))
+                       * lobatto_polynomials[j + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (0))
+                       * lobatto_polynomials[k + 2].value
+                       (this->generalized_support_points[q_point
+                                                         + GeometryInfo<dim>::lines_per_cell
+                                                         * n_edge_points
+                                                         + GeometryInfo<dim>::faces_per_cell
+                                                         * n_face_points]
+                        (1));
+            }
+
+          system_matrix_inv.vmult (solution, system_rhs);
+
+          // Add the computed support_point_values
+          // to the resulting vector
+          // only, if they are not
+          // too small.
+          for (unsigned int i = 0; i <= deg; ++i)
+            for (unsigned int j = 0; j < deg; ++j)
+              for (unsigned int k = 0; k < deg; ++k)
+                if (std::abs (solution ((i * deg + j) * deg + k)) > 1e-14)
+                  nodal_values[i + ((j + 2 * (deg
+                                              + GeometryInfo<dim>::faces_per_cell))
+                                    * deg + k
+                                    + GeometryInfo<dim>::lines_per_cell)
+                               * this->degree]
+                    = solution ((i * deg + j) * deg + k);
+        }
+
+      break;
+    }
+
+    default:
+      Assert (false, ExcNotImplemented ());
+    }
+}
+
+
+
+
+
 // Since this is a vector valued element,
 // we cannot interpolate a scalar function.
 template <int dim>
diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc
index 164864032e..e7314673b9 100644
--- a/source/fe/fe_q_bubbles.cc
+++ b/source/fe/fe_q_bubbles.cc
@@ -328,6 +328,34 @@ FE_Q_Bubbles<dim,spacedim>::clone() const
 
 
 
+template <int dim, int spacedim>
+void
+FE_Q_Bubbles<dim,spacedim>::
+convert_generalized_support_point_values_to_nodal_values(const std::vector<Vector<double> > &support_point_values,
+                                                         std::vector<double>    &nodal_values) const
+{
+  Assert (support_point_values.size() == this->unit_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(),
+                               this->unit_support_points.size()));
+  Assert (nodal_values.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell));
+  Assert (support_point_values[0].size() == this->n_components(),
+          ExcDimensionMismatch(support_point_values[0].size(),this->n_components()));
+
+  for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+    {
+      const std::pair<unsigned int, unsigned int> index
+        = this->system_to_component_index(i);
+      nodal_values[i] = support_point_values[i](index.first);
+    }
+
+  // We don't use the bubble functions for local interpolation
+  for (unsigned int i = 0; i<n_bubbles; ++i)
+    nodal_values[nodal_values.size()-i-1] = 0.;
+}
+
+
+
 template <int dim, int spacedim>
 void
 FE_Q_Bubbles<dim,spacedim>::interpolate(std::vector<double>       &local_dofs,
diff --git a/source/fe/fe_q_dg0.cc b/source/fe/fe_q_dg0.cc
index d708f6b09d..2866943612 100644
--- a/source/fe/fe_q_dg0.cc
+++ b/source/fe/fe_q_dg0.cc
@@ -174,6 +174,33 @@ FE_Q_DG0<dim,spacedim>::clone() const
 
 
 
+template <int dim, int spacedim>
+void
+FE_Q_DG0<dim,spacedim>::
+convert_generalized_support_point_values_to_nodal_values(const std::vector<Vector<double> > &support_point_values,
+                                                         std::vector<double>    &nodal_dofs) const
+{
+  Assert (support_point_values.size() == this->unit_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(),
+                               this->unit_support_points.size()));
+  Assert (nodal_dofs.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_dofs.size(),this->dofs_per_cell));
+  Assert (support_point_values[0].size() == this->n_components(),
+          ExcDimensionMismatch(support_point_values[0].size(),this->n_components()));
+
+  for (unsigned int i=0; i<this->dofs_per_cell-1; ++i)
+    {
+      const std::pair<unsigned int, unsigned int> index
+        = this->system_to_component_index(i);
+      nodal_dofs[i] = support_point_values[i](index.first);
+    }
+
+  // We don't need the discontinuous function for local interpolation
+  nodal_dofs[nodal_dofs.size()-1] = 0.;
+}
+
+
+
 template <int dim, int spacedim>
 void
 FE_Q_DG0<dim,spacedim>::interpolate(std::vector<double>       &local_dofs,
diff --git a/source/fe/fe_rannacher_turek.cc b/source/fe/fe_rannacher_turek.cc
index e158d8783e..4f7624ac4d 100644
--- a/source/fe/fe_rannacher_turek.cc
+++ b/source/fe/fe_rannacher_turek.cc
@@ -101,6 +101,26 @@ void FE_RannacherTurek<dim>::initialize_support_points()
 
 
 
+template <int dim>
+void
+FE_RannacherTurek<dim>::
+convert_generalized_support_point_values_to_nodal_values(const std::vector<Vector<double> > &support_point_values,
+                                                         std::vector<double> &nodal_dofs) const
+{
+  AssertDimension(support_point_values.size(), this->generalized_support_points.size());
+  AssertDimension(nodal_dofs.size(), this->dofs_per_cell);
+
+  // extract component and call scalar version of this function
+  std::vector<double> scalar_values(support_point_values.size());
+  for (unsigned int q = 0; q < support_point_values.size(); ++q)
+    {
+      scalar_values[q] = support_point_values[q][0];
+    }
+  this->interpolate(nodal_dofs, scalar_values);
+}
+
+
+
 template <int dim>
 void FE_RannacherTurek<dim>::interpolate(
   std::vector<double> &local_dofs,
diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index 8c3c4a5188..772bfb5d01 100644
--- a/source/fe/fe_raviart_thomas.cc
+++ b/source/fe/fe_raviart_thomas.cc
@@ -460,6 +460,42 @@ FE_RaviartThomas<dim>::has_support_on_face (
 }
 
 
+
+template <int dim>
+void
+FE_RaviartThomas<dim>::
+convert_generalized_support_point_values_to_nodal_values(const std::vector<Vector<double> > &support_point_values,
+                                                         std::vector<double>    &nodal_values) const
+{
+  Assert (support_point_values.size() == this->generalized_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size()));
+  Assert (nodal_values.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell));
+  Assert (support_point_values[0].size() == this->n_components(),
+          ExcDimensionMismatch(support_point_values[0].size(),this->n_components()));
+
+  std::fill(nodal_values.begin(), nodal_values.end(), 0.);
+
+  const unsigned int n_face_points = boundary_weights.size(0);
+  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+    for (unsigned int k=0; k<n_face_points; ++k)
+      for (unsigned int i=0; i<boundary_weights.size(1); ++i)
+        {
+          nodal_values[i+face*this->dofs_per_face] += boundary_weights(k,i)
+                                                      * support_point_values[face*n_face_points+k](GeometryInfo<dim>::unit_normal_direction[face]);
+        }
+
+  const unsigned int start_cell_dofs = GeometryInfo<dim>::faces_per_cell*this->dofs_per_face;
+  const unsigned int start_cell_points = GeometryInfo<dim>::faces_per_cell*n_face_points;
+
+  for (unsigned int k=0; k<interior_weights.size(0); ++k)
+    for (unsigned int i=0; i<interior_weights.size(1); ++i)
+      for (unsigned int d=0; d<dim; ++d)
+        nodal_values[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * support_point_values[k+start_cell_points](d);
+}
+
+
+
 // Since this is a vector valued element, we cannot interpolate a
 // scalar function
 template <int dim>
diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc
index f208e6e623..56696f83bd 100644
--- a/source/fe/fe_raviart_thomas_nodal.cc
+++ b/source/fe/fe_raviart_thomas_nodal.cc
@@ -280,6 +280,55 @@ FE_RaviartThomasNodal<dim>::has_support_on_face (
 }
 
 
+
+template <int dim>
+void
+FE_RaviartThomasNodal<dim>::
+convert_generalized_support_point_values_to_nodal_values(const std::vector<Vector<double> > &support_point_values,
+                                                         std::vector<double>    &nodal_values) const
+{
+  Assert (support_point_values.size() == this->generalized_support_points.size(),
+          ExcDimensionMismatch(support_point_values.size(), this->generalized_support_points.size()));
+  Assert (nodal_values.size() == this->dofs_per_cell,
+          ExcDimensionMismatch(nodal_values.size(),this->dofs_per_cell));
+  Assert (support_point_values[0].size() == this->n_components(),
+          ExcDimensionMismatch(support_point_values[0].size(),this->n_components()));
+
+  // First do interpolation on
+  // faces. There, the component
+  // evaluated depends on the face
+  // direction and orientation.
+  unsigned int fbase = 0;
+  unsigned int f=0;
+  for (; f<GeometryInfo<dim>::faces_per_cell;
+       ++f, fbase+=this->dofs_per_face)
+    {
+      for (unsigned int i=0; i<this->dofs_per_face; ++i)
+        {
+          nodal_values[fbase+i] = support_point_values[fbase+i](GeometryInfo<dim>::unit_normal_direction[f]);
+        }
+    }
+
+  // The remaining points form dim
+  // chunks, one for each component.
+  const unsigned int istep = (this->dofs_per_cell - fbase) / dim;
+  Assert ((this->dofs_per_cell - fbase) % dim == 0, ExcInternalError());
+
+  f = 0;
+  while (fbase < this->dofs_per_cell)
+    {
+      for (unsigned int i=0; i<istep; ++i)
+        {
+          nodal_values[fbase+i] = support_point_values[fbase+i](f);
+        }
+      fbase+=istep;
+      ++f;
+    }
+  Assert (fbase == this->dofs_per_cell, ExcInternalError());
+}
+
+
+
 template <int dim>
 void
 FE_RaviartThomasNodal<dim>::interpolate(