Remove lex_faces from matrix-free kernels

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index ff05bbc5a0b6e3122738029c404ceede7e3f97a3..bf83df4ea81850a583299634549e0da3b07e6487 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -3536,7 +3536,7 @@ namespace internal
   };
 
 
-  template <int dim, int fe_degree, typename Number, bool lex_faces = false>
+  template <int dim, int fe_degree, typename Number>
   struct FEFaceNormalEvaluationImpl
   {
     template <bool do_evaluate, bool add_into_output>
@@ -3637,20 +3637,17 @@ namespace internal
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
-                                          2,
-                                          lex_faces>(input, output);
+                                          2>(input, output);
               else if (flag & EvaluationFlags::gradients)
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
-                                          1,
-                                          lex_faces>(input, output);
+                                          1>(input, output);
               else
                 evalf.template apply_face<face_direction,
                                           do_evaluate,
                                           add_into_output,
-                                          0,
-                                          lex_faces>(input, output);
+                                          0>(input, output);
               input += in_stride;
               output += out_stride;
             }
@@ -3841,8 +3838,7 @@ namespace internal
       evalf0.template apply_face<face_direction,
                                  do_evaluate,
                                  add_into_output,
-                                 max_derivative,
-                                 lex_faces>(input, output);
+                                 max_derivative>(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;
@@ -3850,8 +3846,7 @@ namespace internal
       evalf1.template apply_face<face_direction,
                                  do_evaluate,
                                  add_into_output,
-                                 max_derivative,
-                                 lex_faces>(input, output);
+                                 max_derivative>(input, output);
 
       if (dim == 3)
         {
@@ -3863,8 +3858,7 @@ namespace internal
           evalf2.template apply_face<face_direction,
                                      do_evaluate,
                                      add_into_output,
-                                     max_derivative,
-                                     lex_faces>(input, output);
+                                     max_derivative>(input, output);
         }
     }
   };

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index 50f31668005e027de460f673158f2e3845ac1132..7a21fcd1859d959a1b50bf101160cfb2c7e15962 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -334,11 +334,6 @@ namespace internal
      *             derivatives, 2 second derivates. Note that all the
      *             derivatives access the data in @p shape_values passed to
      *             the constructor of the class
-     * @tparam lex_faces Sets how the evaluation points on the faces should be
-     *                   sorted: lexicographically or right-hand-system number
-     *                   (special treatment of orientation 1 in 3D). Per default
-     *                   right-hand-system number is enabled, which is only
-     *                   working for dimensions up to 3.
      *
      * @param in address of the input data vector
      * @param out address of the output data vector
@@ -346,8 +341,7 @@ namespace internal
     template <int  face_direction,
               bool contract_onto_face,
               bool add,
-              int  max_derivative,
-              bool lex_faces = false>
+              int  max_derivative>
     void
     apply_face(const Number *DEAL_II_RESTRICT in,
                Number *DEAL_II_RESTRICT       out) const;
@@ -447,8 +441,7 @@ namespace internal
   template <int  face_direction,
             bool contract_onto_face,
             bool add,
-            int  max_derivative,
-            bool lex_faces>
+            int  max_derivative>
   inline void
   EvaluatorTensorProduct<evaluate_general,
                          dim,
@@ -459,21 +452,15 @@ namespace internal
                                               Number *DEAL_II_RESTRICT
                                                 out) const
   {
-    Assert(dim > 0 && (lex_faces || dim < 4),
-           ExcMessage("Only dim=1,2,3 supported"));
+    Assert(dim > 0, ExcMessage("Only dim=1,2,3 supported"));
     static_assert(max_derivative >= 0 && max_derivative < 3,
                   "Only derivative orders 0-2 implemented");
     Assert(shape_values != nullptr,
            ExcMessage(
              "The given array shape_values must not be the null pointer."));
 
-    constexpr int n_blocks1 =
-      lex_faces ? dealii::Utilities::pow<unsigned int>(n_rows, face_direction) :
-                  (dim > 1 ? n_rows : 1);
-    constexpr int n_blocks2 =
-      lex_faces ? dealii::Utilities::pow<unsigned int>(
-                    n_rows, std::max(dim - face_direction - 1, 0)) :
-                  (dim > 2 ? n_rows : 1);
+    constexpr int n_blocks1 = (dim > 1 ? n_rows : 1);
+    constexpr int n_blocks2 = (dim > 2 ? n_rows : 1);
 
     AssertIndexRange(face_direction, dim);
     constexpr int stride     = Utilities::pow(n_rows, face_direction);
@@ -534,55 +521,41 @@ namespace internal
                   }
               }
 
-            if (lex_faces)
+            // 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)
               {
-                ++out;
-                ++in;
+                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 (contract_onto_face)
+                        out += n_rows - 1;
+                      else
+                        in += n_rows - 1;
+                    }
+                  break;
+                case 2:
+                  ++in;
+                  ++out;
+                  break;
+                default:
+                  Assert(false, ExcNotImplemented());
               }
-            else
-              // 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 (contract_onto_face)
-                          out += n_rows - 1;
-                        else
-                          in += n_rows - 1;
-                      }
-                    break;
-                  case 2:
-                    ++in;
-                    ++out;
-                    break;
-                  default:
-                    Assert(false, ExcNotImplemented());
-                }
-          }
-        if (lex_faces)
-          {
-            if (contract_onto_face)
-              in += (dealii::Utilities::pow(n_rows, face_direction + 1) -
-                     n_blocks1);
-            else
-              out += (dealii::Utilities::pow(n_rows, face_direction + 1) -
-                      n_blocks1);
           }
-        else if (face_direction == 1 && dim == 3)
+
+        // adjust for local coordinate system zx
+        if (face_direction == 1 && dim == 3)
           {
-            // adjust for local coordinate system zx
             if (contract_onto_face)
               {
                 in += n_rows * (n_rows - 1);
@@ -734,8 +707,7 @@ namespace internal
     template <int  face_direction,
               bool contract_onto_face,
               bool add,
-              int  max_derivative,
-              bool lex_faces = false>
+              int  max_derivative>
     void
     apply_face(const Number *DEAL_II_RESTRICT in,
                Number *DEAL_II_RESTRICT       out) const;
@@ -899,15 +871,12 @@ namespace internal
   template <int  face_direction,
             bool contract_onto_face,
             bool add,
-            int  max_derivative,
-            bool lex_faces>
+            int  max_derivative>
   inline void
   EvaluatorTensorProduct<evaluate_general, dim, 0, 0, Number, Number2>::
     apply_face(const Number *DEAL_II_RESTRICT in,
                Number *DEAL_II_RESTRICT       out) const
   {
-    static_assert(lex_faces == false, "Not implemented yet.");
-
     Assert(shape_values != nullptr,
            ExcMessage(
              "The given array shape_data must not be the null pointer!"));
@@ -2505,8 +2474,7 @@ namespace internal
     template <int  face_direction,
               bool contract_onto_face,
               bool add,
-              int  max_derivative,
-              bool lex_faces = false>
+              int  max_derivative>
     void
     apply_face(const Number *DEAL_II_RESTRICT in,
                Number *DEAL_II_RESTRICT       out) const;
@@ -2618,8 +2586,7 @@ namespace internal
   template <int  face_direction,
             bool contract_onto_face,
             bool add,
-            int  max_derivative,
-            bool lex_faces>
+            int  max_derivative>
   inline void
   EvaluatorTensorProductAnisotropic<
     evaluate_raviart_thomas,
@@ -2631,12 +2598,9 @@ namespace internal
     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"));
+    Assert(dim > 1 && 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."));