From 382ba57b80899993d29f91ed93f42fdd406528b5 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Sat, 3 Jul 2021 09:51:08 -0400
Subject: [PATCH] Add a more general ctor for FE_SimplexPoly.

This will let us implement FE_SimplexP_Bubbles as an inheriting class in a few
more commits.
---
 include/deal.II/fe/fe_simplex_p.h | 16 +++++++-
 source/fe/fe_simplex_p.cc         | 67 +++++++++++++++++--------------
 2 files changed, 52 insertions(+), 31 deletions(-)

diff --git a/include/deal.II/fe/fe_simplex_p.h b/include/deal.II/fe/fe_simplex_p.h
index fc228114e5..85c05da306 100644
--- a/include/deal.II/fe/fe_simplex_p.h
+++ b/include/deal.II/fe/fe_simplex_p.h
@@ -25,7 +25,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 /**
- * Base class of FE_SimplexP and FE_SimplexDGP.
+ * Base class of FE_SimplexP, FE_SimplexDGP, and FE_SimplexP_Bubbles.
  *
  * @note Only implemented for 2D and 3D.
  *
@@ -42,6 +42,16 @@ public:
                  const std::vector<unsigned int> &                 dpo_vector,
                  const typename FiniteElementData<dim>::Conformity conformity);
 
+  /**
+   * Constructor.
+   */
+  FE_SimplexPoly(
+    const BarycentricPolynomials<dim>              polynomials,
+    const FiniteElementData<dim> &                 fe_data,
+    const std::vector<Point<dim>> &                unit_support_points,
+    const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
+    const FullMatrix<double> &                     interface_constraints);
+
   /**
    * Return a list of constant modes of the element. For this element, the
    * list consists of true arguments for all components.
@@ -103,6 +113,10 @@ public:
     const std::vector<Vector<double>> &support_point_values,
     std::vector<double> &              nodal_values) const override;
 
+protected:
+  /**
+   * Mutex used to guard computation of some internal lookup tables.
+   */
   mutable Threads::Mutex mutex;
 };
 
diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc
index 629623fb54..25b4ff98ce 100644
--- a/source/fe/fe_simplex_p.cc
+++ b/source/fe/fe_simplex_p.cc
@@ -163,8 +163,14 @@ namespace
    */
   template <int dim>
   std::vector<std::vector<Point<dim - 1>>>
-  unit_face_support_points_fe_poly(const unsigned int degree)
+  unit_face_support_points_fe_poly(
+    const unsigned int                          degree,
+    typename FiniteElementData<dim>::Conformity conformity)
   {
+    // Discontinuous elements don't have face support points
+    if (conformity == FiniteElementData<dim>::Conformity::L2)
+      return {};
+
     // this concept doesn't exist in 1D so just return an empty vector
     if (dim == 1)
       return {};
@@ -313,37 +319,38 @@ FE_SimplexPoly<dim, spacedim>::FE_SimplexPoly(
   const unsigned int                                degree,
   const std::vector<unsigned int> &                 dpo_vector,
   const typename FiniteElementData<dim>::Conformity conformity)
+  : FE_SimplexPoly(BarycentricPolynomials<dim>::get_fe_p_basis(degree),
+                   FiniteElementData<dim>(dpo_vector,
+                                          dim == 2 ?
+                                            ReferenceCells::Triangle :
+                                            ReferenceCells::Tetrahedron,
+                                          1,
+                                          degree,
+                                          conformity),
+                   unit_support_points_fe_poly<dim>(degree),
+                   unit_face_support_points_fe_poly<dim>(degree, conformity),
+                   constraints_fe_poly<dim>(degree))
+{}
+
+
+
+template <int dim, int spacedim>
+FE_SimplexPoly<dim, spacedim>::FE_SimplexPoly(
+  const BarycentricPolynomials<dim>              polynomials,
+  const FiniteElementData<dim> &                 fe_data,
+  const std::vector<Point<dim>> &                unit_support_points,
+  const std::vector<std::vector<Point<dim - 1>>> unit_face_support_points,
+  const FullMatrix<double> &                     interface_constraints)
   : dealii::FE_Poly<dim, spacedim>(
-      BarycentricPolynomials<dim>::get_fe_p_basis(degree),
-      FiniteElementData<dim>(dpo_vector,
-                             dim == 2 ? ReferenceCells::Triangle :
-                                        ReferenceCells::Tetrahedron,
-                             1,
-                             degree,
-                             conformity),
-      std::vector<bool>(FiniteElementData<dim>(dpo_vector,
-                                               dim == 2 ?
-                                                 ReferenceCells::Triangle :
-                                                 ReferenceCells::Tetrahedron,
-                                               1,
-                                               degree)
-                          .dofs_per_cell,
-                        true),
-      std::vector<ComponentMask>(
-        FiniteElementData<dim>(dpo_vector,
-                               dim == 2 ? ReferenceCells::Triangle :
-                                          ReferenceCells::Tetrahedron,
-                               1,
-                               degree)
-          .dofs_per_cell,
-        std::vector<bool>(1, true)))
+      polynomials,
+      fe_data,
+      std::vector<bool>(fe_data.dofs_per_cell),
+      std::vector<ComponentMask>(fe_data.dofs_per_cell,
+                                 std::vector<bool>(1, true)))
 {
-  this->unit_support_points = unit_support_points_fe_poly<dim>(degree);
-  // Discontinuous elements don't have face support points
-  if (conformity == FiniteElementData<dim>::Conformity::H1)
-    this->unit_face_support_points =
-      unit_face_support_points_fe_poly<dim>(degree);
-  this->interface_constraints = constraints_fe_poly<dim>(degree);
+  this->unit_support_points      = unit_support_points;
+  this->unit_face_support_points = unit_face_support_points;
+  this->interface_constraints    = interface_constraints;
 }
 
 
-- 
2.39.5