From 5d56330ab09723e677e36f28767b5023f4ee7e25 Mon Sep 17 00:00:00 2001
From: Graham Harper <grahambenharper@gmail.com>
Date: Tue, 3 Jul 2018 21:17:38 +0000
Subject: [PATCH] Add the implementation of the Bernardi-Raugel element.

---
 include/deal.II/fe/fe_bernardi_raugel.h | 103 +++++++++++++
 source/fe/fe_bernardi_raugel.cc         | 183 ++++++++++++++++++++++++
 source/fe/fe_poly_tensor.cc             |   1 +
 source/fe/fe_poly_tensor.inst.in        |   2 +
 4 files changed, 289 insertions(+)
 create mode 100644 include/deal.II/fe/fe_bernardi_raugel.h
 create mode 100644 source/fe/fe_bernardi_raugel.cc

diff --git a/include/deal.II/fe/fe_bernardi_raugel.h b/include/deal.II/fe/fe_bernardi_raugel.h
new file mode 100644
index 0000000000..cce14f340b
--- /dev/null
+++ b/include/deal.II/fe/fe_bernardi_raugel.h
@@ -0,0 +1,103 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_fe_bernardi_raugel_h
+#define dealii_fe_bernardi_raugel_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_bernardi_raugel.h>
+#include <deal.II/base/table.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_poly_tensor.h>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * The Bernardi-Raugel element.
+ *
+ * <h3>Degrees of freedom</h3>
+ * The BR1 element has <i>dim</i> degrees of freedom on each node and 1 on each
+ * face. The shape functions are ordered by the $(Q_1)^d$ shape functions supported on each
+ * vertex, increasing according to vertex ordering on the element in GeometryInfo, then the
+ * bubble functions follow in the ordering given in PolynomialsBernardiRaugel.
+ * 
+ * This element only has 1 degree (degree $p=1$) because it yields an LBB stable pair BR1-P0
+ * for Stokes problems which is lower degree than the Taylor-Hood element. The pair is
+ * sometimes referred to as an enriched P1-P0 element or a reduced P2-P0 element.
+ *
+ * This element does not support hanging nodes or multigrid in the current implementation.
+ *
+ */
+template <int dim>
+class FE_BernardiRaugel
+  : public FE_PolyTensor<PolynomialsBernardiRaugel<dim>, dim>
+{
+public:
+  /**
+   * Constructor for the Bernardi-Raugel element of degree @p p. The only
+   * supported degree is 1.
+   * 
+   * @arg p: The degree of the element $p=1$ for $BR_1$.
+   */
+  FE_BernardiRaugel(const unsigned int p = 1);
+
+  /**
+   * Return a string that uniquely identifies a finite element. This class
+   * returns <tt>FE_BR<dim>(degree)</tt>, with @p dim and @p degree replaced
+   * by appropriate values.
+   */
+  virtual std::string
+  get_name() const;
+
+  virtual std::unique_ptr<FiniteElement<dim, dim>>
+  clone() const;
+
+  // documentation inherited from the base class
+  virtual void
+  convert_generalized_support_point_values_to_dof_values(
+    const std::vector<Vector<double>> &support_point_values,
+    std::vector<double> &              nodal_values) const;
+
+private:
+  /**
+   * Only for internal use. Its full name is @p get_dofs_per_object_vector
+   * function and it creates the @p dofs_per_object vector that is needed
+   * within the constructor to be passed to the constructor of @p
+   * FiniteElementData.
+   */
+  static std::vector<unsigned int>
+  get_dpo_vector();
+
+  /**
+   * Initialize the FiniteElement<dim>::generalized_support_points and
+   * FiniteElement<dim>::generalized_face_support_points fields. Called from
+   * the constructor. See the
+   * @ref GlossGeneralizedSupport "glossary entry on generalized support points"
+   * for more information.
+   */
+  void
+  initialize_support_points();
+};
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/source/fe/fe_bernardi_raugel.cc b/source/fe/fe_bernardi_raugel.cc
new file mode 100644
index 0000000000..bcc27fb2bd
--- /dev/null
+++ b/source/fe/fe_bernardi_raugel.cc
@@ -0,0 +1,183 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/polynomials_p.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/std_cxx14/memory.h>
+#include <deal.II/base/table.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_bernardi_raugel.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping.h>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <iostream>
+#include <sstream>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim>
+FE_BernardiRaugel<dim>::FE_BernardiRaugel(const unsigned int p)
+  : FE_PolyTensor<PolynomialsBernardiRaugel<dim>, dim>(
+      p,
+      FiniteElementData<dim>(get_dpo_vector(),
+                             dim,
+                             2,
+                             FiniteElementData<dim>::Hdiv),
+      std::vector<bool>(PolynomialsBernardiRaugel<dim>::compute_n_pols(p),
+                        true),
+      std::vector<ComponentMask>(PolynomialsBernardiRaugel<dim>::compute_n_pols(
+                                   p),
+                                 std::vector<bool>(dim, true)))
+{
+  Assert(dim == 2 || dim == 3, ExcImpossibleInDim(dim));
+  Assert(p == 1, ExcMessage("Only BR1 elements are available"));
+
+  // const unsigned int n_dofs = this->dofs_per_cell;
+
+  this->mapping_type = mapping_none;
+  // These must be done first, since
+  // they change the evaluation of
+  // basis functions
+
+  // Set up the generalized support
+  // points
+  initialize_support_points();
+}
+
+
+
+template <int dim>
+std::string
+FE_BernardiRaugel<dim>::get_name() const
+{
+  std::ostringstream namebuf;
+  namebuf << "FE_BR<" << dim << ">(" << 1 << ")";
+
+  return namebuf.str();
+}
+
+
+template <int dim>
+std::unique_ptr<FiniteElement<dim, dim>>
+FE_BernardiRaugel<dim>::clone() const
+{
+  return std_cxx14::make_unique<FE_BernardiRaugel<dim>>(*this);
+}
+
+
+
+template <int dim>
+void
+FE_BernardiRaugel<dim>::convert_generalized_support_point_values_to_dof_values(
+  const std::vector<Vector<double>> &support_point_values,
+  std::vector<double> &              nodal_values) const
+{
+  Assert(support_point_values.size() == this->generalized_support_points.size(),
+         ExcDimensionMismatch(support_point_values.size(),
+                              this->generalized_support_points.size()));
+  AssertDimension(support_point_values[0].size(), dim);
+  Assert(nodal_values.size() == this->dofs_per_cell,
+         ExcDimensionMismatch(nodal_values.size(), this->dofs_per_cell));
+
+  static std::vector<Tensor<1, dim>> normals;
+  for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+    {
+      Tensor<1, dim> normal;
+      normal[i / 2] = 1;
+      normals.push_back(normal);
+    }
+
+  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
+    nodal_values[i] = support_point_values[i][i % dim];
+
+  // Compute the support points values for the bubble functions
+  for (unsigned int i = dim * GeometryInfo<dim>::vertices_per_cell;
+       i < dim * GeometryInfo<dim>::vertices_per_cell +
+             GeometryInfo<dim>::faces_per_cell;
+       ++i)
+    {
+      nodal_values[i] = 0;
+      for (unsigned int j = 0; j < dim; ++j)
+        nodal_values[i] +=
+          support_point_values[i][j] *
+          normals[i - dim * GeometryInfo<dim>::vertices_per_cell][j];
+    }
+}
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_BernardiRaugel<dim>::get_dpo_vector()
+{
+  // compute the number of unknowns per cell interior/face/edge
+  //
+  // there are <tt>dim</tt> degrees of freedom per node and there
+  // is 1 degree of freedom per edge in 2D (face in 3D)
+  std::vector<unsigned int> dpo(dim + 1, 0u);
+  dpo[0]       = dim;
+  dpo[dim - 1] = 1u;
+
+  return dpo;
+}
+
+
+template <int dim>
+void
+FE_BernardiRaugel<dim>::initialize_support_points()
+{
+  // The support points for our shape functions are the nodes and
+  // the face midpoints, for a total of #vertices + #faces points
+  this->generalized_support_points.resize(this->dofs_per_cell);
+
+  // We need dim copies of each vertex for the first dim*vertices_per_cell
+  // generalized support points
+  for (unsigned int i = 0; i < dim * GeometryInfo<dim>::vertices_per_cell; ++i)
+    this->generalized_support_points[i] =
+      GeometryInfo<dim>::unit_cell_vertex(i / dim);
+
+  // The remaining 2*dim points are the edge midpoints
+  for (unsigned int i = 0; i < dim; ++i)
+    {
+      for (unsigned int j = 0; j < 2; ++j)
+        {
+          Point<dim> p;
+          p[0] = 0.5;
+          p[1] = 0.5;
+          if (dim == 3)
+            p[2] = 0.5;
+          p[i] = j;
+
+          unsigned int k =
+            dim * GeometryInfo<dim>::vertices_per_cell + i * 2 + j;
+          this->generalized_support_points[k] = p;
+        }
+    }
+}
+
+template class FE_BernardiRaugel<1>;
+template class FE_BernardiRaugel<2>;
+template class FE_BernardiRaugel<3>;
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc
index 98da782066..ae9506d861 100644
--- a/source/fe/fe_poly_tensor.cc
+++ b/source/fe/fe_poly_tensor.cc
@@ -18,6 +18,7 @@
 #include <deal.II/base/derivative_form.h>
 #include <deal.II/base/polynomials_abf.h>
 #include <deal.II/base/polynomials_bdm.h>
+#include <deal.II/base/polynomials_bernardi_raugel.h>
 #include <deal.II/base/polynomials_nedelec.h>
 #include <deal.II/base/polynomials_raviart_thomas.h>
 #include <deal.II/base/polynomials_rt_bubbles.h>
diff --git a/source/fe/fe_poly_tensor.inst.in b/source/fe/fe_poly_tensor.inst.in
index b073e9c8e5..7501ece4f9 100644
--- a/source/fe/fe_poly_tensor.inst.in
+++ b/source/fe/fe_poly_tensor.inst.in
@@ -23,6 +23,8 @@ for (deal_II_dimension : DIMENSIONS)
                                  deal_II_dimension>;
     template class FE_PolyTensor<PolynomialsABF<deal_II_dimension>,
                                  deal_II_dimension>;
+    template class FE_PolyTensor<PolynomialsBernardiRaugel<deal_II_dimension>,
+                                 deal_II_dimension>;
     template class FE_PolyTensor<PolynomialsBDM<deal_II_dimension>,
                                  deal_II_dimension>;
     template class FE_PolyTensor<PolynomialsNedelec<deal_II_dimension>,
-- 
2.39.5