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)
{
}
}
}
-
- /* 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);
- }
- }
};
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,
}
}
- 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;
- }
- }
- }
- }
- }
-
/**
projects / dealii.git / commitdiff