const bool add_into_values_array);
};
- /**
- * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
- * specific sum-factorization kernels and with normal/tangential shape_data
- */
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>
- {
- template <bool integrate>
- static void
- evaluate_or_integrate(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array = false);
-
- private:
- template <typename EvalType, typename Number2>
- static EvalType
- create_evaluator_tensor_product(
- const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape_data)
- {
- return EvalType(shape_data.shape_values,
- shape_data.shape_gradients,
- shape_data.shape_hessians);
- }
-
- template <int normal_dir>
- static void
- evaluate_tensor_product_per_component(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array,
- std::integral_constant<bool, false>);
-
- template <int normal_dir>
- static void
- evaluate_tensor_product_per_component(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array,
- std::integral_constant<bool, true>);
- };
-
template <MatrixFreeFunctions::ElementType type,
Number *values_quad = fe_eval.begin_values();
Number *gradients_quad = fe_eval.begin_gradients();
- Number *hessians_quad = fe_eval.begin_hessians();
switch (dim)
{
if (evaluation_flag & EvaluationFlags::gradients)
eval0.template gradients<0, true, false>(values_dofs,
gradients_quad);
- if (evaluation_flag & EvaluationFlags::hessians)
- eval0.template hessians<0, true, false>(values_dofs,
- hessians_quad);
// advance the next component in 1d array
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += n_q_points;
- hessians_quad += n_q_points;
}
break;
if (evaluation_flag & EvaluationFlags::gradients)
{
eval0.template gradients<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, gradients_quad);
- }
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad xy
- if (!(evaluation_flag & EvaluationFlags::gradients))
- eval0.template gradients<0, true, false>(values_dofs,
- temp1);
- eval1.template gradients<1, true, false>(temp1,
- hessians_quad +
- 2 * n_q_points);
-
- // grad xx
- eval0.template hessians<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, hessians_quad);
+ eval1.template values<1, true, false, 2>(temp1,
+ gradients_quad);
}
// grad y
eval0.template values<0, true, false>(values_dofs, temp1);
if (evaluation_flag & EvaluationFlags::gradients)
- eval1.template gradients<1, true, false>(temp1,
- gradients_quad +
- n_q_points);
-
- // grad yy
- if (evaluation_flag & EvaluationFlags::hessians)
- eval1.template hessians<1, true, false>(temp1,
- hessians_quad +
- n_q_points);
+ eval1.template gradients<1, true, false, 2>(temp1,
+ gradients_quad + 1);
// val: can use values applied in x
if (evaluation_flag & EvaluationFlags::values)
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += 2 * n_q_points;
- hessians_quad += 3 * n_q_points;
}
break;
// grad x
eval0.template gradients<0, true, false>(values_dofs, temp1);
eval1.template values<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2, gradients_quad);
- }
-
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad xz
- if (!(evaluation_flag & EvaluationFlags::gradients))
- {
- eval0.template gradients<0, true, false>(values_dofs,
- temp1);
- eval1.template values<1, true, false>(temp1, temp2);
- }
- eval2.template gradients<2, true, false>(temp2,
- hessians_quad +
- 4 * n_q_points);
-
- // grad xy
- eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- hessians_quad +
- 3 * n_q_points);
-
- // grad xx
- eval0.template hessians<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2, hessians_quad);
+ eval2.template values<2, true, false, 3>(temp2,
+ gradients_quad);
}
// grad y
if (evaluation_flag & EvaluationFlags::gradients)
{
eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- gradients_quad +
- n_q_points);
- }
-
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad yz
- if (!(evaluation_flag & EvaluationFlags::gradients))
- eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template gradients<2, true, false>(temp2,
- hessians_quad +
- 5 * n_q_points);
-
- // grad yy
- eval1.template hessians<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- hessians_quad +
- n_q_points);
+ eval2.template values<2, true, false, 3>(temp2,
+ gradients_quad + 1);
}
// grad z: can use the values applied in x direction stored in
// temp1
eval1.template values<1, true, false>(temp1, temp2);
if (evaluation_flag & EvaluationFlags::gradients)
- eval2.template gradients<2, true, false>(temp2,
- gradients_quad +
- 2 * n_q_points);
-
- // grad zz: can use the values applied in x and y direction stored
- // in temp2
- if (evaluation_flag & EvaluationFlags::hessians)
- eval2.template hessians<2, true, false>(temp2,
- hessians_quad +
- 2 * n_q_points);
+ eval2.template gradients<2, true, false, 3>(temp2,
+ gradients_quad + 2);
// val: can use the values applied in x & y direction stored in
// temp2
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += 3 * n_q_points;
- hessians_quad += 6 * n_q_points;
}
break;
Number *values_quad = fe_eval.begin_values();
Number *gradients_quad = fe_eval.begin_gradients();
- Number *hessians_quad = fe_eval.begin_hessians();
switch (dim)
{
eval0.template gradients<0, false, false>(gradients_quad,
values_dofs);
}
- if ((integration_flag & EvaluationFlags::hessians) != 0u)
- {
- if ((integration_flag & EvaluationFlags::values) != 0u ||
- (integration_flag & EvaluationFlags::gradients) != 0u ||
- add_into_values_array == true)
- eval0.template hessians<0, false, true>(hessians_quad,
- values_dofs);
- else
- eval0.template hessians<0, false, false>(hessians_quad,
- values_dofs);
- }
// advance to the next component in 1d array
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += n_q_points;
- hessians_quad += n_q_points;
}
break;
}
if (integration_flag & EvaluationFlags::gradients)
{
- eval1.template gradients<1, false, false>(gradients_quad +
- n_q_points,
- temp1);
+ eval1.template gradients<1, false, false, 2>(gradients_quad +
+ 1,
+ temp1);
if (integration_flag & EvaluationFlags::values)
eval1.template values<1, false, true>(values_quad, temp1);
if (add_into_values_array == false)
eval0.template values<0, false, false>(temp1, values_dofs);
else
eval0.template values<0, false, true>(temp1, values_dofs);
- eval1.template values<1, false, false>(gradients_quad, temp1);
- eval0.template gradients<0, false, true>(temp1, values_dofs);
- }
- if ((integration_flag & EvaluationFlags::hessians) != 0u)
- {
- // grad xx
- eval1.template values<1, false, false>(hessians_quad, temp1);
-
- if ((integration_flag & EvaluationFlags::values) != 0u ||
- (integration_flag & EvaluationFlags::gradients) != 0u ||
- add_into_values_array == true)
- eval0.template hessians<0, false, true>(temp1, values_dofs);
- else
- eval0.template hessians<0, false, false>(temp1,
- values_dofs);
-
- // grad yy
- eval1.template hessians<1, false, false>(hessians_quad +
- n_q_points,
- temp1);
- eval0.template values<0, false, true>(temp1, values_dofs);
-
- // grad xy
- eval1.template gradients<1, false, false>(hessians_quad +
- 2 * n_q_points,
+ eval1.template values<1, false, false, 2>(gradients_quad,
temp1);
eval0.template gradients<0, false, true>(temp1, values_dofs);
}
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += 2 * n_q_points;
- hessians_quad += 3 * n_q_points;
}
break;
}
if (integration_flag & EvaluationFlags::gradients)
{
- eval2.template gradients<2, false, false>(gradients_quad +
- 2 * n_q_points,
- temp1);
+ eval2.template gradients<2, false, false, 3>(gradients_quad +
+ 2,
+ temp1);
if (integration_flag & EvaluationFlags::values)
eval2.template values<2, false, true>(values_quad, temp1);
eval1.template values<1, false, false>(temp1, temp2);
- eval2.template values<2, false, false>(gradients_quad +
- n_q_points,
- temp1);
+ eval2.template values<2, false, false, 3>(gradients_quad + 1,
+ temp1);
eval1.template gradients<1, false, true>(temp1, temp2);
if (add_into_values_array == false)
eval0.template values<0, false, false>(temp2, values_dofs);
else
eval0.template values<0, false, true>(temp2, values_dofs);
- eval2.template values<2, false, false>(gradients_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval0.template gradients<0, false, true>(temp2, values_dofs);
- }
- if ((integration_flag & EvaluationFlags::hessians) != 0u)
- {
- // grad xx
- eval2.template values<2, false, false>(hessians_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
-
- if ((integration_flag & EvaluationFlags::values) != 0u ||
- (integration_flag & EvaluationFlags::gradients) != 0u ||
- add_into_values_array == true)
- eval0.template hessians<0, false, true>(temp2, values_dofs);
- else
- eval0.template hessians<0, false, false>(temp2,
- values_dofs);
-
- // grad yy
- eval2.template values<2, false, false>(hessians_quad +
- n_q_points,
- temp1);
- eval1.template hessians<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
-
- // grad zz
- eval2.template hessians<2, false, false>(hessians_quad +
- 2 * n_q_points,
- temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
-
- // grad xy
- eval2.template values<2, false, false>(hessians_quad +
- 3 * n_q_points,
- temp1);
- eval1.template gradients<1, false, false>(temp1, temp2);
- eval0.template gradients<0, false, true>(temp2, values_dofs);
-
- // grad xz
- eval2.template gradients<2, false, false>(hessians_quad +
- 4 * n_q_points,
+ eval2.template values<2, false, false, 3>(gradients_quad,
temp1);
eval1.template values<1, false, false>(temp1, temp2);
eval0.template gradients<0, false, true>(temp2, values_dofs);
-
- // grad yz
- eval2.template gradients<2, false, false>(hessians_quad +
- 5 * n_q_points,
- temp1);
- eval1.template gradients<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
}
// advance to the next component in 1d array
values_dofs += dofs_per_comp;
values_quad += n_q_points;
gradients_quad += 3 * n_q_points;
- hessians_quad += 6 * n_q_points;
}
break;
n_dofs,
n_q_points);
- eval.template gradients<0, true, false>(values_dofs_actual_ptr,
- gradients_quad_ptr);
-
- gradients_quad_ptr += n_q_points;
+ eval.template gradients<0, true, false, dim>(
+ values_dofs_actual_ptr, gradients_quad_ptr + d);
}
+ gradients_quad_ptr += n_q_points * dim;
values_dofs_actual_ptr += n_dofs;
}
}
if ((add_into_values_array == false &&
!(integration_flag & EvaluationFlags::values)) &&
d == 0)
- eval.template gradients<0, false, false>(
- gradients_quad_ptr, values_dofs_actual_ptr);
+ eval.template gradients<0, false, false, dim>(
+ gradients_quad_ptr + d, values_dofs_actual_ptr);
else
- eval.template gradients<0, false, true>(
- gradients_quad_ptr, values_dofs_actual_ptr);
-
- gradients_quad_ptr += n_q_points;
+ eval.template gradients<0, false, true, dim>(
+ gradients_quad_ptr + d, values_dofs_actual_ptr);
}
+ gradients_quad_ptr += n_q_points * dim;
values_dofs_actual_ptr += n_dofs;
}
}
}
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <bool integrate>
- inline void
- FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::
- evaluate_or_integrate(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array)
- {
- if (evaluation_flag == EvaluationFlags::nothing)
- return;
-
- AssertDimension(fe_eval.get_shape_info().data.size(), 2);
- // First component:
- evaluate_tensor_product_per_component<0>(
- evaluation_flag,
- values_dofs_actual,
- fe_eval,
- add_into_values_array,
- std::integral_constant<bool, integrate>());
- // Second component :
- evaluate_tensor_product_per_component<1>(
- evaluation_flag,
- values_dofs_actual,
- fe_eval,
- add_into_values_array,
- std::integral_constant<bool, integrate>());
- if (dim == 3)
- {
- // Third component
- evaluate_tensor_product_per_component<2>(
- evaluation_flag,
- values_dofs_actual,
- fe_eval,
- add_into_values_array,
- std::integral_constant<bool, integrate>());
- }
- }
-
- // Helper function that applies the 1d evaluation kernels.
- // std::integral_constant<bool, false> is the interpolation path, and
- // std::integral_constant<bool, true> below is the integration path.
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int normal_dir>
- inline void
- FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::
- evaluate_tensor_product_per_component(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array,
- std::integral_constant<bool, false>)
- {
- (void)add_into_values_array;
-
- using Number2 =
- typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
- using EvalNormal =
- EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
- dim,
- (fe_degree == -1) ? 1 : fe_degree + 1,
- n_q_points_1d,
- normal_dir,
- Number,
- Number2>;
-
- using EvalTangent =
- EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
- dim,
- (fe_degree == -1) ? 1 : fe_degree,
- n_q_points_1d,
- normal_dir,
- Number,
- Number2>;
- using Eval0 = std::conditional_t<normal_dir == 0, EvalNormal, EvalTangent>;
- using Eval1 = std::conditional_t<normal_dir == 1, EvalNormal, EvalTangent>;
- using Eval2 = std::conditional_t<normal_dir == 2, EvalNormal, EvalTangent>;
-
- const auto &shape_info = fe_eval.get_shape_info();
- Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
- ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
- Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
- ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
- Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
- ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
-
- Number *temp1 = fe_eval.get_scratch_data().begin();
- Number *temp2;
-
- temp2 =
- temp1 +
- std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
- Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
-
- const std::size_t n_q_points = shape_info.n_q_points;
- const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
-
- // Initial shift depending on component (normal_dir)
- Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
- Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
- Number *gradients_quad =
- fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
- Number *hessians_quad =
- (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
- fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
-
- switch (dim)
- {
- case 2:
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- eval0.template gradients<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, gradients_quad);
- }
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // The evaluation/integration here *should* work, however
- // the piola transform is not implemented.
- AssertThrow(false, ExcNotImplemented());
- // grad xy
- if (!(evaluation_flag & EvaluationFlags::gradients))
- eval0.template gradients<0, true, false>(values_dofs, temp1);
- eval1.template gradients<1, true, false>(temp1,
- hessians_quad +
- 2 * n_q_points);
-
- // grad xx
- eval0.template hessians<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, hessians_quad);
- }
-
- // grad y
- eval0.template values<0, true, false>(values_dofs, temp1);
- if (evaluation_flag & EvaluationFlags::gradients)
- eval1.template gradients<1, true, false>(temp1,
- gradients_quad +
- n_q_points);
-
- // grad yy
- if (evaluation_flag & EvaluationFlags::hessians)
- eval1.template hessians<1, true, false>(temp1,
- hessians_quad + n_q_points);
-
- // val: can use values applied in x
- if (evaluation_flag & EvaluationFlags::values)
- eval1.template values<1, true, false>(temp1, values_quad);
- break;
- case 3:
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- // grad x
- eval0.template gradients<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2, gradients_quad);
- }
-
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // The evaluation/integration here *should* work, however
- // the piola transform is not implemented.
- AssertThrow(false, ExcNotImplemented());
- // grad xz
- if (!(evaluation_flag & EvaluationFlags::gradients))
- {
- eval0.template gradients<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, temp2);
- }
- eval2.template gradients<2, true, false>(temp2,
- hessians_quad +
- 4 * n_q_points);
-
- // grad xy
- eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- hessians_quad +
- 3 * n_q_points);
-
- // grad xx
- eval0.template hessians<0, true, false>(values_dofs, temp1);
- eval1.template values<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2, hessians_quad);
- }
-
- // grad y
- eval0.template values<0, true, false>(values_dofs, temp1);
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- gradients_quad +
- n_q_points);
- }
-
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad yz
- if (!(evaluation_flag & EvaluationFlags::gradients))
- eval1.template gradients<1, true, false>(temp1, temp2);
- eval2.template gradients<2, true, false>(temp2,
- hessians_quad +
- 5 * n_q_points);
-
- // grad yy
- eval1.template hessians<1, true, false>(temp1, temp2);
- eval2.template values<2, true, false>(temp2,
- hessians_quad + n_q_points);
- }
-
- // grad z: can use the values applied in x direction stored in
- // temp1
- eval1.template values<1, true, false>(temp1, temp2);
- if (evaluation_flag & EvaluationFlags::gradients)
- eval2.template gradients<2, true, false>(temp2,
- gradients_quad +
- 2 * n_q_points);
-
- // grad zz: can use the values applied in x and y direction stored
- // in temp2
- if (evaluation_flag & EvaluationFlags::hessians)
- eval2.template hessians<2, true, false>(temp2,
- hessians_quad +
- 2 * n_q_points);
-
- // val: can use the values applied in x & y direction stored in
- // temp2
- if (evaluation_flag & EvaluationFlags::values)
- eval2.template values<2, true, false>(temp2, values_quad);
- break;
- default:
- AssertThrow(false, ExcNotImplemented());
- }
- }
-
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int normal_dir>
- inline void
- FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::
- evaluate_tensor_product_per_component(
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- Number *values_dofs_actual,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array,
- std::integral_constant<bool, true>)
- {
- using Number2 =
- typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
- using EvalNormal =
- EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
- dim,
- (fe_degree == -1) ? 1 : fe_degree + 1,
- n_q_points_1d,
- normal_dir,
- Number,
- Number2>;
-
- using EvalTangent =
- EvaluatorTensorProductAnisotropic<evaluate_raviart_thomas,
- dim,
- (fe_degree == -1) ? 1 : fe_degree,
- n_q_points_1d,
- normal_dir,
- Number,
- Number2>;
- using Eval0 = std::conditional_t<normal_dir == 0, EvalNormal, EvalTangent>;
- using Eval1 = std::conditional_t<normal_dir == 1, EvalNormal, EvalTangent>;
- using Eval2 = std::conditional_t<normal_dir == 2, EvalNormal, EvalTangent>;
-
- const auto &shape_info = fe_eval.get_shape_info();
- Eval0 eval0 = create_evaluator_tensor_product<Eval0>(
- ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]));
- Eval1 eval1 = create_evaluator_tensor_product<Eval1>(
- ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1]));
- Eval2 eval2 = create_evaluator_tensor_product<Eval2>(
- ((normal_dir == 2) ? shape_info.data[0] : shape_info.data[1]));
-
- Number *temp1 = fe_eval.get_scratch_data().begin();
- Number *temp2;
-
- temp2 =
- temp1 +
- std::max(Utilities::fixed_power<dim>(shape_info.data[0].fe_degree + 1),
- Utilities::fixed_power<dim>(shape_info.data[0].n_q_points_1d));
-
- const std::size_t n_q_points = shape_info.n_q_points;
- const std::size_t dofs_per_comp = shape_info.dofs_per_component_on_cell;
-
- // Initial shift depending on component (normal_dir)
- Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir;
- Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir;
- Number *gradients_quad =
- fe_eval.begin_gradients() + dim * n_q_points * normal_dir;
- Number *hessians_quad =
- (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir :
- fe_eval.begin_hessians() + 6 * n_q_points * normal_dir;
-
- // Integrate path
- switch (dim)
- {
- case 2:
- if ((evaluation_flag & EvaluationFlags::values) &&
- !(evaluation_flag & EvaluationFlags::gradients))
- {
- eval1.template values<1, false, false>(values_quad, temp1);
- if (add_into_values_array == false)
- eval0.template values<0, false, false>(temp1, values_dofs);
- else
- eval0.template values<0, false, true>(temp1, values_dofs);
- }
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- eval1.template gradients<1, false, false>(gradients_quad +
- n_q_points,
- temp1);
- if ((evaluation_flag & EvaluationFlags::values))
- eval1.template values<1, false, true>(values_quad, temp1);
- if (add_into_values_array == false)
- eval0.template values<0, false, false>(temp1, values_dofs);
- else
- eval0.template values<0, false, true>(temp1, values_dofs);
- eval1.template values<1, false, false>(gradients_quad, temp1);
- eval0.template gradients<0, false, true>(temp1, values_dofs);
- }
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad xx
- eval1.template values<1, false, false>(hessians_quad, temp1);
-
- if ((evaluation_flag & EvaluationFlags::values) ||
- (evaluation_flag & EvaluationFlags::gradients) ||
- add_into_values_array == true)
- eval0.template hessians<0, false, true>(temp1, values_dofs);
- else
- eval0.template hessians<0, false, false>(temp1, values_dofs);
-
- // grad yy
- eval1.template hessians<1, false, false>(hessians_quad +
- n_q_points,
- temp1);
- eval0.template values<0, false, true>(temp1, values_dofs);
-
- // grad xy
- eval1.template gradients<1, false, false>(hessians_quad +
- 2 * n_q_points,
- temp1);
- eval0.template gradients<0, false, true>(temp1, values_dofs);
- }
- break;
-
- case 3:
- if ((evaluation_flag & EvaluationFlags::values) &&
- !(evaluation_flag & EvaluationFlags::gradients))
- {
- eval2.template values<2, false, false>(values_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- if (add_into_values_array == false)
- eval0.template values<0, false, false>(temp2, values_dofs);
- else
- eval0.template values<0, false, true>(temp2, values_dofs);
- }
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- eval2.template gradients<2, false, false>(gradients_quad +
- 2 * n_q_points,
- temp1);
- if ((evaluation_flag & EvaluationFlags::values))
- eval2.template values<2, false, true>(values_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval2.template values<2, false, false>(gradients_quad +
- n_q_points,
- temp1);
- eval1.template gradients<1, false, true>(temp1, temp2);
- if (add_into_values_array == false)
- eval0.template values<0, false, false>(temp2, values_dofs);
- else
- eval0.template values<0, false, true>(temp2, values_dofs);
- eval2.template values<2, false, false>(gradients_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval0.template gradients<0, false, true>(temp2, values_dofs);
- }
- if (evaluation_flag & EvaluationFlags::hessians)
- {
- // grad xx
- eval2.template values<2, false, false>(hessians_quad, temp1);
- eval1.template values<1, false, false>(temp1, temp2);
-
- if ((evaluation_flag & EvaluationFlags::values) ||
- (evaluation_flag & EvaluationFlags::gradients) ||
- add_into_values_array == true)
- eval0.template hessians<0, false, true>(temp2, values_dofs);
- else
- eval0.template hessians<0, false, false>(temp2, values_dofs);
-
- // grad yy
- eval2.template values<2, false, false>(hessians_quad + n_q_points,
- temp1);
- eval1.template hessians<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
-
- // grad zz
- eval2.template hessians<2, false, false>(hessians_quad +
- 2 * n_q_points,
- temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
-
- // grad xy
- eval2.template values<2, false, false>(hessians_quad +
- 3 * n_q_points,
- temp1);
- eval1.template gradients<1, false, false>(temp1, temp2);
- eval0.template gradients<0, false, true>(temp2, values_dofs);
-
- // grad xz
- eval2.template gradients<2, false, false>(hessians_quad +
- 4 * n_q_points,
- temp1);
- eval1.template values<1, false, false>(temp1, temp2);
- eval0.template gradients<0, false, true>(temp2, values_dofs);
-
- // grad yz
- eval2.template gradients<2, false, false>(hessians_quad +
- 5 * n_q_points,
- temp1);
- eval1.template gradients<1, false, false>(temp1, temp2);
- eval0.template values<0, false, true>(temp2, values_dofs);
- }
-
- break;
- default:
- AssertThrow(false, ExcNotImplemented());
- }
- }
/**
* This struct implements the change between two different bases. This is an
*
* This class allows for dimension-independent application of the operation,
* implemented by template recursion. It has been tested up to 6d.
- *
- * The last two template arguments in this class are unused. They have a
- * default type set to bool, but any type will compile. These arguments are
- * present for backward compatibility of this (internal) interface with
- * previous versions of deal.II, where the template arguments were used to
- * indicate the number types These are now part of the do_forward() and
- * do_backward() functions, and they will be removed from a future version
- * of deal.II.
*/
template <EvaluatorVariant variant,
EvaluatorQuantity quantity,
int dim,
int basis_size_1,
- int basis_size_2,
- typename = bool,
- typename = bool>
+ int basis_size_2>
struct FEEvaluationImplBasisChange
{
static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
/**
- * This struct performs the evaluation of function values, gradients and
- * Hessians for tensor-product finite elements. This a specialization for
- * elements where the nodal points coincide with the quadrature points like
- * FE_Q shape functions on Gauss-Lobatto elements integrated with
- * Gauss-Lobatto quadrature. The assumption of this class is that the shape
- * 'values' operation is identity, which allows us to write shorter code.
- *
- * In literature, this form of evaluation is often called spectral
- * evaluation, spectral collocation or simply collocation, meaning the same
- * location for shape functions and evaluation space (quadrature points).
+ * Internal function that evaluates the gradients of finite element
+ * functions represented by bases in the collocation space, used by
+ * FEEvaluationImplCollocation and FEEvaluationImplTransformToCollocation.
+ * The evaluation strategy uses sum factorization with the even-odd
+ * optimization and fixed loop bounds.
*/
- template <int dim, int fe_degree, typename Number>
- struct FEEvaluationImplCollocation
+ template <int n_points_1d, int dim, typename Number, typename Number2>
+ inline void
+ evaluate_gradients_collocation(
+ const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
+ const Number *values,
+ Number *gradients)
{
- using Number2 =
- typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
- using Eval = EvaluatorTensorProduct<evaluate_evenodd,
- dim,
- fe_degree + 1,
- fe_degree + 1,
- Number,
- Number2>;
-
- static void
- evaluate(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval);
-
- static void
- do_evaluate(const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- Number *gradients_quad,
- Number *hessians_quad);
-
- static void
- integrate(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array);
+ AssertDimension(shape.shape_gradients_collocation_eo.size(),
+ (n_points_1d + 1) / 2 * n_points_1d);
- static void
- do_integrate(const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
- Number *gradients_quad,
- const Number *hessians_quad,
- const bool add_into_values_array);
- };
+ EvaluatorTensorProduct<evaluate_evenodd,
+ dim,
+ n_points_1d,
+ n_points_1d,
+ Number,
+ Number2>
+ eval({}, shape.shape_gradients_collocation_eo, {});
+ EvaluatorTensorProduct<evaluate_evenodd,
+ 2,
+ n_points_1d,
+ n_points_1d,
+ Number,
+ Number2>
+ eval_2d({}, shape.shape_gradients_collocation_eo, {});
+
+ if (dim == 1)
+ eval.template gradients<0, true, false>(values, gradients);
+ else
+ {
+ if (dim > 2)
+ eval.template gradients<2, true, false, dim>(values, gradients + 2);
+ constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
+ constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
+ const Number *in = values + (loop_bound - 1) * n_points_2d;
+ Number *out = gradients + (loop_bound - 1) * dim * n_points_2d;
+ for (unsigned int l = 0; l < loop_bound; ++l)
+ {
+ eval_2d.template gradients<0, true, false, dim>(in, out);
+ eval_2d.template gradients<1, true, false, dim>(in, out + 1);
+ in -= n_points_2d;
+ out -= dim * n_points_2d;
+ }
+ }
+ }
- template <int dim, int fe_degree, typename Number>
+ /**
+ * Internal function that multiplies by the gradients of test functions and
+ * sums over quadrature points for function representations in the
+ * collocation space, used by FEEvaluationImplCollocation and
+ * FEEvaluationImplTransformToCollocation. The evaluation strategy uses sum
+ * factorization with the even-odd optimization and fixed loop bounds.
+ */
+ template <int n_points_1d, int dim, typename Number, typename Number2>
inline void
- FEEvaluationImplCollocation<dim, fe_degree, Number>::evaluate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval)
+ integrate_gradients_collocation(
+ const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
+ Number *values,
+ const Number *gradients,
+ const bool add_into_values_array)
{
- constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+ AssertDimension(shape.shape_gradients_collocation_eo.size(),
+ (n_points_1d + 1) / 2 * n_points_1d);
- for (unsigned int c = 0; c < n_components; ++c)
+ EvaluatorTensorProduct<evaluate_evenodd,
+ dim,
+ n_points_1d,
+ n_points_1d,
+ Number,
+ Number2>
+ eval({}, shape.shape_gradients_collocation_eo, {});
+ EvaluatorTensorProduct<evaluate_evenodd,
+ 2,
+ n_points_1d,
+ n_points_1d,
+ Number,
+ Number2>
+ eval_2d({}, shape.shape_gradients_collocation_eo, {});
+
+ if (dim == 1)
+ {
+ if (add_into_values_array)
+ eval.template gradients<0, false, true>(gradients, values);
+ else
+ eval.template gradients<0, false, false>(gradients, values);
+ }
+ else
{
- if ((evaluation_flag & EvaluationFlags::values) != 0u)
- for (unsigned int i = 0; i < n_points; ++i)
- fe_eval.begin_values()[n_points * c + i] =
- values_dofs[n_points * c + i];
-
- do_evaluate(fe_eval.get_shape_info().data.front(),
- evaluation_flag,
- values_dofs + c * n_points,
- fe_eval.begin_gradients() + c * dim * n_points,
- fe_eval.begin_hessians() +
- c * dim * (dim + 1) / 2 * n_points);
+ constexpr unsigned int loop_bound = (dim > 2 ? n_points_1d : 1);
+ constexpr unsigned int n_points_2d = n_points_1d * n_points_1d;
+
+ const Number *in = gradients + (loop_bound - 1) * dim * n_points_2d;
+ Number *out = values + (loop_bound - 1) * n_points_2d;
+ for (unsigned int l = 0; l < loop_bound; ++l)
+ {
+ if (add_into_values_array)
+ eval_2d.template gradients<0, false, true, dim>(in, out);
+ else
+ eval_2d.template gradients<0, false, false, dim>(in, out);
+ eval_2d.template gradients<1, false, true, dim>(in + 1, out);
+ in -= dim * n_points_2d;
+ out -= n_points_2d;
+ }
}
+ if (dim > 2)
+ eval.template gradients<2, false, true, dim>(gradients + 2, values);
}
- template <int dim, int fe_degree, typename Number>
+ /**
+ * Internal function that evaluates the Hessians of finite element functions
+ * represented by bases in the collocation space, used by
+ * FEEvaluationImplSelector. The evaluation strategy uses sum
+ * factorization with fixed loop bounds.
+ */
+ template <int n_points_1d, int dim, typename Number>
inline void
- FEEvaluationImplCollocation<dim, fe_degree, Number>::do_evaluate(
- const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- Number *gradients_quad,
- Number *hessians_quad)
+ evaluate_hessians_collocation(const unsigned int n_components,
+ FEEvaluationData<dim, Number, false> &fe_eval)
{
- AssertDimension(shape.shape_gradients_collocation_eo.size(),
- (fe_degree + 2) / 2 * (fe_degree + 1));
- constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
-
- Eval eval({},
- shape.shape_gradients_collocation_eo,
- shape.shape_hessians_collocation_eo);
- if ((evaluation_flag &
- (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
+ using Number2 =
+ typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+
+ // might have non-symmetric quadrature formula, so use the more
+ // conservative 'evaluate_general' scheme rather than 'even_odd' as the
+ // Hessians are not used very often
+ const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+ fe_eval.get_shape_info().data[0];
+ AssertDimension(data.shape_gradients_collocation.size(),
+ data.n_q_points_1d * data.n_q_points_1d);
+ EvaluatorTensorProduct<evaluate_general,
+ dim,
+ n_points_1d,
+ n_points_1d,
+ Number,
+ Number2>
+ eval({},
+ data.shape_gradients_collocation.data(),
+ data.shape_hessians_collocation.data(),
+ data.n_q_points_1d,
+ data.n_q_points_1d);
+
+ const Number *values = fe_eval.begin_values();
+ Number *hessians = fe_eval.begin_hessians();
+ Number *scratch = fe_eval.get_scratch_data().begin();
+ const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
+ for (unsigned int comp = 0; comp < n_components; ++comp)
{
- eval.template gradients<0, true, false>(values_dofs, gradients_quad);
+ // xx derivative
+ eval.template hessians<0, true, false>(values, hessians);
if (dim > 1)
- eval.template gradients<1, true, false>(values_dofs,
- gradients_quad + n_points);
+ {
+ // xy derivative: we might or might not have the gradients already
+ // computed elsewhere, but we recompute them here since it adds
+ // only moderate extra work (at most 25%)
+ eval.template gradients<0, true, false>(values, scratch);
+ eval.template gradients<1, true, false>(scratch,
+ hessians + dim * n_points);
+ // yy derivative
+ eval.template hessians<1, true, false>(values, hessians + n_points);
+ }
if (dim > 2)
- eval.template gradients<2, true, false>(values_dofs,
- gradients_quad +
- 2 * n_points);
+ {
+ // xz derivative
+ eval.template gradients<2, true, false>(scratch,
+ hessians + 4 * n_points);
+ // yz derivative
+ eval.template gradients<1, true, false>(values, scratch);
+ eval.template gradients<2, true, false>(scratch,
+ hessians + 5 * n_points);
+ // zz derivative
+ eval.template hessians<2, true, false>(values,
+ hessians + 2 * n_points);
+ }
+
+ values += n_points;
+ hessians += (dim * (dim + 1)) / 2 * n_points;
}
- if (evaluation_flag & EvaluationFlags::hessians)
+ }
+
+
+
+ /**
+ * Internal function that multiplies by the Hessians of test functions and
+ * sums over quadrature points for function representations in the
+ * collocation space, used by FEEvaluationImplSelector. The evaluation
+ * strategy uses sum factorization with fixed loop bounds.
+ */
+ template <int n_q_points_1d, int dim, typename Number>
+ inline void
+ integrate_hessians_collocation(const unsigned int n_components,
+ FEEvaluationData<dim, Number, false> &fe_eval,
+ const bool add_into_values_array)
+ {
+ using Number2 =
+ typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+
+ const MatrixFreeFunctions::UnivariateShapeData<Number> &data =
+ fe_eval.get_shape_info().data[0];
+ AssertDimension(data.shape_gradients_collocation.size(),
+ data.n_q_points_1d * data.n_q_points_1d);
+ EvaluatorTensorProduct<evaluate_general,
+ dim,
+ n_q_points_1d,
+ n_q_points_1d,
+ Number,
+ Number2>
+ eval({},
+ data.shape_gradients_collocation.data(),
+ data.shape_hessians_collocation.data(),
+ data.n_q_points_1d,
+ data.n_q_points_1d);
+ Number *values = fe_eval.begin_values();
+ const Number *hessians = fe_eval.begin_hessians();
+ Number *scratch = fe_eval.get_scratch_data().begin();
+ const std::size_t n_points = fe_eval.get_shape_info().n_q_points;
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
{
- eval.template hessians<0, true, false>(values_dofs, hessians_quad);
+ // xx derivative
+ if (add_into_values_array == true)
+ eval.template hessians<0, false, true>(hessians, values);
+ else
+ eval.template hessians<0, false, false>(hessians, values);
+
+ // yy derivative
if (dim > 1)
+ eval.template hessians<1, false, true>(hessians + n_points, values);
+ if (dim > 2)
{
- eval.template gradients<1, true, false>(gradients_quad,
- hessians_quad +
- dim * n_points);
- eval.template hessians<1, true, false>(values_dofs,
- hessians_quad + n_points);
+ // zz derivative
+ eval.template hessians<2, false, true>(hessians + 2 * n_points,
+ values);
+ // yz derivative
+ eval.template gradients<2, false, false>(hessians + 5 * n_points,
+ scratch);
+ eval.template gradients<1, false, true>(scratch, values);
+
+ // xz derivative
+ eval.template gradients<2, false, false>(hessians + 4 * n_points,
+ scratch);
}
- if (dim > 2)
+
+ if (dim > 1)
{
- eval.template gradients<2, true, false>(gradients_quad,
- hessians_quad +
- 4 * n_points);
- eval.template gradients<2, true, false>(gradients_quad + n_points,
- hessians_quad +
- 5 * n_points);
- eval.template hessians<2, true, false>(values_dofs,
- hessians_quad +
- 2 * n_points);
+ // xy derivative, combined with xz in 3d
+ eval.template gradients<1, false, (dim > 2)>(hessians +
+ dim * n_points,
+ scratch);
+ eval.template gradients<0, false, true>(scratch, values);
}
+
+ values += n_points;
+ hessians += (dim * (dim + 1)) / 2 * n_points;
}
}
- template <int dim, int fe_degree, typename Number>
- inline void
- FEEvaluationImplCollocation<dim, fe_degree, Number>::integrate(
- const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array)
+ /**
+ * Internal function to evaluate the Hessians of finite element functions in
+ * the non-collocation setting as a fall-back. The evaluation strategy uses
+ * sum factorization with run-time loop bounds and is thus slower than the
+ * collocation case, but it is not as widely used and thus uncritical.
+ */
+ template <int dim, typename Number>
+ void
+ evaluate_hessians_slow(const unsigned int n_components,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, false> &fe_eval)
+ {
+ const auto &univariate_shape_data = fe_eval.get_shape_info().data;
+ using Impl =
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
+ using Eval = typename Impl::Eval;
+ Eval eval0 =
+ Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
+ Eval eval1 = Impl::create_evaluator_tensor_product(
+ &univariate_shape_data[std::min<int>(1,
+ univariate_shape_data.size() - 1)]);
+ Eval eval2 = Impl::create_evaluator_tensor_product(
+ &univariate_shape_data[std::min<int>(2,
+ univariate_shape_data.size() - 1)]);
+
+ const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
+ Number *tmp1 = fe_eval.get_scratch_data().begin();
+ Number *tmp2 =
+ tmp1 + std::max(Utilities::fixed_power<dim>(
+ univariate_shape_data.front().fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ univariate_shape_data.front().n_q_points_1d));
+ Number *hessians = fe_eval.begin_hessians();
+
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp,
+ hessians += n_points * dim * (dim + 1) / 2,
+ values_dofs +=
+ fe_eval.get_shape_info().dofs_per_component_on_cell)
+ switch (dim)
+ {
+ case 1:
+ eval0.template hessians<0, true, false>(values_dofs, hessians);
+ break;
+ case 2:
+ // xx derivative
+ eval0.template hessians<0, true, false>(values_dofs, tmp1);
+ eval1.template values<1, true, false>(tmp1, hessians);
+ // xy derivative
+ eval0.template gradients<0, true, false>(values_dofs, tmp1);
+ eval1.template gradients<1, true, false>(tmp1,
+ hessians + 2 * n_points);
+ // yy derivative
+ eval0.template values<0, true, false>(values_dofs, tmp1);
+ eval1.template hessians<1, true, false>(tmp1, hessians + n_points);
+ break;
+ case 3:
+ // xx derivative
+ eval0.template hessians<0, true, false>(values_dofs, tmp1);
+ eval1.template values<1, true, false>(tmp1, tmp2);
+ eval2.template values<2, true, false>(tmp2, hessians);
+ // xy derivative
+ eval0.template gradients<0, true, false>(values_dofs, tmp1);
+ eval1.template gradients<1, true, false>(tmp1, tmp2);
+ eval2.template values<2, true, false>(tmp2,
+ hessians + 3 * n_points);
+ // xz derivative
+ eval1.template values<1, true, false>(tmp1, tmp2);
+ eval2.template gradients<2, true, false>(tmp2,
+ hessians + 4 * n_points);
+ // yy derivative
+ eval0.template values<0, true, false>(values_dofs, tmp1);
+ eval1.template hessians<1, true, false>(tmp1, tmp2);
+ eval2.template values<2, true, false>(tmp2, hessians + n_points);
+ // yz derivative
+ eval1.template gradients<1, true, false>(tmp1, tmp2);
+ eval2.template gradients<2, true, false>(tmp2,
+ hessians + 5 * n_points);
+ // zz derivative
+ eval1.template values<1, true, false>(tmp1, tmp2);
+ eval2.template hessians<2, true, false>(tmp2,
+ hessians + 2 * n_points);
+ break;
+
+ default:
+ Assert(false,
+ ExcNotImplemented(
+ "Only 1d, 2d and 3d implemented for Hessian"));
+ }
+ }
+
+
+
+ /**
+ * Internal function to multiply by the Hessians of the test functions and
+ * integrate in the non-collocation setting as a fall-back. The evaluation
+ * strategy uses sum factorization with run-time loop bounds and is thus
+ * slower than the collocation case, but it is not as widely used and thus
+ * uncritical.
+ */
+ template <int dim, typename Number>
+ void
+ integrate_hessians_slow(const unsigned int n_components,
+ const FEEvaluationData<dim, Number, false> &fe_eval,
+ Number *values_dofs,
+ const bool add_into_values_array)
{
- constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+ const auto &univariate_shape_data = fe_eval.get_shape_info().data;
+ using Impl =
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_general, dim, -1, 0, Number>;
+ using Eval = typename Impl::Eval;
+ Eval eval0 =
+ Impl::create_evaluator_tensor_product(&univariate_shape_data[0]);
+ Eval eval1 = Impl::create_evaluator_tensor_product(
+ &univariate_shape_data[std::min<int>(1,
+ univariate_shape_data.size() - 1)]);
+ Eval eval2 = Impl::create_evaluator_tensor_product(
+ &univariate_shape_data[std::min<int>(2,
+ univariate_shape_data.size() - 1)]);
+
+ const unsigned int n_points = fe_eval.get_shape_info().n_q_points;
+ Number *tmp1 = fe_eval.get_scratch_data().begin();
+ Number *tmp2 =
+ tmp1 + std::max(Utilities::fixed_power<dim>(
+ univariate_shape_data.front().fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ univariate_shape_data.front().n_q_points_1d));
+ const Number *hessians = fe_eval.begin_hessians();
+
+ for (unsigned int comp = 0; comp < n_components;
+ ++comp,
+ hessians += n_points * dim * (dim + 1) / 2,
+ values_dofs +=
+ fe_eval.get_shape_info().dofs_per_component_on_cell)
+ switch (dim)
+ {
+ case 1:
+ if (add_into_values_array)
+ eval0.template hessians<0, false, true>(hessians, values_dofs);
+ else
+ eval0.template hessians<0, false, false>(hessians, values_dofs);
+ break;
+ case 2:
+ // xx derivative
+ eval1.template values<1, false, false>(hessians, tmp1);
+ if (add_into_values_array)
+ eval0.template hessians<0, false, true>(tmp1, values_dofs);
+ else
+ eval0.template hessians<0, false, false>(tmp1, values_dofs);
+
+ // xy derivative
+ eval1.template gradients<1, false, false>(hessians + 2 * n_points,
+ tmp1);
+ eval0.template gradients<0, false, true>(tmp1, values_dofs);
+ // yy derivative
+ eval1.template hessians<1, false, false>(hessians + n_points, tmp1);
+ eval0.template values<0, false, true>(tmp1, values_dofs);
+ break;
+ case 3:
+ // xx derivative
+ eval2.template values<2, false, false>(hessians, tmp1);
+ eval1.template values<1, false, false>(tmp1, tmp2);
- for (unsigned int c = 0; c < n_components; ++c)
- {
- if ((integration_flag & EvaluationFlags::values) != 0u)
- {
if (add_into_values_array)
- for (unsigned int i = 0; i < n_points; ++i)
- values_dofs[n_points * c + i] +=
- fe_eval.begin_values()[n_points * c + i];
+ eval0.template hessians<0, false, true>(tmp2, values_dofs);
else
- for (unsigned int i = 0; i < n_points; ++i)
- values_dofs[n_points * c + i] =
- fe_eval.begin_values()[n_points * c + i];
- }
+ eval0.template hessians<0, false, false>(tmp2, values_dofs);
+
+ // xy derivative
+ eval2.template values<2, false, false>(hessians + 3 * n_points,
+ tmp1);
+ eval1.template gradients<1, false, false>(tmp1, tmp2);
+ // xz derivative
+ eval2.template gradients<2, false, false>(hessians + 4 * n_points,
+ tmp1);
+ eval1.template values<1, false, true>(tmp1, tmp2);
+ eval1.template values<0, false, true>(tmp2, values_dofs);
+
+ // yy derivative
+ eval2.template values<2, false, false>(hessians + n_points, tmp1);
+ eval1.template hessians<1, false, false>(tmp1, tmp2);
+
+ // yz derivative
+ eval2.template gradients<2, false, false>(hessians + 5 * n_points,
+ tmp1);
+ eval1.template gradients<1, false, true>(tmp1, tmp2);
+
+ // zz derivative
+ eval2.template hessians<2, false, false>(hessians + 2 * n_points,
+ tmp1);
+ eval1.template values<1, false, true>(tmp1, tmp2);
+ eval0.template values<0, false, true>(tmp2, values_dofs);
+ break;
- do_integrate(fe_eval.get_shape_info().data.front(),
- integration_flag,
- values_dofs + c * n_points,
- fe_eval.begin_gradients() + c * dim * n_points,
- fe_eval.begin_hessians() +
- c * dim * (dim + 1) / 2 * n_points,
- add_into_values_array ||
- ((integration_flag & EvaluationFlags::values) != 0u));
- }
+ default:
+ Assert(false,
+ ExcNotImplemented(
+ "Only 1d, 2d and 3d implemented for Hessian"));
+ }
}
+ /**
+ * This struct performs the evaluation of function values, gradients and
+ * Hessians for tensor-product finite elements. This a specialization for
+ * elements where the nodal points coincide with the quadrature points like
+ * FE_Q shape functions on Gauss-Lobatto elements integrated with
+ * Gauss-Lobatto quadrature. The assumption of this class is that the shape
+ * 'values' operation is identity, which allows us to write shorter code.
+ *
+ * In literature, this form of evaluation is often called spectral
+ * evaluation, spectral collocation or simply collocation, meaning the same
+ * location for shape functions and evaluation space (quadrature points).
+ */
template <int dim, int fe_degree, typename Number>
- inline void
- FEEvaluationImplCollocation<dim, fe_degree, Number>::do_integrate(
- const MatrixFreeFunctions::UnivariateShapeData<Number2> &shape,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
- Number *gradients_quad,
- const Number *hessians_quad,
- const bool add_into_values_array)
+ struct FEEvaluationImplCollocation
{
- AssertDimension(shape.shape_gradients_collocation_eo.size(),
- (fe_degree + 2) / 2 * (fe_degree + 1));
+ using Number2 =
+ typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
+ using Eval = EvaluatorTensorProduct<evaluate_evenodd,
+ dim,
+ fe_degree + 1,
+ fe_degree + 1,
+ Number,
+ Number2>;
- Eval eval({},
- shape.shape_gradients_collocation_eo,
- shape.shape_hessians_collocation_eo);
- constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
+ static void
+ evaluate(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, false> &fe_eval)
+ {
+ constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
- if ((integration_flag & EvaluationFlags::hessians) != 0u)
- {
- // diagonal
- // grad xx
- if (add_into_values_array == true)
- eval.template hessians<0, false, true>(hessians_quad, values_dofs);
- else
- eval.template hessians<0, false, false>(hessians_quad, values_dofs);
- // grad yy
- if (dim > 1)
- eval.template hessians<1, false, true>(hessians_quad + n_points,
- values_dofs);
- // grad zz
- if (dim > 2)
- eval.template hessians<2, false, true>(hessians_quad + 2 * n_points,
- values_dofs);
- // off-diagonal
- if (dim == 2)
- {
- // grad xy, queue into gradient
- if (integration_flag & EvaluationFlags::gradients)
- eval.template gradients<1, false, true>(hessians_quad +
- 2 * n_points,
- gradients_quad);
- else
- eval.template gradients<1, false, false>(hessians_quad +
- 2 * n_points,
- gradients_quad);
- }
- if (dim == 3)
- {
- // grad xy, queue into gradient
- if (integration_flag & EvaluationFlags::gradients)
- eval.template gradients<1, false, true>(hessians_quad +
- 3 * n_points,
- gradients_quad);
- else
- eval.template gradients<1, false, false>(hessians_quad +
- 3 * n_points,
- gradients_quad);
-
- // grad xz
- eval.template gradients<2, false, true>(hessians_quad +
- 4 * n_points,
- gradients_quad);
-
- // grad yz
- if (integration_flag & EvaluationFlags::gradients)
- eval.template gradients<2, false, true>(
- hessians_quad + 5 * n_points, gradients_quad + n_points);
- else
- eval.template gradients<2, false, false>(
- hessians_quad + 5 * n_points, gradients_quad + n_points);
- }
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ if ((evaluation_flag & EvaluationFlags::values) != 0u)
+ for (unsigned int i = 0; i < n_points; ++i)
+ fe_eval.begin_values()[n_points * c + i] =
+ values_dofs[n_points * c + i];
+
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ evaluate_gradients_collocation<fe_degree + 1, dim>(
+ fe_eval.get_shape_info().data.front(),
+ values_dofs + c * n_points,
+ fe_eval.begin_gradients() + c * dim * n_points);
+ }
+ }
- // if we did not integrate gradients, set the last slot to zero
- // which was not touched before, in order to avoid the if
- // statement in the gradients loop below
- if ((integration_flag & EvaluationFlags::gradients) == 0u)
- for (unsigned int q = 0; q < n_points; ++q)
- gradients_quad[(dim - 1) * n_points + q] = Number();
- }
+ static void
+ integrate(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, false> &fe_eval,
+ const bool add_into_values_array)
+ {
+ constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim);
- if ((integration_flag &
- (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u)
- {
- if (add_into_values_array ||
- (integration_flag & EvaluationFlags::hessians) != 0u)
- eval.template gradients<0, false, true>(gradients_quad, values_dofs);
- else
- eval.template gradients<0, false, false>(gradients_quad, values_dofs);
- if (dim > 1)
- eval.template gradients<1, false, true>(gradients_quad + n_points,
- values_dofs);
- if (dim > 2)
- eval.template gradients<2, false, true>(gradients_quad + 2 * n_points,
- values_dofs);
- }
- }
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ if ((integration_flag & EvaluationFlags::values) != 0u)
+ {
+ if (add_into_values_array)
+ for (unsigned int i = 0; i < n_points; ++i)
+ values_dofs[n_points * c + i] +=
+ fe_eval.begin_values()[n_points * c + i];
+ else
+ for (unsigned int i = 0; i < n_points; ++i)
+ values_dofs[n_points * c + i] =
+ fe_eval.begin_values()[n_points * c + i];
+ }
+
+ if ((integration_flag & EvaluationFlags::gradients) != 0u)
+ integrate_gradients_collocation<fe_degree + 1, dim>(
+ fe_eval.get_shape_info().data.front(),
+ values_dofs + c * n_points,
+ fe_eval.begin_gradients() + c * dim * n_points,
+ add_into_values_array ||
+ ((integration_flag & EvaluationFlags::values) != 0u));
+ }
+ }
+ };
evaluate(const unsigned int n_components,
const EvaluationFlags::EvaluationFlags evaluation_flag,
const Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval);
+ FEEvaluationData<dim, Number, false> &fe_eval)
+ {
+ const auto &shape_data = fe_eval.get_shape_info().data.front();
+
+ Assert(n_q_points_1d > fe_degree,
+ ExcMessage("You lose information when going to a collocation "
+ "space of lower degree, so the evaluation results "
+ "would be wrong. Thus, this class does not permit "
+ "the chosen operation."));
+ constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ FEEvaluationImplBasisChange<
+ evaluate_evenodd,
+ EvaluatorQuantity::value,
+ dim,
+ (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
+ n_q_points_1d>::do_forward(1,
+ shape_data.shape_values_eo,
+ values_dofs + c * n_dofs,
+ fe_eval.begin_values() + c * n_q_points);
+
+ // apply derivatives in the collocation space
+ if (evaluation_flag & EvaluationFlags::gradients)
+ evaluate_gradients_collocation<n_q_points_1d, dim>(
+ shape_data,
+ fe_eval.begin_values() + c * n_q_points,
+ fe_eval.begin_gradients() + c * dim * n_q_points);
+ }
+ }
static void
integrate(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const EvaluationFlags::EvaluationFlags integration_flag,
Number *values_dofs,
FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array);
+ const bool add_into_values_array)
+ {
+ const auto &shape_data = fe_eval.get_shape_info().data.front();
+
+ Assert(n_q_points_1d > fe_degree,
+ ExcMessage("You lose information when going to a collocation "
+ "space of lower degree, so the evaluation results "
+ "would be wrong. Thus, this class does not permit "
+ "the chosen operation."));
+ constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ // apply derivatives in collocation space
+ if (integration_flag & EvaluationFlags::gradients)
+ integrate_gradients_collocation<n_q_points_1d, dim>(
+ shape_data,
+ fe_eval.begin_values() + c * n_q_points,
+ fe_eval.begin_gradients() + c * dim * n_q_points,
+ /*add_into_values_array=*/
+ integration_flag & EvaluationFlags::values);
+
+ // transform back to the original space
+ FEEvaluationImplBasisChange<
+ evaluate_evenodd,
+ EvaluatorQuantity::value,
+ dim,
+ (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
+ n_q_points_1d>::do_backward(1,
+ shape_data.shape_values_eo,
+ add_into_values_array,
+ fe_eval.begin_values() + c * n_q_points,
+ values_dofs +
+ c *
+ Utilities::pow(fe_degree + 1, dim));
+ }
+ }
};
+ /**
+ * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use
+ * specific sum-factorization kernels and with normal/tangential shape_data
+ */
template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- inline void
- FEEvaluationImplTransformToCollocation<
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::evaluate(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval)
+ struct FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>
{
- const auto &shape_data = fe_eval.get_shape_info().data.front();
+ using Number2 =
+ typename FEEvaluationData<dim, Number, false>::shape_info_number_type;
- Assert(n_q_points_1d > fe_degree,
- ExcMessage("You lose information when going to a collocation space "
- "of lower degree, so the evaluation results would be "
- "wrong. Thus, this class does not permit the desired "
- "operation."));
- constexpr std::size_t n_dofs = Utilities::pow(fe_degree + 1, dim);
- constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+ template <bool integrate>
+ static void
+ evaluate_or_integrate(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number *values_dofs_actual,
+ FEEvaluationData<dim, Number, false> &fe_eval,
+ const bool add_into_values_array = false);
- for (unsigned int c = 0; c < n_components; ++c)
- {
- FEEvaluationImplBasisChange<
- evaluate_evenodd,
- EvaluatorQuantity::value,
- dim,
- (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
- n_q_points_1d>::do_forward(1,
- shape_data.shape_values_eo,
- values_dofs + c * n_dofs,
- fe_eval.begin_values() + c * n_q_points);
-
- // apply derivatives in the collocation space
- if (evaluation_flag &
- (EvaluationFlags::gradients | EvaluationFlags::hessians))
- FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
- do_evaluate(shape_data,
- evaluation_flag & (EvaluationFlags::gradients |
- EvaluationFlags::hessians),
- fe_eval.begin_values() + c * n_q_points,
- fe_eval.begin_gradients() + c * dim * n_q_points,
- fe_eval.begin_hessians() +
- c * dim * (dim + 1) / 2 * n_q_points);
- }
- }
+ private:
+ template <int direction, bool contract_over_rows>
+ static void
+ work_normal(const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+ const Number *in,
+ Number *out,
+ const bool add_into_result = false)
+ {
+ AssertIndexRange(direction, dim);
+ constexpr int n_rows = fe_degree + 1;
+ constexpr int n_columns = n_q_points_1d;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns;
+ constexpr int nn = contract_over_rows ? n_columns : n_rows;
+ const Number2 *shape_data = data.shape_values_eo.data();
+ Assert(shape_data != nullptr, ExcNotInitialized());
+ Assert(contract_over_rows == false || !add_into_result,
+ ExcMessage("Cannot add into result if contract_over_rows = true"));
+
+ constexpr int n_blocks1 = Utilities::pow(fe_degree, direction);
+ constexpr int n_blocks2 = Utilities::pow(fe_degree, dim - direction - 1);
+
+ for (int i2 = 0; i2 < n_blocks2; ++i2)
+ {
+ for (int i1 = 0; i1 < n_blocks1; ++i1)
+ {
+ if (contract_over_rows == false && add_into_result)
+ apply_matrix_vector_product<evaluate_evenodd,
+ EvaluatorQuantity::value,
+ n_rows,
+ n_columns,
+ n_blocks1,
+ n_blocks1,
+ contract_over_rows,
+ true>(shape_data, in, out);
+ else
+ apply_matrix_vector_product<evaluate_evenodd,
+ EvaluatorQuantity::value,
+ n_rows,
+ n_columns,
+ n_blocks1,
+ n_blocks1,
+ contract_over_rows,
+ false>(shape_data, in, out);
+
+ ++in;
+ ++out;
+ }
+ in += n_blocks1 * (mm - 1);
+ out += n_blocks1 * (nn - 1);
+ }
+ }
+
+ template <int direction, int normal_direction, bool contract_over_rows>
+ static void
+ work_tangential(
+ const MatrixFreeFunctions::UnivariateShapeData<Number2> &data,
+ Number *ptr)
+ {
+ AssertIndexRange(direction, dim);
+ static_assert(direction != normal_direction,
+ "Cannot interpolate tangentially in normal direction");
+
+ constexpr int n_rows = fe_degree;
+ constexpr int n_columns = n_q_points_1d;
+ const Number2 *shape_data = data.shape_values_eo.data();
+ Assert(shape_data != nullptr, ExcNotInitialized());
+
+ constexpr int n_blocks1 =
+ (direction > normal_direction) ?
+ Utilities::pow(n_q_points_1d, direction) :
+ (direction > 0 ?
+ (Utilities::pow(fe_degree, direction - 1) * n_q_points_1d) :
+ 1);
+ constexpr int n_blocks2 =
+ (direction > normal_direction) ?
+ Utilities::pow(fe_degree, dim - 1 - direction) :
+ ((direction + 1 < dim) ?
+ (Utilities::pow(fe_degree, dim - 2 - direction) * n_q_points_1d) :
+ 1);
+
+ // Since we perform an in-place interpolation, we must run the step
+ // expanding the size of the basis backward ('contract_over_rows' aka
+ // 'evaluate' case).
+ if (contract_over_rows)
+ {
+ const Number *in =
+ ptr + (n_blocks2 - 1) * n_blocks1 * n_rows + n_blocks1 - 1;
+ Number *out =
+ ptr + (n_blocks2 - 1) * n_blocks1 * n_columns + n_blocks1 - 1;
+ for (int i2 = 0; i2 < n_blocks2; ++i2)
+ {
+ for (int i1 = 0; i1 < n_blocks1; ++i1)
+ {
+ apply_matrix_vector_product<evaluate_evenodd,
+ EvaluatorQuantity::value,
+ n_rows,
+ n_columns,
+ n_blocks1,
+ n_blocks1,
+ true,
+ false>(shape_data, in, out);
+
+ --in;
+ --out;
+ }
+ in -= n_blocks1 * (n_rows - 1);
+ out -= n_blocks1 * (n_columns - 1);
+ }
+ }
+ else
+ {
+ const Number *in = ptr;
+ Number *out = ptr;
+ for (int i2 = 0; i2 < n_blocks2; ++i2)
+ {
+ for (int i1 = 0; i1 < n_blocks1; ++i1)
+ {
+ apply_matrix_vector_product<evaluate_evenodd,
+ EvaluatorQuantity::value,
+ n_rows,
+ n_columns,
+ n_blocks1,
+ n_blocks1,
+ false,
+ false>(shape_data, in, out);
+
+ ++in;
+ ++out;
+ }
+ in += n_blocks1 * (n_columns - 1);
+ out += n_blocks1 * (n_rows - 1);
+ }
+ }
+ }
+ };
template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <bool integrate>
inline void
- FEEvaluationImplTransformToCollocation<
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::integrate(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
- FEEvaluationData<dim, Number, false> &fe_eval,
- const bool add_into_values_array)
+ FEEvaluationImpl<MatrixFreeFunctions::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>::
+ evaluate_or_integrate(
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, false> &fe_eval,
+ const bool add_into_values_array)
{
- const auto &shape_data = fe_eval.get_shape_info().data.front();
+ Assert(dim == 2 || dim == 3,
+ ExcMessage("Only dim = 2,3 implemented for Raviart-Thomas "
+ "evaluation/integration"));
- Assert(n_q_points_1d > fe_degree,
- ExcMessage("You lose information when going to a collocation space "
- "of lower degree, so the evaluation results would be "
- "wrong. Thus, this class does not permit the desired "
- "operation."));
- constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim);
+ if (evaluation_flag == EvaluationFlags::nothing)
+ return;
- for (unsigned int c = 0; c < n_components; ++c)
+ AssertDimension(fe_eval.get_shape_info().data.size(), 2);
+ AssertDimension(n_q_points_1d,
+ fe_eval.get_shape_info().data[0].n_q_points_1d);
+ AssertDimension(n_q_points_1d,
+ fe_eval.get_shape_info().data[1].n_q_points_1d);
+ AssertDimension(fe_degree, fe_eval.get_shape_info().data[0].fe_degree);
+ AssertDimension(fe_degree, fe_eval.get_shape_info().data[1].fe_degree + 1);
+
+ const auto &shape_data = fe_eval.get_shape_info().data;
+ const unsigned int dofs_per_component =
+ Utilities::pow(fe_degree, dim - 1) * (fe_degree + 1);
+ const unsigned int n_points = Utilities::pow(n_q_points_1d, dim);
+ Number *gradients = fe_eval.begin_gradients();
+ Number *values = fe_eval.begin_values();
+
+ if (integrate)
{
- // apply derivatives in collocation space
- if (integration_flag &
- (EvaluationFlags::gradients | EvaluationFlags::hessians))
- FEEvaluationImplCollocation<dim, n_q_points_1d - 1, Number>::
- do_integrate(shape_data,
- integration_flag & (EvaluationFlags::gradients |
- EvaluationFlags::hessians),
- fe_eval.begin_values() + c * n_q_points,
- fe_eval.begin_gradients() + c * dim * n_q_points,
- fe_eval.begin_hessians() +
- c * dim * (dim + 1) / 2 * n_q_points,
- /*add_into_values_array=*/
- integration_flag & EvaluationFlags::values);
-
- // transform back to the original space
- FEEvaluationImplBasisChange<
- evaluate_evenodd,
- EvaluatorQuantity::value,
- dim,
- (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1),
- n_q_points_1d>::do_backward(1,
- shape_data.shape_values_eo,
- add_into_values_array,
- fe_eval.begin_values() + c * n_q_points,
- values_dofs +
- c * Utilities::pow(fe_degree + 1, dim));
+ const bool do_values = evaluation_flag & EvaluationFlags::values;
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients,
+ do_values);
+ if constexpr (dim > 2)
+ work_tangential<2, 0, false>(shape_data[1], values);
+ work_tangential<1, 0, false>(shape_data[1], values);
+ work_normal<0, false>(shape_data[0],
+ values,
+ values_dofs,
+ add_into_values_array);
+
+ values += n_points;
+ gradients += n_points * dim;
+ values_dofs += dofs_per_component;
+
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients,
+ do_values);
+ if constexpr (dim > 2)
+ work_tangential<2, 1, false>(shape_data[1], values);
+ work_tangential<0, 1, false>(shape_data[1], values);
+ work_normal<1, false>(shape_data[0],
+ values,
+ values_dofs,
+ add_into_values_array);
+
+ if constexpr (dim > 2)
+ {
+ values += n_points;
+ gradients += n_points * dim;
+ values_dofs += dofs_per_component;
+
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ integrate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients,
+ do_values);
+ work_tangential<1, 2, false>(shape_data[1], values);
+ work_tangential<0, 2, false>(shape_data[1], values);
+ work_normal<2, false>(shape_data[0],
+ values,
+ values_dofs,
+ add_into_values_array);
+ }
+ }
+ else
+ {
+ work_normal<0, true>(shape_data[0], values_dofs, values);
+ work_tangential<1, 0, true>(shape_data[1], values);
+ if constexpr (dim > 2)
+ work_tangential<2, 0, true>(shape_data[1], values);
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients);
+
+ values += n_points;
+ gradients += n_points * dim;
+ values_dofs += dofs_per_component;
+
+ work_normal<1, true>(shape_data[0], values_dofs, values);
+ work_tangential<0, 1, true>(shape_data[1], values);
+ if constexpr (dim > 2)
+ work_tangential<2, 1, true>(shape_data[1], values);
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients);
+
+ if constexpr (dim > 2)
+ {
+ values += n_points;
+ gradients += n_points * dim;
+ values_dofs += dofs_per_component;
+
+ work_normal<2, true>(shape_data[0], values_dofs, values);
+ work_tangential<0, 2, true>(shape_data[1], values);
+ work_tangential<1, 2, true>(shape_data[1], values);
+ if ((evaluation_flag & EvaluationFlags::gradients) != 0u)
+ evaluate_gradients_collocation<n_q_points_1d, dim>(shape_data[0],
+ values,
+ gradients);
+ }
}
}
const EvaluationFlags::EvaluationFlags evaluation_flag,
OtherNumber *values_dofs,
FEEvaluationData<dim, Number, false> &fe_eval,
- const bool sum_into_values_array = false)
+ const bool sum_into_values_array_in = false)
{
// `OtherNumber` is either `const Number` (evaluate()) or `Number`
// (integrate())
element_type == ElementType::tensor_raviart_thomas,
ExcNotImplemented());
+ EvaluationFlags::EvaluationFlags actual_flag = evaluation_flag;
+ bool sum_into_values_array = sum_into_values_array_in;
+ if (evaluation_flag & EvaluationFlags::hessians)
+ {
+ actual_flag |= EvaluationFlags::values;
+ Assert(element_type != MatrixFreeFunctions::tensor_none,
+ ExcNotImplemented());
+ if constexpr (do_integrate)
+ {
+ if (fe_eval.get_shape_info().data[0].fe_degree <
+ fe_eval.get_shape_info().data[0].n_q_points_1d)
+ integrate_hessians_collocation<n_q_points_1d>(
+ n_components,
+ fe_eval,
+ evaluation_flag & EvaluationFlags::values);
+ else
+ {
+ integrate_hessians_slow(n_components,
+ fe_eval,
+ values_dofs,
+ sum_into_values_array);
+ sum_into_values_array = true;
+ }
+ }
+ }
+
if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
element_type == ElementType::tensor_symmetric_collocation)
{
evaluate_or_integrate<
FEEvaluationImplCollocation<dim, fe_degree, Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
n_q_points_1d,
Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
n_q_points_1d,
Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
fe_degree,
n_q_points_1d,
Number>>(n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
n_q_points_1d,
Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
n_q_points_1d,
Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
}
else if (element_type == ElementType::tensor_raviart_thomas)
{
- FEEvaluationImpl<ElementType::tensor_raviart_thomas,
- dim,
- (fe_degree == -1) ? 1 : fe_degree,
- (n_q_points_1d < 1) ? 1 : n_q_points_1d,
- Number>::
- template evaluate_or_integrate<do_integrate>(evaluation_flag,
- const_cast<Number *>(
- values_dofs),
- fe_eval,
- sum_into_values_array);
+ if constexpr (fe_degree > 0 && n_q_points_1d > 0 && dim > 1)
+ {
+ FEEvaluationImpl<ElementType::tensor_raviart_thomas,
+ dim,
+ fe_degree,
+ n_q_points_1d,
+ Number>::
+ template evaluate_or_integrate<do_integrate>(
+ actual_flag,
+ const_cast<Number *>(values_dofs),
+ fe_eval,
+ sum_into_values_array);
+ }
+ else
+ {
+ Assert(false,
+ ExcNotImplemented("Raviart-Thomas currently only possible "
+ "in 2d/3d and with templated degree"));
+ }
}
else
{
n_q_points_1d,
Number>>(
n_components,
- evaluation_flag,
+ actual_flag,
values_dofs,
fe_eval,
sum_into_values_array);
}
+ if ((evaluation_flag & EvaluationFlags::hessians) && !do_integrate)
+ {
+ Assert(element_type != MatrixFreeFunctions::tensor_none,
+ ExcNotImplemented());
+ if (fe_eval.get_shape_info().data[0].fe_degree <
+ fe_eval.get_shape_info().data[0].n_q_points_1d)
+ evaluate_hessians_collocation<n_q_points_1d>(n_components, fe_eval);
+ else
+ evaluate_hessians_slow(n_components, values_dofs, fe_eval);
+ }
+
return false;
}
Number,
Number2>
eval_grad({}, data.shape_gradients_collocation_eo, {});
- eval_grad.template gradients<0, true, false>(
+ eval_grad.template gradients<0, true, false, 3>(
values_quad, gradients_quad);
- eval_grad.template gradients<1, true, false>(
- values_quad, gradients_quad + n_q_points);
+ eval_grad.template gradients<1, true, false, 3>(
+ values_quad, gradients_quad + 1);
}
else
{
// grad x
eval0.template gradients<0, true, false>(values_dofs,
scratch_data);
- eval1.template values<1, true, false>(scratch_data,
- gradients_quad);
+ eval1.template values<1, true, false, 3>(scratch_data,
+ gradients_quad);
// grad y
eval0.template values<0, true, false>(values_dofs,
scratch_data);
- eval1.template gradients<1, true, false>(scratch_data,
- gradients_quad +
- n_q_points);
+ eval1.template gradients<1, true, false, 3>(
+ scratch_data, gradients_quad + 1);
if ((evaluation_flag & EvaluationFlags::values) != 0u)
eval1.template values<1, true, false>(scratch_data,
// grad z
eval0.template values<0, true, false>(values_dofs + n_dofs,
scratch_data);
- eval1.template values<1, true, false>(
- scratch_data, gradients_quad + (dim - 1) * n_q_points);
+ eval1.template values<1, true, false, 3>(scratch_data,
+ gradients_quad + 2);
break;
case 2:
- eval0.template values<0, true, false>(values_dofs + n_dofs,
- gradients_quad +
- n_q_points);
- eval0.template gradients<0, true, false>(values_dofs,
- gradients_quad);
+ eval0.template values<0, true, false, 2>(values_dofs + n_dofs,
+ gradients_quad + 1);
+ eval0.template gradients<0, true, false, 2>(values_dofs,
+ gradients_quad);
if ((evaluation_flag & EvaluationFlags::values) != 0u)
eval0.template values<0, true, false>(values_dofs,
values_quad);
{
case 3:
// grad z
- eval1.template values<1, false, false>(gradients_quad +
- 2 * n_q_points,
- gradients_quad +
- 2 * n_q_points);
- eval0.template values<0, false, false>(gradients_quad +
- 2 * n_q_points,
+ eval1.template values<1, false, false, 3>(gradients_quad + 2,
+ scratch_data);
+ eval0.template values<0, false, false>(scratch_data,
values_dofs + n_dofs);
if (symmetric_evaluate &&
use_collocation_evaluation(fe_degree, n_q_points_1d))
Number2>
eval_grad({}, data.shape_gradients_collocation_eo, {});
if ((integration_flag & EvaluationFlags::values) != 0u)
- eval_grad.template gradients<1, false, true>(
- gradients_quad + n_q_points, values_quad);
+ eval_grad.template gradients<1, false, true, 3>(
+ gradients_quad + 1, values_quad);
else
- eval_grad.template gradients<1, false, false>(
- gradients_quad + n_q_points, values_quad);
- eval_grad.template gradients<0, false, true>(
+ eval_grad.template gradients<1, false, false, 3>(
+ gradients_quad + 1, values_quad);
+ eval_grad.template gradients<0, false, true, 3>(
gradients_quad, values_quad);
eval0.template values<1, false, false>(values_quad,
values_quad);
{
eval1.template values<1, false, false>(values_quad,
scratch_data);
- eval1.template gradients<1, false, true>(
- gradients_quad + n_q_points, scratch_data);
+ eval1.template gradients<1, false, true, 3>(
+ gradients_quad + 1, scratch_data);
}
else
- eval1.template gradients<1, false, false>(
- gradients_quad + n_q_points, scratch_data);
+ eval1.template gradients<1, false, false, 3>(
+ gradients_quad + 1, scratch_data);
// grad y
eval0.template values<0, false, false>(scratch_data,
values_dofs);
// grad x
- eval1.template values<1, false, false>(gradients_quad,
- scratch_data);
+ eval1.template values<1, false, false, 3>(gradients_quad,
+ scratch_data);
eval0.template gradients<0, false, true>(scratch_data,
values_dofs);
}
break;
case 2:
- eval0.template values<0, false, false>(gradients_quad +
- n_q_points,
- values_dofs + n_dofs);
- eval0.template gradients<0, false, false>(gradients_quad,
- values_dofs);
+ eval0.template values<0, false, false, 2>(gradients_quad + 1,
+ values_dofs +
+ n_dofs);
+ eval0.template gradients<0, false, false, 2>(gradients_quad,
+ values_dofs);
if ((integration_flag & EvaluationFlags::values) != 0u)
eval0.template values<0, false, true>(values_quad,
values_dofs);
}
};
+
+
template <int dim, int fe_degree, int n_q_points_1d, typename Number>
struct FEFaceEvaluationImplRaviartThomas
{
{
if (face_direction == face_no / 2)
{
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree + 1,
- 0,
- Number,
- Number2>
- evalf(shape_data[face_no % 2].begin(),
- nullptr,
- nullptr,
- n_points_1d,
- 0);
-
- 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;
+ constexpr int stride_ = Utilities::pow(fe_degree + 1, face_direction);
+
+ const int n_rows = fe_degree != -1 ? fe_degree + 1 : n_points_1d;
+ const int stride = Utilities::pow(n_rows, face_direction);
+ const std::array<int, 2> n_blocks{
+ {(dim > 1 ? n_rows : 1), (dim > 2 ? n_rows : 1)}};
+ 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}};
for (unsigned int c = 0; c < n_components; ++c)
{
if (flag & EvaluationFlags::hessians)
- evalf.template apply_face<face_direction,
- do_evaluate,
- add_into_output,
- 2>(input, output);
+ interpolate_to_face<fe_degree + 1,
+ stride_,
+ do_evaluate,
+ add_into_output,
+ 2>(shape_data[face_no % 2].begin(),
+ n_blocks,
+ steps,
+ input,
+ output,
+ n_rows,
+ stride);
else if (flag & EvaluationFlags::gradients)
- evalf.template apply_face<face_direction,
- do_evaluate,
- add_into_output,
- 1>(input, output);
+ interpolate_to_face<fe_degree + 1,
+ stride_,
+ do_evaluate,
+ add_into_output,
+ 1>(shape_data[face_no % 2].begin(),
+ n_blocks,
+ steps,
+ input,
+ output,
+ n_rows,
+ stride);
else
- evalf.template apply_face<face_direction,
- do_evaluate,
- add_into_output,
- 0>(input, output);
- input += in_stride;
- output += out_stride;
+ interpolate_to_face<fe_degree + 1,
+ stride_,
+ do_evaluate,
+ add_into_output,
+ 0>(shape_data[face_no % 2].begin(),
+ n_blocks,
+ steps,
+ input,
+ output,
+ n_rows,
+ stride);
+ if (do_evaluate)
+ {
+ input += dofs_per_component_on_cell;
+ output += dofs_per_component_on_face;
+ }
+ else
+ {
+ output += dofs_per_component_on_cell;
+ input += dofs_per_component_on_face;
+ }
}
}
else if (face_direction < dim)
if (integrate)
for (unsigned int q = 0; q < n_q_points; ++q)
tmp_values[q] =
- gradients_quad[(c * dim + d) * n_q_points + orientation[q]];
+ gradients_quad[(c * n_q_points + orientation[q]) * dim + d];
else
for (unsigned int q = 0; q < n_q_points; ++q)
tmp_values[orientation[q]] =
- gradients_quad[(c * dim + d) * n_q_points + q];
+ gradients_quad[(c * n_q_points + q) * dim + d];
for (unsigned int q = 0; q < n_q_points; ++q)
- gradients_quad[(c * dim + d) * n_q_points + q] = tmp_values[q];
+ gradients_quad[(c * n_q_points + q) * dim + d] = tmp_values[q];
}
if (flag & EvaluationFlags::hessians)
{
Assert(gradients_quad != nullptr, ExcInternalError());
if (integrate)
for (unsigned int q = 0; q < n_q_points; ++q)
- tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
- orientation[q]][v];
+ tmp_values[q] =
+ gradients_quad[(c * n_q_points + orientation[q]) * dim + d]
+ [v];
else
for (unsigned int q = 0; q < n_q_points; ++q)
tmp_values[orientation[q]] =
- gradients_quad[(c * dim + d) * n_q_points + q][v];
+ gradients_quad[(c * n_q_points + q) * dim + d][v];
for (unsigned int q = 0; q < n_q_points; ++q)
- gradients_quad[(c * dim + d) * n_q_points + q][v] =
+ gradients_quad[(c * n_q_points + q) * dim + d][v] =
tmp_values[q];
}
if (flag & EvaluationFlags::hessians)
n_dofs,
n_q_points);
- eval.template gradients<0, true, false>(
- values_dofs_actual_ptr, gradients_quad_ptr);
-
- gradients_quad_ptr += n_q_points;
+ eval.template gradients<0, true, false, dim>(
+ values_dofs_actual_ptr, gradients_quad_ptr + d);
}
+ gradients_quad_ptr += n_q_points * dim;
values_dofs_actual_ptr += n_dofs;
}
}
if (!(integration_flag & EvaluationFlags::values) &&
d == 0)
- eval.template gradients<0, false, false>(
- gradients_quad_ptr, values_dofs_actual_ptr);
+ eval.template gradients<0, false, false, dim>(
+ gradients_quad_ptr + d, values_dofs_actual_ptr);
else
- eval.template gradients<0, false, true>(
- gradients_quad_ptr, values_dofs_actual_ptr);
-
- gradients_quad_ptr += n_q_points;
+ eval.template gradients<0, false, true, dim>(
+ gradients_quad_ptr + d, values_dofs_actual_ptr);
}
+ gradients_quad_ptr += n_q_points * dim;
values_dofs_actual_ptr += n_dofs;
}
}
projects / dealii.git / commitdiff