]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Matrix-free evaluation of FE_RaviartThomasNodal 13591/head
authorNiklas Wik <niiklaswiik@gmail.com>
Tue, 1 Mar 2022 17:30:34 +0000 (18:30 +0100)
committerNiklas Wik <niiklaswiik@gmail.com>
Mon, 25 Apr 2022 23:06:20 +0000 (01:06 +0200)
Test added for matrix-free RT
Split of evaluation_template_factory.templates.h
Decrease FE_EVAL_FACTORY_DEGREE_MAX to 2 for windows vm

30 files changed:
.github/workflows/windows.yml
include/deal.II/fe/fe_raviart_thomas.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_template_face_factory.templates.h [new file with mode: 0644]
include/deal.II/matrix_free/evaluation_template_factory.templates.h
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
include/deal.II/matrix_free/tensor_product_kernels.h
source/fe/fe_raviart_thomas_nodal.cc
source/matrix_free/CMakeLists.txt
source/matrix_free/evaluation_template_face_factory.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory.inst.in [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory_inst2.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory_inst3.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory_inst4.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory_inst5.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_face_factory_inst6.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory.cc
source/matrix_free/evaluation_template_factory.inst.in
source/matrix_free/evaluation_template_factory_inst10.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_inst7.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_inst8.cc [new file with mode: 0644]
source/matrix_free/evaluation_template_factory_inst9.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_01.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_01.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_common.h [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_face_01.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_face_01.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_rt_face_common.h [new file with mode: 0644]

index b46066f10563f5d64fe5641488e698b929e6d5f0..cd90aafcaf9cd5a26e2d41ab1cd0c7257aaed8fc 100644 (file)
@@ -29,7 +29,7 @@ jobs:
         mkdir build
         mkdir c:/project
         cd build
-        cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX" -T host=x64 -A x64 ..
+        cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX /D FE_EVAL_FACTORY_DEGREE_MAX=2" -T host=x64 -A x64 ..
     - name: archive logs
       uses: actions/upload-artifact@v1
       with:
index ab12990c5e8410454f85e232d20f091d2090b957..fcd9e93e4d0097709b222d7d4e4c927cc453cddf 100644 (file)
@@ -327,6 +327,14 @@ private:
  * FiniteElementData<dim>::degree is higher by one than the constructor
  * argument!
  */
+
+namespace internal
+{
+  template <int dim>
+  std::vector<unsigned int>
+  get_lexicographic_numbering_rt_nodal(const unsigned int degree);
+} // namespace internal
+
 template <int dim>
 class FE_RaviartThomasNodal : public FE_PolyTensor<dim>
 {
index e2dc5a51472648fcbaf53083570c8e7a4297d4c4..1faeb7a335e6823ba2628a3dae6dfd2826e6cad8 100644 (file)
@@ -86,6 +86,12 @@ namespace internal
     static const EvaluatorVariant variant = evaluate_evenodd;
   };
 
+  template <bool is_long>
+  struct EvaluatorSelector<MatrixFreeFunctions::tensor_raviart_thomas, is_long>
+  {
+    static const EvaluatorVariant variant = evaluate_raviart_thomas;
+  };
+
 
 
   /**
@@ -180,6 +186,79 @@ namespace internal
               const bool                             add_into_values_array);
   };
 
+  /**
+   * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
+   * specific sum-factorization kernels and with normal/tangential shape_data
+   */
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                          dim,
+                          fe_degree,
+                          n_q_points_1d,
+                          Number>
+  {
+    using EvalNormal =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim,
+                             (fe_degree == -1) ? 1 : fe_degree + 1,
+                             n_q_points_1d,
+                             Number>;
+
+    using EvalTangent =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim,
+                             (fe_degree == -1) ? 1 : fe_degree,
+                             n_q_points_1d,
+                             Number>;
+    template <bool integrate>
+    static void
+    evaluate_or_integrate(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array = false);
+
+  private:
+    template <typename EvalType>
+    static EvalType
+    create_evaluator_tensor_product(
+      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+    {
+      // Shape_data is of size 2 with normal direction in [0] and tangential
+      // direction in [1]. See shape_info.template.h
+      const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
+      if (std::is_same<EvalType, EvalNormal>::value)
+        {
+          shape_data = &shape_info.data[0];
+        }
+      else
+        {
+          shape_data = &shape_info.data[1];
+        }
+      return EvalType(shape_data->shape_values,
+                      shape_data->shape_gradients,
+                      shape_data->shape_hessians);
+    }
+
+    template <int normal_dir>
+    static void
+    evaluate_tensor_product_per_component(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array,
+      std::integral_constant<bool, false>);
+
+    template <int normal_dir>
+    static void
+    evaluate_tensor_product_per_component(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array,
+      std::integral_constant<bool, true>);
+  };
+
 
 
   template <MatrixFreeFunctions::ElementType type,
@@ -918,6 +997,454 @@ namespace internal
   }
 
 
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <bool integrate>
+  inline void
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::
+    evaluate_or_integrate(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array)
+  {
+    if (evaluation_flag == EvaluationFlags::nothing)
+      return;
+
+    AssertDimension(fe_eval.get_shape_info().data.size(), 2);
+    // First component:
+    evaluate_tensor_product_per_component<0>(
+      evaluation_flag,
+      values_dofs_actual,
+      fe_eval,
+      add_into_values_array,
+      std::integral_constant<bool, integrate>());
+    // Second component :
+    evaluate_tensor_product_per_component<1>(
+      evaluation_flag,
+      values_dofs_actual,
+      fe_eval,
+      add_into_values_array,
+      std::integral_constant<bool, integrate>());
+    if (dim == 3)
+      {
+        // Third component
+        evaluate_tensor_product_per_component<2>(
+          evaluation_flag,
+          values_dofs_actual,
+          fe_eval,
+          add_into_values_array,
+          std::integral_constant<bool, integrate>());
+      }
+  }
+
+  // Helper function that applies the 1d evaluation kernals.
+  // std::integral_constant<bool, false> is the interpolation path, and
+  // std::integral_constant<bool, true> bellow is the integration path.
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int normal_dir>
+  inline void
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::
+    evaluate_tensor_product_per_component(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array,
+      std::integral_constant<bool, false>)
+  {
+    (void)add_into_values_array;
+    using Eval0 =
+      typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+    using Eval1 =
+      typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
+    using Eval2 =
+      typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
+
+    Eval0 eval0 =
+      create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
+    Eval1 eval1 =
+      create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
+    Eval2 eval2 =
+      create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+
+    Number *temp1 = fe_eval.get_scratch_data().begin();
+    Number *temp2;
+
+    temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+                               fe_eval.get_shape_info().data[0].fe_degree + 1),
+                             Utilities::fixed_power<dim>(
+                               fe_eval.get_shape_info().data[0].n_q_points_1d));
+
+    const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
+    const std::size_t dofs_per_comp =
+      fe_eval.get_shape_info().dofs_per_component_on_cell;
+
+    // Initial shift depending on component (normal_dir)
+    Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
+    Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
+    Number *gradients_quad =
+      fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
+    Number *hessians_quad =
+      (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
+                   fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
+
+    switch (dim)
+      {
+        case 2:
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              eval0.template gradients<0, true, false, normal_dir>(values_dofs,
+                                                                   temp1);
+              eval1.template values<1, true, false, normal_dir>(temp1,
+                                                                gradients_quad);
+            }
+          if (evaluation_flag & EvaluationFlags::hessians)
+            {
+              // The evaluation/integration here *should* work, however
+              // the piola transform is not implemented.
+              AssertThrow(false, ExcNotImplemented());
+              // grad xy
+              if (!(evaluation_flag & EvaluationFlags::gradients))
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs, temp1);
+              eval1.template gradients<1, true, false, normal_dir>(
+                temp1, hessians_quad + 2 * n_q_points);
+
+              // grad xx
+              eval0.template hessians<0, true, false, normal_dir>(values_dofs,
+                                                                  temp1);
+              eval1.template values<1, true, false, normal_dir>(temp1,
+                                                                hessians_quad);
+            }
+
+          // grad y
+          eval0.template values<0, true, false, normal_dir>(values_dofs, temp1);
+          if (evaluation_flag & EvaluationFlags::gradients)
+            eval1.template gradients<1, true, false, normal_dir>(
+              temp1, gradients_quad + n_q_points);
+
+          // grad yy
+          if (evaluation_flag & EvaluationFlags::hessians)
+            eval1.template hessians<1, true, false, normal_dir>(temp1,
+                                                                hessians_quad +
+                                                                  n_q_points);
+
+          // val: can use values applied in x
+          if (evaluation_flag & EvaluationFlags::values)
+            eval1.template values<1, true, false, normal_dir>(temp1,
+                                                              values_quad);
+          break;
+        case 3:
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              // grad x
+              eval0.template gradients<0, true, false, normal_dir>(values_dofs,
+                                                                   temp1);
+              eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+              eval2.template values<2, true, false, normal_dir>(temp2,
+                                                                gradients_quad);
+            }
+
+          if (evaluation_flag & EvaluationFlags::hessians)
+            {
+              // The evaluation/integration here *should* work, however
+              // the piola transform is not implemented.
+              AssertThrow(false, ExcNotImplemented());
+              // grad xz
+              if (!(evaluation_flag & EvaluationFlags::gradients))
+                {
+                  eval0.template gradients<0, true, false, normal_dir>(
+                    values_dofs, temp1);
+                  eval1.template values<1, true, false, normal_dir>(temp1,
+                                                                    temp2);
+                }
+              eval2.template gradients<2, true, false, normal_dir>(
+                temp2, hessians_quad + 4 * n_q_points);
+
+              // grad xy
+              eval1.template gradients<1, true, false, normal_dir>(temp1,
+                                                                   temp2);
+              eval2.template values<2, true, false, normal_dir>(
+                temp2, hessians_quad + 3 * n_q_points);
+
+              // grad xx
+              eval0.template hessians<0, true, false, normal_dir>(values_dofs,
+                                                                  temp1);
+              eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+              eval2.template values<2, true, false, normal_dir>(temp2,
+                                                                hessians_quad);
+            }
+
+          // grad y
+          eval0.template values<0, true, false, normal_dir>(values_dofs, temp1);
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              eval1.template gradients<1, true, false, normal_dir>(temp1,
+                                                                   temp2);
+              eval2.template values<2, true, false, normal_dir>(temp2,
+                                                                gradients_quad +
+                                                                  n_q_points);
+            }
+
+          if (evaluation_flag & EvaluationFlags::hessians)
+            {
+              // grad yz
+              if (!(evaluation_flag & EvaluationFlags::gradients))
+                eval1.template gradients<1, true, false, normal_dir>(temp1,
+                                                                     temp2);
+              eval2.template gradients<2, true, false, normal_dir>(
+                temp2, hessians_quad + 5 * n_q_points);
+
+              // grad yy
+              eval1.template hessians<1, true, false, normal_dir>(temp1, temp2);
+              eval2.template values<2, true, false, normal_dir>(temp2,
+                                                                hessians_quad +
+                                                                  n_q_points);
+            }
+
+          // grad z: can use the values applied in x direction stored in
+          // temp1
+          eval1.template values<1, true, false, normal_dir>(temp1, temp2);
+          if (evaluation_flag & EvaluationFlags::gradients)
+            eval2.template gradients<2, true, false, normal_dir>(
+              temp2, gradients_quad + 2 * n_q_points);
+
+          // grad zz: can use the values applied in x and y direction stored
+          // in temp2
+          if (evaluation_flag & EvaluationFlags::hessians)
+            eval2.template hessians<2, true, false, normal_dir>(
+              temp2, hessians_quad + 2 * n_q_points);
+
+          // val: can use the values applied in x & y direction stored in
+          // temp2
+          if (evaluation_flag & EvaluationFlags::values)
+            eval2.template values<2, true, false, normal_dir>(temp2,
+                                                              values_quad);
+          break;
+        default:
+          AssertThrow(false, ExcNotImplemented());
+      }
+  }
+
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  template <int normal_dir>
+  inline void
+  FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+                   dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   Number>::
+    evaluate_tensor_product_per_component(
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      Number *                               values_dofs_actual,
+      FEEvaluationData<dim, Number, false> & fe_eval,
+      const bool                             add_into_values_array,
+      std::integral_constant<bool, true>)
+  {
+    using Eval0 =
+      typename std::conditional<normal_dir == 0, EvalNormal, EvalTangent>::type;
+    using Eval1 =
+      typename std::conditional<normal_dir == 1, EvalNormal, EvalTangent>::type;
+    using Eval2 =
+      typename std::conditional<normal_dir == 2, EvalNormal, EvalTangent>::type;
+
+    Eval0 eval0 =
+      create_evaluator_tensor_product<Eval0>(fe_eval.get_shape_info());
+    Eval1 eval1 =
+      create_evaluator_tensor_product<Eval1>(fe_eval.get_shape_info());
+    Eval2 eval2 =
+      create_evaluator_tensor_product<Eval2>(fe_eval.get_shape_info());
+
+    Number *temp1 = fe_eval.get_scratch_data().begin();
+    Number *temp2;
+
+    temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+                               fe_eval.get_shape_info().data[0].fe_degree + 1),
+                             Utilities::fixed_power<dim>(
+                               fe_eval.get_shape_info().data[0].n_q_points_1d));
+
+    const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points;
+    const std::size_t dofs_per_comp =
+      fe_eval.get_shape_info().dofs_per_component_on_cell;
+
+    // Initial shift depending on component (normal_dir)
+    Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
+    Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
+    Number *gradients_quad =
+      fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
+    Number *hessians_quad =
+      (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
+                   fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
+
+    // Integrate path
+    switch (dim)
+      {
+        case 2:
+          if ((evaluation_flag & EvaluationFlags::values) &&
+              !(evaluation_flag & EvaluationFlags::gradients))
+            {
+              eval1.template values<1, false, false, normal_dir>(values_quad,
+                                                                 temp1);
+              if (add_into_values_array == false)
+                eval0.template values<0, false, false, normal_dir>(temp1,
+                                                                   values_dofs);
+              else
+                eval0.template values<0, false, true, normal_dir>(temp1,
+                                                                  values_dofs);
+            }
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              eval1.template gradients<1, false, false, normal_dir>(
+                gradients_quad + n_q_points, temp1);
+              if ((evaluation_flag & EvaluationFlags::values))
+                eval1.template values<1, false, true, normal_dir>(values_quad,
+                                                                  temp1);
+              if (add_into_values_array == false)
+                eval0.template values<0, false, false, normal_dir>(temp1,
+                                                                   values_dofs);
+              else
+                eval0.template values<0, false, true, normal_dir>(temp1,
+                                                                  values_dofs);
+              eval1.template values<1, false, false, normal_dir>(gradients_quad,
+                                                                 temp1);
+              eval0.template gradients<0, false, true, normal_dir>(temp1,
+                                                                   values_dofs);
+            }
+          if (evaluation_flag & EvaluationFlags::hessians)
+            {
+              // grad xx
+              eval1.template values<1, false, false, normal_dir>(hessians_quad,
+                                                                 temp1);
+
+              if ((evaluation_flag & EvaluationFlags::values) ||
+                  (evaluation_flag & EvaluationFlags::gradients) ||
+                  add_into_values_array == true)
+                eval0.template hessians<0, false, true, normal_dir>(
+                  temp1, values_dofs);
+              else
+                eval0.template hessians<0, false, false, normal_dir>(
+                  temp1, values_dofs);
+
+              // grad yy
+              eval1.template hessians<1, false, false, normal_dir>(
+                hessians_quad + n_q_points, temp1);
+              eval0.template values<0, false, true, normal_dir>(temp1,
+                                                                values_dofs);
+
+              // grad xy
+              eval1.template gradients<1, false, false, normal_dir>(
+                hessians_quad + 2 * n_q_points, temp1);
+              eval0.template gradients<0, false, true, normal_dir>(temp1,
+                                                                   values_dofs);
+            }
+          break;
+
+        case 3:
+          if ((evaluation_flag & EvaluationFlags::values) &&
+              !(evaluation_flag & EvaluationFlags::gradients))
+            {
+              eval2.template values<2, false, false, normal_dir>(values_quad,
+                                                                 temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              if (add_into_values_array == false)
+                eval0.template values<0, false, false, normal_dir>(temp2,
+                                                                   values_dofs);
+              else
+                eval0.template values<0, false, true, normal_dir>(temp2,
+                                                                  values_dofs);
+            }
+          if (evaluation_flag & EvaluationFlags::gradients)
+            {
+              eval2.template gradients<2, false, false, normal_dir>(
+                gradients_quad + 2 * n_q_points, temp1);
+              if ((evaluation_flag & EvaluationFlags::values))
+                eval2.template values<2, false, true, normal_dir>(values_quad,
+                                                                  temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              eval2.template values<2, false, false, normal_dir>(
+                gradients_quad + n_q_points, temp1);
+              eval1.template gradients<1, false, true, normal_dir>(temp1,
+                                                                   temp2);
+              if (add_into_values_array == false)
+                eval0.template values<0, false, false, normal_dir>(temp2,
+                                                                   values_dofs);
+              else
+                eval0.template values<0, false, true, normal_dir>(temp2,
+                                                                  values_dofs);
+              eval2.template values<2, false, false, normal_dir>(gradients_quad,
+                                                                 temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              eval0.template gradients<0, false, true, normal_dir>(temp2,
+                                                                   values_dofs);
+            }
+          if (evaluation_flag & EvaluationFlags::hessians)
+            {
+              // grad xx
+              eval2.template values<2, false, false, normal_dir>(hessians_quad,
+                                                                 temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+
+              if ((evaluation_flag & EvaluationFlags::values) ||
+                  (evaluation_flag & EvaluationFlags::gradients) ||
+                  add_into_values_array == true)
+                eval0.template hessians<0, false, true, normal_dir>(
+                  temp2, values_dofs);
+              else
+                eval0.template hessians<0, false, false, normal_dir>(
+                  temp2, values_dofs);
+
+              // grad yy
+              eval2.template values<2, false, false, normal_dir>(hessians_quad +
+                                                                   n_q_points,
+                                                                 temp1);
+              eval1.template hessians<1, false, false, normal_dir>(temp1,
+                                                                   temp2);
+              eval0.template values<0, false, true, normal_dir>(temp2,
+                                                                values_dofs);
+
+              // grad zz
+              eval2.template hessians<2, false, false, normal_dir>(
+                hessians_quad + 2 * n_q_points, temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              eval0.template values<0, false, true, normal_dir>(temp2,
+                                                                values_dofs);
+
+              // grad xy
+              eval2.template values<2, false, false, normal_dir>(
+                hessians_quad + 3 * n_q_points, temp1);
+              eval1.template gradients<1, false, false, normal_dir>(temp1,
+                                                                    temp2);
+              eval0.template gradients<0, false, true, normal_dir>(temp2,
+                                                                   values_dofs);
+
+              // grad xz
+              eval2.template gradients<2, false, false, normal_dir>(
+                hessians_quad + 4 * n_q_points, temp1);
+              eval1.template values<1, false, false, normal_dir>(temp1, temp2);
+              eval0.template gradients<0, false, true, normal_dir>(temp2,
+                                                                   values_dofs);
+
+              // grad yz
+              eval2.template gradients<2, false, false, normal_dir>(
+                hessians_quad + 5 * n_q_points, temp1);
+              eval1.template gradients<1, false, false, normal_dir>(temp1,
+                                                                    temp2);
+              eval0.template values<0, false, true, normal_dir>(temp2,
+                                                                values_dofs);
+            }
+
+          break;
+        default:
+          AssertThrow(false, ExcNotImplemented());
+      }
+  }
 
   /**
    * This struct implements the change between two different bases. This is an
@@ -1779,7 +2306,8 @@ namespace internal
 
       Assert(fe_eval.get_shape_info().data.size() == 1 ||
                (fe_eval.get_shape_info().data.size() == dim &&
-                element_type == ElementType::tensor_general),
+                element_type == ElementType::tensor_general) ||
+               element_type == ElementType::tensor_raviart_thomas,
              ExcNotImplemented());
 
       if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
@@ -1847,6 +2375,18 @@ namespace internal
                                              values_dofs,
                                              fe_eval);
         }
+      else if (element_type == ElementType::tensor_raviart_thomas)
+        {
+          FEEvaluationImpl<
+            ElementType::tensor_raviart_thomas,
+            dim,
+            (fe_degree == -1) ? 1 : fe_degree,
+            (n_q_points_1d < 1) ? 1 : n_q_points_1d,
+            Number>::template evaluate_or_integrate<false>(evaluation_flag,
+                                                           const_cast<Number *>(
+                                                             values_dofs),
+                                                           fe_eval);
+        }
       else
         {
           FEEvaluationImpl<ElementType::tensor_general,
@@ -1896,7 +2436,8 @@ namespace internal
 
       Assert(fe_eval.get_shape_info().data.size() == 1 ||
                (fe_eval.get_shape_info().data.size() == dim &&
-                element_type == ElementType::tensor_general),
+                element_type == ElementType::tensor_general) ||
+               element_type == ElementType::tensor_raviart_thomas,
              ExcNotImplemented());
 
       if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
@@ -1973,6 +2514,19 @@ namespace internal
                                               fe_eval,
                                               sum_into_values_array);
         }
+      else if (element_type == ElementType::tensor_raviart_thomas)
+        {
+          FEEvaluationImpl<ElementType::tensor_raviart_thomas,
+                           dim,
+                           (fe_degree == -1) ? 1 : fe_degree,
+                           (n_q_points_1d < 1) ? 1 : n_q_points_1d,
+                           Number>::
+            template evaluate_or_integrate<true>(integration_flag,
+                                                 const_cast<Number *>(
+                                                   values_dofs),
+                                                 fe_eval,
+                                                 sum_into_values_array);
+        }
       else
         {
           FEEvaluationImpl<ElementType::tensor_general,
@@ -2477,11 +3031,552 @@ namespace internal
     }
   };
 
+  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+  struct FEFaceEvaluationImplRaviartThomas
+  {
+    using EvalNormal =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim - 1,
+                             (fe_degree == -1) ? 1 : fe_degree + 1,
+                             n_q_points_1d,
+                             Number>;
+    using EvalTangent =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim - 1,
+                             (fe_degree == -1) ? 1 : fe_degree,
+                             n_q_points_1d,
+                             Number>;
+
+    using EvalGeneral = EvaluatorTensorProduct<evaluate_general,
+                                               dim - 1,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               Number>;
+    template <typename EvalType>
+    static EvalType
+    create_evaluator_tensor_product(
+      const MatrixFreeFunctions::UnivariateShapeData<Number> &data,
+      const unsigned int                                      subface_index,
+      const unsigned int                                      direction)
+    {
+      if (subface_index >= GeometryInfo<dim>::max_children_per_cell)
+        return EvalType(data.shape_values,
+                        data.shape_gradients,
+                        data.shape_hessians);
+      else
+        {
+          Assert(false, ExcNotImplemented());
+
+          const unsigned int index =
+            direction == 0 ? subface_index % 2 : subface_index / 2;
+          return EvalType(data.values_within_subface[index],
+                          data.gradients_within_subface[index],
+                          data.hessians_within_subface[index]);
+        }
+    }
+
+    template <bool integrate>
+    static void
+    evaluate_or_integrate_in_face(
+      const EvaluationFlags::EvaluationFlags        evaluation_flag,
+      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+      Number *                                      values_dofs,
+      FEEvaluationData<dim, Number, true> &         fe_eval,
+      Number *                                      scratch_data,
+      const unsigned int                            subface_index,
+      const unsigned int                            face_no)
+    {
+      // TODO. Make sure hanging nodes also are supported.
+      // The following part probably needs a rethink.
+      EvalNormal eval_normal =
+        create_evaluator_tensor_product<EvalNormal>(shape_info.data.front(),
+                                                    subface_index,
+                                                    0);
+      EvalTangent eval_tangent =
+        create_evaluator_tensor_product<EvalTangent>(shape_info.data.back(),
+                                                     subface_index,
+                                                     1);
+
+      // Used for normal faces which are isotropic
+      EvalGeneral eval_general =
+        create_evaluator_tensor_product<EvalGeneral>(shape_info.data.back(),
+                                                     subface_index,
+                                                     0);
+
+      // Note, n_dofs on tangent face
+      const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face;
+      const std::size_t n_dofs_normal =
+        n_dofs_tangent - Utilities::pow(fe_degree, dim - 2);
+
+      const unsigned int face_direction = face_no / 2;
+
+      if (face_direction == 0)
+        {
+          evaluate_in_face_apply<-1, 0>(
+            eval_general,
+            eval_general,
+            values_dofs,
+            fe_eval,
+            scratch_data,
+            evaluation_flag,
+            n_dofs_normal,
+            std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_normal;
+
+          evaluate_in_face_apply<0, 1>(
+            eval_normal,
+            eval_tangent,
+            values_dofs,
+            fe_eval,
+            scratch_data,
+            evaluation_flag,
+            n_dofs_tangent,
+            std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_tangent;
+
+          if (dim == 3)
+            {
+              evaluate_in_face_apply<1, 2>(
+                eval_tangent,
+                eval_normal,
+                values_dofs,
+                fe_eval,
+                scratch_data,
+                evaluation_flag,
+                n_dofs_tangent,
+                std::integral_constant<bool, integrate>());
+            }
+        }
+      else if (face_direction == 1)
+        {
+          // NOTE. Take zx-coordinates into account for dim == 3
+          if (dim == 3)
+            evaluate_in_face_apply<1, 0>(
+              eval_tangent,
+              eval_normal,
+              values_dofs,
+              fe_eval,
+              scratch_data,
+              evaluation_flag,
+              n_dofs_tangent,
+              std::integral_constant<bool, integrate>());
+          else
+            evaluate_in_face_apply<0, 0>(
+              eval_normal,
+              eval_tangent,
+              values_dofs,
+              fe_eval,
+              scratch_data,
+              evaluation_flag,
+              n_dofs_tangent,
+              std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_tangent;
+
+          evaluate_in_face_apply<-1, 1>(
+            eval_general,
+            eval_general,
+            values_dofs,
+            fe_eval,
+            scratch_data,
+            evaluation_flag,
+            n_dofs_normal,
+            std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_normal;
+
+          if (dim == 3)
+            {
+              // NOTE. Take zx-coordinates into account
+              evaluate_in_face_apply<0, 2>(
+                eval_normal,
+                eval_tangent,
+                values_dofs,
+                fe_eval,
+                scratch_data,
+                evaluation_flag,
+                n_dofs_tangent,
+                std::integral_constant<bool, integrate>());
+            }
+        }
+      else
+        {
+          evaluate_in_face_apply<0, 0>(
+            eval_normal,
+            eval_tangent,
+            values_dofs,
+            fe_eval,
+            scratch_data,
+            evaluation_flag,
+            n_dofs_tangent,
+            std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_tangent;
+
+          evaluate_in_face_apply<1, 1>(
+            eval_tangent,
+            eval_normal,
+            values_dofs,
+            fe_eval,
+            scratch_data,
+            evaluation_flag,
+            n_dofs_tangent,
+            std::integral_constant<bool, integrate>());
+
+          values_dofs += 3 * n_dofs_tangent;
+
+          if (dim == 3)
+            {
+              evaluate_in_face_apply<-1, 2>(
+                eval_general,
+                eval_general,
+                values_dofs,
+                fe_eval,
+                scratch_data,
+                evaluation_flag,
+                n_dofs_normal,
+                std::integral_constant<bool, integrate>());
+            }
+        }
+    }
+
+    /*
+     * Helper function which applies the 1D kernels for on one
+     * component in a face. normal_dir indicates the direction of the continuous
+     * component of the RT space. std::integral_constant<bool, false> is the
+     * evaluation path, and std::integral_constant<bool, true> bellow is the
+     * integration path.
+     */
+    template <int normal_dir, int component, typename Eval0, typename Eval1>
+    static inline void
+    evaluate_in_face_apply(
+      const Eval0 &                          eval0,
+      const Eval1 &                          eval1,
+      Number *                               values_dofs,
+      FEEvaluationData<dim, Number, true> &  fe_eval,
+      Number *                               scratch_data,
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      const std::size_t                      dofs_stride,
+      std::integral_constant<bool, false>)
+    {
+      constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+
+      Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+      Number *gradients_quad =
+        fe_eval.begin_gradients() + dim * n_q_points * component;
+      Number *hessians_quad =
+        fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+
+      // Evaluation path
+
+      if ((evaluation_flag & EvaluationFlags::values) &&
+          !(evaluation_flag & EvaluationFlags::gradients))
+        {
+          switch (dim)
+            {
+              case 3:
+                eval0.template values<0, true, false, normal_dir>(values_dofs,
+                                                                  values_quad);
+                eval1.template values<1, true, false, normal_dir>(values_quad,
+                                                                  values_quad);
+                break;
+              case 2:
+                eval0.template values<0, true, false, normal_dir>(values_dofs,
+                                                                  values_quad);
+                break;
+              default:
+                Assert(false, ExcNotImplemented());
+            }
+        }
+      else if (evaluation_flag & EvaluationFlags::gradients)
+        {
+          switch (dim)
+            {
+              case 3:
+                // grad x
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs, scratch_data);
+                eval1.template values<1, true, false, normal_dir>(
+                  scratch_data, gradients_quad);
+
+                // grad y
+                eval0.template values<0, true, false, normal_dir>(values_dofs,
+                                                                  scratch_data);
+                eval1.template gradients<1, true, false, normal_dir>(
+                  scratch_data, gradients_quad + n_q_points);
+
+                if (evaluation_flag & EvaluationFlags::values)
+                  eval1.template values<1, true, false, normal_dir>(
+                    scratch_data, values_quad);
+
+                // grad z
+                eval0.template values<0, true, false, normal_dir>(values_dofs +
+                                                                    dofs_stride,
+                                                                  scratch_data);
+                eval1.template values<1, true, false, normal_dir>(
+                  scratch_data, gradients_quad + 2 * n_q_points);
+
+                break;
+              case 2:
+                eval0.template values<0, true, false, normal_dir>(
+                  values_dofs + dofs_stride, gradients_quad + n_q_points);
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs, gradients_quad);
+                if ((evaluation_flag & EvaluationFlags::values))
+                  eval0.template values<0, true, false, normal_dir>(
+                    values_dofs, values_quad);
+                break;
+              default:
+                AssertThrow(false, ExcNotImplemented());
+            }
+        }
+
+      if (evaluation_flag & EvaluationFlags::hessians)
+        {
+          switch (dim)
+            {
+              case 3:
+                // grad xx
+                eval0.template hessians<0, true, false, normal_dir>(
+                  values_dofs, scratch_data);
+                eval1.template values<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad);
+
+                // grad yy
+                eval0.template values<0, true, false, normal_dir>(values_dofs,
+                                                                  scratch_data);
+                eval1.template hessians<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad + n_q_points);
+
+                // grad zz
+                eval0.template values<0, true, false, normal_dir>(
+                  values_dofs + 2 * dofs_stride, scratch_data);
+                eval1.template values<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad + 2 * n_q_points);
+
+                // grad xy
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs, scratch_data);
+                eval1.template gradients<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad + 3 * n_q_points);
+
+                // grad xz
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs + dofs_stride, scratch_data);
+                eval1.template values<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad + 4 * n_q_points);
+
+                // grad yz
+                eval0.template values<0, true, false, normal_dir>(values_dofs +
+                                                                    dofs_stride,
+                                                                  scratch_data);
+                eval1.template gradients<1, true, false, normal_dir>(
+                  scratch_data, hessians_quad + 5 * n_q_points);
+
+                break;
+              case 2:
+                // grad xx
+                eval0.template hessians<0, true, false, normal_dir>(
+                  values_dofs, hessians_quad);
+                // grad yy
+                eval0.template values<0, true, false, normal_dir>(
+                  values_dofs + 2 * dofs_stride, hessians_quad + n_q_points);
+                // grad xy
+                eval0.template gradients<0, true, false, normal_dir>(
+                  values_dofs + dofs_stride, hessians_quad + 2 * n_q_points);
+                break;
+              default:
+                AssertThrow(false, ExcNotImplemented());
+            }
+        }
+    }
+
+    template <int normal_dir, int component, typename Eval0, typename Eval1>
+    static inline void
+    evaluate_in_face_apply(
+      const Eval0 &                          eval0,
+      const Eval1 &                          eval1,
+      Number *                               values_dofs,
+      FEEvaluationData<dim, Number, true> &  fe_eval,
+      Number *                               scratch_data,
+      const EvaluationFlags::EvaluationFlags evaluation_flag,
+      const std::size_t                      dofs_stride,
+      std::integral_constant<bool, true>)
+    {
+      constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1);
+
+      Number *values_quad = fe_eval.begin_values() + n_q_points * component;
+      Number *gradients_quad =
+        fe_eval.begin_gradients() + dim * n_q_points * component;
+      Number *hessians_quad =
+        fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component;
+
+      // Integration path
+      if ((evaluation_flag & EvaluationFlags::values) &&
+          !(evaluation_flag & EvaluationFlags::gradients))
+        {
+          switch (dim)
+            {
+              case 3:
+                eval1.template values<1, false, false, normal_dir>(values_quad,
+                                                                   values_quad);
+                eval0.template values<0, false, false, normal_dir>(values_quad,
+                                                                   values_dofs);
+                break;
+              case 2:
+                eval0.template values<0, false, false, normal_dir>(values_quad,
+                                                                   values_dofs);
+                break;
+              default:
+                Assert(false, ExcNotImplemented());
+            }
+        }
+      else if (evaluation_flag & EvaluationFlags::gradients)
+        {
+          switch (dim)
+            {
+              case 3:
+                // grad z
+                eval1.template values<1, false, false, normal_dir>(
+                  gradients_quad + 2 * n_q_points,
+                  gradients_quad + 2 * n_q_points);
+                eval0.template values<0, false, false, normal_dir>(
+                  gradients_quad + 2 * n_q_points, values_dofs + dofs_stride);
+
+                if (evaluation_flag & EvaluationFlags::values)
+                  {
+                    eval1.template values<1, false, false, normal_dir>(
+                      values_quad, scratch_data);
+                    eval1.template gradients<1, false, true, normal_dir>(
+                      gradients_quad + n_q_points, scratch_data);
+                  }
+                else
+                  eval1.template gradients<1, false, false, normal_dir>(
+                    gradients_quad + n_q_points, scratch_data);
+
+                // grad y
+                eval0.template values<0, false, false, normal_dir>(scratch_data,
+                                                                   values_dofs);
+
+                // grad x
+                eval1.template values<1, false, false, normal_dir>(
+                  gradients_quad, scratch_data);
+                eval0.template gradients<0, false, true, normal_dir>(
+                  scratch_data, values_dofs);
+
+                break;
+              case 2:
+                eval0.template values<0, false, false, normal_dir>(
+                  gradients_quad + n_q_points, values_dofs + dofs_stride);
+                eval0.template gradients<0, false, false, normal_dir>(
+                  gradients_quad, values_dofs);
+                if (evaluation_flag & EvaluationFlags::values)
+                  eval0.template values<0, false, true, normal_dir>(
+                    values_quad, values_dofs);
+                break;
+              default:
+                AssertThrow(false, ExcNotImplemented());
+            }
+        }
+
+      if (evaluation_flag & EvaluationFlags::hessians)
+        {
+          switch (dim)
+            {
+              case 3:
+                // grad xx
+                eval1.template values<1, false, false, normal_dir>(
+                  hessians_quad, scratch_data);
+                if ((evaluation_flag &
+                     (EvaluationFlags::values | EvaluationFlags::gradients)))
+                  eval0.template hessians<0, false, true, normal_dir>(
+                    scratch_data, values_dofs);
+                else
+                  eval0.template hessians<0, false, false, normal_dir>(
+                    scratch_data, values_dofs);
+
+                // grad yy
+                eval1.template hessians<1, false, false, normal_dir>(
+                  hessians_quad + n_q_points, scratch_data);
+                eval0.template values<0, false, true, normal_dir>(scratch_data,
+                                                                  values_dofs);
+
+                // grad zz
+                eval1.template values<1, false, false, normal_dir>(
+                  hessians_quad + 2 * n_q_points, scratch_data);
+                eval0.template values<0, false, false, normal_dir>(
+                  scratch_data, values_dofs + 2 * dofs_stride);
+
+                // grad xy
+                eval1.template gradients<1, false, false, normal_dir>(
+                  hessians_quad + 3 * n_q_points, scratch_data);
+                eval0.template gradients<0, false, true, normal_dir>(
+                  scratch_data, values_dofs);
+
+                // grad xz
+                eval1.template values<1, false, false, normal_dir>(
+                  hessians_quad + 4 * n_q_points, scratch_data);
+                if ((evaluation_flag & EvaluationFlags::gradients))
+                  eval0.template gradients<0, false, true, normal_dir>(
+                    scratch_data, values_dofs + dofs_stride);
+                else
+                  eval0.template gradients<0, false, false, normal_dir>(
+                    scratch_data, values_dofs + dofs_stride);
+
+                // grad yz
+                eval1.template gradients<1, false, false, normal_dir>(
+                  hessians_quad + 5 * n_q_points, scratch_data);
+                eval0.template values<0, false, true, normal_dir>(
+                  scratch_data, values_dofs + dofs_stride);
+
+                break;
+              case 2:
+                // grad xx
+                if (evaluation_flag &
+                    (EvaluationFlags::values | EvaluationFlags::gradients))
+                  eval0.template hessians<0, false, true, normal_dir>(
+                    hessians_quad, values_dofs);
+                else
+                  eval0.template hessians<0, false, false, normal_dir>(
+                    hessians_quad, values_dofs);
+
+                // grad yy
+                eval0.template values<0, false, false, normal_dir>(
+                  hessians_quad + n_q_points, values_dofs + 2 * dofs_stride);
+                // grad xy
+                if ((evaluation_flag & EvaluationFlags::gradients))
+                  eval0.template gradients<0, false, true, normal_dir>(
+                    hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
+                else
+                  eval0.template gradients<0, false, false, normal_dir>(
+                    hessians_quad + 2 * n_q_points, values_dofs + dofs_stride);
+                break;
+              default:
+                AssertThrow(false, ExcNotImplemented());
+            }
+        }
+    }
+  };
 
 
   template <int dim, int fe_degree, typename Number, bool lex_faces = false>
   struct FEFaceNormalEvaluationImpl
   {
+    using EvalNormal =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim,
+                             (fe_degree == -1) ? 1 : fe_degree + 1,
+                             0,
+                             Number>;
+    using EvalTangent =
+      EvaluatorTensorProduct<evaluate_raviart_thomas,
+                             dim,
+                             (fe_degree == -1) ? 1 : fe_degree,
+                             0,
+                             Number>;
+
     template <bool do_evaluate, bool add_into_output>
     static void
     interpolate(const unsigned int                            n_components,
@@ -2495,17 +3590,20 @@ namespace internal
                  shape_info.data.front().fe_degree ||
                fe_degree == -1,
              ExcInternalError());
-
-      interpolate_generic<do_evaluate, add_into_output>(
-        n_components,
-        input,
-        output,
-        flags,
-        face_no,
-        shape_info.data.front().fe_degree + 1,
-        shape_info.data.front().shape_data_on_face,
-        shape_info.dofs_per_component_on_cell,
-        3 * shape_info.dofs_per_component_on_face);
+      if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        interpolate_generic_raviart_thomas<do_evaluate, add_into_output>(
+          n_components, input, output, flags, face_no, shape_info);
+      else
+        interpolate_generic<do_evaluate, add_into_output>(
+          n_components,
+          input,
+          output,
+          flags,
+          face_no,
+          shape_info.data.front().fe_degree + 1,
+          shape_info.data.front().shape_data_on_face,
+          shape_info.dofs_per_component_on_cell,
+          3 * shape_info.dofs_per_component_on_face);
     }
 
     /**
@@ -2611,6 +3709,168 @@ namespace internal
             dofs_per_component_on_face);
         }
     }
+
+  private:
+    template <typename EvalType>
+    static EvalType
+    create_evaluator_tensor_product(
+      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+      const unsigned int                            face_no)
+    {
+      // Shape_data is of size 2 with normal direction in [0] and tangential
+      // direction in [1]. See shape_info.template.h
+      const MatrixFreeFunctions::UnivariateShapeData<Number> *shape_data;
+      if (std::is_same<EvalType, EvalNormal>::value)
+        {
+          shape_data = &shape_info.data[0];
+        }
+      else
+        {
+          shape_data = &shape_info.data[1];
+        }
+      return EvalType(shape_data->shape_data_on_face[face_no % 2],
+                      AlignedVector<Number>(),
+                      AlignedVector<Number>());
+    }
+
+    template <bool do_evaluate,
+              bool add_into_output,
+              int  face_direction = 0,
+              int  max_derivative = 0>
+    static void
+    interpolate_generic_raviart_thomas(
+      const unsigned int                            n_components,
+      const Number *                                input,
+      Number *                                      output,
+      const EvaluationFlags::EvaluationFlags        flag,
+      const unsigned int                            face_no,
+      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info)
+    {
+      bool increase_max_der = false;
+      if ((flag & EvaluationFlags::hessians && max_derivative < 2) ||
+          (flag & EvaluationFlags::gradients && max_derivative < 1))
+        increase_max_der = true;
+
+      if (face_direction == face_no / 2 && !increase_max_der)
+        {
+          interpolate_generic_raviart_thomas_apply_face<do_evaluate,
+                                                        add_into_output,
+                                                        face_direction,
+                                                        max_derivative>(
+            shape_info, face_no, input, output);
+        }
+      else if (face_direction == face_no / 2)
+        {
+          // Only increase max_derivative
+          interpolate_generic_raviart_thomas<do_evaluate,
+                                             add_into_output,
+                                             face_direction,
+                                             std::min(max_derivative + 1, 2)>(
+            n_components, input, output, flag, face_no, shape_info);
+        }
+      else if (face_direction < dim)
+        {
+          if (increase_max_der)
+            {
+              interpolate_generic_raviart_thomas<
+                do_evaluate,
+                add_into_output,
+                std::min(face_direction + 1, dim - 1),
+                std::min(max_derivative + 1, 2)>(
+                n_components, input, output, flag, face_no, shape_info);
+            }
+          else
+            {
+              interpolate_generic_raviart_thomas<do_evaluate,
+                                                 add_into_output,
+                                                 std::min(face_direction + 1,
+                                                          dim - 1),
+                                                 max_derivative>(
+                n_components, input, output, flag, face_no, shape_info);
+            }
+        }
+    }
+
+    /* Help function for interpolate_generic_raviart_thomas */
+    template <bool do_evaluate,
+              bool add_into_output,
+              int  face_direction,
+              int  max_derivative>
+    static inline void
+    interpolate_generic_raviart_thomas_apply_face(
+      const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+      const unsigned int                            face_no,
+      const Number *                                input,
+      Number *                                      output)
+    {
+      using Evalf0 = typename std::
+        conditional<face_direction == 0, EvalNormal, EvalTangent>::type;
+      using Evalf1 = typename std::
+        conditional<face_direction == 1, EvalNormal, EvalTangent>::type;
+      using Evalf2 = typename std::
+        conditional<face_direction == 2, EvalNormal, EvalTangent>::type;
+
+      Evalf0 evalf0 =
+        create_evaluator_tensor_product<Evalf0>(shape_info, face_no);
+      Evalf1 evalf1 =
+        create_evaluator_tensor_product<Evalf1>(shape_info, face_no);
+      Evalf2 evalf2 =
+        create_evaluator_tensor_product<Evalf2>(shape_info, face_no);
+
+      const unsigned int dofs_per_component_on_cell =
+        shape_info.dofs_per_component_on_cell;
+      const unsigned int dofs_per_component_on_face =
+        3 * shape_info.dofs_per_component_on_face;
+
+      // NOTE! dofs_per_component_on_face is in the tangent direction,
+      // i.e (fe.degree+1)*fe.degree. Normal faces are only
+      // fe.degree*fe.degree
+      const unsigned int in_stride =
+        do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face;
+      const unsigned int out_stride =
+        do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell;
+
+      const unsigned int in_stride_after_normal =
+        do_evaluate ?
+          dofs_per_component_on_cell :
+          dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2);
+      const unsigned int out_stride_after_normal =
+        do_evaluate ?
+          dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) :
+          dofs_per_component_on_cell;
+
+      evalf0.template apply_face<face_direction,
+                                 do_evaluate,
+                                 add_into_output,
+                                 max_derivative,
+                                 lex_faces,
+                                 0>(input, output);
+      // stride to next component
+      input += (face_direction == 0) ? in_stride_after_normal : in_stride;
+      output += (face_direction == 0) ? out_stride_after_normal : out_stride;
+
+      evalf1.template apply_face<face_direction,
+                                 do_evaluate,
+                                 add_into_output,
+                                 max_derivative,
+                                 lex_faces,
+                                 1>(input, output);
+
+      if (dim == 3)
+        {
+          // stride to next component
+          input += (face_direction == 1) ? in_stride_after_normal : in_stride;
+          output +=
+            (face_direction == 1) ? out_stride_after_normal : out_stride;
+
+          evalf2.template apply_face<face_direction,
+                                     do_evaluate,
+                                     add_into_output,
+                                     max_derivative,
+                                     lex_faces,
+                                     2>(input, output);
+        }
+    }
   };
 
 
@@ -2949,9 +4209,25 @@ namespace internal
       constexpr unsigned int n_q_points_1d_actual =
         fe_degree > -1 ? n_q_points_1d : 0;
 
-      if (fe_degree > -1 &&
-          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-          shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+      if (fe_degree >= 1 &&
+          shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        {
+          FEFaceEvaluationImplRaviartThomas<dim,
+                                            (fe_degree == -1) ? 1 : fe_degree,
+                                            (n_q_points_1d < 1) ? 1 :
+                                                                  n_q_points_1d,
+                                            Number>::
+            template evaluate_or_integrate_in_face<false>(evaluation_flag,
+                                                          shape_info,
+                                                          temp,
+                                                          fe_eval,
+                                                          scratch_data,
+                                                          subface_index,
+                                                          face_no);
+        }
+      else if (fe_degree > -1 &&
+               subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
         FEFaceEvaluationImpl<true,
                              dim,
                              fe_degree,
@@ -3107,10 +4383,26 @@ namespace internal
         fe_degree > -1 ? n_q_points_1d : 0;
       const unsigned int subface_index = fe_eval.get_subface_index();
 
-      if (fe_degree > -1 &&
-          fe_eval.get_subface_index() >=
-            GeometryInfo<dim - 1>::max_children_per_cell &&
-          shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
+      if (fe_degree >= 1 &&
+          shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas)
+        {
+          FEFaceEvaluationImplRaviartThomas<dim,
+                                            (fe_degree == -1) ? 1 : fe_degree,
+                                            (n_q_points_1d < 1) ? 1 :
+                                                                  n_q_points_1d,
+                                            Number>::
+            template evaluate_or_integrate_in_face<true>(integration_flag,
+                                                         shape_info,
+                                                         temp,
+                                                         fe_eval,
+                                                         scratch_data,
+                                                         subface_index,
+                                                         face_no);
+        }
+      else if (fe_degree > -1 &&
+               fe_eval.get_subface_index() >=
+                 GeometryInfo<dim - 1>::max_children_per_cell &&
+               shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric)
         FEFaceEvaluationImpl<
           true,
           dim,
diff --git a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h
new file mode 100644 (file)
index 0000000..fa06896
--- /dev/null
@@ -0,0 +1,125 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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_matrix_free_evaluation_template_factory_templates_h
+#define dealii_matrix_free_evaluation_template_factory_templates_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/evaluation_selector.h>
+#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/evaluation_template_factory_internal.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  template <int dim, typename Number>
+  void
+  FEFaceEvaluationFactory<dim, Number>::evaluate(
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags evaluation_flag,
+    const Number *                         values_dofs,
+    FEEvaluationData<dim, Number, true> &  fe_eval)
+  {
+    instantiation_helper_run<1,
+                             FEFaceEvaluationImplEvaluateSelector<dim, Number>>(
+      fe_eval.get_shape_info().data[0].fe_degree,
+      fe_eval.get_shape_info().data[0].n_q_points_1d,
+      n_components,
+      evaluation_flag,
+      values_dofs,
+      fe_eval);
+  }
+
+
+
+  template <int dim, typename Number>
+  void
+  FEFaceEvaluationFactory<dim, Number>::integrate(
+    const unsigned int                     n_components,
+    const EvaluationFlags::EvaluationFlags integration_flag,
+    Number *                               values_dofs,
+    FEEvaluationData<dim, Number, true> &  fe_eval)
+  {
+    instantiation_helper_run<
+      1,
+      FEFaceEvaluationImplIntegrateSelector<dim, Number>>(
+      fe_eval.get_shape_info().data[0].fe_degree,
+      fe_eval.get_shape_info().data[0].n_q_points_1d,
+      n_components,
+      integration_flag,
+      values_dofs,
+      fe_eval);
+  }
+
+
+
+  template <int dim, typename Number, typename VectorizedArrayType>
+  void
+  FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::evaluate(
+    const unsigned int                                n_components,
+    const EvaluationFlags::EvaluationFlags            evaluation_flag,
+    const Number *                                    src_ptr,
+    const std::vector<ArrayView<const Number>> *      sm_ptr,
+    FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+  {
+    instantiation_helper_run<
+      1,
+      FEFaceEvaluationImplGatherEvaluateSelector<dim,
+                                                 Number,
+                                                 VectorizedArrayType>>(
+      fe_eval.get_shape_info().data[0].fe_degree,
+      fe_eval.get_shape_info().data[0].n_q_points_1d,
+      n_components,
+      evaluation_flag,
+      src_ptr,
+      sm_ptr,
+      fe_eval);
+  }
+
+
+
+  template <int dim, typename Number, typename VectorizedArrayType>
+  void
+  FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::integrate(
+    const unsigned int                                n_components,
+    const EvaluationFlags::EvaluationFlags            integration_flag,
+    Number *                                          dst_ptr,
+    const std::vector<ArrayView<const Number>> *      sm_ptr,
+    FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
+  {
+    instantiation_helper_run<
+      1,
+      FEFaceEvaluationImplIntegrateScatterSelector<dim,
+                                                   Number,
+                                                   VectorizedArrayType>>(
+      fe_eval.get_shape_info().data[0].fe_degree,
+      fe_eval.get_shape_info().data[0].n_q_points_1d,
+      n_components,
+      integration_flag,
+      dst_ptr,
+      sm_ptr,
+      fe_eval);
+  }
+
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index cdaca44458f197220c9d40615a2a016c2183d7e2..d70ba5758ded07c70a831bed311b7eb88c6cd754 100644 (file)
@@ -93,97 +93,6 @@ namespace internal
 
 
 
-  template <int dim, typename Number>
-  void
-  FEFaceEvaluationFactory<dim, Number>::evaluate(
-    const unsigned int                     n_components,
-    const EvaluationFlags::EvaluationFlags evaluation_flag,
-    const Number *                         values_dofs,
-    FEEvaluationData<dim, Number, true> &  fe_eval)
-  {
-    instantiation_helper_run<1,
-                             FEFaceEvaluationImplEvaluateSelector<dim, Number>>(
-      fe_eval.get_shape_info().data[0].fe_degree,
-      fe_eval.get_shape_info().data[0].n_q_points_1d,
-      n_components,
-      evaluation_flag,
-      values_dofs,
-      fe_eval);
-  }
-
-
-
-  template <int dim, typename Number>
-  void
-  FEFaceEvaluationFactory<dim, Number>::integrate(
-    const unsigned int                     n_components,
-    const EvaluationFlags::EvaluationFlags integration_flag,
-    Number *                               values_dofs,
-    FEEvaluationData<dim, Number, true> &  fe_eval)
-  {
-    instantiation_helper_run<
-      1,
-      FEFaceEvaluationImplIntegrateSelector<dim, Number>>(
-      fe_eval.get_shape_info().data[0].fe_degree,
-      fe_eval.get_shape_info().data[0].n_q_points_1d,
-      n_components,
-      integration_flag,
-      values_dofs,
-      fe_eval);
-  }
-
-
-
-  template <int dim, typename Number, typename VectorizedArrayType>
-  void
-  FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::evaluate(
-    const unsigned int                                n_components,
-    const EvaluationFlags::EvaluationFlags            evaluation_flag,
-    const Number *                                    src_ptr,
-    const std::vector<ArrayView<const Number>> *      sm_ptr,
-    FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
-  {
-    instantiation_helper_run<
-      1,
-      FEFaceEvaluationImplGatherEvaluateSelector<dim,
-                                                 Number,
-                                                 VectorizedArrayType>>(
-      fe_eval.get_shape_info().data[0].fe_degree,
-      fe_eval.get_shape_info().data[0].n_q_points_1d,
-      n_components,
-      evaluation_flag,
-      src_ptr,
-      sm_ptr,
-      fe_eval);
-  }
-
-
-
-  template <int dim, typename Number, typename VectorizedArrayType>
-  void
-  FEFaceEvaluationGatherFactory<dim, Number, VectorizedArrayType>::integrate(
-    const unsigned int                                n_components,
-    const EvaluationFlags::EvaluationFlags            integration_flag,
-    Number *                                          dst_ptr,
-    const std::vector<ArrayView<const Number>> *      sm_ptr,
-    FEEvaluationData<dim, VectorizedArrayType, true> &fe_eval)
-  {
-    instantiation_helper_run<
-      1,
-      FEFaceEvaluationImplIntegrateScatterSelector<dim,
-                                                   Number,
-                                                   VectorizedArrayType>>(
-      fe_eval.get_shape_info().data[0].fe_degree,
-      fe_eval.get_shape_info().data[0].n_q_points_1d,
-      n_components,
-      integration_flag,
-      dst_ptr,
-      sm_ptr,
-      fe_eval);
-  }
-
-
-
   template <int dim, typename Number>
   void
   CellwiseInverseMassFactory<dim, Number>::apply(
index f406690258fad0aa0cab8c01117bd61c6ec9e669..4f1ce43dd94c9928492280dbae575f667ba0cc31 100644 (file)
@@ -1160,16 +1160,21 @@ namespace internal
             dof_info[no].fe_index_conversion[fe_index].clear();
             for (unsigned int c = 0; c < dof_info[no].n_base_elements; ++c)
               {
-                dof_info[no].n_components[c] = fe.element_multiplicity(c);
+                dof_info[no].n_components[c] =
+                  fe.element_multiplicity(c) *
+                  fe.base_element(c).n_components();
                 for (unsigned int l = 0; l < dof_info[no].n_components[c]; ++l)
                   {
                     dof_info[no].component_to_base_index.push_back(c);
                     dof_info[no]
                       .component_dof_indices_offset[fe_index]
-                      .push_back(dof_info[no]
-                                   .component_dof_indices_offset[fe_index]
-                                   .back() +
-                                 fe.base_element(c).n_dofs_per_cell());
+                      .push_back(
+                        dof_info[no]
+                          .component_dof_indices_offset[fe_index]
+                          .back() +
+                        shape_infos(dof_info[no].global_base_element_offset + c,
+                                    fe_index)
+                          .dofs_per_component_on_cell);
                     dof_info[no].fe_index_conversion[fe_index].push_back(
                       fe.base_element(c).degree);
                   }
index 1187d0708c3359224d955599ece5b9f73c0f1bf0..ee7394ee587526ab4048245cad96944d5358aafe 100644 (file)
@@ -100,9 +100,18 @@ namespace internal
       tensor_symmetric_plus_dg0 = 5,
 
       /**
-       * Shape functions without an tensor product properties.
+       * Special case of the FE_RaviartThomasNodal element with anisotropic
+       * tensor product shape functions, i.e degree (k + 1) in normal direction,
+       * and k in tangential direction.
        */
-      tensor_none = 6
+      tensor_raviart_thomas = 6,
+
+      /**
+       * Shape functions without a tensor product properties.
+       */
+      tensor_none = 7
+
+
     };
 
 
index 9f1a084a81b8a6d8e942e4be31971571c1962c43..63a2d6e92a30fc6be9379478c57ed4b8d3c239e3 100644 (file)
@@ -33,6 +33,7 @@
 #include <deal.II/fe/fe_pyramid_p.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
 #include <deal.II/fe/fe_simplex_p.h>
 #include <deal.II/fe/fe_simplex_p_bubbles.h>
 #include <deal.II/fe/fe_tools.h>
@@ -170,8 +171,6 @@ namespace internal
         }
     }
 
-
-
     template <int dim_to, int dim, int spacedim>
     std::unique_ptr<FiniteElement<dim_to, dim_to>>
     create_fe(const FiniteElement<dim, spacedim> &fe)
@@ -194,7 +193,6 @@ namespace internal
     }
 
 
-
     // ----------------- actual ShapeInfo implementation --------------------
 
     template <typename Number>
@@ -238,15 +236,171 @@ namespace internal
     {
       static_assert(dim == spacedim,
                     "Currently, only the case dim=spacedim is implemented");
-      if (quad_in.is_tensor_product() == false ||
-          dynamic_cast<const FE_SimplexP<dim> *>(
-            &fe_in.base_element(base_element_number)) ||
-          dynamic_cast<const FE_SimplexDGP<dim> *>(
-            &fe_in.base_element(base_element_number)) ||
-          dynamic_cast<const FE_WedgeP<dim> *>(
-            &fe_in.base_element(base_element_number)) ||
-          dynamic_cast<const FE_PyramidP<dim> *>(
+
+      // ShapeInfo for RT elements. Here, data is of size 2 instead of 1.
+      // data[0] is univariate_shape_data in normal direction and
+      // data[1] is univariate_shape_data in tangential direction
+      //
+      if (dynamic_cast<const FE_RaviartThomasNodal<dim> *>(
             &fe_in.base_element(base_element_number)))
+        {
+          element_type = tensor_raviart_thomas;
+
+          const auto quad = quad_in.get_tensor_basis()[0];
+
+          const FiniteElement<dim> &fe =
+            fe_in.base_element(base_element_number);
+          n_dimensions = dim;
+          n_components = fe_in.n_components();
+
+          data.resize(2);
+          const unsigned int n_q_points_1d = quad.size();
+
+          n_q_points      = Utilities::fixed_power<dim>(n_q_points_1d);
+          n_q_points_face = Utilities::fixed_power<dim - 1>(n_q_points_1d);
+
+          dofs_per_component_on_cell = fe_in.n_dofs_per_cell() / n_components;
+
+          // NOTE dofs_per_component_on_face is in tangential direction!
+          dofs_per_component_on_face =
+            fe_in.n_dofs_per_face() + Utilities::pow(fe_in.degree, dim - 2);
+          const unsigned int dofs_per_face_normal = fe_in.n_dofs_per_face();
+
+          lexicographic_numbering =
+            get_lexicographic_numbering_rt_nodal<dim>(fe_in.degree);
+
+          // To get the right shape_values of the RT element
+          std::vector<unsigned int> lex_normal, lex_tangent;
+          for (unsigned int i = 0; i < fe.degree; ++i)
+            lex_tangent.push_back(i);
+
+          lex_normal.push_back(0);
+          for (unsigned int i = dofs_per_face_normal * 2 * dim;
+               i < dofs_per_face_normal * 2 * dim + fe.degree - 1;
+               ++i)
+            lex_normal.push_back(i);
+          lex_normal.push_back(dofs_per_face_normal);
+
+          // 'direction' distingusishes between normal and tangential direction
+          for (unsigned int direction = 0; direction < 2; ++direction)
+            {
+              UnivariateShapeData<Number> &univariate_shape_data =
+                (direction == 0) ? data.front() : data.back();
+
+              univariate_shape_data.element_type  = tensor_raviart_thomas;
+              univariate_shape_data.quadrature    = quad;
+              univariate_shape_data.n_q_points_1d = n_q_points_1d;
+              univariate_shape_data.fe_degree     = fe.degree - direction;
+
+              // grant write access to common univariate shape data
+              auto &shape_values    = univariate_shape_data.shape_values;
+              auto &shape_gradients = univariate_shape_data.shape_gradients;
+              auto &shape_hessians  = univariate_shape_data.shape_hessians;
+
+              auto &values_within_subface =
+                univariate_shape_data.values_within_subface;
+              auto &gradients_within_subface =
+                univariate_shape_data.gradients_within_subface;
+              auto &hessians_within_subface =
+                univariate_shape_data.hessians_within_subface;
+
+              auto &shape_data_on_face =
+                univariate_shape_data.shape_data_on_face;
+
+              const unsigned int n_dofs_1d  = fe.degree + 1 - direction;
+              const unsigned int array_size = n_dofs_1d * n_q_points_1d;
+
+              shape_gradients.resize_fast(array_size);
+              shape_values.resize_fast(array_size);
+              shape_hessians.resize_fast(array_size);
+
+              values_within_subface[0].resize(array_size);
+              values_within_subface[1].resize(array_size);
+              gradients_within_subface[0].resize(array_size);
+              gradients_within_subface[1].resize(array_size);
+              hessians_within_subface[0].resize(array_size);
+              hessians_within_subface[1].resize(array_size);
+
+              shape_data_on_face[0].resize(3 * n_dofs_1d);
+              shape_data_on_face[1].resize(3 * n_dofs_1d);
+
+              Point<dim> unit_point;
+              for (unsigned int i = 0; i < n_dofs_1d; ++i)
+                {
+                  // need to reorder from hierarchical to lexicographic to get
+                  // the DoFs correct
+                  const unsigned int my_i =
+                    (direction == 0) ? lex_normal[i] : lex_tangent[i];
+                  for (unsigned int q = 0; q < n_q_points_1d; ++q)
+                    {
+                      Point<dim> q_point = unit_point;
+                      q_point[direction] = quad.get_points()[q][0];
+
+                      shape_values[i * n_q_points_1d + q] =
+                        fe.shape_value_component(my_i, q_point, 0);
+                      shape_gradients[i * n_q_points_1d + q] =
+                        fe.shape_grad_component(my_i, q_point, 0)[direction];
+                      shape_hessians[i * n_q_points_1d + q] =
+                        fe.shape_grad_grad_component(my_i,
+                                                     q_point,
+                                                     0)[direction][direction];
+
+                      // evaluate basis functions on the two 1D subfaces (i.e.,
+                      // at the positions divided by one half and shifted by one
+                      // half, respectively) for hanging nodes
+                      q_point[direction] *= 0.5;
+                      values_within_subface[0][i * n_q_points_1d + q] =
+                        fe.shape_value_component(my_i, q_point, 0);
+                      gradients_within_subface[0][i * n_q_points_1d + q] =
+                        fe.shape_grad_component(my_i, q_point, 0)[direction];
+                      hessians_within_subface[0][i * n_q_points_1d + q] =
+                        fe.shape_grad_grad_component(my_i,
+                                                     q_point,
+                                                     0)[direction][direction];
+                      q_point[direction] += 0.5;
+                      values_within_subface[1][i * n_q_points_1d + q] =
+                        fe.shape_value_component(my_i, q_point, 0);
+                      gradients_within_subface[1][i * n_q_points_1d + q] =
+                        fe.shape_grad_component(my_i, q_point, 0)[direction];
+                      hessians_within_subface[1][i * n_q_points_1d + q] =
+                        fe.shape_grad_grad_component(my_i,
+                                                     q_point,
+                                                     0)[direction][direction];
+                    }
+                  // evaluate basis functions on the 1D faces, i.e., in zero and
+                  // one
+                  Point<dim> q_point = unit_point;
+                  q_point[direction] = 0;
+                  shape_data_on_face[0][i] =
+                    fe.shape_value_component(my_i, q_point, 0);
+                  shape_data_on_face[0][i + n_dofs_1d] =
+                    fe.shape_grad_component(my_i, q_point, 0)[direction];
+                  shape_data_on_face[0][i + 2 * n_dofs_1d] =
+                    fe.shape_grad_grad_component(my_i,
+                                                 q_point,
+                                                 0)[direction][direction];
+                  q_point[direction] = 1;
+                  shape_data_on_face[1][i] =
+                    fe.shape_value_component(my_i, q_point, 0);
+                  shape_data_on_face[1][i + n_dofs_1d] =
+                    fe.shape_grad_component(my_i, q_point, 0)[direction];
+                  shape_data_on_face[1][i + 2 * n_dofs_1d] =
+                    fe.shape_grad_grad_component(my_i,
+                                                 q_point,
+                                                 0)[direction][direction];
+                }
+            }
+          return;
+        }
+      else if (quad_in.is_tensor_product() == false ||
+               dynamic_cast<const FE_SimplexP<dim> *>(
+                 &fe_in.base_element(base_element_number)) ||
+               dynamic_cast<const FE_SimplexDGP<dim> *>(
+                 &fe_in.base_element(base_element_number)) ||
+               dynamic_cast<const FE_WedgeP<dim> *>(
+                 &fe_in.base_element(base_element_number)) ||
+               dynamic_cast<const FE_PyramidP<dim> *>(
+                 &fe_in.base_element(base_element_number)))
         {
           // specialization for arbitrary finite elements and quadrature rules
           // as needed in the context, e.g., of simplices
index 608367b3b209743bfbcb900d8c384fac7a812e6e..2e8885f01fa6c8b3a354ab1aa900dec5d889e938 100644 (file)
@@ -63,7 +63,11 @@ namespace internal
      * coefficient arrays. See the documentation of the EvaluatorTensorProduct
      * specialization for more information.
      */
-    evaluate_symmetric_hierarchical
+    evaluate_symmetric_hierarchical,
+    /**
+     * Raviart-Thomas elements with anisotropic polynomials.
+     */
+    evaluate_raviart_thomas
   };
 
 
@@ -194,21 +198,30 @@ namespace internal
       (void)dummy2;
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int  direction,
+              bool contract_over_rows,
+              bool add,
+              int  normal_dir = 0>
     void
     values(const Number in[], Number out[]) const
     {
       apply<direction, contract_over_rows, add>(shape_values, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int  direction,
+              bool contract_over_rows,
+              bool add,
+              int  normal_dir = 0>
     void
     gradients(const Number in[], Number out[]) const
     {
       apply<direction, contract_over_rows, add>(shape_gradients, in, out);
     }
 
-    template <int direction, bool contract_over_rows, bool add>
+    template <int  direction,
+              bool contract_over_rows,
+              bool add,
+              int  normal_dir = 0>
     void
     hessians(const Number in[], Number out[]) const
     {
@@ -2328,6 +2341,482 @@ namespace internal
 
 
 
+  /**
+   * Internal evaluator for shape function in 2D and 3D using the
+   * tensor product form of the anisotropic basis functions of the
+   * raviart-thomas element, with degree k+1 in normal direction and
+   * k in tangential direction.
+   *
+   * @tparam dim Space dimension in which this class is applied
+   * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+   *                to the number of 1d shape functions in the usual tensor
+   *                contraction setting
+   * @tparam n_columns Number of columns in the transformation matrix, which
+   *                   corresponds to the number of 1d shape functions in the
+   *                   usual tensor contraction setting
+   * @tparam Number Abstract number type for input and output arrays
+   * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+   *                 same type as the input/output arrays); must implement
+   *                 operator* with Number and produce Number as an output to
+   *                 be a valid type
+   */
+  template <int dim,
+            int n_rows,
+            int n_columns,
+            typename Number,
+            typename Number2>
+  struct EvaluatorTensorProduct<evaluate_raviart_thomas,
+                                dim,
+                                n_rows,
+                                n_columns,
+                                Number,
+                                Number2>
+  {
+    static constexpr unsigned int n_rows_of_product =
+      numbers::invalid_unsigned_int;
+    static constexpr unsigned int n_columns_of_product =
+      numbers::invalid_unsigned_int;
+
+    /**
+     * Empty constructor. Does nothing. Be careful when using 'values' and
+     * related methods because they need to be filled with the other pointer
+     */
+    EvaluatorTensorProduct()
+      : shape_values(nullptr)
+      , shape_gradients(nullptr)
+      , shape_hessians(nullptr)
+    {}
+
+    /**
+     * Constructor, taking the data from ShapeInfo
+     */
+    EvaluatorTensorProduct(const AlignedVector<Number2> &shape_values,
+                           const AlignedVector<Number2> &shape_gradients,
+                           const AlignedVector<Number2> &shape_hessians,
+                           const unsigned int            dummy1 = 0,
+                           const unsigned int            dummy2 = 0)
+      : shape_values(shape_values.begin())
+      , shape_gradients(shape_gradients.begin())
+      , shape_hessians(shape_hessians.begin())
+    {
+      // We can enter this function either for the apply() path that has
+      // n_rows * n_columns entries or for the apply_face() path that only has
+      // n_rows * 3 entries in the array. Since we cannot decide about the use
+      // we must allow for both here.
+      Assert(shape_values.size() == 0 ||
+               shape_values.size() == n_rows * n_columns ||
+               shape_values.size() == 3 * n_rows,
+             ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
+      Assert(shape_gradients.size() == 0 ||
+               shape_gradients.size() == n_rows * n_columns,
+             ExcDimensionMismatch(shape_gradients.size(), n_rows * n_columns));
+      Assert(shape_hessians.size() == 0 ||
+               shape_hessians.size() == n_rows * n_columns,
+             ExcDimensionMismatch(shape_hessians.size(), n_rows * n_columns));
+      (void)dummy1;
+      (void)dummy2;
+    }
+
+    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    void
+    values(const Number in[], Number out[]) const
+    {
+      apply<direction, contract_over_rows, add, normal_dir>(shape_values,
+                                                            in,
+                                                            out);
+    }
+
+    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    void
+    gradients(const Number in[], Number out[]) const
+    {
+      apply<direction, contract_over_rows, add, normal_dir>(shape_gradients,
+                                                            in,
+                                                            out);
+    }
+
+    template <int direction, bool contract_over_rows, bool add, int normal_dir>
+    void
+    hessians(const Number in[], Number out[]) const
+    {
+      apply<direction, contract_over_rows, add, normal_dir>(shape_hessians,
+                                                            in,
+                                                            out);
+    }
+
+    /**
+     * This function applies the tensor product kernel, corresponding to a
+     * multiplication of 1D stripes, along the given @p direction of the tensor
+     * data in the input array. This function allows the @p in and @p out
+     * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
+     * perform the contraction in place where @p in and @p out point to the
+     * same address. For the case n_rows != n_columns, the output is only
+     * correct if @p one_line is set to true.
+     *
+     * @tparam direction Direction that is evaluated
+     * @tparam contract_over_rows If true, the tensor contraction sums
+     *                            over the rows in the given @p shape_data
+     *                            array, otherwise it sums over the columns
+     * @tparam add If true, the result is added to the output vector, else
+     *             the computed values overwrite the content in the output
+     * @tparam normal_dir Indicates the direction of the continuous component of the
+     *                    RT space in terms of the normal onto the face, e.g
+     *                    0 if the  is in x-direction, 1 if in y-direction
+     *                    etc.
+     * @tparam one_line If true, the kernel is only applied along a single 1D
+     *                  stripe within a dim-dimensional tensor, not the full
+     *                  n_rows^dim points as in the @p false case.
+     *
+     * @param shape_data Transformation matrix with @p n_rows rows and
+     *                   @p n_columns columns, stored in row-major format
+     * @param in Pointer to the start of the input data vector
+     * @param out Pointer to the start of the output data vector
+     */
+    template <int  direction,
+              bool contract_over_rows,
+              bool add,
+              int  normal_dir,
+              bool one_line = false>
+    static void
+    apply(const Number2 *DEAL_II_RESTRICT shape_data,
+          const Number *                  in,
+          Number *                        out);
+
+    template <int  face_direction,
+              bool contract_onto_face,
+              bool add,
+              int  max_derivative,
+              bool lex_faces = false,
+              int  normal_direction>
+    void
+    apply_face(const Number *DEAL_II_RESTRICT in,
+               Number *DEAL_II_RESTRICT       out) const;
+
+    const Number2 *shape_values;
+    const Number2 *shape_gradients;
+    const Number2 *shape_hessians;
+  };
+
+  template <int dim,
+            int n_rows,
+            int n_columns,
+            typename Number,
+            typename Number2>
+  template <int  direction,
+            bool contract_over_rows,
+            bool add,
+            int  normal_dir,
+            bool one_line>
+  inline void
+  EvaluatorTensorProduct<evaluate_raviart_thomas,
+                         dim,
+                         n_rows,
+                         n_columns,
+                         Number,
+                         Number2>::apply(const Number2 *DEAL_II_RESTRICT
+                                                       shape_data,
+                                         const Number *in,
+                                         Number *      out)
+  {
+    static_assert(one_line == false || direction == dim - 1,
+                  "Single-line evaluation only works for direction=dim-1.");
+    Assert(shape_data != nullptr,
+           ExcMessage(
+             "The given array shape_data must not be the null pointer!"));
+    Assert(dim == direction + 1 || one_line == true || n_rows == n_columns ||
+             in != out,
+           ExcMessage("In-place operation only supported for "
+                      "n_rows==n_columns or single-line interpolation"));
+    AssertIndexRange(direction, dim);
+    constexpr int mm = contract_over_rows ? n_rows : n_columns,
+                  nn = contract_over_rows ? n_columns : n_rows;
+
+    constexpr int stride    = Utilities::pow(n_columns, direction);
+    constexpr int n_blocks1 = one_line ? 1 : stride;
+
+    // The number of blocks depend on both direction and dimension.
+    constexpr int n_blocks2 =
+      (dim - direction - 1 == 0) ?
+        1 :
+        ((direction == normal_dir) ?
+           Utilities::pow((n_rows - 1),
+                          (direction >= dim) ? 0 : dim - direction - 1) :
+           (((direction < normal_dir) ? (n_rows + 1) : n_rows) *
+            ((dim - direction == 3) ? n_rows : 1)));
+
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            Number x[mm];
+            for (int i = 0; i < mm; ++i)
+              x[i] = in[stride * i];
+
+            for (int col = 0; col < nn; ++col)
+              {
+                Number2 val0;
+
+                if (contract_over_rows)
+                  val0 = shape_data[col];
+                else
+                  val0 = shape_data[col * n_columns];
+
+                Number res0 = val0 * x[0];
+                for (int i = 1; i < mm; ++i)
+                  {
+                    if (contract_over_rows)
+                      val0 = shape_data[i * n_columns + col];
+                    else
+                      val0 = shape_data[col * n_columns + i];
+
+                    res0 += val0 * x[i];
+                  }
+                if (add)
+                  out[stride * col] += res0;
+
+                else
+                  out[stride * col] = res0;
+              }
+
+            if (one_line == false)
+              {
+                ++in;
+                ++out;
+              }
+          }
+        if (one_line == false)
+          {
+            in += stride * (mm - 1);
+            out += stride * (nn - 1);
+          }
+      }
+  }
+
+  template <int dim,
+            int n_rows,
+            int n_columns,
+            typename Number,
+            typename Number2>
+  template <int  face_direction,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative,
+            bool lex_faces,
+            int  normal_direction>
+  inline void
+  EvaluatorTensorProduct<evaluate_raviart_thomas,
+                         dim,
+                         n_rows,
+                         n_columns,
+                         Number,
+                         Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
+                                              Number *DEAL_II_RESTRICT
+                                                out) const
+  {
+    Assert(dim > 1 && (lex_faces || dim < 4),
+           ExcMessage("Only dim=2,3 supported"));
+    static_assert(max_derivative >= 0 && max_derivative < 3,
+                  "Only derivative orders 0-2 implemented");
+    static_assert(!lex_faces,
+                  "lex_faces = True is not implemented for Raviart-Thomas");
+    Assert(shape_values != nullptr,
+           ExcMessage(
+             "The given array shape_values must not be the null pointer."));
+
+    // Determine the number of blocks depending on the face and normaldirection,
+    // as well as dimension.
+    constexpr int n_blocks1 = (face_direction == normal_direction) ?
+                                (n_rows - 1) :
+                              ((face_direction == 0 && normal_direction == 2) ||
+                               (face_direction == 1 && normal_direction == 2) ||
+                               (face_direction == 2 && normal_direction == 1)) ?
+                                n_rows :
+                                (n_rows + 1);
+    constexpr int n_blocks2 =
+      (dim == 2) ? 1 :
+                   ((face_direction == normal_direction) ?
+                      (n_rows - 1) :
+                      (((face_direction == 0 && normal_direction == 1) ||
+                        (face_direction == 1 && normal_direction == 0) ||
+                        (face_direction == 2 && normal_direction == 0)) ?
+                         n_rows :
+                         (n_rows + 1)));
+
+    AssertIndexRange(face_direction, dim);
+
+    constexpr int stride =
+      (face_direction == normal_direction) ?
+        Utilities::pow(n_rows - 1, face_direction) :
+        ((face_direction == 0) ?
+           1 :
+           ((face_direction == 2) ?
+              n_rows * (n_rows + 1) :
+              ((face_direction == 1 && normal_direction == 0) ? (n_rows + 1) :
+                                                                n_rows)));
+    constexpr int out_stride = n_blocks1 * n_blocks2;
+
+    const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
+
+    for (int i2 = 0; i2 < n_blocks2; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks1; ++i1)
+          {
+            if (contract_onto_face == true)
+              {
+                Number res0 = shape_values[0] * in[0];
+                Number res1, res2;
+
+                if (max_derivative > 0)
+                  res1 = shape_values[n_rows] * in[0];
+
+                if (max_derivative > 1)
+                  res2 = shape_values[2 * n_rows] * in[0];
+
+                for (int ind = 1; ind < n_rows; ++ind)
+                  {
+                    res0 += shape_values[ind] * in[stride * ind];
+                    if (max_derivative > 0)
+                      res1 += shape_values[ind + n_rows] * in[stride * ind];
+
+                    if (max_derivative > 1)
+                      res2 += shape_values[ind + 2 * n_rows] * in[stride * ind];
+                  }
+                if (add)
+                  {
+                    out[0] += res0;
+
+                    if (max_derivative > 0)
+                      out[out_stride] += res1;
+
+                    if (max_derivative > 1)
+                      out[2 * out_stride] += res2;
+                  }
+                else
+                  {
+                    out[0] = res0;
+
+                    if (max_derivative > 0)
+                      out[out_stride] = res1;
+
+                    if (max_derivative > 1)
+                      out[2 * out_stride] = res2;
+                  }
+              }
+            else
+              {
+                for (int col = 0; col < n_rows; ++col)
+                  {
+                    if (add)
+                      out[col * stride] += shape_values[col] * in[0];
+                    else
+                      out[col * stride] = shape_values[col] * in[0];
+
+                    if (max_derivative > 0)
+                      out[col * stride] +=
+                        shape_values[col + n_rows] * in[out_stride];
+
+                    if (max_derivative > 1)
+                      out[col * stride] +=
+                        shape_values[col + 2 * n_rows] * in[2 * out_stride];
+                  }
+              }
+
+            // increment: in regular case, just go to the next point in
+            // x-direction. If we are at the end of one chunk in x-dir, need
+            // to jump over to the next layer in z-direction
+            switch (face_direction)
+              {
+                case 0:
+                  in += contract_onto_face ? n_rows : 1;
+                  out += contract_onto_face ? 1 : n_rows;
+                  break;
+
+                case 1:
+                  ++in;
+                  ++out;
+                  // faces 2 and 3 in 3D use local coordinate system zx, which
+                  // is the other way around compared to the tensor
+                  // product. Need to take that into account.
+                  if (dim == 3)
+                    {
+                      if (normal_direction == 0)
+                        {
+                          if (contract_onto_face)
+                            out += n_rows - 1;
+                          else
+                            in += n_rows - 1;
+                        }
+                      if (normal_direction == 1)
+                        {
+                          if (contract_onto_face)
+                            out += n_rows - 2;
+                          else
+                            in += n_rows - 2;
+                        }
+                      if (normal_direction == 2)
+                        {
+                          if (contract_onto_face)
+                            out += n_rows;
+                          else
+                            in += n_rows;
+                        }
+                    }
+                  break;
+
+                case 2:
+                  ++in;
+                  ++out;
+                  break;
+
+                default:
+                  Assert(false, ExcNotImplemented());
+              }
+          }
+        if (face_direction == 1 && dim == 3)
+          {
+            // adjust for local coordinate system zx
+            if (contract_onto_face)
+              {
+                if (normal_direction == 0)
+                  {
+                    in += (n_rows + 1) * (n_rows - 1);
+                    out -= n_rows * (n_rows + 1) - 1;
+                  }
+                if (normal_direction == 1)
+                  {
+                    in += (n_rows - 1) * (n_rows - 1);
+                    out -= (n_rows - 1) * (n_rows - 1) - 1;
+                  }
+                if (normal_direction == 2)
+                  {
+                    in += (n_rows - 1) * (n_rows);
+                    out -= (n_rows) * (n_rows + 1) - 1;
+                  }
+              }
+            else
+              {
+                if (normal_direction == 0)
+                  {
+                    out += (n_rows + 1) * (n_rows - 1);
+                    in -= n_rows * (n_rows + 1) - 1;
+                  }
+                if (normal_direction == 1)
+                  {
+                    out += (n_rows - 1) * (n_rows - 1);
+                    in -= (n_rows - 1) * (n_rows - 1) - 1;
+                  }
+                if (normal_direction == 2)
+                  {
+                    out += (n_rows - 1) * (n_rows);
+                    in -= (n_rows) * (n_rows + 1) - 1;
+                  }
+              }
+          }
+      }
+  }
+
+
+
   /**
    * Struct to avoid using Tensor<1, dim, Point<dim2>> in
    * evaluate_tensor_product_value_and_gradient because a Point cannot be used
index 08ecfed20f342e57c65f6115caaa5d014afd2aca..8397f71be1172bbd45fa6e1b9a4b16bbee33d785 100644 (file)
@@ -147,67 +147,6 @@ namespace
   }
 
 
-
-  // 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)
@@ -219,11 +158,11 @@ namespace
     // 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);
+      internal::get_lexicographic_numbering_rt_nodal<dim>(degree + 1);
+
+    hierarchic_to_lexicographic =
+      Utilities::invert_permutation(lexicographic_to_hierarchic);
 
     // since we only store an anisotropic polynomial for the first component,
     // we set up a second numbering to switch out the actual coordinate
@@ -426,7 +365,61 @@ namespace
   }
 } // namespace
 
+namespace internal
+{
+  template <int dim>
+  std::vector<unsigned int>
+  get_lexicographic_numbering_rt_nodal(const unsigned int degree)
+  {
+    const unsigned int        n_dofs_face = Utilities::pow(degree, dim - 1);
+    std::vector<unsigned int> lexicographic_numbering;
+    // component 1
+    for (unsigned int j = 0; j < n_dofs_face; j++)
+      {
+        lexicographic_numbering.push_back(j);
+        for (unsigned int i = n_dofs_face * 2 * dim;
+             i < n_dofs_face * 2 * dim + degree - 1;
+             i++)
+          lexicographic_numbering.push_back(i + j * (degree - 1));
+        lexicographic_numbering.push_back(n_dofs_face + j);
+      }
+
+    // component 2
+    unsigned int layers = (dim == 3) ? degree : 1;
+    for (unsigned int k = 0; k < layers; k++)
+      {
+        unsigned int k_add = k * degree;
+        for (unsigned int j = n_dofs_face * 2; j < n_dofs_face * 2 + degree;
+             j++)
+          lexicographic_numbering.push_back(j + k_add);
+
+        for (unsigned int i = n_dofs_face * (2 * dim + (degree - 1));
+             i < n_dofs_face * (2 * dim + (degree - 1)) + (degree - 1) * degree;
+             i++)
+          {
+            lexicographic_numbering.push_back(i + k_add * (degree - 1));
+          }
+        for (unsigned int j = n_dofs_face * 3; j < n_dofs_face * 3 + degree;
+             j++)
+          lexicographic_numbering.push_back(j + k_add);
+      }
 
+    // component 3
+    if (dim == 3)
+      {
+        for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; i++)
+          lexicographic_numbering.push_back(i);
+        for (unsigned int i = 6 * n_dofs_face + n_dofs_face * 2 * (degree - 1);
+             i < 6 * n_dofs_face + n_dofs_face * 3 * (degree - 1);
+             i++)
+          lexicographic_numbering.push_back(i);
+        for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; i++)
+          lexicographic_numbering.push_back(i);
+      }
+
+    return lexicographic_numbering;
+  }
+} // namespace internal
 
 // --------------------- actual implementation of element --------------------
 
index bdaa5a92208c7d6b80d65bbb2cedc3ac8abd3cc8..ecc6b35ab37bad44d0b3e58e7ae448dda2184cfc 100644 (file)
@@ -23,7 +23,17 @@ SET(_src
   evaluation_template_factory_inst4.cc
   evaluation_template_factory_inst5.cc
   evaluation_template_factory_inst6.cc
+  evaluation_template_factory_inst7.cc
+  evaluation_template_factory_inst8.cc
+  evaluation_template_factory_inst9.cc
+  evaluation_template_factory_inst10.cc
   evaluation_template_factory_hanging_nodes.cc
+  evaluation_template_face_factory.cc
+  evaluation_template_face_factory_inst2.cc
+  evaluation_template_face_factory_inst3.cc
+  evaluation_template_face_factory_inst4.cc
+  evaluation_template_face_factory_inst5.cc
+  evaluation_template_face_factory_inst6.cc
   fe_point_evaluation.cc
   mapping_info.cc
   mapping_info_inst2.cc
@@ -36,6 +46,7 @@ SET(_src
 
 SET(_inst
   evaluation_template_factory.inst.in
+  evaluation_template_face_factory.inst.in
   evaluation_template_factory_hanging_nodes.inst.in
   fe_point_evaluation.inst.in
   mapping_info.inst.in
diff --git a/source/matrix_free/evaluation_template_face_factory.cc b/source/matrix_free/evaluation_template_face_factory.cc
new file mode 100644 (file)
index 0000000..9a2df3a
--- /dev/null
@@ -0,0 +1,28 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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/matrix_free/evaluation_template_face_factory.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+#define SPLIT_INSTANTIATIONS_COUNT 6
+#ifndef SPLIT_INSTANTIATIONS_INDEX
+#  define SPLIT_INSTANTIATIONS_INDEX 0
+#endif
+
+#include "evaluation_template_face_factory.inst"
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/matrix_free/evaluation_template_face_factory.inst.in b/source/matrix_free/evaluation_template_face_factory.inst.in
new file mode 100644 (file)
index 0000000..d98cac1
--- /dev/null
@@ -0,0 +1,32 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS;
+     deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
+  {
+    template struct dealii::internal::
+      FEFaceEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
+
+    template struct dealii::internal::FEFaceEvaluationGatherFactory<
+      deal_II_dimension,
+      double,
+      deal_II_scalar_vectorized>;
+
+    template struct dealii::internal::FEFaceEvaluationGatherFactory<
+      deal_II_dimension,
+      float,
+      deal_II_scalar_vectorized>;
+  }
diff --git a/source/matrix_free/evaluation_template_face_factory_inst2.cc b/source/matrix_free/evaluation_template_face_factory_inst2.cc
new file mode 100644 (file)
index 0000000..5b3b840
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 1
+#include "evaluation_template_face_factory.cc"
diff --git a/source/matrix_free/evaluation_template_face_factory_inst3.cc b/source/matrix_free/evaluation_template_face_factory_inst3.cc
new file mode 100644 (file)
index 0000000..979c7e4
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 2
+#include "evaluation_template_face_factory.cc"
diff --git a/source/matrix_free/evaluation_template_face_factory_inst4.cc b/source/matrix_free/evaluation_template_face_factory_inst4.cc
new file mode 100644 (file)
index 0000000..b444b11
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 3
+#include "evaluation_template_face_factory.cc"
diff --git a/source/matrix_free/evaluation_template_face_factory_inst5.cc b/source/matrix_free/evaluation_template_face_factory_inst5.cc
new file mode 100644 (file)
index 0000000..bcfd448
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 4
+#include "evaluation_template_face_factory.cc"
diff --git a/source/matrix_free/evaluation_template_face_factory_inst6.cc b/source/matrix_free/evaluation_template_face_factory_inst6.cc
new file mode 100644 (file)
index 0000000..8aaeede
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 5
+#include "evaluation_template_face_factory.cc"
index cab955b323d39a041a31ee7fe932afbc7da373a3..9d2bf6931984391de507a5a6523fb2340ee7f7e2 100644 (file)
@@ -18,7 +18,7 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-#define SPLIT_INSTANTIATIONS_COUNT 6
+#define SPLIT_INSTANTIATIONS_COUNT 10
 #ifndef SPLIT_INSTANTIATIONS_INDEX
 #  define SPLIT_INSTANTIATIONS_INDEX 0
 #endif
index 03b1de0a40c356c7ab51b057c36ceff6c1162776..441855d054fda70ab27303de92f561846b4a0539 100644 (file)
@@ -20,19 +20,6 @@ for (deal_II_dimension : DIMENSIONS;
     template struct dealii::internal::
       FEEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
 
-    template struct dealii::internal::
-      FEFaceEvaluationFactory<deal_II_dimension, deal_II_scalar_vectorized>;
-
-    template struct dealii::internal::FEFaceEvaluationGatherFactory<
-      deal_II_dimension,
-      double,
-      deal_II_scalar_vectorized>;
-
-    template struct dealii::internal::FEFaceEvaluationGatherFactory<
-      deal_II_dimension,
-      float,
-      deal_II_scalar_vectorized>;
-
     // inverse mass
     template struct dealii::internal::
       CellwiseInverseMassFactory<deal_II_dimension, deal_II_scalar_vectorized>;
diff --git a/source/matrix_free/evaluation_template_factory_inst10.cc b/source/matrix_free/evaluation_template_factory_inst10.cc
new file mode 100644 (file)
index 0000000..b93c12d
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 9
+#include "evaluation_template_factory.cc"
diff --git a/source/matrix_free/evaluation_template_factory_inst7.cc b/source/matrix_free/evaluation_template_factory_inst7.cc
new file mode 100644 (file)
index 0000000..d71730b
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 6
+#include "evaluation_template_factory.cc"
diff --git a/source/matrix_free/evaluation_template_factory_inst8.cc b/source/matrix_free/evaluation_template_factory_inst8.cc
new file mode 100644 (file)
index 0000000..47f1f7a
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 7
+#include "evaluation_template_factory.cc"
diff --git a/source/matrix_free/evaluation_template_factory_inst9.cc b/source/matrix_free/evaluation_template_factory_inst9.cc
new file mode 100644 (file)
index 0000000..94dc699
--- /dev/null
@@ -0,0 +1,17 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+#define SPLIT_INSTANTIATIONS_INDEX 8
+#include "evaluation_template_factory.cc"
diff --git a/tests/matrix_free/matrix_vector_rt_01.cc b/tests/matrix_free/matrix_vector_rt_01.cc
new file mode 100644 (file)
index 0000000..d358717
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+// This function tests the correctness of the matrix-free implementation
+// of the FE_RaviartThomasNodal element by evaluating a simple fe operator
+// and comparing the result with FEVaules which is considered the
+// reference. The mesh is a hypercube mesh with no hanging nodes and no other
+// constraints
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+
+  DoFHandler<dim> dof(tria);
+
+  dof.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+  do_test<dim, fe_degree, double>(dof, constraints);
+}
diff --git a/tests/matrix_free/matrix_vector_rt_01.output b/tests/matrix_free/matrix_vector_rt_01.output
new file mode 100644 (file)
index 0000000..0e53c41
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 144
+DEAL:2d::
+DEAL:2d::Norm of difference: 6.20418e-16
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 16
+DEAL:2d::Number of degrees of freedom: 312
+DEAL:2d::
+DEAL:2d::Norm of difference: 7.03559e-16
+DEAL:2d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 64
+DEAL:3d::Number of degrees of freedom: 1728
+DEAL:3d::
+DEAL:3d::Norm of difference: 1.04798e-15
+DEAL:3d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Number of cells: 64
+DEAL:3d::Number of degrees of freedom: 5616
+DEAL:3d::
+DEAL:3d::Norm of difference: 1.49180e-15
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_rt_common.h b/tests/matrix_free/matrix_vector_rt_common.h
new file mode 100644 (file)
index 0000000..e2d0ae2
--- /dev/null
@@ -0,0 +1,241 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+// This test template evaluates a simple operator on FE_PolyTensor
+// elements using FEEvaluation and compares the result with the output
+// of FEValues (which is considered to be the reference) on cartesian
+// meshes without hanging nodes. (It will be extended to also handle general
+// meshes and hanging nodes in the future.) The test do not include
+// multithreading because FEValues is not thread-safe.
+// See matrix_vector_rt_01.cc for an example that uses this template.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/matrix_free.templates.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+// forward declare this function. will be implemented in .cc files
+template <int dim, int fe_degree>
+void
+test();
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  virtual ~MatrixFreeTest(){};
+
+  void
+  operator()(const MatrixFree<dim, Number> &              data,
+             Vector<Number> &                             dst,
+             const Vector<Number> &                       src,
+             const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data);
+
+    // OBS! This will need to be modified once the Piola transform is
+    // implemented
+    unsigned int n_cells =
+      data.get_dof_handler().get_triangulation().n_active_cells();
+    Number piola =
+      (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        fe_eval.reinit(cell);
+        fe_eval.gather_evaluate(src,
+                                EvaluationFlags::values |
+                                  EvaluationFlags::gradients);
+
+        for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+          {
+            fe_eval.submit_value(Number(10 * piola) * fe_eval.get_value(q), q);
+            fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
+          }
+
+        fe_eval.integrate_scatter(EvaluationFlags::values |
+                                    EvaluationFlags::gradients,
+                                  dst);
+      }
+  };
+
+
+  void
+  test_functions(Vector<Number> &dst, const Vector<Number> &src) const
+  {
+    data.cell_loop(&MatrixFreeTest::operator(), this, dst, src);
+  };
+
+protected:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+template <int dim, int fe_degree, typename Number>
+void
+do_test(const DoFHandler<dim> &          dof,
+        const AffineConstraints<double> &constraints)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+          << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+          << std::endl;
+
+
+  //   constraints.distribute(solution);
+  MatrixFree<dim, Number> mf_data;
+  {
+    const QGaussLobatto<1>                           quad(fe_degree + 2);
+    typename MatrixFree<dim, Number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
+    data.mapping_update_flags  = update_gradients | update_hessians;
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  // create vector with random entries
+  Vector<Number> solution, initial_condition;
+
+  mf_data.initialize_dof_vector(solution);
+  mf_data.initialize_dof_vector(initial_condition);
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      initial_condition[i] = random_value<Number>();
+    }
+
+  // MatrixFree solution
+  MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+  mf.test_functions(solution, initial_condition);
+
+
+  SparsityPattern        sp;
+  SparseMatrix<double>   system_matrix;
+  DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
+  DoFTools::make_sparsity_pattern(dof, dsp);
+  sp.copy_from(dsp);
+  system_matrix.reinit(sp);
+
+  FEValues<dim> fe_val(dof.get_fe(),
+                       Quadrature<dim>(mf_data.get_quadrature(0)),
+                       update_values | update_gradients | update_JxW_values);
+
+
+  const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+  FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+
+
+  const FEValuesExtractors::Vector velocities(0);
+
+  // Assemble matrix
+  for (const auto &cell : dof.active_cell_iterators())
+    {
+      fe_val.reinit(cell);
+      local_matrix = 0;
+
+      for (const auto q : fe_val.quadrature_point_indices())
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.;
+            const Tensor<2, dim> grad_phi_i = fe_val[velocities].gradient(i, q);
+
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              {
+                const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q);
+                const Tensor<2, dim> grad_phi_j =
+                  fe_val[velocities].gradient(j, q);
+
+                local_matrix(i, j) +=
+                  (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) *
+                  fe_val.JxW(q);
+              }
+          }
+      cell->get_dof_indices(local_dof_indices);
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+          system_matrix.add(local_dof_indices[i],
+                            local_dof_indices[j],
+                            local_matrix(i, j));
+    }
+
+  Vector<Number> ref(solution.size());
+
+  // Compute reference
+  system_matrix.vmult(ref, initial_condition);
+
+  ref -= solution;
+
+  const double diff_norm = ref.linfty_norm() / solution.linfty_norm();
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("2d");
+    test<2, 2>();
+    test<2, 3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 2>();
+    test<3, 3>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_rt_face_01.cc b/tests/matrix_free/matrix_vector_rt_face_01.cc
new file mode 100644 (file)
index 0000000..efb1e82
--- /dev/null
@@ -0,0 +1,44 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+// This function tests the correctness of the matrix-free implementation
+// of the FE_RaviartThomasNodal element by evaluating a face operator
+// and comparing the result with FEVaules which is considered the
+// reference. The mesh is a hypercube mesh with no hanging nodes and no other
+// constraints
+
+#include "../tests.h"
+
+#include "matrix_vector_rt_face_common.h"
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+
+  FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+
+  DoFHandler<dim> dof(tria);
+
+  dof.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+  do_test<dim, fe_degree, double>(dof, constraints);
+}
diff --git a/tests/matrix_free/matrix_vector_rt_face_01.output b/tests/matrix_free/matrix_vector_rt_face_01.output
new file mode 100644 (file)
index 0000000..99a9ae9
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 40
+DEAL:2d::
+DEAL:2d::Norm of difference: 3.44596e-16
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Number of cells: 4
+DEAL:2d::Number of degrees of freedom: 84
+DEAL:2d::
+DEAL:2d::Norm of difference: 3.47432e-15
+DEAL:2d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Number of cells: 8
+DEAL:3d::Number of degrees of freedom: 240
+DEAL:3d::
+DEAL:3d::Norm of difference: 2.28531e-15
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_rt_face_common.h b/tests/matrix_free/matrix_vector_rt_face_common.h
new file mode 100644 (file)
index 0000000..922ab00
--- /dev/null
@@ -0,0 +1,311 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+// This test template evaluates a simple operator on FE_PolyTensor
+// elements using FEFaceEvaluation and compares the result with the output
+// of FEFaceValues (which is considered to be the reference) on cartesian
+// meshes without hanging nodes. (It will be extended to also handle general
+// meshes and hanging nodes in the future.) The test do not include
+// multithreading because FEFaceValues is not thread-safe.
+// See matrix_vector_rt_face_01.cc for an example that uses this template.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/matrix_free.templates.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+// forward declare this function. will be implemented in .cc files
+template <int dim, int fe_degree>
+void
+test();
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  virtual ~MatrixFreeTest(){};
+
+  void
+  dummy(const MatrixFree<dim, Number> &,
+        Vector<Number> &,
+        const Vector<Number> &,
+        const std::pair<unsigned int, unsigned int> &) const
+  {}
+
+  void
+  operator_face(const MatrixFree<dim, Number> &              data,
+                Vector<Number> &                             dst,
+                const Vector<Number> &                       src,
+                const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data,
+                                                                         true);
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval_n(
+      data, false);
+
+    // Note that this will need to be modified once the Piola transform is
+    // implemented
+    const unsigned int n_cells =
+      data.get_dof_handler().get_triangulation().n_active_cells();
+    const Number piola =
+      (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        fe_eval.reinit(face);
+        fe_eval_n.reinit(face);
+
+        fe_eval.gather_evaluate(src,
+                                EvaluationFlags::values |
+                                  EvaluationFlags::gradients);
+        fe_eval_n.gather_evaluate(src,
+                                  EvaluationFlags::values |
+                                    EvaluationFlags::gradients);
+
+        for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+          {
+            fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q);
+            fe_eval_n.submit_value(Number(10. * piola) * fe_eval_n.get_value(q),
+                                   q);
+
+            fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
+            fe_eval_n.submit_gradient(Number(piola) * fe_eval_n.get_gradient(q),
+                                      q);
+          }
+
+        fe_eval.integrate_scatter(EvaluationFlags::values |
+                                    EvaluationFlags::gradients,
+                                  dst);
+        fe_eval_n.integrate_scatter(EvaluationFlags::values |
+                                      EvaluationFlags::gradients,
+                                    dst);
+      }
+  };
+
+  void
+  operator_boundary(
+    const MatrixFree<dim, Number> &              data,
+    Vector<Number> &                             dst,
+    const Vector<Number> &                       src,
+    const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data,
+                                                                         true);
+
+    // Note that this will need to be modified once the Piola transform is
+    // implemented
+    const unsigned int n_cells =
+      data.get_dof_handler().get_triangulation().n_active_cells();
+    const Number piola =
+      (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4);
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        fe_eval.reinit(face);
+        fe_eval.gather_evaluate(src,
+                                EvaluationFlags::values |
+                                  EvaluationFlags::gradients);
+
+        for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+          {
+            fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q);
+            fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q);
+          }
+
+        fe_eval.integrate_scatter(EvaluationFlags::values |
+                                    EvaluationFlags::gradients,
+                                  dst);
+      }
+  };
+
+
+  void
+  test_functions(Vector<Number> &dst, const Vector<Number> &src) const
+  {
+    data.loop(&MatrixFreeTest::dummy,
+              &MatrixFreeTest::operator_face,
+              &MatrixFreeTest::operator_boundary,
+              this,
+              dst,
+              src);
+  };
+
+protected:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename Number>
+void
+do_test(const DoFHandler<dim> &          dof,
+        const AffineConstraints<double> &constraints)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  deallog << "Number of cells: " << dof.get_triangulation().n_active_cells()
+          << std::endl;
+  deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl
+          << std::endl;
+
+
+  //   constraints.distribute(solution);
+  MatrixFree<dim, Number> mf_data;
+  {
+    const QGaussLobatto<1>                           quad(fe_degree + 2);
+    const MappingQ<dim>                              mapping(fe_degree);
+    typename MatrixFree<dim, Number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, Number>::AdditionalData::none;
+    data.mapping_update_flags  = update_gradients | update_JxW_values;
+    data.mapping_update_flags_inner_faces =
+      (update_gradients | update_JxW_values);
+    data.mapping_update_flags_boundary_faces =
+      (update_gradients | update_JxW_values);
+    mf_data.reinit(mapping, dof, constraints, quad, data);
+  }
+
+  Vector<Number> solution, initial_condition;
+
+  // create vector with random entries
+  mf_data.initialize_dof_vector(solution);
+  mf_data.initialize_dof_vector(initial_condition);
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      initial_condition[i] = random_value<Number>();
+    }
+
+  MatrixFreeTest<dim, fe_degree, fe_degree + 2, Number> mf(mf_data);
+  mf.test_functions(solution, initial_condition);
+
+
+  SparsityPattern        sp;
+  SparseMatrix<double>   system_matrix;
+  DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs());
+  DoFTools::make_sparsity_pattern(dof, dsp);
+  sp.copy_from(dsp);
+  system_matrix.reinit(sp);
+
+  FEFaceValues<dim> fe_val(
+    dof.get_fe(),
+    QGaussLobatto<dim - 1>(mf_data.get_quadrature(0).size()),
+    update_values | update_gradients | update_JxW_values);
+
+
+  const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell;
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+  FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
+
+
+  const FEValuesExtractors::Vector velocities(0);
+  // Assemble matrix
+  for (const auto &cell : dof.active_cell_iterators())
+    for (const auto &face : cell->face_iterators())
+      {
+        fe_val.reinit(cell, face);
+        local_matrix = 0;
+
+        for (const auto q : fe_val.quadrature_point_indices())
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.;
+              const Tensor<2, dim> grad_phi_i =
+                fe_val[velocities].gradient(i, q);
+
+              for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                {
+                  const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q);
+                  const Tensor<2, dim> grad_phi_j =
+                    fe_val[velocities].gradient(j, q);
+                  local_matrix(i, j) +=
+                    (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) *
+                    fe_val.JxW(q);
+                }
+            }
+        cell->get_dof_indices(local_dof_indices);
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          for (unsigned int j = 0; j < dofs_per_cell; ++j)
+            system_matrix.add(local_dof_indices[i],
+                              local_dof_indices[j],
+                              local_matrix(i, j));
+      }
+
+  Vector<Number> ref(solution.size());
+
+  // Compute reference
+  system_matrix.vmult(ref, initial_condition);
+
+  ref -= solution;
+
+  const double diff_norm = ref.linfty_norm() / solution.linfty_norm();
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("2d");
+    test<2, 2>();
+    test<2, 3>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 2>();
+    // test<3, 3>(); //Takes way too long
+    deallog.pop();
+  }
+}

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.