From 25599fba0e0667c1a261b1926ae12ccb1bc8f520 Mon Sep 17 00:00:00 2001
From: kinnewig <Sebastian@Kinnewig.org>
Date: Fri, 12 Apr 2024 17:54:56 +0200
Subject: [PATCH] Add get_prolongation_matrix() to FE_NedelecSZ.

---
 include/deal.II/fe/fe_nedelec_sz.h | 33 +++++++++++++++++++++
 source/fe/fe_nedelec_sz.cc         | 47 ++++++++++++++++++++++++++++++
 2 files changed, 80 insertions(+)

diff --git a/include/deal.II/fe/fe_nedelec_sz.h b/include/deal.II/fe/fe_nedelec_sz.h
index ee2852cece..6fa29c2f49 100644
--- a/include/deal.II/fe/fe_nedelec_sz.h
+++ b/include/deal.II/fe/fe_nedelec_sz.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/derivative_form.h>
 #include <deal.II/base/polynomials_integrated_legendre_sz.h>
 #include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
 
 #include <deal.II/fe/fe.h>
 #include <deal.II/fe/fe_values.h>
@@ -156,6 +157,32 @@ public:
                             const Point<dim>  &p,
                             const unsigned int component) const override;
 
+  /**
+   * Embedding matrix between grids.
+   *
+   * The identity operator from a coarse grid space into a fine grid space is
+   * associated with a matrix @p P. The restriction of this matrix @p P_i to a
+   * single child cell is returned here.
+   *
+   * The matrix @p P is the concatenation, not the sum of the cell matrices @p
+   * P_i. That is, if the same non-zero entry <tt>j,k</tt> exists in two
+   * different child matrices @p P_i, the value should be the same in both
+   * matrices and it is copied into the matrix @p P only once.
+   *
+   * Row and column indices are related to fine grid and coarse grid spaces,
+   * respectively, consistent with the definition of the associated operator.
+   *
+   * These matrices are used by routines assembling the prolongation matrix
+   * for multi-level methods.  Upon assembling the transfer matrix between
+   * cells using this matrix array, zero elements in the prolongation matrix
+   * are discarded and will not fill up the transfer matrix.
+   */
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix(
+    const unsigned int         child,
+    const RefinementCase<dim> &refinement_case =
+      RefinementCase<dim>::isotropic_refinement) const override;
+
 protected:
   /**
    * The mapping kind to be used to map shape functions from the reference
@@ -464,6 +491,12 @@ private:
   fill_face_values(const typename Triangulation<dim, dim>::cell_iterator &cell,
                    const Quadrature<dim> &quadrature,
                    const InternalData    &fedata) const;
+
+  /**
+   * Mutex variables used for protecting the initialization of restriction
+   * and embedding matrices.
+   */
+  mutable Threads::Mutex prolongation_matrix_mutex;
 };
 
 
diff --git a/source/fe/fe_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc
index d891c9631f..ff533e331b 100644
--- a/source/fe/fe_nedelec_sz.cc
+++ b/source/fe/fe_nedelec_sz.cc
@@ -3417,6 +3417,53 @@ FE_NedelecSZ<dim, spacedim>::create_polynomials(const unsigned int degree)
 
 
 
+template <int dim, int spacedim>
+const FullMatrix<double> &
+FE_NedelecSZ<dim, spacedim>::get_prolongation_matrix(
+  const unsigned int         child,
+  const RefinementCase<dim> &refinement_case) const
+{
+  AssertIndexRange(refinement_case,
+                   RefinementCase<dim>::isotropic_refinement + 1);
+  Assert(refinement_case != RefinementCase<dim>::no_refinement,
+         ExcMessage(
+           "Prolongation matrices are only available for refined cells!"));
+  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
+
+  // initialization upon first request
+  if (this->prolongation[refinement_case - 1][child].n() == 0)
+    {
+      std::lock_guard<std::mutex> lock(prolongation_matrix_mutex);
+
+      // if matrix got updated while waiting for the lock
+      if (this->prolongation[refinement_case - 1][child].n() ==
+          this->n_dofs_per_cell())
+        return this->prolongation[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
+      FE_NedelecSZ<dim> &this_nonconst = const_cast<FE_NedelecSZ<dim> &>(*this);
+
+      // Reinit the vectors of
+      // restriction and prolongation
+      // matrices to the right sizes.
+      // Restriction only for isotropic
+      // refinement
+      this_nonconst.reinit_restriction_and_prolongation_matrices();
+      // Fill prolongation matrices with embedding operators
+      FETools::compute_embedding_matrices(
+        this_nonconst,
+        this_nonconst.prolongation,
+        true,
+        1.e-15 * std::exp(std::pow(this->degree, 1.075)));
+    }
+
+  // we use refinement_case-1 here. the -1 takes care of the origin of the
+  // vector, as for RefinementCase<dim>::no_refinement (=0) there is no data
+  // available and so the vector indices are shifted
+  return this->prolongation[refinement_case - 1][child];
+}
+
 // explicit instantiations
 #include "fe_nedelec_sz.inst"
 
-- 
2.39.5