Prevent FE_Bernstein from implementing nonsensical interpolation

diff --git a/include/deal.II/fe/fe_bernstein.h b/include/deal.II/fe/fe_bernstein.h
index f55fd70789c5cc8a6c0d94dc1c2513103343a336..cb9052473325db0ad134e9206a63691ab464d272 100644 (file)
--- a/include/deal.II/fe/fe_bernstein.h
+++ b/include/deal.II/fe/fe_bernstein.h
@@ -70,6 +70,39 @@ public:
    */
   FE_Bernstein (const unsigned int p);
 
+  /**
+   * FE_Bernstein is not interpolatory in the element interior, which prevents
+   * this element from defining an interpolation matrix. An exception will be
+   * thrown.
+   *
+   * Overrides the implementation from FE_Q_Base.
+   */
+  virtual void
+  get_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+                            FullMatrix<double>       &matrix) const;
+
+  /**
+   * FE_Bernstein is not interpolatory in the element interior, which prevents
+   * this element from defining a restriction matrix. An exception will be
+   * thrown.
+   *
+   * Overrides the implementation from FE_Q_Base.
+   */
+  virtual const FullMatrix<double> &
+  get_restriction_matrix (const unsigned int child,
+                          const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
+  /**
+   * FE_Bernstein is not interpolatory in the element interior, which prevents
+   * this element from defining a prolongation matrix. An exception will be
+   * thrown.
+   *
+   * Overrides the implementation from FE_Q_Base.
+   */
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix (const unsigned int child,
+                           const RefinementCase<dim> &refinement_case=RefinementCase<dim>::isotropic_refinement) const;
+
   /**
    * Return the matrix interpolating from a face of one element to the face of
    * the neighboring element.  The size of the matrix is then

diff --git a/source/fe/fe_bernstein.cc b/source/fe/fe_bernstein.cc
index 1e828e207a9769218217ecf16b2bbb58dbe67c28..6fd68ab662d2d2cf1512b7a8e78b6ad072faa985 100644 (file)
--- a/source/fe/fe_bernstein.cc
+++ b/source/fe/fe_bernstein.cc
@@ -43,6 +43,51 @@ FE_Bernstein<dim,spacedim>::FE_Bernstein (const unsigned int degree)
 {}
 
 
+
+template <int dim, int spacedim>
+void
+FE_Bernstein<dim,spacedim>::
+get_interpolation_matrix (const FiniteElement<dim,spacedim> &,
+                          FullMatrix<double> &) const
+{
+  // no interpolation possible. throw exception, as documentation says
+  typedef FiniteElement<dim,spacedim> FEE;
+  AssertThrow (false,
+               typename FEE::ExcInterpolationNotImplemented());
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_Bernstein<dim,spacedim>::get_restriction_matrix (const unsigned int,
+                                                    const RefinementCase<dim> &) const
+{
+  typedef FiniteElement<dim,spacedim> FEE;
+  AssertThrow (false,
+               typename FEE::ExcProjectionVoid());
+  // return dummy, nothing will happen because the base class FE_Q_Base
+  // implements lazy evaluation of those matrices
+  return this->restriction[0][0];
+}
+
+
+
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_Bernstein<dim,spacedim>::get_prolongation_matrix (const unsigned int,
+                                                     const RefinementCase<dim> &) const
+{
+  typedef FiniteElement<dim,spacedim> FEE;
+  AssertThrow (false,
+               typename FEE::ExcEmbeddingVoid());
+  // return dummy, nothing will happen because the base class FE_Q_Base
+  // implements lazy evaluation of those matrices
+  return this->prolongation[0][0];
+}
+
+
+
 template <int dim, int spacedim>
 void
 FE_Bernstein<dim,spacedim>::