}
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,
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;
- 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
}
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,
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;
- 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];
+ }
}
projects / dealii.git / commitdiff