]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Delete obsolete code
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Wed, 30 Aug 2023 08:38:57 +0000 (10:38 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Mon, 11 Sep 2023 19:58:51 +0000 (21:58 +0200)
include/deal.II/matrix_free/tensor_product_kernels.h

index 11aa7fefc738d5eeea65d63a7ae93e229004909c..fc9c8429458f43eca3fd8e7a53fb88dbd12714d3 100644 (file)
@@ -1130,13 +1130,8 @@ namespace internal
         }
       else
         {
-          // 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.empty() ||
-                   shape_values.size() == n_rows * n_columns ||
-                   shape_values.size() == 3 * n_rows,
+                   shape_values.size() == n_rows * n_columns,
                  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
           Assert(shape_gradients.empty() ||
                    shape_gradients.size() == n_rows * n_columns,
@@ -1268,43 +1263,6 @@ namespace internal
           const Number                   *in,
           Number                         *out);
 
-    /**
-     * This function applies the tensor product operation to produce face values
-     * from cell values. As opposed to the apply method, this method assumes
-     * that the directions orthogonal to the face have n_rows degrees of
-     * freedom per direction and not n_columns for those directions lower than
-     * the one currently applied. In other words, apply_face() must be called
-     * before calling any interpolation within the face.
-     *
-     * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc)
-     * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
-     *                            and interpolation into n_rows^(dim-1) points
-     *                            is performed. This is a typical scenario in
-     *                            FEFaceEvaluation::evaluate() calls. If false,
-     *                            data from n_rows^(dim-1) points is expanded
-     *                            into the n_rows^dim points of the higher-
-     *                            dimensional data array. Derivatives in the
-     *                            case contract_onto_face==false are summed
-     *                            together
-     * @tparam add If true, the result is added to the output vector, else
-     *             the computed values overwrite the content in the output
-     * @tparam max_derivative Sets the number of derivatives that should be
-     *             computed. 0 means only values, 1 means values and first
-     *             derivatives, 2 second derivates. Note that all the
-     *             derivatives access the data in @p shape_values passed to
-     *             the constructor of the class
-     *
-     * @param in address of the input data vector
-     * @param out address of the output data vector
-     */
-    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;
@@ -1380,174 +1338,6 @@ namespace internal
 
 
 
-  template <int  n_rows_template,
-            int  stride_template,
-            bool contract_onto_face,
-            bool add,
-            int  max_derivative,
-            typename Number,
-            typename Number2>
-  inline std::enable_if_t<contract_onto_face, void>
-  interpolate_to_face(const Number2            *shape_values,
-                      const std::array<int, 2> &n_blocks,
-                      const std::array<int, 2> &steps,
-                      const Number             *input,
-                      Number *DEAL_II_RESTRICT  output,
-                      const int                 n_rows_runtime = 0,
-                      const int                 stride_runtime = 1)
-  {
-    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
-    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
-
-    Number *output1 = output + n_blocks[0] * n_blocks[1];
-    Number *output2 = output1 + n_blocks[0] * n_blocks[1];
-    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
-          {
-            Number res0 = shape_values[0] * input[0];
-            Number res1, res2;
-            if (max_derivative > 0)
-              res1 = shape_values[n_rows] * input[0];
-            if (max_derivative > 1)
-              res2 = shape_values[2 * n_rows] * input[0];
-            for (int ind = 1; ind < n_rows; ++ind)
-              {
-                res0 += shape_values[ind] * input[stride * ind];
-                if (max_derivative > 0)
-                  res1 += shape_values[ind + n_rows] * input[stride * ind];
-                if (max_derivative > 1)
-                  res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
-              }
-            if (add)
-              {
-                output[i1] += res0;
-                if (max_derivative > 0)
-                  output1[i1] += res1;
-                if (max_derivative > 1)
-                  output2[i2] += res2;
-              }
-            else
-              {
-                output[i1] = res0;
-                if (max_derivative > 0)
-                  output1[i1] = res1;
-                if (max_derivative > 1)
-                  output2[i1] = res2;
-              }
-            input += steps[0];
-          }
-        output += n_blocks[0];
-        if (max_derivative > 0)
-          output1 += n_blocks[0];
-        if (max_derivative > 1)
-          output2 += n_blocks[0];
-        input += steps[1];
-      }
-  }
-
-
-
-  template <int  n_rows_template,
-            int  stride_template,
-            bool contract_onto_face,
-            bool add,
-            int  max_derivative,
-            typename Number,
-            typename Number2>
-  inline std::enable_if_t<!contract_onto_face, void>
-  interpolate_to_face(const Number2            *shape_values,
-                      const std::array<int, 2> &n_blocks,
-                      const std::array<int, 2> &steps,
-                      const Number             *input,
-                      Number *DEAL_II_RESTRICT  output,
-                      const int                 n_rows_runtime = 0,
-                      const int                 stride_runtime = 1)
-  {
-    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
-    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
-
-    const Number *input1 = input + n_blocks[0] * n_blocks[1];
-    const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
-    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
-      {
-        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
-          {
-            const Number in = input[i1];
-            Number       in1, in2;
-            if (max_derivative > 0)
-              in1 = input1[i1];
-            if (max_derivative > 1)
-              in2 = input2[i1];
-            for (int col = 0; col < n_rows; ++col)
-              {
-                Number result =
-                  add ? (output[col * stride] + shape_values[col] * in) :
-                        (shape_values[col] * in);
-                if (max_derivative > 0)
-                  result += shape_values[col + n_rows] * in1;
-                if (max_derivative > 1)
-                  result += shape_values[col + 2 * n_rows] * in2;
-
-                output[col * stride] = result;
-              }
-            output += steps[0];
-          }
-        input += n_blocks[0];
-        if (max_derivative > 0)
-          input1 += n_blocks[0];
-        if (max_derivative > 1)
-          input2 += n_blocks[0];
-        output += steps[1];
-      }
-  }
-
-
-
-  template <EvaluatorVariant variant,
-            int              dim,
-            int              n_rows,
-            int              n_columns,
-            typename Number,
-            typename Number2>
-  template <int  face_direction,
-            bool contract_to_face,
-            bool add,
-            int  max_derivative>
-  inline void
-  EvaluatorTensorProduct<variant, dim, n_rows, n_columns, Number, Number2>::
-    apply_face(const Number *DEAL_II_RESTRICT in,
-               Number *DEAL_II_RESTRICT       out) const
-  {
-    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      stride = Utilities::pow(n_rows, face_direction);
-    std::array<int, 2> steps;
-    if constexpr (face_direction == 0)
-      steps = {{n_rows, 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_rows, -n_rows * n_rows * n_rows + 1}};
-    else if constexpr (face_direction == 2)
-      steps = {{1, 0}};
-
-    interpolate_to_face<n_rows, stride, contract_to_face, add, max_derivative>(
-      this->shape_values,
-      {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}},
-      steps,
-      in,
-      out);
-  }
-
-
-
   /**
    * Internal evaluator for shape function using the tensor product form
    * of the basis functions. The same as the other templated class but
@@ -1612,13 +1402,8 @@ namespace internal
         }
       else
         {
-          // 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.empty() ||
-                   shape_values.size() == n_rows * n_columns ||
-                   shape_values.size() == n_rows * 3,
+                   shape_values.size() == n_rows * n_columns,
                  ExcDimensionMismatch(shape_values.size(), n_rows * n_columns));
           Assert(shape_gradients.empty() ||
                    shape_gradients.size() == n_rows * n_columns,
@@ -1721,14 +1506,6 @@ namespace internal
           const Number                   *in,
           Number                         *out) const;
 
-    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;
-
     const Number2     *shape_values;
     const Number2     *shape_gradients;
     const Number2     *shape_hessians;
@@ -1812,45 +1589,165 @@ namespace internal
 
 
 
-  template <EvaluatorVariant variant,
-            int              dim,
+  /**
+   * This function applies the tensor product operation to produce face
+   * values from cell values. The algorithm involved here can be interpreted
+   * the first sweep in sum factorization, reducing the dimensionality of
+   * the data set from dim-dimensional cell values to (dim-1)-dimensional
+   * face values. This step is always done before we evaluate within the
+   * face, as it reduces the dimensionality.
+   *
+   * @tparam n_rows_template The number of entries within the interpolation,
+   *             typically equal to the polynomial degree plus one, if known
+   *             at compile time, otherwise n_rows_runtime is used
+   * @tparam stride_template The stride between successive entries in the
+   *             one-dimensional operation of sum factorization, if known at
+   *             compile time, otherwise stride_runtime is used
+   * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
+   *                            and interpolation into n_rows^(dim-1) points
+   *                            is performed. This is a typical scenario in
+   *                            FEFaceEvaluation::evaluate() calls. If false,
+   *                            data from n_rows^(dim-1) points is expanded
+   *                            into the n_rows^dim points of the higher-
+   *                            dimensional data array. Derivatives in the
+   *                            case contract_onto_face==false are summed
+   *                            together
+   * @tparam add If true, the result is added to the output vector, else
+   *             the computed values overwrite the content in the output
+   * @tparam max_derivative Sets the number of derivatives that should be
+   *             computed. 0 means only values, 1 means values and first
+   *             derivatives, 2 second derivates. Note that all the
+   *             derivatives access the data in @p shape_values passed to
+   *             the constructor of the class
+   *
+   * @param shape_values address of the interpolation matrix
+   * @param n_blocks Number of interpolation layer used along the two other
+   *             dimensions tangential to the interpolation direction
+   * @param steps Increments in the input array from one step to the next,
+   *             varied in conjunction with the @p stride variable.
+   * @param input Address of the input data vector
+   * @param output Address of the output data vector
+   */
+  template <int  n_rows_template,
+            int  stride_template,
+            bool contract_onto_face,
+            bool add,
+            int  max_derivative,
             typename Number,
             typename Number2>
-  template <int  face_direction,
-            bool contract_to_face,
+  inline std::enable_if_t<contract_onto_face, void>
+  interpolate_to_face(const Number2            *shape_values,
+                      const std::array<int, 2> &n_blocks,
+                      const std::array<int, 2> &steps,
+                      const Number             *input,
+                      Number *DEAL_II_RESTRICT  output,
+                      const int                 n_rows_runtime = 0,
+                      const int                 stride_runtime = 1)
+  {
+    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
+    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
+
+    Number *output1 = output + n_blocks[0] * n_blocks[1];
+    Number *output2 = output1 + n_blocks[0] * n_blocks[1];
+    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
+          {
+            Number res0 = shape_values[0] * input[0];
+            Number res1, res2;
+            if (max_derivative > 0)
+              res1 = shape_values[n_rows] * input[0];
+            if (max_derivative > 1)
+              res2 = shape_values[2 * n_rows] * input[0];
+            for (int ind = 1; ind < n_rows; ++ind)
+              {
+                res0 += shape_values[ind] * input[stride * ind];
+                if (max_derivative > 0)
+                  res1 += shape_values[ind + n_rows] * input[stride * ind];
+                if (max_derivative > 1)
+                  res2 += shape_values[ind + 2 * n_rows] * input[stride * ind];
+              }
+            if (add)
+              {
+                output[i1] += res0;
+                if (max_derivative > 0)
+                  output1[i1] += res1;
+                if (max_derivative > 1)
+                  output2[i2] += res2;
+              }
+            else
+              {
+                output[i1] = res0;
+                if (max_derivative > 0)
+                  output1[i1] = res1;
+                if (max_derivative > 1)
+                  output2[i1] = res2;
+              }
+            input += steps[0];
+          }
+        output += n_blocks[0];
+        if (max_derivative > 0)
+          output1 += n_blocks[0];
+        if (max_derivative > 1)
+          output2 += n_blocks[0];
+        input += steps[1];
+      }
+  }
+
+
+
+  template <int  n_rows_template,
+            int  stride_template,
+            bool contract_onto_face,
             bool add,
-            int  max_derivative>
-  inline void
-  EvaluatorTensorProduct<variant, dim, 0, 0, Number, Number2>::apply_face(
-    const Number *DEAL_II_RESTRICT in,
-    Number *DEAL_II_RESTRICT       out) const
+            int  max_derivative,
+            typename Number,
+            typename Number2>
+  inline std::enable_if_t<!contract_onto_face, void>
+  interpolate_to_face(const Number2            *shape_values,
+                      const std::array<int, 2> &n_blocks,
+                      const std::array<int, 2> &steps,
+                      const Number             *input,
+                      Number *DEAL_II_RESTRICT  output,
+                      const int                 n_rows_runtime = 0,
+                      const int                 stride_runtime = 1)
   {
-    Assert(shape_values != nullptr,
-           ExcMessage(
-             "The given array shape_data must not be the null pointer!"));
-    static_assert(dim > 0 && dim < 4, "Only dim=1,2,3 supported");
-
-    const int          stride = Utilities::pow(n_rows, face_direction);
-    const int          n_rows = this->n_rows;
-    std::array<int, 2> steps;
-    if constexpr (face_direction == 0)
-      steps = {{n_rows, 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_rows, -n_rows * n_rows * n_rows + 1}};
-    else if constexpr (face_direction == 2)
-      steps = {{1, 0}};
-
-    interpolate_to_face<0, 0, contract_to_face, add, max_derivative>(
-      this->shape_values,
-      {{(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}},
-      steps,
-      in,
-      out,
-      n_rows,
-      stride);
+    const int n_rows = n_rows_template > 0 ? n_rows_template : n_rows_runtime;
+    const int stride = n_rows_template > 0 ? stride_template : stride_runtime;
+
+    const Number *input1 = input + n_blocks[0] * n_blocks[1];
+    const Number *input2 = input1 + n_blocks[0] * n_blocks[1];
+    for (int i2 = 0; i2 < n_blocks[1]; ++i2)
+      {
+        for (int i1 = 0; i1 < n_blocks[0]; ++i1)
+          {
+            const Number in = input[i1];
+            Number       in1, in2;
+            if (max_derivative > 0)
+              in1 = input1[i1];
+            if (max_derivative > 1)
+              in2 = input2[i1];
+            for (int col = 0; col < n_rows; ++col)
+              {
+                Number result =
+                  add ? (output[col * stride] + shape_values[col] * in) :
+                        (shape_values[col] * in);
+                if (max_derivative > 0)
+                  result += shape_values[col + n_rows] * in1;
+                if (max_derivative > 1)
+                  result += shape_values[col + 2 * n_rows] * in2;
+
+                output[col * stride] = result;
+              }
+            output += steps[0];
+          }
+        input += n_blocks[0];
+        if (max_derivative > 0)
+          input1 += n_blocks[0];
+        if (max_derivative > 1)
+          input2 += n_blocks[0];
+        output += steps[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.