From d6bed8055741f9a00ab245b2a1c72d8286db3d56 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 12 Jul 2016 14:19:01 +0200
Subject: [PATCH] Prevent FE_Bernstein from implementing nonsensical
 interpolation

---
 include/deal.II/fe/fe_bernstein.h | 33 +++++++++++++++++++++++
 source/fe/fe_bernstein.cc         | 45 +++++++++++++++++++++++++++++++
 2 files changed, 78 insertions(+)

diff --git a/include/deal.II/fe/fe_bernstein.h b/include/deal.II/fe/fe_bernstein.h
index f55fd70789..cb90524733 100644
--- 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 1e828e207a..6fd68ab662 100644
--- 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>::
-- 
2.39.5