Merge pull request #13591 from NiklasWik/RT_matrix_free

Matrixfree evaluation with `FE_RaviartThomasNodal`
dealii · Apr 26, 2022 · e99aafd · e99aafd
2 parents d10093c + 6918f72
commit e99aafd
Show file tree

Hide file tree

Showing 30 changed files with 3,092 additions and 213 deletions.
diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml
@@ -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:

diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h
@@ -327,6 +327,14 @@ class FE_RaviartThomas : public FE_PolyTensor<dim>
  * 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>
 {

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
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
@@ -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
diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h
@@ -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(

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -1159,16 +1159,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);
                   }

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
@@ -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
+
+
     };