From: Simon Sticko <simon@sticko.se>
Date: Tue, 1 Jun 2021 07:01:01 +0000 (+0200)
Subject: Implement face interpolation between FE_Bernstein and elements with support points.
X-Git-Tag: v9.4.0-rc1~1293^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2e453c2c7563544a09c923dfeb3724afe3fc2368;p=dealii.git

Implement face interpolation between FE_Bernstein and elements with support points.

By using the implementation of get_face_interpolation_matrix in the
base class (FE_Q_Base) when the incoming element isn't a Bernstein
element.
---

diff --git a/include/deal.II/fe/fe_bernstein.h b/include/deal.II/fe/fe_bernstein.h
index 362d0d7ec4..4d471eef78 100644
--- a/include/deal.II/fe/fe_bernstein.h
+++ b/include/deal.II/fe/fe_bernstein.h
@@ -110,8 +110,8 @@ public:
    * the neighboring element.  The size of the matrix is then
    * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
    * FE_Bernstein element family only provides interpolation matrices for
-   * elements of the same type and FE_Nothing. For all other elements, an
-   * exception of type
+   * elements of the same type, for elements that have support points, and
+   * FE_Nothing. For all other elements, an exception of type
    * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
    */
   virtual void
@@ -124,8 +124,8 @@ public:
    * the neighboring element.  The size of the matrix is then
    * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. The
    * FE_Bernstein element family only provides interpolation matrices for
-   * elements of the same type and FE_Nothing. For all other elements, an
-   * exception of type
+   * elements of the same type, for elements that have support points, and
+   * FE_Nothing. For all other elements, an exception of type
    * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
    */
   virtual void
diff --git a/source/fe/fe_bernstein.cc b/source/fe/fe_bernstein.cc
index cff7465e48..b25f172e44 100644
--- a/source/fe/fe_bernstein.cc
+++ b/source/fe/fe_bernstein.cc
@@ -195,15 +195,13 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
           Assert(std::fabs(sum - 1) < eps, ExcInternalError());
         }
     }
-  else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
+  else
     {
-      // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+      // When the incoming element is not FE_Bernstein we can just delegate to
+      // the base class to create the interpolation matrix.
+      FE_Q_Base<dim, spacedim>::get_subface_interpolation_matrix(
+        x_source_fe, subface, interpolation_matrix, face_no);
     }
-  else
-    AssertThrow(
-      false,
-      (typename FiniteElement<dim,
-                              spacedim>::ExcInterpolationNotImplemented()));
 }
 
 
diff --git a/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.cc b/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.cc
new file mode 100644
index 0000000000..8dc698a658
--- /dev/null
+++ b/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.cc
@@ -0,0 +1,65 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test that we can use FE_Bernstein::get_face_interpolation_matrix to create
+// an interpolation matrix between FE_Bernstein and an element that has support
+// points.
+
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+
+// Set up an FE_Q element and an FE_Bernstein element with the incoming order.
+// Call get_face_interpolation_matrix to create the interpolation matrix
+// between these and write it to deallog.
+template <int dim>
+void
+create_and_print_face_interpolation_matrix(const unsigned int order)
+{
+  deallog << "dim = " << dim << ", order = " << order << std::endl;
+
+  const FE_Q<dim>         fe_q(order);
+  const FE_Bernstein<dim> fe_bernstein(order);
+
+  FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(),
+                                          fe_bernstein.n_dofs_per_face());
+  const unsigned int face_index = 0;
+
+
+  fe_bernstein.get_face_interpolation_matrix(fe_q,
+                                             interpolation_matrix,
+                                             face_index);
+
+  interpolation_matrix.print(deallog);
+  deallog << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  const int dim = 2;
+
+  const std::vector<unsigned int> orders = {1, 2};
+  for (const unsigned int order : orders)
+    create_and_print_face_interpolation_matrix<dim>(order);
+}
diff --git a/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.output b/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.output
new file mode 100644
index 0000000000..9096def154
--- /dev/null
+++ b/tests/fe/face_interpolation_fe_Bernstein_2_fe_q.output
@@ -0,0 +1,10 @@
+
+DEAL::dim = 2, order = 1
+DEAL::1.0  0.0  
+DEAL::0.0  1.0  
+DEAL::
+DEAL::dim = 2, order = 2
+DEAL::1.0  0.0  0.0  
+DEAL::0.0  1.0  0.0  
+DEAL::0.25 0.25 0.50 
+DEAL::