]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove the deprecated FiniteElement::interpolate() function. 4204/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 5 Apr 2017 17:15:06 +0000 (11:15 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Thu, 6 Apr 2017 15:24:16 +0000 (09:24 -0600)
Also remove all of its implementations.

21 files changed:
doc/news/changes/incompatibilities/20170405Bangerth [new file with mode: 0644]
include/deal.II/fe/fe.h
include/deal.II/fe/fe_abf.h
include/deal.II/fe/fe_bdm.h
include/deal.II/fe/fe_dg_vector.h
include/deal.II/fe/fe_dg_vector.templates.h
include/deal.II/fe/fe_nedelec.h
include/deal.II/fe/fe_q_bubbles.h
include/deal.II/fe/fe_q_hierarchical.h
include/deal.II/fe/fe_rannacher_turek.h
include/deal.II/fe/fe_raviart_thomas.h
source/fe/fe.cc
source/fe/fe_abf.cc
source/fe/fe_bdm.cc
source/fe/fe_nedelec.cc
source/fe/fe_q_base.cc
source/fe/fe_q_bubbles.cc
source/fe/fe_q_hierarchical.cc
source/fe/fe_rannacher_turek.cc
source/fe/fe_raviart_thomas.cc
source/fe/fe_raviart_thomas_nodal.cc

diff --git a/doc/news/changes/incompatibilities/20170405Bangerth b/doc/news/changes/incompatibilities/20170405Bangerth
new file mode 100644 (file)
index 0000000..2fb1492
--- /dev/null
@@ -0,0 +1,7 @@
+Removed: The FiniteElement::interpolate() function and all of its
+implementations in derived classes has been removed. It was
+previously already deprecated. Use
+FiniteElement::convert_generalized_support_point_values_to_nodal_values()
+instead.
+<br>
+(Wolfgang Bangerth, 2017/04/05)
index 1d7c58cbe08be06cb8198dbb5c47efc00690e1bf..dd0663a6a7d309e3e705542567449ef88a41c928 100644 (file)
@@ -2008,50 +2008,6 @@ public:
   convert_generalized_support_point_values_to_nodal_values (const std::vector<Vector<double> > &support_point_values,
                                                             std::vector<double>                &nodal_values) const;
 
-  /**
-   * Interpolate a set of scalar values, computed in the generalized support
-   * points.
-   *
-   * @note This function is implemented in FiniteElement for the case that the
-   * element has support points. In this case, the resulting coefficients are
-   * just the values in the support points. All other elements must
-   * reimplement it.
-   *
-   * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead.
-   */
-  virtual
-  void
-  interpolate(std::vector<double>       &local_dofs,
-              const std::vector<double> &values) const DEAL_II_DEPRECATED;
-
-  /**
-   * Interpolate a set of vector values, computed in the generalized support
-   * points.
-   *
-   * Since a finite element often only interpolates part of a vector,
-   * <tt>offset</tt> is used to determine the first component of the vector to
-   * be interpolated. Maybe consider changing your data structures to use the
-   * next function.
-   *
-   * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead.
-   */
-  virtual
-  void
-  interpolate(std::vector<double>                &local_dofs,
-              const std::vector<Vector<double> > &values,
-              const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  /**
-   * Interpolate a set of vector values, computed in the generalized support
-   * points.
-   *
-   * @deprecated Use convert_generalized_support_point_values_to_nodal_values() instead.
-   */
-  virtual
-  void
-  interpolate(std::vector<double>                                         &local_dofs,
-              const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   //@}
 
   /**
index 8a6f0a7c160618352040aec2076e58ac46276944..eef630f3886578f3c8d7594397585d16ec576af8 100644 (file)
@@ -129,16 +129,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   virtual std::size_t memory_consumption () const;
   virtual FiniteElement<dim> *clone() const;
 
index b6a1d384dcc40c0cd802e581e2fdb09f2eb410a5..c7cc26e1ea64e57ea388f100384c6a69bf19f279 100644 (file)
@@ -78,16 +78,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
 private:
   /**
    * Only for internal use. Its full name is @p get_dofs_per_object_vector
index 516e14c03c9b29f1b127773bab9e8e70249a2edb..fdafa6a7039841c88c04c916db8747308de1175b 100644 (file)
@@ -80,16 +80,6 @@ public:
   virtual bool has_support_on_face (const unsigned int shape_index,
                                     const unsigned int face_index) const;
 
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<double>          &values) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   virtual std::size_t memory_consumption () const;
 
 private:
index ebaf40e9fbe6cc04104d5c8b5bc590c9963d1679..4a67fc7b8fa1bc01a3f5eb4dd979b60ef53d9d39 100644 (file)
@@ -90,35 +90,6 @@ FE_DGVector<PolynomialType,dim,spacedim>::has_support_on_face
 }
 
 
-template <class PolynomialType, int dim, int spacedim>
-void
-FE_DGVector<PolynomialType,dim,spacedim>::interpolate
-(std::vector<double> &,
- const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-template <class PolynomialType, int dim, int spacedim>
-void
-FE_DGVector<PolynomialType,dim,spacedim>::interpolate
-(std::vector<double> & /*local_dofs*/,
- const std::vector<Vector<double> > & /*values*/,
- const unsigned int /*offset*/) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-template <class PolynomialType, int dim, int spacedim>
-void
-FE_DGVector<PolynomialType,dim,spacedim>::interpolate
-(std::vector<double> & /*local_dofs*/,
- const VectorSlice<const std::vector<std::vector<double> > > & /*values*/) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
 
 template <class PolynomialType, int dim, int spacedim>
 std::size_t
index 9ea468c1221283e16fa9ff40fdc4ca21d1d73559..480b52e43eb6732dbb28150a105ed5af899a6a54 100644 (file)
@@ -261,16 +261,6 @@ public:
   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;
-
-  virtual void interpolate (std::vector<double> &local_dofs,
-                            const std::vector<Vector<double> > &values,
-                            const unsigned int offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate (std::vector<double> &local_dofs,
-                            const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   /**
    * Return a list of constant modes of the element.
    */
index 27ca79a0b17ce8e501bf4027240e2072b508843b..20ac7c0720f17e9e4f35010717e0ab227c5c90af 100644 (file)
@@ -110,16 +110,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double>          &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   /**
    * Return the matrix interpolating from the given finite element to the
    * present one.  The size of the matrix is then @p dofs_per_cell times
index b11550778d18b7d101e98cf43d3249729a038d5b..6b956843e593201c6a26c237614c819dcf20b6aa 100644 (file)
@@ -688,31 +688,6 @@ public:
   virtual std::pair<Table<2,bool>, std::vector<unsigned int> >
   get_constant_modes () const;
 
-  /**
-   * This function is not implemented and throws an exception if called.
-   */
-  virtual
-  void interpolate(std::vector<double>       &local_dofs,
-                   const std::vector<double> &values) const DEAL_II_DEPRECATED;
-
-  /**
-   * This function is not implemented and throws an exception if called.
-   */
-  virtual
-  void
-  interpolate(std::vector<double>                &local_dofs,
-              const std::vector<Vector<double> > &values,
-              const unsigned int offset = 0) const DEAL_II_DEPRECATED;
-
-  /**
-   * This function is not implemented and throws an exception if called.
-   */
-  virtual
-  void
-  interpolate(std::vector<double> &local_dofs,
-              const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
-
 protected:
   /**
    * @p clone function instead of a copy constructor.
index 80c5b7a2ae9e7e1be50deffdb712ed9ce64bc4cd..c7739b9952ae88e807271ec12856499460254863 100644 (file)
@@ -73,15 +73,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
 private:
   /**
    * Order of this element.
index bb7151844ff9907a2fc54a03f6bab5ca589b1461..ebcab3ef3c69dcffda30989ffca26ff09894860f 100644 (file)
@@ -134,16 +134,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
   /**
    * Return a list of constant modes of the element. This method is currently
    * not correctly implemented because it returns ones for all components.
@@ -271,16 +261,6 @@ public:
   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;
-
-  virtual void interpolate(std::vector<double>                &local_dofs,
-                           const std::vector<Vector<double> > &values,
-                           const unsigned int                  offset = 0) const DEAL_II_DEPRECATED;
-
-  virtual void interpolate(std::vector<double> &local_dofs,
-                           const VectorSlice<const std::vector<std::vector<double> > > &values) const DEAL_II_DEPRECATED;
-
 
   virtual void get_face_interpolation_matrix (const FiniteElement<dim> &source,
                                               FullMatrix<double>       &matrix) const;
index 6b2bdd88a5dfca3b27d0cfcb32b0ff00d61d5556..dbcd3c97fff0007576904a4d50d0096ed37728ba 100644 (file)
@@ -1113,75 +1113,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector<Vect
 }
 
 
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::interpolate(
-  std::vector<double>       &local_dofs,
-  const std::vector<double> &values) const
-{
-  Assert (has_support_points(), ExcFEHasNoSupportPoints());
-  Assert (values.size() == unit_support_points.size(),
-          ExcDimensionMismatch(values.size(), unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (this->n_components() == 1,
-          ExcDimensionMismatch(this->n_components(), 1));
-
-  std::copy(values.begin(), values.end(), local_dofs.begin());
-}
-
-
-
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::interpolate(
-  std::vector<double>                &local_dofs,
-  const std::vector<Vector<double> > &values,
-  const unsigned int                  offset) const
-{
-  Assert (has_support_points(), ExcFEHasNoSupportPoints());
-  Assert (values.size() == unit_support_points.size(),
-          ExcDimensionMismatch(values.size(), unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values[0].size() >= offset+this->n_components(),
-          ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
-
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    {
-      const std::pair<unsigned int, unsigned int> index
-        = this->system_to_component_index(i);
-      local_dofs[i] = values[i](offset+index.first);
-    }
-}
-
-
-
-
-template <int dim, int spacedim>
-void
-FiniteElement<dim,spacedim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  Assert (has_support_points(), ExcFEHasNoSupportPoints());
-  Assert (values[0].size() == unit_support_points.size(),
-          ExcDimensionMismatch(values.size(), unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values.size() == this->n_components(),
-          ExcDimensionMismatch(values.size(), this->n_components()));
-
-  for (unsigned int i=0; i<this->dofs_per_cell; ++i)
-    {
-      const std::pair<unsigned int, unsigned int> index
-        = this->system_to_component_index(i);
-      local_dofs[i] = values[index.first][i];
-    }
-}
-
-
 
 template <int dim, int spacedim>
 std::size_t
index 4ad9967d18e0acc6563efcd571782a064fc32580..03fa7b26209cf9cb14d50bb8bc3d2256f08eb009 100644 (file)
@@ -570,142 +570,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector<Vect
 
 
 
-template <int dim>
-void
-FE_ABF<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <int dim>
-void
-FE_ABF<dim>::interpolate(
-  std::vector<double>    &local_dofs,
-  const std::vector<Vector<double> > &values,
-  const unsigned int offset) const
-{
-  Assert (values.size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values[0].size() >= offset+this->n_components(),
-          ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
-
-  std::fill(local_dofs.begin(), local_dofs.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)
-        {
-          local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
-                                                    * values[face*n_face_points+k](GeometryInfo<dim>::unit_normal_direction[face]+offset);
-        }
-
-  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)
-        local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[k+start_cell_points](d+offset);
-
-  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)
-        local_dofs[start_abf_dofs+i] += interior_weights_abf(k,i,d) * values[k+start_cell_points](d+offset);
-
-  // 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)
-            local_dofs[start_abf_dofs+i] += n_orient * boundary_weights_abf(k + fp, i)
-                                            * values[face*n_face_points+fp](GeometryInfo<dim>::unit_normal_direction[face]+offset);
-        }
-    }
-
-  // TODO: Check if this "correction" can be removed.
-  for (unsigned int i=0; i<boundary_weights_abf.size(1); ++i)
-    if (std::fabs (local_dofs[start_abf_dofs+i]) < 1.0e-16)
-      local_dofs[start_abf_dofs+i] = 0.0;
-}
-
-
-
-template <int dim>
-void
-FE_ABF<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  Assert (values.size() == this->n_components(),
-          ExcDimensionMismatch(values.size(), this->n_components()));
-  Assert (values[0].size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values[0].size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-
-  std::fill(local_dofs.begin(), local_dofs.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)
-        {
-          local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
-                                                    * values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+k];
-        }
-
-  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)
-        local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[d][k+start_cell_points];
-
-  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)
-        local_dofs[start_abf_dofs+i] += interior_weights_abf(k,i,d) * values[d][k+start_cell_points];
-
-  // 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)
-            local_dofs[start_abf_dofs+i] += n_orient * boundary_weights_abf(k + fp, i)
-                                            * values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+fp];
-        }
-    }
-
-  // TODO: Check if this "correction" can be removed.
-  for (unsigned int i=0; i<boundary_weights_abf.size(1); ++i)
-    if (std::fabs (local_dofs[start_abf_dofs+i]) < 1.0e-16)
-      local_dofs[start_abf_dofs+i] = 0.0;
-}
-
-
-
 template <int dim>
 std::size_t
 FE_ABF<dim>::memory_consumption () const
index 557c8fc7c3848a85d14cc0fa5dbd25eaa4446a70..3b1909aadf918cc529c9ca97947d253a5bec51fe 100644 (file)
@@ -193,100 +193,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector<Vect
 
 
 
-template <int dim>
-void
-FE_BDM<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-template <int dim>
-void
-FE_BDM<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<Vector<double> > &,
-  unsigned int) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-
-template <int dim>
-void
-FE_BDM<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  AssertDimension (values.size(), dim);
-  Assert (values[0].size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.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)
-            local_dofs[dbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][pbase+i];
-          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 += values[GeometryInfo<dim>::unit_normal_direction[f]][pbase+k] * test_values_face[k][i];
-              local_dofs[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 += values[d][pbase+k] * test_values_cell[k][i];
-          local_dofs[dbase+i] = s;
-        }
-    }
-
-  Assert (dbase == this->dofs_per_cell, ExcInternalError());
-}
-
-
-
-
 template <int dim>
 std::vector<unsigned int>
 FE_BDM<dim>::get_dpo_vector (const unsigned int deg)
index 39fef899742a0c8d707786ff383060f10cfbcc38..e4aedd7233ebf685a6821b237af68436717b8892 100644 (file)
@@ -4012,2505 +4012,6 @@ convert_generalized_support_point_values_to_nodal_values (const std::vector<Vect
 
 
 
-
-
-// Since this is a vector valued element,
-// we cannot interpolate a scalar function.
-template <int dim>
-void FE_Nedelec<dim>::interpolate (std::vector<double> &, const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented ());
-}
-
-
-// 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>::interpolate (std::vector<double> &local_dofs,
-                              const std::vector<Vector<double> > &values,
-                              const unsigned int offset) const
-{
-  const unsigned int deg = this->degree-1;
-
-  Assert (values.size () == this->generalized_support_points.size (),
-          ExcDimensionMismatch (values.size (),
-                                this->generalized_support_points.size ()));
-  Assert (local_dofs.size () == this->dofs_per_cell,
-          ExcDimensionMismatch (local_dofs.size (),this->dofs_per_cell));
-  Assert (values[0].size () >= offset + this->n_components (),
-          ExcDimensionMismatch (values[0].size (),
-                                offset + this->n_components ()));
-  std::fill (local_dofs.begin (), local_dofs.end (), 0.);
-
-  if (offset < dim)
-    switch (dim)
-      {
-      case 2:
-      {
-        const QGauss<1> reference_edge_quadrature (this->degree);
-        const unsigned int &n_edge_points
-          = reference_edge_quadrature.size ();
-
-        // Let us begin with the
-        // interpolation part.
-        for (unsigned int i = 0; i < 2; ++i)
-          {
-            for (unsigned int q_point = 0; q_point < n_edge_points;
-                 ++q_point)
-              local_dofs[i * this->degree]
-              += reference_edge_quadrature.weight (q_point)
-                 * values[q_point + i * n_edge_points] (1);
-
-            // Add the computed values
-            // to the resulting vector
-            // only, if they are not
-            // too small.
-            if (std::abs (local_dofs[i * this->degree]) < 1e-14)
-              local_dofs[i * this->degree] = 0.0;
-          }
-
-        if (offset == 0)
-          for (unsigned int i = 0; i < 2; ++i)
-            {
-              for (unsigned int q_point = 0; q_point < n_edge_points;
-                   ++q_point)
-                local_dofs[(i + 2) * this->degree]
-                += reference_edge_quadrature.weight (q_point)
-                   * values[q_point + (i + 2) * n_edge_points] (0);
-
-              if (std::abs (local_dofs[(i + 2) * this->degree]) < 1e-14)
-                local_dofs[(i + 2) * 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 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);
-            const unsigned int
-            line_coordinate[GeometryInfo<2>::lines_per_cell]
-              = {1, 1, 0, 0};
-            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.
-            FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
-
-            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);
-
-            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)
-              if ((line < 2) || (offset == 0))
-                {
-                  // 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
-                        = values[line * n_edge_points + q_point]
-                          (line_coordinate[line])
-                          - local_dofs[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 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)
-                      local_dofs[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 std::vector<Polynomials::Polynomial<double> > &
-            legendre_polynomials
-              = Polynomials::Legendre::generate_complete_basis (this->degree-1);
-            const unsigned int &n_interior_points
-              = reference_quadrature.size ();
-
-            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);
-            solution.reinit (system_matrix_inv.m ());
-            system_rhs.reinit (system_matrix.m ());
-
-            if (offset == 0)
-              {
-                // Set up the right hand side
-                // for the horizontal shape
-                // functions.
-                system_rhs = 0;
-
-                for (unsigned int q_point = 0;
-                     q_point < n_interior_points; ++q_point)
-                  {
-                    double tmp
-                      = 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 < this->degree; ++j)
-                        tmp -= local_dofs[(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 < this->degree; ++i)
-                      for (unsigned int j = 0; j < this->degree-1; ++j)
-                        system_rhs (i * (this->degree-1) + 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));
-                  }
-
-                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 < this->degree; ++i)
-                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                    if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                      local_dofs[(i + GeometryInfo<dim>::lines_per_cell)
-                                 * (this->degree-1) + j
-                                 + GeometryInfo<dim>::lines_per_cell]
-                        = solution (i * (this->degree-1) + j);
-              }
-
-            // Set up the right hand side
-            // for the vertical shape
-            // functions.
-            system_rhs = 0;
-
-            for (unsigned int q_point = 0; q_point < n_interior_points;
-                 ++q_point)
-              {
-                double tmp
-                  = 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 < this->degree; ++j)
-                    tmp -= local_dofs[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 < this->degree; ++i)
-                  for (unsigned int j = 0; j < this->degree-1; ++j)
-                    system_rhs (i * (this->degree-1) + 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 values
-            // to the resulting vector
-            // only, if they are not
-            // too small.
-            for (unsigned int i = 0; i < this->degree; ++i)
-              for (unsigned int j = 0; j < this->degree-1; ++j)
-                if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                  local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                  + this->degree-1) * this->degree]
-                    = solution (i * (this->degree-1) + j);
-          }
-
-        break;
-      }
-
-      case 3:
-      {
-        const QGauss<1>
-        reference_edge_quadrature (this->degree);
-        const unsigned int &
-        n_edge_points = reference_edge_quadrature.size ();
-
-        // Let us begin with the
-        // interpolation part.
-        for (unsigned int i = 0; i < 4; ++i)
-          {
-            for (unsigned int q_point = 0; q_point < n_edge_points;
-                 ++q_point)
-              local_dofs[(i + 8) * this->degree]
-              += reference_edge_quadrature.weight (q_point)
-                 * values[q_point + (i + 8) * n_edge_points] (2);
-
-            // Add the computed values
-            // to the resulting vector
-            // only, if they are not
-            // too small.
-            if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
-              local_dofs[(i + 8) * this->degree] = 0.0;
-          }
-
-        if (offset + 1 < dim)
-          {
-            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)
-                    local_dofs[(i + 4 * j) * this->degree]
-                    += reference_edge_quadrature.weight (q_point)
-                       * values[q_point + (i + 4 * j) * n_edge_points]
-                       (1);
-
-                  // Add the computed values
-                  // to the resulting vector
-                  // only, if they are not
-                  // too small.
-                  if (std::abs (local_dofs[(i + 4 * j) * this->degree])
-                      < 1e-14)
-                    local_dofs[(i + 4 * j) * this->degree] = 0.0;
-                }
-
-            if (offset == 0)
-              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)
-                      local_dofs[(i + 4 * j + 2) * this->degree]
-                      += reference_edge_quadrature.weight (q_point)
-                         * values[q_point + (i + 4 * j + 2)
-                                  * n_edge_points] (0);
-
-                    // Add the computed values
-                    // to the resulting vector
-                    // only, if they are not
-                    // too small.
-                    if (std::abs (local_dofs[(i + 4 * j + 2)
-                                             * this->degree]) < 1e-14)
-                      local_dofs[(i + 4 * j + 2) * 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 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);
-            const unsigned int
-            line_coordinate[GeometryInfo<3>::lines_per_cell]
-              = {1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
-            FullMatrix<double> system_matrix (this->degree-1, this->degree-1);
-            FullMatrix<double> system_matrix_inv (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 ();
-
-            Vector<double> system_rhs (system_matrix.m ());
-            Vector<double> solution (system_rhs.size ());
-
-            // 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));
-
-            system_matrix_inv.invert (system_matrix);
-
-            for (unsigned int line = 0;
-                 line < GeometryInfo<dim>::lines_per_cell; ++line)
-              {
-                // Set up the right hand side.
-                system_rhs = 0;
-
-                if ((((line == 0) || (line == 1) || (line == 4) ||
-                      (line == 5)) && (offset + 1 < dim)) ||
-                    (((line == 2) || (line == 3) || (line == 6) ||
-                      (line == 7)) && (offset == 0)) || (line > 7))
-                  {
-                    for (unsigned int q_point = 0; q_point < n_edge_points;
-                         ++q_point)
-                      {
-                        double tmp
-                          = values[line * n_edge_points + q_point]
-                            (line_coordinate[line])
-                            - local_dofs[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)
-                        local_dofs[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.n ());
-            system_matrix_inv.invert (system_matrix);
-            solution.reinit (system_matrix.m ());
-            system_rhs.reinit (system_matrix.m ());
-
-            for (unsigned int face = 0;
-                 face < GeometryInfo<dim>::faces_per_cell; ++face)
-              {
-                switch (face)
-                  {
-                  case 0:
-                  {
-                    if (offset + 1 < dim)
-                      {
-                        // 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
-                              = 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 < this->degree; ++j)
-                                tmp
-                                -= local_dofs[4 * i * this->degree
-                                              + j]
-                                   * this->shape_value_component
-                                   (4 * 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 < this->degree; ++i)
-                              for (unsigned int j = 0; j < this->degree-1; ++j)
-                                system_rhs (i * (this->degree-1) + j)
-                                += boundary_weights
-                                   (q_point + n_edge_points,
-                                    2 * (i * (this->degree-1) + j)) * 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 < this->degree; ++i)
-                          for (unsigned int j = 0; j < this->degree-1; ++j)
-                            if (std::abs (solution (i * (this->degree-1) + j))
-                                > 1e-14)
-                              local_dofs[(i
-                                          + GeometryInfo<dim>::lines_per_cell)
-                                         * (this->degree-1) + j
-                                         + GeometryInfo<dim>::lines_per_cell]
-                                = solution (i * (this->degree-1) + 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
-                          = values[q_point
-                                   + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points] (2);
-
-                        for (unsigned int i = 0; i < 2; ++i)
-                          for (unsigned int j = 0; j < this->degree; ++j)
-                            tmp -= local_dofs[2 * (i + 4)
-                                              * this->degree + j]
-                                   * this->shape_value_component
-                                   (2 * (i + 4) * this->degree + j,
-                                    this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points],
-                                    2);
-
-                        for (unsigned int i = 0; i < this->degree; ++i)
-                          for (unsigned int j = 0; j < this->degree-1; ++j)
-                            system_rhs (i * (this->degree-1) + j)
-                            += boundary_weights
-                               (q_point + n_edge_points,
-                                2 * (i * (this->degree-1) + j) + 1)
-                               * 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 < this->degree; ++i)
-                      for (unsigned int j = 0; j < this->degree-1; ++j)
-                        if (std::abs (solution (i * (this->degree-1) + j)) > 1e-14)
-                          local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                          + this->degree-1)
-                                     * this->degree]
-                            = solution (i * (this->degree-1) + j);
-
-                    break;
-                  }
-
-                  case 1:
-                  {
-                    if (offset + 1 < dim)
-                      {
-                        // 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
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points
-                                       + n_face_points] (1);
-
-                            for (unsigned int i = 0; i < 2; ++i)
-                              for (unsigned int j = 0; j <= deg; ++j)
-                                tmp -= local_dofs[(4 * i + 1)
-                                                  * this->degree + j]
-                                       * this->shape_value_component
-                                       ((4 * i + 1) * this->degree
-                                        + j,
-                                        this->generalized_support_points[q_point
-                                                                         + GeometryInfo<dim>::lines_per_cell
-                                                                         * n_edge_points
-                                                                         + n_face_points],
-                                        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)) * 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 <= deg; ++i)
-                          for (unsigned int j = 0; j < deg; ++j)
-                            if (std::abs (solution (i * deg + j))
-                                > 1e-14)
-                              local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                          + 2 * this->degree) * 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
-                          = values[q_point
-                                   + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points + n_face_points]
-                            (2);
-
-                        for (unsigned int i = 0; i < 2; ++i)
-                          for (unsigned int j = 0; j <= deg; ++j)
-                            tmp -= local_dofs[(2 * (i + 4) + 1)
-                                              * this->degree + j]
-                                   * this->shape_value_component
-                                   ((2 * (i + 4) + 1) * this->degree
-                                    + j,
-                                    this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points
-                                                                     + n_face_points],
-                                    2);
-
-                        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 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)
-                          local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                          + 3 * deg)
-                                     * this->degree]
-                            = solution (i * deg + j);
-
-                    break;
-                  }
-
-                  case 2:
-                  {
-                    if (offset == 0)
-                      {
-                        // 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
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points + 2 * n_face_points]
-                                (2);
-
-                            for (unsigned int i = 0; i < 2; ++i)
-                              for (unsigned int j = 0; j <= deg; ++j)
-                                tmp -= local_dofs[(i + 8) * this->degree
-                                                  + j]
-                                       * this->shape_value_component
-                                       ((i + 8) * this->degree + j,
-                                        this->generalized_support_points[q_point
-                                                                         + GeometryInfo<dim>::lines_per_cell
-                                                                         * n_edge_points
-                                                                         + 2
-                                                                         * n_face_points],
-                                        2);
-
-                            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 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)
-                              local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                          + 4 * this->degree) * 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
-                          = values[q_point
-                                   + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points
-                                   + 2 * n_face_points] (0);
-
-                        for (unsigned int i = 0; i < 2; ++i)
-                          for (unsigned int j = 0; j <= deg; ++j)
-                            tmp -= local_dofs[(4 * i + 2)
-                                              * this->degree + j]
-                                   * this->shape_value_component
-                                   ((4 * i + 2) * this->degree
-                                    + j,
-                                    this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points
-                                                                     + 2
-                                                                     * n_face_points],
-                                    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) + 1) * 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 <= deg; ++i)
-                      for (unsigned int j = 0; j < deg; ++j)
-                        if (std::abs (solution (i * deg + j)) > 1e-14)
-                          local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                          + 5 * deg) * this->degree]
-                            = solution (i * deg + j);
-
-                    break;
-                  }
-
-                  case 3:
-                  {
-                    if (offset == 0)
-                      {
-                        // 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
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points + 3 * n_face_points]
-                                (2);
-
-                            for (unsigned int i = 0; i < 2; ++i)
-                              for (unsigned int j = 0; j <= deg; ++j)
-                                tmp -= local_dofs[(i + 10) * this->degree
-                                                  + j]
-                                       * this->shape_value_component
-                                       ((i + 10) * this->degree + j,
-                                        this->generalized_support_points[q_point
-                                                                         + GeometryInfo<dim>::lines_per_cell
-                                                                         * n_edge_points
-                                                                         + 3
-                                                                         * n_face_points],
-                                        2);
-
-                            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 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)
-                              local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                          + 6 * this->degree) * 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
-                          = values[q_point
-                                   + GeometryInfo<dim>::lines_per_cell
-                                   * n_edge_points + 3
-                                   * n_face_points] (0);
-
-                        for (unsigned int i = 0; i < 2; ++i)
-                          for (unsigned int j = 0; j <= deg; ++j)
-                            tmp -= local_dofs[(4 * i + 3)
-                                              * this->degree + j]
-                                   * this->shape_value_component
-                                   ((4 * i + 3) * this->degree
-                                    + j,
-                                    this->generalized_support_points[q_point
-                                                                     + GeometryInfo<dim>::lines_per_cell
-                                                                     * n_edge_points
-                                                                     + 3
-                                                                     * n_face_points],
-                                    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) + 1) * 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 <= deg; ++i)
-                      for (unsigned int j = 0; j < deg; ++j)
-                        if (std::abs (solution (i * deg + j)) > 1e-14)
-                          local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                          + 7 * deg) * this->degree]
-                            = solution (i * deg + j);
-
-                    break;
-                  }
-
-                  case 4:
-                  {
-                    if (offset + 1 < dim)
-                      {
-                        // Set up the right hand side
-                        // for the horizontal shape
-                        // functions.
-                        if (offset == 0)
-                          {
-                            system_rhs = 0;
-
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points + 4
-                                           * n_face_points] (0);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j <= deg; ++j)
-                                    tmp -= local_dofs[(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
-                                                                             + 4
-                                                                             * n_face_points],
-                                            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 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)
-                                  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                              + 8 * this->degree) * 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
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points + 4
-                                       * n_face_points] (1);
-
-                            for (unsigned int i = 0; i < 2; ++i)
-                              for (unsigned int j = 0; j <= deg; ++j)
-                                tmp -= local_dofs[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
-                                                                         + 4
-                                                                         * n_face_points],
-                                        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 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)
-                              local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                              + 9 * deg)
-                                         * this->degree]
-                                = solution (i * deg + j);
-                      }
-
-                    break;
-                  }
-
-                  default:
-                    if (offset + 1 < dim)
-                      {
-                        // Set up the right hand side
-                        // for the horizontal shape
-                        // functions.
-                        if (offset == 0)
-                          {
-                            system_rhs = 0;
-
-                            for (unsigned int q_point = 0;
-                                 q_point < n_face_points; ++q_point)
-                              {
-                                double tmp
-                                  = values[q_point
-                                           + GeometryInfo<dim>::lines_per_cell
-                                           * n_edge_points
-                                           + 5 * n_face_points] (0);
-
-                                for (unsigned int i = 0; i < 2; ++i)
-                                  for (unsigned int j = 0; j <= deg; ++j)
-                                    tmp -= local_dofs[(i + 6)
-                                                      * this->degree + j]
-                                           * this->shape_value_component
-                                           ((i + 6) * this->degree + j,
-                                            this->generalized_support_points[q_point
-                                                                             + GeometryInfo<dim>::lines_per_cell
-                                                                             * n_edge_points
-                                                                             + 5
-                                                                             * n_face_points],
-                                            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 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)
-                                  local_dofs[(i + GeometryInfo<dim>::lines_per_cell
-                                              + 10 * this->degree)
-                                             * 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
-                              = values[q_point
-                                       + GeometryInfo<dim>::lines_per_cell
-                                       * n_edge_points + 5
-                                       * n_face_points] (1);
-
-                            for (unsigned int i = 0; i < 2; ++i)
-                              for (unsigned int j = 0; j <= deg; ++j)
-                                tmp -= local_dofs[(i + 4)
-                                                  * this->degree + j]
-                                       * this->shape_value_component
-                                       ((i + 4) * this->degree + j,
-                                        this->generalized_support_points[q_point
-                                                                         + GeometryInfo<dim>::lines_per_cell
-                                                                         * n_edge_points
-                                                                         + 5
-                                                                         * n_face_points],
-                                        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 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)
-                              local_dofs[i + (j + GeometryInfo<dim>::lines_per_cell
-                                              + 11 * deg) * this->degree]
-                                = 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);
-            system_rhs.reinit (system_matrix_inv.m ());
-            solution.reinit (system_matrix.m ());
-
-            if (offset + 1 < dim)
-              {
-                if (offset == 0)
-                  {
-                    // Set up the right hand side.
-                    system_rhs = 0;
-
-                    for (unsigned int q_point = 0;
-                         q_point < n_interior_points; ++q_point)
-                      {
-                        double tmp
-                          = 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 -= local_dofs[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 -= local_dofs[(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));
-                      }
-
-                    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)
-                            local_dofs[((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
-                      = 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 -= local_dofs[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 -= local_dofs[(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)
-                                   + local_dofs[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 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)
-                        local_dofs[((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
-                  = 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 -= local_dofs[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 -= local_dofs[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 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)
-                    local_dofs[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 ());
-      }
-}
-
-
-// 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>::interpolate (std::vector<double> &local_dofs,
-                              const VectorSlice<const std::vector<std::vector<double> > > &values)
-const
-{
-  const unsigned int deg = this->degree-1;
-  Assert (values.size () == this->n_components (),
-          ExcDimensionMismatch (values.size (), this->n_components ()));
-  Assert (values[0].size () == this->generalized_support_points.size (),
-          ExcDimensionMismatch (values[0].size (),
-                                this->generalized_support_points.size ()));
-  Assert (local_dofs.size () == this->dofs_per_cell,
-          ExcDimensionMismatch (local_dofs.size (), this->dofs_per_cell));
-  std::fill (local_dofs.begin (), local_dofs.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)
-              local_dofs[(i + 2 * j) * this->degree]
-              += reference_edge_quadrature.weight (q_point)
-                 * values[1 - j][q_point + (i + 2 * j) * n_edge_points];
-
-            // Add the computed values
-            // to the resulting vector
-            // only, if they are not
-            // too small.
-            if (std::abs (local_dofs[(i + 2 * j) * this->degree]) < 1e-14)
-              local_dofs[(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 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
-                    = values[line_coordinate[line]][line * n_edge_points
-                                                    + q_point]
-                      - local_dofs[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 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)
-                  local_dofs[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
-                = values[0][q_point + GeometryInfo<dim>::lines_per_cell
-                            * n_edge_points];
-
-              for (unsigned int i = 0; i < 2; ++i)
-                for (unsigned int j = 0; j <= deg; ++j)
-                  tmp -= local_dofs[(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 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)
-                local_dofs[(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
-                = values[1][q_point + GeometryInfo<dim>::lines_per_cell
-                            * n_edge_points];
-
-              for (unsigned int i = 0; i < 2; ++i)
-                for (unsigned int j = 0; j <= deg; ++j)
-                  tmp -= local_dofs[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 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)
-                local_dofs[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)
-            local_dofs[(i + 8) * this->degree]
-            += reference_edge_quadrature.weight (q_point)
-               * values[2][q_point + (i + 8) * n_edge_points];
-
-          for (unsigned int i = 0; i < 2; ++i)
-            for (unsigned int j = 0; j < 2; ++j)
-              for (unsigned int k = 0; k < 2; ++k)
-                local_dofs[(i + 2 * (2 * j + k)) * this->degree]
-                += reference_edge_quadrature.weight (q_point)
-                   * values[1 - k][q_point + (i + 2 * (2 * j + k))
-                                   * n_edge_points];
-        }
-
-      // Add the computed values
-      // to the resulting vector
-      // only, if they are not
-      // too small.
-      for (unsigned int i = 0; i < 4; ++i)
-        if (std::abs (local_dofs[(i + 8) * this->degree]) < 1e-14)
-          local_dofs[(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 (local_dofs[(i + 2 * (2 * j + k)) * this->degree])
-                < 1e-14)
-              local_dofs[(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 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
-                    = values[line_coordinate[line]][line * this->degree
-                                                    + q_point]
-                      - local_dofs[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)
-                  local_dofs[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
-                    = values[face_coordinates[face][0]][q_point
-                                                        + GeometryInfo<dim>::lines_per_cell
-                                                        * n_edge_points];
-
-                  for (unsigned int i = 0; i < 2; ++i)
-                    for (unsigned int j = 0; j <= deg; ++j)
-                      tmp -= local_dofs[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 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)
-                    local_dofs[(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
-                    = values[face_coordinates[face][1]][q_point
-                                                        + GeometryInfo<dim>::lines_per_cell
-                                                        * n_edge_points];
-
-                  for (int i = 2; i < (int) GeometryInfo<dim>::lines_per_face; ++i)
-                    for (unsigned int j = 0; j <= deg; ++j)
-                      tmp -= local_dofs[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 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)
-                    local_dofs[((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
-                = values[0][q_point + GeometryInfo<dim>::lines_per_cell
-                            * n_edge_points
-                            + GeometryInfo<dim>::faces_per_cell
-                            * n_face_points];
-
-              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 -= local_dofs[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 -= local_dofs[(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)
-                  local_dofs[((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
-                = values[1][q_point + GeometryInfo<dim>::lines_per_cell
-                            * n_edge_points
-                            + GeometryInfo<dim>::faces_per_cell
-                            * n_face_points];
-
-              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 -= local_dofs[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 -= local_dofs[(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)
-                             + local_dofs[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 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)
-                  local_dofs[((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
-                = values[2][q_point + GeometryInfo<dim>::lines_per_cell
-                            * n_edge_points
-                            + GeometryInfo<dim>::faces_per_cell
-                            * n_face_points];
-
-              for (unsigned int i = 0; i <= deg; ++i)
-                for (unsigned int j = 0; j < 4; ++j)
-                  {
-                    tmp -= local_dofs[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 -= local_dofs[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 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)
-                  local_dofs[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 ());
-    }
-}
-
-
-
 template <int dim>
 std::pair<Table<2,bool>, std::vector<unsigned int> >
 FE_Nedelec<dim>::get_constant_modes() const
index cb0b1d848bcee19704f8042aaf91e115d1f88f32..097a6e5c232aadbee21fbb04572ccb9fea822383 100644 (file)
@@ -1601,6 +1601,7 @@ FE_Q_Base<PolynomialType,dim,spacedim>::get_constant_modes () const
 }
 
 
+
 // explicit instantiations
 #include "fe_q_base.inst"
 
index da85d736f25fb706aff7d2f3855ddfca64347e21..eccf22f2ff48dcba949f270cb66de983665b3a2d 100644 (file)
@@ -356,84 +356,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector<Vecto
 
 
 
-template <int dim, int spacedim>
-void
-FE_Q_Bubbles<dim,spacedim>::interpolate(std::vector<double>       &local_dofs,
-                                        const std::vector<double> &values) const
-{
-  Assert (values.size() == this->unit_support_points.size(),
-          ExcDimensionMismatch(values.size(),
-                               this->unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (this->n_components() == 1,
-          ExcDimensionMismatch(this->n_components(), 1));
-
-  std::copy(values.begin(), values.end(), local_dofs.begin());
-
-  // We don't use the bubble functions for local interpolation
-  for (unsigned int i = 0; i<n_bubbles; ++i)
-    local_dofs[local_dofs.size()-i-1] = 0.;
-}
-
-
-
-template <int dim, int spacedim>
-void
-FE_Q_Bubbles<dim,spacedim>::interpolate(std::vector<double>    &local_dofs,
-                                        const std::vector<Vector<double> > &values,
-                                        const unsigned int offset) const
-{
-  Assert (values.size() == this->unit_support_points.size(),
-          ExcDimensionMismatch(values.size(),
-                               this->unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values[0].size() >= offset+this->n_components(),
-          ExcDimensionMismatch(values[0].size(),offset+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);
-      local_dofs[i] = values[i](offset+index.first);
-    }
-
-  // We don't use the bubble functions for local interpolation
-  for (unsigned int i = 0; i<n_bubbles; ++i)
-    local_dofs[local_dofs.size()-i-1] = 0.;
-}
-
-
-
-template <int dim, int spacedim>
-void
-FE_Q_Bubbles<dim,spacedim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  Assert (values[0].size() == this->unit_support_points.size(),
-          ExcDimensionMismatch(values.size(),
-                               this->unit_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values.size() == this->n_components(),
-          ExcDimensionMismatch(values.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);
-      local_dofs[i] = values[index.first][i];
-    }
-
-  // We don't use the bubble functions for local interpolation
-  for (unsigned int i = 0; i<n_bubbles; ++i)
-    local_dofs[local_dofs.size()-i-1] = 0.;
-}
-
-
-
 template <int dim, int spacedim>
 void
 FE_Q_Bubbles<dim,spacedim>::
index 519b6009946620f95859fa86309ab3fe5123cddc..0e4e9b90887d330a6c12fc811c95f3b55423ec6b 100644 (file)
@@ -123,39 +123,6 @@ FE_Q_Hierarchical<dim>::get_name () const
 }
 
 
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<double> &) const
-{
-  // The default implementation assumes that the FE has a delta property,
-  // i.e., values at the support points equal the corresponding DoFs. This
-  // is obviously not the case here.
-  Assert (false, ExcNotImplemented());
-}
-
-
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<Vector<double> > &,
-  unsigned int) const
-{
-  Assert (false, ExcNotImplemented());
-}
-
-
-template <int dim>
-void
-FE_Q_Hierarchical<dim>::interpolate(
-  std::vector<double> &,
-  const VectorSlice<const std::vector<std::vector<double> > > &) const
-{
-  Assert (false, ExcNotImplemented());
-}
-
 
 template <int dim>
 FiniteElement<dim> *
@@ -164,6 +131,8 @@ FE_Q_Hierarchical<dim>::clone() const
   return new FE_Q_Hierarchical<dim>(*this);
 }
 
+
+
 template <int dim>
 void
 FE_Q_Hierarchical<dim>::get_interpolation_matrix(const FiniteElement< dim> &source,
index 4f7624ac4dac705ddf5697203ed13246c808e779..8b1b8e7bfcc446f22de1f7e52d82f45daa4a1026 100644 (file)
@@ -105,34 +105,15 @@ 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
+                                                         std::vector<double> &nodal_values) 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,
-  const std::vector<double> &values) const
-{
-  AssertDimension(values.size(), this->generalized_support_points.size());
-  AssertDimension(local_dofs.size(), this->dofs_per_cell);
+  AssertDimension(nodal_values.size(), this->dofs_per_cell);
 
   const unsigned int q_points_per_face = this->weights.size();
-  std::fill(local_dofs.begin(), local_dofs.end(), 0.0);
+  std::fill(nodal_values.begin(), nodal_values.end(), 0.0);
 
-  std::vector<double>::const_iterator value = values.begin();
+  std::vector<Vector<double> >::const_iterator value = support_point_values.begin();
   for (unsigned int face = 0;
        face < dealii::GeometryInfo<dim>::faces_per_cell;
        ++face)
@@ -141,7 +122,7 @@ void FE_RannacherTurek<dim>::interpolate(
            q < q_points_per_face;
            ++q)
         {
-          local_dofs[face] += (*value) * this->weights[q];
+          nodal_values[face] += (*value)[0] * this->weights[q];
           ++value;
         }
     }
@@ -149,46 +130,6 @@ void FE_RannacherTurek<dim>::interpolate(
 
 
 
-template <int dim>
-void FE_RannacherTurek<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const std::vector<Vector<double> > &values,
-  const unsigned int offset) const
-{
-  AssertDimension(values.size(), this->generalized_support_points.size());
-  AssertDimension(local_dofs.size(), this->dofs_per_cell);
-
-  // extract component at offset and call scalar version of this function
-  std::vector<double> scalar_values(values.size());
-  for (unsigned int q = 0; q < values.size(); ++q)
-    {
-      scalar_values[q] = values[q][offset];
-    }
-  this->interpolate(local_dofs, scalar_values);
-}
-
-
-
-template <int dim>
-void FE_RannacherTurek<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  AssertDimension(values.size(), 1);
-  AssertDimension(values[0].size(), this->generalized_support_points.size());
-  AssertDimension(local_dofs.size(), this->dofs_per_cell);
-
-  // convert data structure to use scalar version of this function
-  std::vector<double> scalar_values(values[0].size());
-  for (unsigned int q = 0; q < values[0].size(); ++q)
-    {
-      scalar_values[q] = values[0][q];
-    }
-  this->interpolate(local_dofs, scalar_values);
-}
-
-
-
 // explicit instantiations
 #include "fe_rannacher_turek.inst"
 
index 772bfb5d01fb359462df0e988a000c66471ce140..7964b5a6ada9e54a80c99e43d046af53376c8e2b 100644 (file)
@@ -496,88 +496,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector<Vecto
 
 
 
-// Since this is a vector valued element, we cannot interpolate a
-// scalar function
-template <int dim>
-void
-FE_RaviartThomas<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::interpolate(
-  std::vector<double>    &local_dofs,
-  const std::vector<Vector<double> > &values,
-  const unsigned int offset) const
-{
-  Assert (values.size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values[0].size() >= offset+this->n_components(),
-          ExcDimensionMismatch(values[0].size(),offset+this->n_components()));
-
-  std::fill(local_dofs.begin(), local_dofs.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)
-        {
-          local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
-                                                    * values[face*n_face_points+k](GeometryInfo<dim>::unit_normal_direction[face]+offset);
-        }
-
-  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)
-        local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[k+start_cell_points](d+offset);
-}
-
-
-template <int dim>
-void
-FE_RaviartThomas<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  Assert (values.size() == this->n_components(),
-          ExcDimensionMismatch(values.size(), this->n_components()));
-  Assert (values[0].size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values[0].size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-
-  std::fill(local_dofs.begin(), local_dofs.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)
-        {
-          local_dofs[i+face*this->dofs_per_face] += boundary_weights(k,i)
-                                                    * values[GeometryInfo<dim>::unit_normal_direction[face]][face*n_face_points+k];
-        }
-
-  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)
-        local_dofs[start_cell_dofs+i*dim+d] += interior_weights(k,i,d) * values[d][k+start_cell_points];
-}
-
-
-
 template <int dim>
 std::size_t
 FE_RaviartThomas<dim>::memory_consumption () const
index 56696f83bd574b20e3859bfc87b4432d90339afa..b5d7f269e8645aac57b30f64b482ac5bc5536a34 100644 (file)
@@ -329,109 +329,6 @@ convert_generalized_support_point_values_to_nodal_values(const std::vector<Vecto
 
 
 
-template <int dim>
-void
-FE_RaviartThomasNodal<dim>::interpolate(
-  std::vector<double> &,
-  const std::vector<double> &) const
-{
-  Assert(false, ExcNotImplemented());
-}
-
-
-template <int dim>
-void
-FE_RaviartThomasNodal<dim>::interpolate(
-  std::vector<double>    &local_dofs,
-  const std::vector<Vector<double> > &values,
-  const unsigned int offset) const
-{
-  Assert (values.size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  Assert (values[0].size() >= offset+this->n_components(),
-          ExcDimensionMismatch(values[0].size(),offset+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)
-        {
-          local_dofs[fbase+i] = values[fbase+i](offset+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)
-        {
-          local_dofs[fbase+i] = values[fbase+i](offset+f);
-        }
-      fbase+=istep;
-      ++f;
-    }
-  Assert (fbase == this->dofs_per_cell, ExcInternalError());
-}
-
-
-template <int dim>
-void
-FE_RaviartThomasNodal<dim>::interpolate(
-  std::vector<double> &local_dofs,
-  const VectorSlice<const std::vector<std::vector<double> > > &values) const
-{
-  Assert (values.size() == this->n_components(),
-          ExcDimensionMismatch(values.size(), this->n_components()));
-  Assert (values[0].size() == this->generalized_support_points.size(),
-          ExcDimensionMismatch(values.size(), this->generalized_support_points.size()));
-  Assert (local_dofs.size() == this->dofs_per_cell,
-          ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell));
-  // 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)
-        {
-          local_dofs[fbase+i] = values[GeometryInfo<dim>::unit_normal_direction[f]][fbase+i];
-        }
-    }
-  // 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)
-        {
-          local_dofs[fbase+i] = values[f][fbase+i];
-        }
-      fbase+=istep;
-      ++f;
-    }
-  Assert (fbase == this->dofs_per_cell, ExcInternalError());
-}
-
-
 //TODO: There are tests that check that the following few functions don't produce assertion failures, but none that actually check whether they do the right thing. one example for such a test would be to project a function onto an hp space and make sure that the convergence order is correct with regard to the lowest used polynomial degree
 
 template <int dim>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.