From dc07aedaa737d90d8ce625cfc4fb2cc4e0280e4a Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 28 Feb 2021 07:13:53 +0100
Subject: [PATCH] Implement FE_SimplexPoly::get_restriction_matrix()

---
 include/deal.II/fe/fe_simplex_p.h       | 11 +++++++
 include/deal.II/fe/fe_tools.templates.h | 25 ++++++++++-----
 source/fe/fe_simplex_p.cc               | 42 +++++++++++++++++++++++++
 3 files changed, 70 insertions(+), 8 deletions(-)

diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h
index 43fb6b08a3..2d240489c6 100644
--- a/include/deal.II/fe/fe_simplex_p.h
+++ b/include/deal.II/fe/fe_simplex_p.h
@@ -60,6 +60,17 @@ public:
     const RefinementCase<dim> &refinement_case =
       RefinementCase<dim>::isotropic_refinement) const override;
 
+  /**
+   * @copydoc dealii::FiniteElement::get_restriction_matrix()
+   *
+   * @note Only implemented for RefinementCase::isotropic_refinement.
+   */
+  virtual const FullMatrix<double> &
+  get_restriction_matrix(
+    const unsigned int         child,
+    const RefinementCase<dim> &refinement_case =
+      RefinementCase<dim>::isotropic_refinement) const override;
+
   /**
    * @copydoc dealii::FiniteElement::get_face_interpolation_matrix()
    */
diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h
index f96914254c..48512f06fc 100644
--- a/include/deal.II/fe/fe_tools.templates.h
+++ b/include/deal.II/fe/fe_tools.templates.h
@@ -2164,9 +2164,15 @@ namespace FETools
     const unsigned int nd     = fe.n_components();
     const unsigned int degree = fe.degree;
 
+    const ReferenceCell reference_cell = fe.reference_cell();
+
+    const auto &mapping =
+      reference_cell.template get_default_linear_mapping<dim, spacedim>();
+    const auto &q_fine =
+      reference_cell.get_gauss_type_quadrature<dim>(degree + 1);
+
     // prepare FEValues, quadrature etc on
     // coarse cell
-    QGauss<dim>        q_fine(degree + 1);
     const unsigned int nq = q_fine.size();
 
     // create mass matrix on coarse cell.
@@ -2174,9 +2180,10 @@ namespace FETools
     {
       // set up a triangulation for coarse cell
       Triangulation<dim, spacedim> tr;
-      GridGenerator::hyper_cube(tr, 0, 1);
+      GridGenerator::reference_cell(reference_cell, tr);
 
-      FEValues<dim, spacedim> coarse(fe,
+      FEValues<dim, spacedim> coarse(mapping,
+                                     fe,
                                      q_fine,
                                      update_JxW_values | update_values);
 
@@ -2212,9 +2219,10 @@ namespace FETools
 
 
     const auto compute_one_case =
-      [&fe, &q_fine, n, nd, nq](const unsigned int        ref_case,
-                                const FullMatrix<double> &inverse_mass_matrix,
-                                std::vector<FullMatrix<double>> &matrices) {
+      [&reference_cell, &mapping, &fe, &q_fine, n, nd, nq](
+        const unsigned int               ref_case,
+        const FullMatrix<double> &       inverse_mass_matrix,
+        std::vector<FullMatrix<double>> &matrices) {
         const unsigned int nc =
           GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
 
@@ -2228,11 +2236,12 @@ namespace FETools
 
         // create a respective refinement on the triangulation
         Triangulation<dim, spacedim> tr;
-        GridGenerator::hyper_cube(tr, 0, 1);
+        GridGenerator::reference_cell(reference_cell, tr);
         tr.begin_active()->set_refine_flag(RefinementCase<dim>(ref_case));
         tr.execute_coarsening_and_refinement();
 
-        FEValues<dim, spacedim> fine(fe,
+        FEValues<dim, spacedim> fine(mapping,
+                                     fe,
                                      q_fine,
                                      update_quadrature_points |
                                        update_JxW_values | update_values);
diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc
index 9f32a2f293..fe9b39cab5 100644
--- a/source/fe/fe_simplex_p.cc
+++ b/source/fe/fe_simplex_p.cc
@@ -354,6 +354,48 @@ FE_SimplexPoly<dim, spacedim>::get_prolongation_matrix(
 
 
 
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_SimplexPoly<dim, spacedim>::get_restriction_matrix(
+  const unsigned int         child,
+  const RefinementCase<dim> &refinement_case) const
+{
+  Assert(refinement_case == RefinementCase<dim>::isotropic_refinement,
+         ExcNotImplemented());
+  AssertDimension(dim, spacedim);
+
+  // initialization upon first request
+  if (this->restriction[refinement_case - 1][child].n() == 0)
+    {
+      std::lock_guard<std::mutex> lock(this->mutex);
+
+      // if matrix got updated while waiting for the lock
+      if (this->restriction[refinement_case - 1][child].n() ==
+          this->n_dofs_per_cell())
+        return this->restriction[refinement_case - 1][child];
+
+      // now do the work. need to get a non-const version of data in order to
+      // be able to modify them inside a const function
+      auto &this_nonconst = const_cast<FE_SimplexPoly<dim, spacedim> &>(*this);
+
+      std::vector<std::vector<FullMatrix<double>>> isotropic_matrices(
+        RefinementCase<dim>::isotropic_refinement);
+      isotropic_matrices.back().resize(
+        GeometryInfo<dim>::n_children(RefinementCase<dim>(refinement_case)),
+        FullMatrix<double>(this->n_dofs_per_cell(), this->n_dofs_per_cell()));
+
+      FETools::compute_projection_matrices(*this, isotropic_matrices, true);
+
+      this_nonconst.restriction[refinement_case - 1].swap(
+        isotropic_matrices.back());
+    }
+
+  // finally return the matrix
+  return this->restriction[refinement_case - 1][child];
+}
+
+
+
 template <int dim, int spacedim>
 void
 FE_SimplexPoly<dim, spacedim>::get_face_interpolation_matrix(
-- 
2.39.5