]> https://gitweb.dealii.org/ - dealii.git/commitdiff
More straight-forward version of FE_RaviartThomasNodal
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 12 Jun 2021 19:06:39 +0000 (21:06 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 5 Jul 2021 17:42:41 +0000 (19:42 +0200)
include/deal.II/fe/fe_raviart_thomas.h
source/fe/fe_raviart_thomas.cc
source/fe/fe_raviart_thomas_nodal.cc
tests/fe/fe_conformity_dim_3_fe_rt_nodal.cc [new file with mode: 0644]
tests/fe/fe_conformity_dim_3_fe_rt_nodal.output [new file with mode: 0644]

index 0679aa3429bdf6bd74ba716f1803610d34923447..ab12990c5e8410454f85e232d20f091d2090b957 100644 (file)
 
 #include <deal.II/base/config.h>
 
-#include <deal.II/base/geometry_info.h>
-#include <deal.II/base/polynomial.h>
-#include <deal.II/base/polynomials_raviart_thomas.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>
@@ -309,29 +305,23 @@ private:
  * For this Raviart-Thomas element, the node values are not cell and face
  * moments with respect to certain polynomials, but the values in quadrature
  * points. Following the general scheme for numbering degrees of freedom, the
- * node values on edges are first, edge by edge, according to the natural
- * ordering of the edges of a cell. The interior degrees of freedom are last.
+ * node values on faces (edges in 2D, quads in 3D) are first, face by face,
+ * according to the natural ordering of the faces of a cell. The interior
+ * degrees of freedom are last.
  *
  * For an RT-element of degree <i>k</i>, we choose <i>(k+1)<sup>d-1</sup></i>
  * Gauss points on each face. These points are ordered lexicographically with
  * respect to the orientation of the face. This way, the normal component
- * which is in <i>Q<sub>k</sub></i> is uniquely determined. Furthermore, since
- * this Gauss-formula is exact on <i>Q<sub>2k+1</sub></i>, these node values
- * correspond to the exact integration of the moments of the RT-space.
+ * which is in <i>Q<sub>k</sub></i>, is uniquely determined. Furthermore,
+ * since this Gauss-formula is exact for polynomials of degree <i>2k+1</i>,
+ * these node values correspond to the exact integration of the moments of the
+ * RT-space.
  *
- * In the interior of the cells, the moments are with respect to an
- * anisotropic <i>Q<sub>k</sub></i> space, where the test functions are one
- * degree lower in the direction corresponding to the vector component under
- * consideration. This is emulated by using an anisotropic Gauss formula for
- * integration.
- *
- * @todo The current implementation is for Cartesian meshes only. You must use
- * MappingCartesian.
- *
- * @todo Even if this element is implemented for two and three space
- * dimensions, the definition of the node values relies on consistently
- * oriented faces in 3D. Therefore, care should be taken on complicated
- * meshes.
+ * These face polynomials are extended into the interior by the means of a
+ * QGaussLobatto formula for the normal direction. In other words, the
+ * polynomials are the tensor product of Lagrange polynomials on the points of
+ * a QGaussLobatto formula in the normal direction with Lagrange polynomials
+ * on the points of a QGauss quadrature formula.
  *
  * @note The degree stored in the member variable
  * FiniteElementData<dim>::degree is higher by one than the constructor
@@ -358,11 +348,6 @@ public:
   virtual std::unique_ptr<FiniteElement<dim, dim>>
   clone() const override;
 
-  virtual void
-  convert_generalized_support_point_values_to_dof_values(
-    const std::vector<Vector<double>> &support_point_values,
-    std::vector<double> &              nodal_values) const override;
-
   virtual void
   get_face_interpolation_matrix(const FiniteElement<dim> &source,
                                 FullMatrix<double> &      matrix,
@@ -374,6 +359,12 @@ public:
     const unsigned int        subface,
     FullMatrix<double> &      matrix,
     const unsigned int        face_no = 0) const override;
+
+  virtual void
+  convert_generalized_support_point_values_to_dof_values(
+    const std::vector<Vector<double>> &support_point_values,
+    std::vector<double> &              nodal_values) const override;
+
   virtual bool
   hp_constraints_are_implemented() const override;
 
@@ -394,51 +385,37 @@ public:
   compare_for_domination(const FiniteElement<dim> &fe_other,
                          const unsigned int codim = 0) const override final;
 
-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(const unsigned int degree);
+  virtual const FullMatrix<double> &
+  get_restriction_matrix(
+    const unsigned int         child,
+    const RefinementCase<dim> &refinement_case =
+      RefinementCase<dim>::isotropic_refinement) const override;
 
-  /**
-   * Compute the vector used for the @p restriction_is_additive field passed
-   * to the base class's constructor.
-   */
-  static std::vector<bool>
-  get_ria_vector(const unsigned int degree);
+  virtual const FullMatrix<double> &
+  get_prolongation_matrix(
+    const unsigned int         child,
+    const RefinementCase<dim> &refinement_case =
+      RefinementCase<dim>::isotropic_refinement) const override;
 
+private:
   /**
    * This function returns @p true, if the shape function @p shape_index has
    * non-zero function values somewhere on the face @p face_index.
-   *
-   * Right now, this is only implemented for RT0 in 1D. Otherwise, returns
-   * always @p true.
    */
   virtual bool
   has_support_on_face(const unsigned int shape_index,
                       const unsigned int face_index) const override;
 
-  /**
-   * 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(const unsigned int rt_degree);
-
   /**
    * Initialize the permutation pattern and the pattern of sign change.
    */
   void
   initialize_quad_dof_index_permutation_and_sign_change();
+
+  /*
+   * Mutex for protecting initialization of restriction and embedding matrix.
+   */
+  mutable Threads::Mutex mutex;
 };
 
 
index bcf3513c5f2daf866a21d7a1f054ef04a40067ce..85c78506b83c0eb0e26b061e3373b02088a71b7c 100644 (file)
@@ -14,6 +14,8 @@
 // ---------------------------------------------------------------------
 
 
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/polynomials_raviart_thomas.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/quadrature_lib.h>
index 086679cb463d058d1f6505124fe78f3be618fa16..aa4a1b9d73d8a54133914077b0c2885a06931694 100644 (file)
 // ---------------------------------------------------------------------
 
 
+#include <deal.II/base/polynomial.h>
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor_product_polynomials.h>
 
 #include <deal.II/dofs/dof_accessor.h>
 
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_raviart_thomas.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 <memory>
 #include <sstream>
 
 
 DEAL_II_NAMESPACE_OPEN
 
-// TODO: implement the adjust_quad_dof_index_for_face_orientation_table and
-// adjust_line_dof_index_for_line_orientation_table fields, and write tests
-// similar to bits/face_orientation_and_fe_q_*
+
+// ---------------- polynomial class for FE_RaviartThomasNodal ---------------
+
+namespace
+{
+  template <int dim>
+  class PolynomialsRaviartThomasNodal : public TensorPolynomialsBase<dim>
+  {
+  public:
+    PolynomialsRaviartThomasNodal(const unsigned int degree);
+
+    /**
+     * Compute the value and derivatives of each Raviart-Thomas polynomial at
+     * @p unit_point.
+     *
+     * The size of the vectors must either be zero or equal <tt>n()</tt>. In
+     * the first case, the function will not compute these values.
+     */
+    void
+    evaluate(const Point<dim> &           unit_point,
+             std::vector<Tensor<1, dim>> &values,
+             std::vector<Tensor<2, dim>> &grads,
+             std::vector<Tensor<3, dim>> &grad_grads,
+             std::vector<Tensor<4, dim>> &third_derivatives,
+             std::vector<Tensor<5, dim>> &fourth_derivatives) const override;
+
+    /**
+     * Return the name of the space, which is <tt>PolynomialsRaviartThomas</tt>.
+     */
+    std::string
+    name() const override;
+
+    /**
+     * Return the number of polynomials in the space without requiring to
+     * build an object of PolynomialsRaviartThomas. This is required by the
+     * FiniteElement classes.
+     */
+    static unsigned int
+    n_polynomials(const unsigned int degree);
+
+    const std::vector<unsigned int> &
+    get_renumbering() const;
+
+    /**
+     * @copydoc TensorPolynomialsBase<dim>::clone()
+     */
+    virtual std::unique_ptr<TensorPolynomialsBase<dim>>
+    clone() const override;
+
+    /**
+     * Compute the generalized support points of the associated element in the
+     * ordering of the element. Note that they are not support points in the
+     * classical sense as polynomials of the different components have
+     * different points, which need to be combined in terms of Piola
+     * transforms.
+     */
+    std::vector<Point<dim>>
+    get_polynomial_support_points() const;
+
+  private:
+    /**
+     * The degree variable passed to the constructor.
+     */
+    const unsigned int degree;
+
+    /**
+     * An object representing the polynomial space for a single component. We
+     * can re-use it by rotating the coordinates of the evaluation point.
+     */
+    const AnisotropicPolynomials<dim> polynomial_space;
+
+    /**
+     * Renumbering from lexicographic to hierarchic order.
+     */
+    std::vector<unsigned int> lexicographic_to_hierarchic;
+
+    /**
+     * Renumbering from hierarchic to lexicographic order. Inverse of
+     * lexicographic_to_hierarchic.
+     */
+    std::vector<unsigned int> hierarchic_to_lexicographic;
+
+    /**
+     * Renumbering from shifted polynomial spaces to lexicographic one
+     */
+    std::array<std::vector<unsigned int>, dim> renumber_aniso;
+  };
+
+
+
+  // Create nodal Raviart-Thomas polynomials as the tensor product of Lagrange
+  // polynomials on Gauss-Lobatto points of degree + 2 points in the
+  // continuous direction and degree + 1 points in the discontinuous
+  // directions (we could also choose Lagrange polynomials on Gauss points but
+  // those are slightly more expensive to handle in classes).
+  std::vector<std::vector<Polynomials::Polynomial<double>>>
+  create_rt_polynomials(const unsigned int dim, const unsigned int degree)
+  {
+    std::vector<std::vector<Polynomials::Polynomial<double>>> pols(dim);
+    pols[0] = Polynomials::generate_complete_Lagrange_basis(
+      QGaussLobatto<1>(degree + 2).get_points());
+    if (degree > 0)
+      for (unsigned int d = 1; d < dim; ++d)
+        pols[d] = Polynomials::generate_complete_Lagrange_basis(
+          QGaussLobatto<1>(degree + 1).get_points());
+    else
+      for (unsigned int d = 1; d < dim; ++d)
+        pols[d] = Polynomials::generate_complete_Lagrange_basis(
+          QMidpoint<1>().get_points());
+
+    return pols;
+  }
+
+
+
+  // set up the numbering of the rt polynomials
+  std::vector<unsigned int>
+  compute_rt_hierarchic_to_lexicographic(const unsigned int dim,
+                                         const unsigned int degree)
+  {
+    const unsigned int n_pols =
+      (degree + 2) * Utilities::pow(degree + 1, dim - 1);
+
+    std::vector<unsigned int> hierarchic_to_lexicographic;
+
+    // dofs on faces
+    for (unsigned int face_no = 0; face_no < 2 * dim; ++face_no)
+      {
+        const unsigned int stride_x = face_no < 2 ? degree + 2 : 1;
+        const unsigned int stride_y =
+          face_no < 4 ? (degree + 2) * (degree + 1) : degree + 1;
+        const unsigned int offset =
+          (face_no % 2) * Utilities::pow(degree + 1, 1 + face_no / 2);
+        for (unsigned int j = 0; j < (dim > 2 ? degree + 1 : 1); ++j)
+          for (unsigned int i = 0; i < degree + 1; ++i)
+            hierarchic_to_lexicographic.push_back(
+              (face_no / 2) * n_pols + offset + i * stride_x + j * stride_y);
+      }
+    // dofs on cells, starting with x component...
+    for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k)
+      for (unsigned int j = 0; j < (dim > 1 ? degree + 1 : 1); ++j)
+        for (unsigned int i = 1; i < degree + 1; ++i)
+          hierarchic_to_lexicographic.push_back(
+            k * (degree + 1) * (degree + 2) + j * (degree + 2) + i);
+    // ... then y component ...
+    if (dim > 1)
+      for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k)
+        for (unsigned int j = 1; j < degree + 1; ++j)
+          for (unsigned int i = 0; i < degree + 1; ++i)
+            hierarchic_to_lexicographic.push_back(
+              n_pols + k * (degree + 1) * (degree + 2) + j * (degree + 1) + i);
+    // ... and finally z component
+    if (dim > 2)
+      for (unsigned int k = 1; k < degree + 1; ++k)
+        for (unsigned int j = 0; j < degree + 1; ++j)
+          for (unsigned int i = 0; i < degree + 1; ++i)
+            hierarchic_to_lexicographic.push_back(
+              2 * n_pols + k * (degree + 1) * (degree + 1) + j * (degree + 1) +
+              i);
+
+    AssertDimension(hierarchic_to_lexicographic.size(), n_pols * dim);
+
+#ifdef DEBUG
+    // assert that we have a valid permutation
+    std::vector<unsigned int> copy(hierarchic_to_lexicographic);
+    std::sort(copy.begin(), copy.end());
+    for (unsigned int i = 0; i < copy.size(); ++i)
+      AssertDimension(i, copy[i]);
+#endif
+
+    return hierarchic_to_lexicographic;
+  }
+
+
+
+  template <int dim>
+  PolynomialsRaviartThomasNodal<dim>::PolynomialsRaviartThomasNodal(
+    const unsigned int degree)
+    : TensorPolynomialsBase<dim>(degree, n_polynomials(degree))
+    , degree(degree)
+    , polynomial_space(create_rt_polynomials(dim, degree))
+  {
+    // create renumbering of the unknowns from the lexicographic order to the
+    // actual order required by the finite element class with unknowns on
+    // faces placed first
+    const unsigned int n_pols = polynomial_space.n();
+    hierarchic_to_lexicographic =
+      compute_rt_hierarchic_to_lexicographic(dim, degree);
+
+    lexicographic_to_hierarchic =
+      Utilities::invert_permutation(hierarchic_to_lexicographic);
+
+    // since we only store an anisotropic polynomial for the first component,
+    // we set up a second numbering to switch out the actual coordinate
+    // direction
+    renumber_aniso[0].resize(n_pols);
+    for (unsigned int i = 0; i < n_pols; ++i)
+      renumber_aniso[0][i] = i;
+    if (dim > 1)
+      {
+        // switch x and y component (i and j loops)
+        renumber_aniso[1].resize(n_pols);
+        for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k)
+          for (unsigned int j = 0; j < degree + 2; ++j)
+            for (unsigned int i = 0; i < degree + 1; ++i)
+              renumber_aniso[1][(k * (degree + 2) + j) * (degree + 1) + i] =
+                j + i * (degree + 2) + k * (degree + 2) * (degree + 1);
+      }
+    if (dim > 2)
+      {
+        // switch x and z component (i and k loops)
+        renumber_aniso[2].resize(n_pols);
+        for (unsigned int k = 0; k < degree + 2; ++k)
+          for (unsigned int j = 0; j < degree + 1; ++j)
+            for (unsigned int i = 0; i < degree + 1; ++i)
+              renumber_aniso[2][(k * (degree + 1) + j) * (degree + 1) + i] =
+                k + i * (degree + 2) + j * (degree + 2) * (degree + 1);
+      }
+  }
+
+
+
+  template <int dim>
+  void
+  PolynomialsRaviartThomasNodal<dim>::evaluate(
+    const Point<dim> &           unit_point,
+    std::vector<Tensor<1, dim>> &values,
+    std::vector<Tensor<2, dim>> &grads,
+    std::vector<Tensor<3, dim>> &grad_grads,
+    std::vector<Tensor<4, dim>> &third_derivatives,
+    std::vector<Tensor<5, dim>> &fourth_derivatives) const
+  {
+    Assert(values.size() == this->n() || values.size() == 0,
+           ExcDimensionMismatch(values.size(), this->n()));
+    Assert(grads.size() == this->n() || grads.size() == 0,
+           ExcDimensionMismatch(grads.size(), this->n()));
+    Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
+           ExcDimensionMismatch(grad_grads.size(), this->n()));
+    Assert(third_derivatives.size() == this->n() ||
+             third_derivatives.size() == 0,
+           ExcDimensionMismatch(third_derivatives.size(), this->n()));
+    Assert(fourth_derivatives.size() == this->n() ||
+             fourth_derivatives.size() == 0,
+           ExcDimensionMismatch(fourth_derivatives.size(), this->n()));
+
+    std::vector<double>         p_values;
+    std::vector<Tensor<1, dim>> p_grads;
+    std::vector<Tensor<2, dim>> p_grad_grads;
+    std::vector<Tensor<3, dim>> p_third_derivatives;
+    std::vector<Tensor<4, dim>> p_fourth_derivatives;
+
+    const unsigned int n_sub = polynomial_space.n();
+    p_values.resize((values.size() == 0) ? 0 : n_sub);
+    p_grads.resize((grads.size() == 0) ? 0 : n_sub);
+    p_grad_grads.resize((grad_grads.size() == 0) ? 0 : n_sub);
+    p_third_derivatives.resize((third_derivatives.size() == 0) ? 0 : n_sub);
+    p_fourth_derivatives.resize((fourth_derivatives.size() == 0) ? 0 : n_sub);
+
+    for (unsigned int d = 0; d < dim; ++d)
+      {
+        // First we copy the point. The polynomial space for component d
+        // consists of polynomials of degree k in x_d and degree k+1 in the
+        // other variables. in order to simplify this, we use the same
+        // AnisotropicPolynomial space and simply rotate the coordinates
+        // through all directions.
+        Point<dim> p;
+        for (unsigned int c = 0; c < dim; ++c)
+          p(c) = unit_point((c + d) % dim);
+
+        polynomial_space.evaluate(p,
+                                  p_values,
+                                  p_grads,
+                                  p_grad_grads,
+                                  p_third_derivatives,
+                                  p_fourth_derivatives);
+
+        for (unsigned int i = 0; i < p_values.size(); ++i)
+          values[lexicographic_to_hierarchic[i + d * n_sub]][d] =
+            p_values[renumber_aniso[d][i]];
+
+        for (unsigned int i = 0; i < p_grads.size(); ++i)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            grads[lexicographic_to_hierarchic[i + d * n_sub]][d]
+                 [(d1 + d) % dim] = p_grads[renumber_aniso[d][i]][d1];
+
+        for (unsigned int i = 0; i < p_grad_grads.size(); ++i)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              grad_grads[lexicographic_to_hierarchic[i + d * n_sub]][d]
+                        [(d1 + d) % dim][(d2 + d) % dim] =
+                          p_grad_grads[renumber_aniso[d][i]][d1][d2];
+
+        for (unsigned int i = 0; i < p_third_derivatives.size(); ++i)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              for (unsigned int d3 = 0; d3 < dim; ++d3)
+                third_derivatives[lexicographic_to_hierarchic[i + d * n_sub]][d]
+                                 [(d1 + d) % dim][(d2 + d) % dim]
+                                 [(d3 + d) % dim] =
+                                   p_third_derivatives[renumber_aniso[d][i]][d1]
+                                                      [d2][d3];
+
+        for (unsigned int i = 0; i < p_fourth_derivatives.size(); ++i)
+          for (unsigned int d1 = 0; d1 < dim; ++d1)
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              for (unsigned int d3 = 0; d3 < dim; ++d3)
+                for (unsigned int d4 = 0; d4 < dim; ++d4)
+                  fourth_derivatives[lexicographic_to_hierarchic[i + d * n_sub]]
+                                    [d][(d1 + d) % dim][(d2 + d) % dim]
+                                    [(d3 + d) % dim][(d4 + d) % dim] =
+                                      p_fourth_derivatives[renumber_aniso[d][i]]
+                                                          [d1][d2][d3][d4];
+      }
+  }
+
+
+
+  template <int dim>
+  std::string
+  PolynomialsRaviartThomasNodal<dim>::name() const
+  {
+    return "PolynomialsRaviartThomasNodal";
+  }
+
+
+
+  template <int dim>
+  unsigned int
+  PolynomialsRaviartThomasNodal<dim>::n_polynomials(unsigned int degree)
+  {
+    return dim * (degree + 2) * Utilities::pow(degree + 1, dim - 1);
+  }
+
+
+
+  template <int dim>
+  std::unique_ptr<TensorPolynomialsBase<dim>>
+  PolynomialsRaviartThomasNodal<dim>::clone() const
+  {
+    return std::make_unique<PolynomialsRaviartThomasNodal<dim>>(*this);
+  }
+
+
+
+  template <int dim>
+  std::vector<Point<dim>>
+  PolynomialsRaviartThomasNodal<dim>::get_polynomial_support_points() const
+  {
+    Assert(dim > 0 && dim <= 3, ExcImpossibleInDim(dim));
+    const Quadrature<1> low(
+      degree == 0 ? static_cast<Quadrature<1>>(QMidpoint<1>()) :
+                    static_cast<Quadrature<1>>(QGaussLobatto<1>(degree + 1)));
+    const QGaussLobatto<1>  high(degree + 2);
+    const QAnisotropic<dim> quad =
+      (dim == 1 ? QAnisotropic<dim>(high) :
+                  (dim == 2 ? QAnisotropic<dim>(high, low) :
+                              QAnisotropic<dim>(high, low, low)));
+
+    const unsigned int      n_sub = polynomial_space.n();
+    std::vector<Point<dim>> points(dim * n_sub);
+    points.resize(n_polynomials(degree));
+    for (unsigned int d = 0; d < dim; ++d)
+      for (unsigned int i = 0; i < n_sub; ++i)
+        points[lexicographic_to_hierarchic[i + d * n_sub]] =
+          quad.point(renumber_aniso[d][i]);
+    return points;
+  }
+
+
+
+  // Return a vector of "dofs per object" where the components of the returned
+  // vector refer to:
+  // 0 = vertex
+  // 1 = edge
+  // 2 = face (which is a cell in 2D)
+  // 3 = cell
+  std::vector<unsigned int>
+  get_rt_dpo_vector(const unsigned int dim, const unsigned int degree)
+  {
+    std::vector<unsigned int> dpo(dim + 1);
+    dpo[0]                     = 0;
+    dpo[1]                     = 0;
+    unsigned int dofs_per_face = 1;
+    for (unsigned int d = 1; d < dim; ++d)
+      dofs_per_face *= (degree + 1);
+
+    dpo[dim - 1] = dofs_per_face;
+    dpo[dim]     = dim * degree * dofs_per_face;
+
+    return dpo;
+  }
+} // namespace
+
+
+
+// --------------------- actual implementation of element --------------------
 
 template <int dim>
-FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int deg)
-  : FE_PolyTensor<dim>(PolynomialsRaviartThomas<dim>(deg),
-                       FiniteElementData<dim>(get_dpo_vector(deg),
+FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int degree)
+  : FE_PolyTensor<dim>(PolynomialsRaviartThomasNodal<dim>(degree),
+                       FiniteElementData<dim>(get_rt_dpo_vector(dim, degree),
                                               dim,
-                                              deg + 1,
+                                              degree + 1,
                                               FiniteElementData<dim>::Hdiv),
-                       get_ria_vector(deg),
+                       std::vector<bool>(1, false),
                        std::vector<ComponentMask>(
-                         PolynomialsRaviartThomas<dim>::n_polynomials(deg),
+                         PolynomialsRaviartThomasNodal<dim>::n_polynomials(
+                           degree),
                          std::vector<bool>(dim, true)))
 {
   Assert(dim >= 2, ExcImpossibleInDim(dim));
-  const unsigned int n_dofs = this->n_dofs_per_cell();
 
   this->mapping_kind = {mapping_raviart_thomas};
-  // First, initialize the
-  // generalized support points and
-  // quadrature weights, since they
-  // are required for interpolation.
-  initialize_support_points(deg);
-
-  // Now compute the inverse node matrix, generating the correct
-  // basis functions from the raw ones. For a discussion of what
-  // exactly happens here, see FETools::compute_node_matrix.
-  const FullMatrix<double> M = FETools::compute_node_matrix(*this);
-  this->inverse_node_matrix.reinit(n_dofs, n_dofs);
-  this->inverse_node_matrix.invert(M);
-  // From now on, the shape functions provided by FiniteElement::shape_value
-  // and similar functions will be the correct ones, not
-  // the raw shape functions from the polynomial space anymore.
-
-  // Reinit the vectors of
-  // prolongation matrices to the
-  // right sizes. There are no
-  // restriction matrices implemented
-  for (unsigned int ref_case = RefinementCase<dim>::cut_x;
-       ref_case < RefinementCase<dim>::isotropic_refinement + 1;
-       ++ref_case)
-    {
-      const unsigned int nc =
-        GeometryInfo<dim>::n_children(RefinementCase<dim>(ref_case));
 
-      for (unsigned int i = 0; i < nc; ++i)
-        this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
-    }
+  // First, initialize the generalized support points and quadrature weights,
+  // since they are required for interpolation.
+  this->generalized_support_points =
+    PolynomialsRaviartThomasNodal<dim>(degree).get_polynomial_support_points();
+  AssertDimension(this->generalized_support_points.size(),
+                  this->n_dofs_per_cell());
 
-  // TODO: the implementation makes the assumption that all faces have the
-  // same number of dofs
-  AssertDimension(this->n_unique_faces(), 1);
   const unsigned int face_no = 0;
+  if (dim > 1)
+    this->generalized_face_support_points[face_no] =
+      degree == 0 ? QGauss<dim - 1>(1).get_points() :
+                    QGaussLobatto<dim - 1>(degree + 1).get_points();
 
-  // Fill prolongation matrices with embedding operators
-  FETools::compute_embedding_matrices(*this, this->prolongation);
-  // TODO[TL]: for anisotropic refinement we will probably need a table of
-  // submatrices with an array for each refine case
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
     face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
@@ -103,7 +468,7 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int deg)
                                                         face_embeddings,
                                                         0,
                                                         0);
-  this->interface_constraints.reinit((1 << (dim - 1)) *
+  this->interface_constraints.reinit(GeometryInfo<dim>::max_children_per_face *
                                        this->n_dofs_per_face(face_no),
                                      this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
@@ -126,21 +491,14 @@ template <int dim>
 std::string
 FE_RaviartThomasNodal<dim>::get_name() const
 {
-  // note that the
-  // FETools::get_fe_by_name
-  // function depends on the
-  // particular format of the string
-  // this function returns, so they
-  // have to be kept in synch
-
-  // note that this->degree is the maximal
-  // polynomial degree and is thus one higher
-  // than the argument given to the
-  // constructor
-  std::ostringstream namebuf;
-  namebuf << "FE_RaviartThomasNodal<" << dim << ">(" << this->degree - 1 << ")";
-
-  return namebuf.str();
+  // note that the FETools::get_fe_by_name function depends on the particular
+  // format of the string this function returns, so they have to be kept in
+  // synch
+
+  // note that this->degree is the maximal polynomial degree and is thus one
+  // higher than the argument given to the constructor
+  return "FE_RaviartThomasNodal<" + std::to_string(dim) + ">(" +
+         std::to_string(this->degree - 1) + ")";
 }
 
 
@@ -158,94 +516,6 @@ FE_RaviartThomasNodal<dim>::clone() const
 
 
 
-template <int dim>
-void
-FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
-{
-  // TODO: the implementation makes the assumption that all faces have the
-  // same number of dofs
-  AssertDimension(this->n_unique_faces(), 1);
-  const unsigned int face_no = 0;
-
-  this->generalized_support_points.resize(this->n_dofs_per_cell());
-  this->generalized_face_support_points[face_no].resize(
-    this->n_dofs_per_face(face_no));
-
-  // Number of the point being entered
-  unsigned int current = 0;
-
-  // On the faces, we choose as many
-  // Gauss points as necessary to
-  // determine the normal component
-  // uniquely. This is the deg of
-  // the Raviart-Thomas element plus
-  // one.
-  if (dim > 1)
-    {
-      QGauss<dim - 1> face_points(deg + 1);
-      Assert(face_points.size() == this->n_dofs_per_face(face_no),
-             ExcInternalError());
-      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
-        this->generalized_face_support_points[face_no][k] =
-          face_points.point(k);
-      Quadrature<dim> faces =
-        QProjector<dim>::project_to_all_faces(this->reference_cell(),
-                                              face_points);
-      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
-                                     GeometryInfo<dim>::faces_per_cell;
-           ++k)
-        this->generalized_support_points[k] = faces.point(
-          k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell(),
-                                                       0,
-                                                       true,
-                                                       false,
-                                                       false,
-                                                       this->n_dofs_per_face(
-                                                         face_no)));
-
-      current =
-        this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
-    }
-
-  if (deg == 0)
-    return;
-  // In the interior, we need
-  // anisotropic Gauss quadratures,
-  // different for each direction.
-  QGauss<1> high(deg + 1);
-  QGauss<1> low(deg);
-
-  for (unsigned int d = 0; d < dim; ++d)
-    {
-      std::unique_ptr<QAnisotropic<dim>> quadrature;
-      switch (dim)
-        {
-          case 1:
-            quadrature = std::make_unique<QAnisotropic<dim>>(high);
-            break;
-          case 2:
-            quadrature =
-              std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
-                                                  ((d == 1) ? low : high));
-            break;
-          case 3:
-            quadrature =
-              std::make_unique<QAnisotropic<dim>>(((d == 0) ? low : high),
-                                                  ((d == 1) ? low : high),
-                                                  ((d == 2) ? low : high));
-            break;
-          default:
-            Assert(false, ExcNotImplemented());
-        }
-
-      for (unsigned int k = 0; k < quadrature->size(); ++k)
-        this->generalized_support_points[current++] = quadrature->point(k);
-    }
-  Assert(current == this->n_dofs_per_cell(), ExcInternalError());
-}
-
-
-
 template <int dim>
 void
 FE_RaviartThomasNodal<
@@ -255,68 +525,56 @@ FE_RaviartThomasNodal<
   if (dim < 3)
     return;
 
-  // TODO: Implement this for this class
-  return;
-}
-
-
-
-template <int dim>
-std::vector<unsigned int>
-FE_RaviartThomasNodal<dim>::get_dpo_vector(const unsigned int deg)
-{
-  // the element is face-based and we have
-  // (deg+1)^(dim-1) DoFs per face
-  unsigned int dofs_per_face = 1;
-  for (unsigned int d = 1; d < dim; ++d)
-    dofs_per_face *= deg + 1;
-
-  // and then there are interior dofs
-  const unsigned int interior_dofs = dim * deg * dofs_per_face;
-
-  std::vector<unsigned int> dpo(dim + 1);
-  dpo[dim - 1] = dofs_per_face;
-  dpo[dim]     = interior_dofs;
-
-  return dpo;
-}
-
-
-
-template <>
-std::vector<bool>
-FE_RaviartThomasNodal<1>::get_ria_vector(const unsigned int)
-{
-  Assert(false, ExcImpossibleInDim(1));
-  return std::vector<bool>();
+  const unsigned int n       = this->degree;
+  const unsigned int face_no = 0;
+  Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
+  for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
+    // face support points are in lexicographic ordering with x running
+    // fastest. invert that (y running fastest)
+    {
+      unsigned int i = local % n, j = local / n;
+
+      // face_orientation=false, face_flip=false, face_rotation=false
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      0) =
+        j + i * n - local;
+      // face_orientation=false, face_flip=false, face_rotation=true
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      1) =
+        i + (n - 1 - j) * n - local;
+      // face_orientation=false, face_flip=true,  face_rotation=false
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      2) =
+        (n - 1 - j) + (n - 1 - i) * n - local;
+      // face_orientation=false, face_flip=true,  face_rotation=true
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      3) =
+        (n - 1 - i) + j * n - local;
+      // face_orientation=true,  face_flip=false, face_rotation=false
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      4) = 0;
+      // face_orientation=true,  face_flip=false, face_rotation=true
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      5) =
+        j + (n - 1 - i) * n - local;
+      // face_orientation=true,  face_flip=true,  face_rotation=false
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      6) =
+        (n - 1 - i) + (n - 1 - j) * n - local;
+      // face_orientation=true,  face_flip=true,  face_rotation=true
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      7) =
+        (n - 1 - j) + i * n - local;
+
+      // for face_orientation == false, we need to switch the sign
+      for (unsigned int i = 0; i < 4; ++i)
+        this->adjust_quad_dof_sign_for_face_orientation_table[face_no](local,
+                                                                       i) = 1;
+    }
 }
 
 
 
-template <int dim>
-std::vector<bool>
-FE_RaviartThomasNodal<dim>::get_ria_vector(const unsigned int deg)
-{
-  const unsigned int dofs_per_cell =
-    PolynomialsRaviartThomas<dim>::n_polynomials(deg);
-  unsigned int dofs_per_face = deg + 1;
-  for (unsigned int d = 2; d < dim; ++d)
-    dofs_per_face *= deg + 1;
-  // all face dofs need to be
-  // non-additive, since they have
-  // continuity requirements.
-  // however, the interior dofs are
-  // made additive
-  std::vector<bool> ret_val(dofs_per_cell, false);
-  for (unsigned int i = GeometryInfo<dim>::faces_per_cell * dofs_per_face;
-       i < dofs_per_cell;
-       ++i)
-    ret_val[i] = true;
-
-  return ret_val;
-}
-
-
 template <int dim>
 bool
 FE_RaviartThomasNodal<dim>::has_support_on_face(
@@ -326,20 +584,16 @@ FE_RaviartThomasNodal<dim>::has_support_on_face(
   AssertIndexRange(shape_index, this->n_dofs_per_cell());
   AssertIndexRange(face_index, GeometryInfo<dim>::faces_per_cell);
 
-  // The first degrees of freedom are
-  // on the faces and each face has
-  // degree degrees.
-  const unsigned int support_face = shape_index / this->degree;
+  // The first degrees of freedom are on the faces and each face has degree
+  // degrees.
+  const unsigned int support_face = shape_index / this->n_dofs_per_face();
 
-  // The only thing we know for sure
-  // is that shape functions with
-  // support on one face are zero on
-  // the opposite face.
+  // The only thing we know for sure is that shape functions with support on
+  // one face are zero on the opposite face.
   if (support_face < GeometryInfo<dim>::faces_per_cell)
     return (face_index != GeometryInfo<dim>::opposite_face[support_face]);
 
-  // In all other cases, return true,
-  // which is safe
+  // In all other cases, return true, which is safe
   return true;
 }
 
@@ -361,10 +615,8 @@ FE_RaviartThomasNodal<dim>::
          ExcDimensionMismatch(support_point_values[0].size(),
                               this->n_components()));
 
-  // First do interpolation on
-  // faces. There, the component
-  // evaluated depends on the face
-  // direction and orientation.
+  // First do interpolation on faces. There, the component evaluated depends
+  // on the face direction and orientation.
   unsigned int fbase = 0;
   unsigned int f     = 0;
   for (; f < GeometryInfo<dim>::faces_per_cell;
@@ -377,8 +629,7 @@ FE_RaviartThomasNodal<dim>::
         }
     }
 
-  // The remaining points form dim
-  // chunks, one for each component.
+  // The remaining points form dim chunks, one for each component
   const unsigned int istep = (this->n_dofs_per_cell() - fbase) / dim;
   Assert((this->n_dofs_per_cell() - fbase) % dim == 0, ExcInternalError());
 
@@ -416,11 +667,9 @@ std::vector<std::pair<unsigned int, unsigned int>>
 FE_RaviartThomasNodal<dim>::hp_vertex_dof_identities(
   const FiniteElement<dim> &fe_other) const
 {
-  // we can presently only compute these
-  // identities if both FEs are
-  // FE_RaviartThomasNodals or the other is FE_Nothing.
-  // In either case, no dofs are assigned on the vertex,
-  // so we shouldn't be getting here at all.
+  // we can presently only compute these identities if both FEs are
+  // FE_RaviartThomasNodals or the other is FE_Nothing.  In either case, no
+  // dofs are assigned on the vertex, so we shouldn't be getting here at all.
   if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other) != nullptr)
     return std::vector<std::pair<unsigned int, unsigned int>>();
   else if (dynamic_cast<const FE_Nothing<dim> *>(&fe_other) != nullptr)
@@ -439,36 +688,27 @@ std::vector<std::pair<unsigned int, unsigned int>>
 FE_RaviartThomasNodal<dim>::hp_line_dof_identities(
   const FiniteElement<dim> &fe_other) const
 {
-  // we can presently only compute
-  // these identities if both FEs are
-  // FE_RaviartThomasNodals or if the other
-  // one is FE_Nothing
+  // we can presently only compute these identities if both FEs are
+  // FE_RaviartThomasNodals or if the other one is FE_Nothing
   if (const FE_RaviartThomasNodal<dim> *fe_q_other =
         dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
     {
-      // dofs are located on faces; these are
-      // only lines in 2d
+      // dofs are located on faces; these are only lines in 2d
       if (dim != 2)
         return std::vector<std::pair<unsigned int, unsigned int>>();
 
-      // dofs are located along lines, so two
-      // dofs are identical only if in the
-      // following two cases (remember that
-      // the face support points are Gauss
-      // points):
+      // dofs are located along lines, so two dofs are identical only if in
+      // the following two cases (remember that the face support points are
+      // Gauss points):
       // 1. this->degree = fe_q_other->degree,
-      //   in the case, all the dofs on
-      //   the line are identical
+      //   in the case, all the dofs on the line are identical
       // 2. this->degree-1 and fe_q_other->degree-1
-      //   are both even, i.e. this->dof_per_line
-      //   and fe_q_other->dof_per_line are both odd,
-      //   there exists only one point (the middle one)
-      //   such that dofs are identical on this point
+      //   are both even, i.e. this->dof_per_line and fe_q_other->dof_per_line
+      //   are both odd, there exists only one point (the middle one) such
+      //   that dofs are identical on this point
       //
-      // to understand this, note that
-      // this->degree is the *maximal*
-      // polynomial degree, and is thus one
-      // higher than the argument given to
+      // to understand this, note that this->degree is the *maximal*
+      // polynomial degree, and is thus one higher than the argument given to
       // the constructor
       const unsigned int p = this->degree - 1;
       const unsigned int q = fe_q_other->degree - 1;
@@ -504,20 +744,16 @@ FE_RaviartThomasNodal<dim>::hp_quad_dof_identities(
   const FiniteElement<dim> &fe_other,
   const unsigned int        face_no) const
 {
-  // we can presently only compute
-  // these identities if both FEs are
-  // FE_RaviartThomasNodals or if the other
-  // one is FE_Nothing
+  // we can presently only compute these identities if both FEs are
+  // FE_RaviartThomasNodals or if the other one is FE_Nothing
   if (const FE_RaviartThomasNodal<dim> *fe_q_other =
         dynamic_cast<const FE_RaviartThomasNodal<dim> *>(&fe_other))
     {
-      // dofs are located on faces; these are
-      // only quads in 3d
+      // dofs are located on faces; these are only quads in 3d
       if (dim != 3)
         return std::vector<std::pair<unsigned int, unsigned int>>();
 
-      // this works exactly like the line
-      // case above
+      // this works exactly like the line case above
       const unsigned int p = this->n_dofs_per_quad(face_no);
 
       AssertDimension(fe_q_other->n_unique_faces(), 1);
@@ -633,44 +869,31 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
          ExcDimensionMismatch(interpolation_matrix.m(),
                               x_source_fe.n_dofs_per_face(face_no)));
 
-  // ok, source is a RaviartThomasNodal element, so
-  // we will be able to do the work
+  // ok, source is a RaviartThomasNodal element, so we will be able to do the
+  // work
   const FE_RaviartThomasNodal<dim> &source_fe =
     dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
 
-  // Make sure, that the element,
-  // for which the DoFs should be
-  // constrained is the one with
-  // the higher polynomial degree.
-  // Actually the procedure will work
-  // also if this assertion is not
-  // satisfied. But the matrices
-  // produced in that case might
-  // lead to problems in the
-  // hp-procedures, which use this
+  // Make sure that the element for which the DoFs should be constrained is
+  // the one with the higher polynomial degree.  Actually the procedure will
+  // work also if this assertion is not satisfied. But the matrices produced
+  // in that case might lead to problems in the hp-procedures, which use this
   // method.
   Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
-  // generate a quadrature
-  // with the generalized support points.
-  // This is later based as a
-  // basis for the QProjector,
-  // which returns the support
-  // points on the face.
+  // generate a quadrature with the generalized support points.  This is later
+  // based as a basis for the QProjector, which returns the support points on
+  // the face.
   Quadrature<dim - 1> quad_face_support(
     source_fe.generalized_face_support_points[face_no]);
 
-  // Rule of thumb for FP accuracy,
-  // that can be expected for a
-  // given polynomial degree.
-  // This value is used to cut
-  // off values close to zero.
+  // Rule of thumb for FP accuracy, that can be expected for a given
+  // polynomial degree.  This value is used to cut off values close to zero.
   double eps = 2e-13 * this->degree * (dim - 1);
 
-  // compute the interpolation
-  // matrix by simply taking the
-  // value at the support points.
+  // compute the interpolation matrix by simply taking the value at the
+  // support points.
   const Quadrature<dim> face_projection =
     QProjector<dim>::project_to_face(this->reference_cell(),
                                      quad_face_support,
@@ -685,12 +908,9 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
           double matrix_entry =
             this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
 
-          // Correct the interpolated
-          // value. I.e. if it is close
-          // to 1 or 0, make it exactly
-          // 1 or 0. Unfortunately, this
-          // is required to avoid problems
-          // with higher order elements.
+          // Correct the interpolated value. I.e. if it is close to 1 or 0,
+          // make it exactly 1 or 0. Unfortunately, this is required to avoid
+          // problems with higher order elements.
           if (std::fabs(matrix_entry - 1.0) < eps)
             matrix_entry = 1.0;
           if (std::fabs(matrix_entry) < eps)
@@ -700,11 +920,8 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
         }
     }
 
-  // make sure that the row sum of
-  // each of the matrices is 1 at
-  // this point. this must be so
-  // since the shape functions sum up
-  // to 1
+  // make sure that the row sum of each of the matrices is 1 at this
+  // point. this must be so since the shape functions sum up to 1
   for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
     {
       double sum = 0.;
@@ -726,9 +943,8 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
   FullMatrix<double> &      interpolation_matrix,
   const unsigned int        face_no) const
 {
-  // this is only implemented, if the
-  // source FE is also a
-  // RaviartThomasNodal element
+  // this is only implemented, if the source FE is also a RaviartThomasNodal
+  // element
   AssertThrow((x_source_fe.get_name().find("FE_RaviartThomasNodal<") == 0) ||
                 (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
                    &x_source_fe) != nullptr),
@@ -741,45 +957,31 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
          ExcDimensionMismatch(interpolation_matrix.m(),
                               x_source_fe.n_dofs_per_face(face_no)));
 
-  // ok, source is a RaviartThomasNodal element, so
-  // we will be able to do the work
+  // ok, source is a RaviartThomasNodal element, so we will be able to do the
+  // work
   const FE_RaviartThomasNodal<dim> &source_fe =
     dynamic_cast<const FE_RaviartThomasNodal<dim> &>(x_source_fe);
 
-  // Make sure, that the element,
-  // for which the DoFs should be
-  // constrained is the one with
-  // the higher polynomial degree.
-  // Actually the procedure will work
-  // also if this assertion is not
-  // satisfied. But the matrices
-  // produced in that case might
-  // lead to problems in the
-  // hp-procedures, which use this
+  // Make sure that the element for which the DoFs should be constrained is
+  // the one with the higher polynomial degree. Actually the procedure will
+  // work also if this assertion is not satisfied. But the matrices produced
+  // in that case might lead to problems in the hp-procedures, which use this
   // method.
   Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
-  // generate a quadrature
-  // with the generalized support points.
-  // This is later based as a
-  // basis for the QProjector,
-  // which returns the support
-  // points on the face.
+  // generate a quadrature with the generalized support points.  This is later
+  // based as a basis for the QProjector, which returns the support points on
+  // the face.
   Quadrature<dim - 1> quad_face_support(
     source_fe.generalized_face_support_points[face_no]);
 
-  // Rule of thumb for FP accuracy,
-  // that can be expected for a
-  // given polynomial degree.
-  // This value is used to cut
-  // off values close to zero.
+  // Rule of thumb for FP accuracy, that can be expected for a given
+  // polynomial degree. This value is used to cut off values close to zero.
   double eps = 2e-13 * this->degree * (dim - 1);
 
-  // compute the interpolation
-  // matrix by simply taking the
-  // value at the support points.
-
+  // compute the interpolation matrix by simply taking the value at the
+  // support points.
   const Quadrature<dim> subface_projection =
     QProjector<dim>::project_to_subface(this->reference_cell(),
                                         quad_face_support,
@@ -795,12 +997,9 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
           double matrix_entry =
             this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
 
-          // Correct the interpolated
-          // value. I.e. if it is close
-          // to 1 or 0, make it exactly
-          // 1 or 0. Unfortunately, this
-          // is required to avoid problems
-          // with higher order elements.
+          // Correct the interpolated value. I.e. if it is close to 1 or 0,
+          // make it exactly 1 or 0. Unfortunately, this is required to avoid
+          // problems with higher order elements.
           if (std::fabs(matrix_entry - 1.0) < eps)
             matrix_entry = 1.0;
           if (std::fabs(matrix_entry) < eps)
@@ -810,11 +1009,8 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
         }
     }
 
-  // make sure that the row sum of
-  // each of the matrices is 1 at
-  // this point. this must be so
-  // since the shape functions sum up
-  // to 1
+  // make sure that the row sum of each of the matrices is 1 at this
+  // point. this must be so since the shape functions sum up to 1
   for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
     {
       double sum = 0.;
@@ -829,6 +1025,122 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
 
 
 
+template <int dim>
+const FullMatrix<double> &
+FE_RaviartThomasNodal<dim>::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(this->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_RaviartThomasNodal<dim> &this_nonconst =
+        const_cast<FE_RaviartThomasNodal<dim> &>(*this);
+      if (refinement_case == RefinementCase<dim>::isotropic_refinement)
+        {
+          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_embedding_matrices(*this, isotropic_matrices, true);
+          this_nonconst.prolongation[refinement_case - 1].swap(
+            isotropic_matrices.back());
+        }
+      else
+        {
+          // must compute both restriction and prolongation matrices because
+          // we only check for their size and the reinit call initializes them
+          // all
+          this_nonconst.reinit_restriction_and_prolongation_matrices();
+          FETools::compute_embedding_matrices(*this,
+                                              this_nonconst.prolongation);
+          FETools::compute_projection_matrices(*this,
+                                               this_nonconst.restriction);
+        }
+    }
+
+  // finally return the matrix
+  return this->prolongation[refinement_case - 1][child];
+}
+
+
+
+template <int dim>
+const FullMatrix<double> &
+FE_RaviartThomasNodal<dim>::get_restriction_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(
+           "Restriction matrices are only available for refined cells!"));
+  AssertIndexRange(child, GeometryInfo<dim>::n_children(refinement_case));
+
+  // 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
+      FE_RaviartThomasNodal<dim> &this_nonconst =
+        const_cast<FE_RaviartThomasNodal<dim> &>(*this);
+      if (refinement_case == RefinementCase<dim>::isotropic_refinement)
+        {
+          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());
+        }
+      else
+        {
+          // must compute both restriction and prolongation matrices because
+          // we only check for their size and the reinit call initializes them
+          // all
+          this_nonconst.reinit_restriction_and_prolongation_matrices();
+          FETools::compute_embedding_matrices(*this,
+                                              this_nonconst.prolongation);
+          FETools::compute_projection_matrices(*this,
+                                               this_nonconst.restriction);
+        }
+    }
+
+  // finally return the matrix
+  return this->restriction[refinement_case - 1][child];
+}
+
+
+
 // explicit instantiations
 #include "fe_raviart_thomas_nodal.inst"
 
diff --git a/tests/fe/fe_conformity_dim_3_fe_rt_nodal.cc b/tests/fe/fe_conformity_dim_3_fe_rt_nodal.cc
new file mode 100644 (file)
index 0000000..4df987b
--- /dev/null
@@ -0,0 +1,93 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+
+#include "../tests.h"
+
+// STL
+#include <fstream>
+#include <iostream>
+
+// My test headers
+#include "fe_conformity_test.h"
+
+#define PRECISION 4
+
+int
+main(int, char **)
+{
+  std::ofstream logfile("output");
+  dealii::deallog << std::setprecision(PRECISION);
+  dealii::deallog << std::fixed;
+  logfile << std::setprecision(PRECISION);
+  logfile << std::fixed;
+  dealii::deallog.attach(logfile);
+
+  try
+    {
+      using namespace FEConforimityTest;
+
+      constexpr int dim = 3;
+
+      for (unsigned int fe_degree = 0; fe_degree < 3; ++fe_degree)
+        {
+          // H(div) conformal
+          FE_RaviartThomasNodal<dim> fe(fe_degree);
+
+          {
+            for (unsigned int this_switch = 0; this_switch < (dim == 2 ? 4 : 8);
+                 ++this_switch)
+              {
+                deallog << std::endl
+                        << "*******   degree " << fe_degree
+                        << "   *******   orientation case " << this_switch
+                        << "   *******" << std::endl;
+
+                FEConformityTest<dim> fe_conformity_tester(fe, this_switch);
+                fe_conformity_tester.run();
+              }
+          }
+        } // ++fe_degree
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/fe/fe_conformity_dim_3_fe_rt_nodal.output b/tests/fe/fe_conformity_dim_3_fe_rt_nodal.output
new file mode 100644 (file)
index 0000000..c8f0cab
--- /dev/null
@@ -0,0 +1,577 @@
+
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 0   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 1   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 2   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 3   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 4   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 5   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 6   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 0   *******   orientation case 7   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(0)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 6
+DEAL::   n_dofs_per_face      : 1
+DEAL::   n_dofs_per_quad      : 1
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 0   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 1   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 2   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 3   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 4   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 5   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 6   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 1   *******   orientation case 7   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(1)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 36
+DEAL::   n_dofs_per_face      : 4
+DEAL::   n_dofs_per_quad      : 4
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 0   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 1   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 2   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 3   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 4   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 5   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 6   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::
+DEAL::*******   degree 2   *******   orientation case 7   *******
+DEAL::Element Info:  
+DEAL::   name                 : FE_RaviartThomasNodal<3>(2)
+DEAL::   is_primitive         : 0
+DEAL::   n_dofs_per_cell      : 108
+DEAL::   n_dofs_per_face      : 9
+DEAL::   n_dofs_per_quad      : 9
+DEAL::   n_dofs_per_line      : 0
+DEAL::   n_dofs_per_vertex    : 0
+DEAL::
+DEAL::   first_line_index     : 0
+DEAL::   first_quad_index     : 0
+DEAL::   first_face_line_index: 0
+DEAL::   first_face_quad_index: 0
+DEAL::
+DEAL::   n_components         : 3
+DEAL::   n_blocks             : 1
+DEAL::   n_base_elements      : 1
+DEAL::
+DEAL::Normal jumps (at quad points) in cell   0_0:   at face   1
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   
+DEAL::Normal jumps (at quad points) in cell   1_0:   at face   0
+DEAL::   interface jumps:   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   0.0000   

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.