]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify normal evaluation for Raviart-Thomas
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Wed, 30 Aug 2023 11:59:15 +0000 (13:59 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 11 Sep 2023 19:58:51 +0000 (21:58 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/tensor_product_kernels.h

index 1401a041807bde48d97f4846dee00c68e2f686be..9e7dbc1a7522458e82aa9bdec66d097fac84c956 100644 (file)
@@ -3586,11 +3586,134 @@ namespace internal
 
       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);
+          constexpr int stride1 = Utilities::pow(fe_degree + 1, face_direction);
+          constexpr int stride0 = Utilities::pow(fe_degree, face_direction);
+          constexpr int stride2 = fe_degree * (fe_degree + 1);
+
+          const int degree =
+            fe_degree != -1 ? fe_degree : shape_info.data[0].fe_degree;
+          const int n_rows_n = degree + 1;
+          const int n_rows_t = degree;
+
+          std::array<int, 3> strides{{1, 1, 1}};
+          if (face_direction > 0)
+            {
+              strides[0] =
+                n_rows_n * Utilities::pow(n_rows_t, face_direction - 1);
+              strides[1] = n_rows_t * (face_direction == 3 ? n_rows_n : 1);
+              strides[2] = Utilities::pow(n_rows_t, face_direction);
+            }
+          const dealii::ndarray<int, 3, 3> dofs_per_direction{
+            {{{n_rows_n, n_rows_t, n_rows_t}},
+             {{n_rows_t, n_rows_n, n_rows_t}},
+             {{n_rows_t, n_rows_t, n_rows_n}}}};
+
+          std::array<int, 2> steps, n_blocks;
+
+          if constexpr (face_direction == 0)
+            steps = {{degree + (face_direction == 0), 0}};
+          else if constexpr (face_direction == 1 && dim == 2)
+            steps = {{1, 0}};
+          else if constexpr (face_direction == 1)
+            // in 3d, the coordinate system is zx, not xz -> switch indices
+            steps = {
+              {n_rows_n * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
+          else if constexpr (face_direction == 2)
+            steps = {{1, 0}};
+
+          n_blocks[0] = dofs_per_direction[0][(face_direction + 1) % dim];
+          n_blocks[1] =
+            dim > 2 ? dofs_per_direction[0][(face_direction + 2) % dim] : 1;
+
+          interpolate_to_face<
+            (fe_degree != -1 ? (fe_degree + (face_direction == 0)) : 0),
+            ((face_direction < 2) ? stride1 : stride2),
+            do_evaluate,
+            add_into_output,
+            max_derivative>(shape_info.data[face_direction != 0]
+                              .shape_data_on_face[face_no % 2]
+                              .begin(),
+                            n_blocks,
+                            steps,
+                            input,
+                            output,
+                            degree + (face_direction == 0),
+                            strides[0]);
+
+          if (do_evaluate)
+            {
+              input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+              output += 3 * n_blocks[0] * n_blocks[1];
+            }
+          else
+            {
+              output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+              input += 3 * n_blocks[0] * n_blocks[1];
+            }
+
+          // must only change steps only for face direction 0
+          if constexpr (face_direction == 0)
+            steps = {{degree, 0}};
+
+          n_blocks[0] = dofs_per_direction[1][(face_direction + 1) % dim];
+          n_blocks[1] =
+            dim > 2 ? dofs_per_direction[1][(face_direction + 2) % dim] : 1;
+
+          interpolate_to_face<
+            (fe_degree != -1 ? (fe_degree + (face_direction == 1)) : 0),
+            ((face_direction < 2) ? stride0 : stride2),
+            do_evaluate,
+            add_into_output,
+            max_derivative>(shape_info.data[face_direction != 1]
+                              .shape_data_on_face[face_no % 2]
+                              .begin(),
+                            n_blocks,
+                            steps,
+                            input,
+                            output,
+                            degree + (face_direction == 1),
+                            strides[1]);
+
+          if constexpr (dim > 2)
+            {
+              if (do_evaluate)
+                {
+                  input += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+                  output += 3 * n_blocks[0] * n_blocks[1];
+                }
+              else
+                {
+                  output += n_rows_n * Utilities::pow(n_rows_t, dim - 1);
+                  input += 3 * n_blocks[0] * n_blocks[1];
+                }
+
+              if constexpr (face_direction == 0)
+                steps = {{degree, 0}};
+              else if constexpr (face_direction == 1)
+                // in 3d, the coordinate system is zx, not xz -> switch indices
+                steps = {
+                  {n_rows_t * n_rows_t, -n_rows_n * n_rows_t * n_rows_t + 1}};
+              else if constexpr (face_direction == 2)
+                steps = {{1, 0}};
+
+              n_blocks[0] = dofs_per_direction[2][(face_direction + 1) % dim];
+              n_blocks[1] = dofs_per_direction[2][(face_direction + 2) % dim];
+
+              interpolate_to_face<
+                (fe_degree != -1 ? (fe_degree + (face_direction == 2)) : 0),
+                stride0,
+                do_evaluate,
+                add_into_output,
+                max_derivative>(shape_info.data[face_direction != 2]
+                                  .shape_data_on_face[face_no % 2]
+                                  .begin(),
+                                n_blocks,
+                                steps,
+                                input,
+                                output,
+                                degree + (face_direction == 2),
+                                strides[2]);
+            }
         }
       else if (face_direction == face_no / 2)
         {
@@ -3623,135 +3746,6 @@ namespace internal
             }
         }
     }
-
-    /* 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<Number2> &shape_info,
-      const unsigned int                             face_no,
-      const Number                                  *input,
-      Number                                        *output)
-    {
-      // These types are evaluators in either normal or tangential direction
-      // depending on the face direction, with different normal directions for
-      // the different components.
-      using Evalf0 = typename std::conditional<
-        face_direction == 0,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree + 1,
-                                          0,
-                                          0,
-                                          Number,
-                                          Number2>,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree,
-                                          0,
-                                          0,
-                                          Number,
-                                          Number2>>::type;
-      using Evalf1 = typename std::conditional<
-        face_direction == 1,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree + 1,
-                                          0,
-                                          1,
-                                          Number,
-                                          Number2>,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree,
-                                          0,
-                                          1,
-                                          Number,
-                                          Number2>>::type;
-      using Evalf2 = typename std::conditional<
-        face_direction == 2,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree + 1,
-                                          0,
-                                          2,
-                                          Number,
-                                          Number2>,
-        EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
-                                          dim,
-                                          (fe_degree == -1) ? 1 : fe_degree,
-                                          0,
-                                          2,
-                                          Number,
-                                          Number2>>::type;
-
-      Evalf0 evalf0 =
-        create_evaluator_tensor_product<Evalf0>((face_direction == 0) ?
-                                                  shape_info.data[0] :
-                                                  shape_info.data[1],
-                                                face_no);
-      Evalf1 evalf1 =
-        create_evaluator_tensor_product<Evalf1>((face_direction == 1) ?
-                                                  shape_info.data[0] :
-                                                  shape_info.data[1],
-                                                face_no);
-      Evalf2 evalf2 =
-        create_evaluator_tensor_product<Evalf2>((face_direction == 2) ?
-                                                  shape_info.data[0] :
-                                                  shape_info.data[1],
-                                                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>(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>(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>(input, output);
-        }
-    }
   };
 
 
index fc9c8429458f43eca3fd8e7a53fb88dbd12714d3..c9717aab71bea4218b157d8e293a2cb1b9b1c2eb 100644 (file)
@@ -1928,20 +1928,14 @@ namespace internal
           const Number                   *in,
           Number                         *out);
 
-    template <int  face_direction,
-              bool contract_onto_face,
-              bool add,
-              int  max_derivative>
-    void
-    apply_face(const Number *DEAL_II_RESTRICT in,
-               Number *DEAL_II_RESTRICT       out) const;
-
   private:
     const Number2 *shape_values;
     const Number2 *shape_gradients;
     const Number2 *shape_hessians;
   };
 
+
+
   template <int dim,
             int n_rows,
             int n_columns,
@@ -2014,226 +2008,6 @@ namespace internal
       }
   }
 
-  template <int dim,
-            int n_rows,
-            int n_columns,
-            int normal_dir,
-            typename Number,
-            typename Number2>
-  template <int  face_direction,
-            bool contract_onto_face,
-            bool add,
-            int  max_derivative>
-  inline void
-  EvaluatorTensorProductAnisotropic<
-    evaluate_raviart_thomas,
-    dim,
-    n_rows,
-    n_columns,
-    normal_dir,
-    Number,
-    Number2>::apply_face(const Number *DEAL_II_RESTRICT in,
-                         Number *DEAL_II_RESTRICT       out) const
-  {
-    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");
-    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_dir) ? (n_rows - 1) :
-                              ((face_direction == 0 && normal_dir == 2) ||
-                               (face_direction == 1 && normal_dir == 2) ||
-                               (face_direction == 2 && normal_dir == 1)) ?
-                                                               n_rows :
-                                                               (n_rows + 1);
-    constexpr int n_blocks2 = (dim == 2) ?
-                                1 :
-                                ((face_direction == normal_dir) ?
-                                   (n_rows - 1) :
-                                   (((face_direction == 0 && normal_dir == 1) ||
-                                     (face_direction == 1 && normal_dir == 0) ||
-                                     (face_direction == 2 && normal_dir == 0)) ?
-                                      n_rows :
-                                      (n_rows + 1)));
-
-    AssertIndexRange(face_direction, dim);
-
-    constexpr int in_stride =
-      (face_direction == normal_dir) ?
-        Utilities::pow(n_rows - 1, face_direction) :
-        ((face_direction == 0) ?
-           1 :
-           ((face_direction == 2) ?
-              n_rows * (n_rows + 1) :
-              ((face_direction == 1 && normal_dir == 0) ? (n_rows + 1) :
-                                                          n_rows)));
-    constexpr int out_stride = n_blocks1 * n_blocks2;
-
-    const Number2 *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[in_stride * ind];
-                    if (max_derivative > 0)
-                      res1 += shape_values[ind + n_rows] * in[in_stride * ind];
-
-                    if (max_derivative > 1)
-                      res2 +=
-                        shape_values[ind + 2 * n_rows] * in[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 * in_stride] += shape_values[col] * in[0];
-                    else
-                      out[col * in_stride] = shape_values[col] * in[0];
-
-                    if (max_derivative > 0)
-                      out[col * in_stride] +=
-                        shape_values[col + n_rows] * in[out_stride];
-
-                    if (max_derivative > 1)
-                      out[col * in_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_dir == 0)
-                        {
-                          if (contract_onto_face)
-                            out += n_rows - 1;
-                          else
-                            in += n_rows - 1;
-                        }
-                      if (normal_dir == 1)
-                        {
-                          if (contract_onto_face)
-                            out += n_rows - 2;
-                          else
-                            in += n_rows - 2;
-                        }
-                      if (normal_dir == 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_dir == 0)
-                  {
-                    in += (n_rows + 1) * (n_rows - 1);
-                    out -= n_rows * (n_rows + 1) - 1;
-                  }
-                if (normal_dir == 1)
-                  {
-                    in += (n_rows - 1) * (n_rows - 1);
-                    out -= (n_rows - 1) * (n_rows - 1) - 1;
-                  }
-                if (normal_dir == 2)
-                  {
-                    in += (n_rows - 1) * (n_rows);
-                    out -= (n_rows) * (n_rows + 1) - 1;
-                  }
-              }
-            else
-              {
-                if (normal_dir == 0)
-                  {
-                    out += (n_rows + 1) * (n_rows - 1);
-                    in -= n_rows * (n_rows + 1) - 1;
-                  }
-                if (normal_dir == 1)
-                  {
-                    out += (n_rows - 1) * (n_rows - 1);
-                    in -= (n_rows - 1) * (n_rows - 1) - 1;
-                  }
-                if (normal_dir == 2)
-                  {
-                    out += (n_rows - 1) * (n_rows);
-                    in -= (n_rows) * (n_rows + 1) - 1;
-                  }
-              }
-          }
-      }
-  }
-
 
 
   /**

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.