template <int dim, typename Number>
struct FEFaceEvaluationImplEvaluateSelector
{
- template <int fe_degree, int n_q_points_1d>
static bool
- run(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags evaluation_flag,
- const Number *values_dofs,
- FEEvaluationData<dim, Number, true> &fe_eval)
+ evaluate_tensor_none(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval)
{
const auto &shape_info = fe_eval.get_shape_info();
const auto &shape_data = shape_info.data.front();
using Number2 =
typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
- if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
- {
- Assert((fe_eval.get_dof_access_index() ==
- MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- fe_eval.is_interior_face() == false) == false,
- ExcNotImplemented());
+ Assert((fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) == false,
+ ExcNotImplemented());
- const unsigned int face_no = fe_eval.get_face_no();
- const unsigned int face_orientation = fe_eval.get_face_orientation();
- const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
- const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
+ const unsigned int face_no = fe_eval.get_face_no();
+ const unsigned int face_orientation = fe_eval.get_face_orientation();
+ const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
+ const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
- using Eval =
- EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
+ using Eval =
+ EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
- if (evaluation_flag & EvaluationFlags::values)
- {
- const auto *const shape_values =
- &shape_data.shape_values_face(face_no, face_orientation, 0);
+ if (evaluation_flag & EvaluationFlags::values)
+ {
+ const auto *const shape_values =
+ &shape_data.shape_values_face(face_no, face_orientation, 0);
- auto *values_quad_ptr = fe_eval.begin_values();
- auto *values_dofs_actual_ptr = values_dofs;
+ auto *values_quad_ptr = fe_eval.begin_values();
+ auto *values_dofs_actual_ptr = values_dofs;
- Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
- for (unsigned int c = 0; c < n_components; ++c)
- {
- eval.template values<0, true, false>(values_dofs_actual_ptr,
- values_quad_ptr);
+ Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ eval.template values<0, true, false>(values_dofs_actual_ptr,
+ values_quad_ptr);
- values_quad_ptr += n_q_points;
- values_dofs_actual_ptr += n_dofs;
- }
+ values_quad_ptr += n_q_points;
+ values_dofs_actual_ptr += n_dofs;
}
+ }
- if (evaluation_flag & EvaluationFlags::gradients)
- {
- auto *gradients_quad_ptr = fe_eval.begin_gradients();
- const auto *values_dofs_actual_ptr = values_dofs;
+ if (evaluation_flag & EvaluationFlags::gradients)
+ {
+ auto *gradients_quad_ptr = fe_eval.begin_gradients();
+ const auto *values_dofs_actual_ptr = values_dofs;
- std::array<const Number2 *, dim> shape_gradients;
- for (unsigned int d = 0; d < dim; ++d)
- shape_gradients[d] = &shape_data.shape_gradients_face(
- face_no, face_orientation, d, 0);
+ std::array<const Number2 *, dim> shape_gradients;
+ for (unsigned int d = 0; d < dim; ++d)
+ shape_gradients[d] =
+ &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
- for (unsigned int c = 0; c < n_components; ++c)
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ for (unsigned int d = 0; d < dim; ++d)
{
- for (unsigned int d = 0; d < dim; ++d)
- {
- Eval eval(nullptr,
- shape_gradients[d],
- nullptr,
- n_dofs,
- 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;
+ Eval eval(
+ nullptr, shape_gradients[d], nullptr, n_dofs, 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;
}
-
- Assert(!(evaluation_flag & EvaluationFlags::hessians),
- ExcNotImplemented());
-
- return true;
}
- const unsigned int dofs_per_face =
- fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
- Utilities::pow(shape_data.fe_degree + 1, dim - 1);
+ Assert(!(evaluation_flag & EvaluationFlags::hessians),
+ ExcNotImplemented());
- // Note: we always keep storage of values, 1st and 2nd derivatives in an
- // array, so reserve space for all three here
- Number *temp = fe_eval.get_scratch_data().begin();
- Number *scratch_data = temp + 3 * n_components * dofs_per_face;
-
- bool use_vectorization = true;
+ return true;
+ }
- if (fe_eval.get_dof_access_index() ==
- MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
- for (unsigned int v = 0; v < Number::size(); ++v)
- if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
- fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
- use_vectorization = false;
+ template <int fe_degree>
+ static void
+ project_to_face(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool use_vectorization,
+ Number *temp,
+ Number *scratch_data)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
if (use_vectorization == false)
{
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+ const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
+
for (unsigned int v = 0; v < Number::size(); ++v)
{
// the loop breaks once an invalid_unsigned_int is hit for
// some faces might be at the boundaries but others not)
if (fe_eval.get_cell_ids()[v] == numbers::invalid_unsigned_int)
{
- for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
- ++i)
+ for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
temp[i][v] = 0;
continue;
}
scratch_data,
fe_eval.get_face_no(v));
- for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
- ++i)
+ for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
temp[i][v] = scratch_data[i][v];
}
}
values_dofs,
temp,
fe_eval.get_face_no());
+ }
+
+
+ template <int fe_degree, int n_q_points_1d>
+ static void
+ evaluate_in_face(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ Number *temp,
+ Number *scratch_data)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
const unsigned int subface_index = fe_eval.get_subface_index();
constexpr unsigned int n_q_points_1d_actual =
fe_eval.begin_hessians(),
scratch_data,
subface_index);
+ }
+
+ static void
+ adjust_quadrature_for_face_orientation(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool use_vectorization,
+ Number *temp)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
if (use_vectorization == false)
{
n_components,
v,
evaluation_flag,
- &fe_eval.get_shape_info().face_orientations_quad(
+ &shape_info.face_orientations_quad(
fe_eval.get_face_orientation(v), 0),
false,
shape_info.n_q_points_face,
dim,
n_components,
evaluation_flag,
- &fe_eval.get_shape_info().face_orientations_quad(
- fe_eval.get_face_orientation(), 0),
+ &shape_info.face_orientations_quad(fe_eval.get_face_orientation(), 0),
false,
shape_info.n_q_points_face,
temp,
fe_eval.begin_values(),
fe_eval.begin_gradients(),
fe_eval.begin_hessians());
+ }
+
+
+
+ template <int fe_degree, int n_q_points_1d>
+ static bool
+ evaluate_tensor(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ // Note: we always keep storage of values, 1st and 2nd derivatives in an
+ // array, so reserve space for all three here
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ bool use_vectorization = true;
+ if (fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+ for (unsigned int v = 0; v < Number::size(); ++v)
+ if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
+ fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
+ use_vectorization = false;
+
+ project_to_face<fe_degree>(n_components,
+ evaluation_flag,
+ values_dofs,
+ fe_eval,
+ use_vectorization,
+ temp,
+ scratch_data);
+
+ evaluate_in_face<fe_degree, n_q_points_1d>(
+ n_components, evaluation_flag, fe_eval, temp, scratch_data);
+
+ if (dim == 3)
+ adjust_quadrature_for_face_orientation(
+ n_components, evaluation_flag, fe_eval, use_vectorization, temp);
return false;
}
+
+ template <int fe_degree, int n_q_points_1d>
+ static bool
+ run(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+
+ if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
+ return evaluate_tensor_none(n_components,
+ evaluation_flag,
+ values_dofs,
+ fe_eval);
+ else
+ return evaluate_tensor<fe_degree, n_q_points_1d>(n_components,
+ evaluation_flag,
+ values_dofs,
+ fe_eval);
+ }
};
template <int dim, typename Number>
- struct FEFaceEvaluationImplIntegrateSelector
+ struct FEFaceEvaluationImplProjectToFaceSelector
+ {
+ template <int fe_degree>
+ static bool
+ run(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ // Note: we always keep storage of values, 1st and 2nd derivatives in an
+ // array, so reserve space for all three here
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ bool use_vectorization = true;
+ if (fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+ for (unsigned int v = 0; v < Number::size(); ++v)
+ if (fe_eval.get_cell_ids()[v] != numbers::invalid_unsigned_int &&
+ fe_eval.get_face_no(v) != fe_eval.get_face_no(0))
+ use_vectorization = false;
+
+ FEFaceEvaluationImplEvaluateSelector<dim, Number>::
+ template project_to_face<fe_degree>(n_components,
+ evaluation_flag,
+ values_dofs,
+ fe_eval,
+ use_vectorization,
+ temp,
+ scratch_data);
+
+ return false;
+ }
+ };
+
+
+
+ template <int dim, typename Number>
+ struct FEFaceEvaluationImplEvaluateInFaceSelector
{
template <int fe_degree, int n_q_points_1d>
static bool
run(const unsigned int n_components,
- const EvaluationFlags::EvaluationFlags integration_flag,
- Number *values_dofs,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
FEEvaluationData<dim, Number, true> &fe_eval)
{
const auto &shape_info = fe_eval.get_shape_info();
const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ // Note: we always keep storage of values, 1st and 2nd derivatives in an
+ // array, so reserve space for all three here
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ FEFaceEvaluationImplEvaluateSelector<dim, Number>::
+ template evaluate_in_face<fe_degree, n_q_points_1d>(
+ n_components, evaluation_flag, fe_eval, temp, scratch_data);
+
+ return false;
+ }
+ };
+
+
+
+ template <int dim, typename Number>
+ struct FEFaceEvaluationImplIntegrateSelector
+ {
+ static bool
+ integrate_tensor_none(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool sum_into_values)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
using Number2 =
typename FEEvaluationData<dim, Number, true>::shape_info_number_type;
- if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
- {
- Assert((fe_eval.get_dof_access_index() ==
- MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- fe_eval.is_interior_face() == false) == false,
- ExcNotImplemented());
-
- const unsigned int face_no = fe_eval.get_face_no();
- const unsigned int face_orientation = fe_eval.get_face_orientation();
- const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
- const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
+ Assert((fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) == false,
+ ExcNotImplemented());
- using Eval =
- EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
+ const unsigned int face_no = fe_eval.get_face_no();
+ const unsigned int face_orientation = fe_eval.get_face_orientation();
+ const std::size_t n_dofs = shape_info.dofs_per_component_on_cell;
+ const std::size_t n_q_points = shape_info.n_q_points_faces[face_no];
- if (integration_flag & EvaluationFlags::values)
- {
- const auto *const shape_values =
- &shape_data.shape_values_face(face_no, face_orientation, 0);
+ using Eval =
+ EvaluatorTensorProduct<evaluate_general, 1, 0, 0, Number, Number2>;
- auto *values_quad_ptr = fe_eval.begin_values();
- auto *values_dofs_actual_ptr = values_dofs;
+ if (integration_flag & EvaluationFlags::values)
+ {
+ const auto *const shape_values =
+ &shape_data.shape_values_face(face_no, face_orientation, 0);
- Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
- for (unsigned int c = 0; c < n_components; ++c)
- {
- eval.template values<0, false, false>(values_quad_ptr,
- values_dofs_actual_ptr);
+ auto *values_quad_ptr = fe_eval.begin_values();
+ auto *values_dofs_actual_ptr = values_dofs;
- values_quad_ptr += n_q_points;
- values_dofs_actual_ptr += n_dofs;
- }
+ Eval eval(shape_values, nullptr, nullptr, n_dofs, n_q_points);
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ if (sum_into_values)
+ eval.template values<0, false, true>(values_quad_ptr,
+ values_dofs_actual_ptr);
+ else
+ eval.template values<0, false, false>(values_quad_ptr,
+ values_dofs_actual_ptr);
+ values_quad_ptr += n_q_points;
+ values_dofs_actual_ptr += n_dofs;
}
+ }
- if (integration_flag & EvaluationFlags::gradients)
- {
- auto *gradients_quad_ptr = fe_eval.begin_gradients();
- auto *values_dofs_actual_ptr = values_dofs;
+ if (integration_flag & EvaluationFlags::gradients)
+ {
+ auto *gradients_quad_ptr = fe_eval.begin_gradients();
+ auto *values_dofs_actual_ptr = values_dofs;
- std::array<const Number2 *, dim> shape_gradients;
- for (unsigned int d = 0; d < dim; ++d)
- shape_gradients[d] = &shape_data.shape_gradients_face(
- face_no, face_orientation, d, 0);
+ std::array<const Number2 *, dim> shape_gradients;
+ for (unsigned int d = 0; d < dim; ++d)
+ shape_gradients[d] =
+ &shape_data.shape_gradients_face(face_no, face_orientation, d, 0);
- for (unsigned int c = 0; c < n_components; ++c)
+ for (unsigned int c = 0; c < n_components; ++c)
+ {
+ for (unsigned int d = 0; d < dim; ++d)
{
- for (unsigned int d = 0; d < dim; ++d)
- {
- Eval eval(nullptr,
- shape_gradients[d],
- nullptr,
- n_dofs,
- n_q_points);
-
- if (!(integration_flag & EvaluationFlags::values) &&
- d == 0)
- eval.template gradients<0, false, false, dim>(
- gradients_quad_ptr + d, values_dofs_actual_ptr);
- else
- 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;
+ Eval eval(
+ nullptr, shape_gradients[d], nullptr, n_dofs, n_q_points);
+
+ if (!sum_into_values &&
+ !(integration_flag & EvaluationFlags::values) && d == 0)
+ eval.template gradients<0, false, false, dim>(
+ gradients_quad_ptr + d, values_dofs_actual_ptr);
+ else
+ 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;
}
-
- Assert(!(integration_flag & EvaluationFlags::hessians),
- ExcNotImplemented());
-
- return true;
}
- const unsigned int dofs_per_face =
- fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
- Utilities::pow(shape_data.fe_degree + 1, dim - 1);
-
- Number *temp = fe_eval.get_scratch_data().begin();
- Number *scratch_data = temp + 3 * n_components * dofs_per_face;
+ Assert(!(integration_flag & EvaluationFlags::hessians),
+ ExcNotImplemented());
- bool use_vectorization = true;
+ return true;
+ }
- if (fe_eval.get_dof_access_index() ==
- MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
- use_vectorization =
- fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
- std::all_of(fe_eval.get_cell_ids().begin() + 1,
- fe_eval.get_cell_ids().end(),
- [&](const auto &v) {
- return v == fe_eval.get_cell_ids()[0] ||
- v == numbers::invalid_unsigned_int;
- });
+ static void
+ adjust_quadrature_for_face_orientation(
+ const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool use_vectorization,
+ Number *temp)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
if (use_vectorization == false)
{
fe_eval.begin_values(),
fe_eval.begin_gradients(),
fe_eval.begin_hessians());
+ }
+
+ template <int fe_degree, int n_q_points_1d>
+ static void
+ integrate_in_face(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ Number *temp,
+ Number *scratch_data)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
const unsigned int n_q_points_1d_actual =
fe_degree > -1 ? n_q_points_1d : 0;
fe_eval.begin_hessians(),
scratch_data,
subface_index);
+ }
+
+ template <int fe_degree>
+ static void
+ collect_from_face(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool use_vectorization,
+ const Number *temp,
+ Number *scratch_data,
+ const bool sum_into_values)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+ const unsigned int dofs_per_face = n_components * dofs_per_comp_face;
if (use_vectorization == false)
{
template interpolate<false, false>(n_components,
integration_flag,
shape_info,
- values_dofs,
+ temp,
scratch_data,
fe_eval.get_face_no(v));
- for (unsigned int i = 0; i < 3 * n_components * dofs_per_face;
- ++i)
- temp[i][v] = scratch_data[i][v];
+ if (sum_into_values)
+ for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
+ values_dofs[i][v] += scratch_data[i][v];
+ else
+ for (unsigned int i = 0; i < 3 * dofs_per_face; ++i)
+ values_dofs[i][v] = scratch_data[i][v];
}
}
else
- FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
- template interpolate<false, false>(n_components,
- integration_flag,
- shape_info,
- temp,
- values_dofs,
- fe_eval.get_face_no());
+ {
+ if (sum_into_values)
+ FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+ template interpolate<false, true>(n_components,
+ integration_flag,
+ shape_info,
+ temp,
+ values_dofs,
+ fe_eval.get_face_no());
+ else
+ FEFaceNormalEvaluationImpl<dim, fe_degree, Number>::
+ template interpolate<false, false>(n_components,
+ integration_flag,
+ shape_info,
+ temp,
+ values_dofs,
+ fe_eval.get_face_no());
+ }
+ }
+
+ template <int fe_degree, int n_q_points_1d>
+ static bool
+ integrate_tensor(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool sum_into_values)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ bool use_vectorization = true;
+
+ if (fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+ use_vectorization =
+ fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
+ std::all_of(fe_eval.get_cell_ids().begin() + 1,
+ fe_eval.get_cell_ids().end(),
+ [&](const auto &v) {
+ return v == fe_eval.get_cell_ids()[0] ||
+ v == numbers::invalid_unsigned_int;
+ });
+
+ if (dim == 3)
+ adjust_quadrature_for_face_orientation(
+ n_components, integration_flag, fe_eval, use_vectorization, temp);
+
+ integrate_in_face<fe_degree, n_q_points_1d>(
+ n_components, integration_flag, fe_eval, temp, scratch_data);
+
+ collect_from_face<fe_degree>(n_components,
+ integration_flag,
+ values_dofs,
+ fe_eval,
+ use_vectorization,
+ temp,
+ scratch_data,
+ sum_into_values);
+
+ return false;
+ }
+
+ template <int fe_degree, int n_q_points_1d>
+ static bool
+ run(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool sum_into_values)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+
+ if (shape_info.element_type == MatrixFreeFunctions::tensor_none)
+ return integrate_tensor_none(n_components,
+ integration_flag,
+ values_dofs,
+ fe_eval,
+ sum_into_values);
+ else
+ return integrate_tensor<fe_degree, n_q_points_1d>(n_components,
+ integration_flag,
+ values_dofs,
+ fe_eval,
+ sum_into_values);
+ }
+ };
+
+
+
+ template <int dim, typename Number>
+ struct FEFaceEvaluationImplCollectFromFaceSelector
+ {
+ template <int fe_degree>
+ static bool
+ run(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+ Number *values_dofs,
+ FEEvaluationData<dim, Number, true> &fe_eval,
+ const bool sum_into_values)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ bool use_vectorization = true;
+
+ if (fe_eval.get_dof_access_index() ==
+ MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ fe_eval.is_interior_face() == false) // exterior faces in the ECL loop
+ use_vectorization =
+ fe_eval.get_cell_ids()[0] != numbers::invalid_unsigned_int &&
+ std::all_of(fe_eval.get_cell_ids().begin() + 1,
+ fe_eval.get_cell_ids().end(),
+ [&](const auto &v) {
+ return v == fe_eval.get_cell_ids()[0] ||
+ v == numbers::invalid_unsigned_int;
+ });
+
+ FEFaceEvaluationImplIntegrateSelector<dim, Number>::
+ template collect_from_face<fe_degree>(n_components,
+ integration_flag,
+ values_dofs,
+ fe_eval,
+ use_vectorization,
+ temp,
+ scratch_data,
+ sum_into_values);
+
+ return false;
+ }
+ };
+
+
+
+ template <int dim, typename Number>
+ struct FEFaceEvaluationImplIntegrateInFaceSelector
+ {
+ template <int fe_degree, int n_q_points_1d>
+ static bool
+ run(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags integration_flag,
+
+ FEEvaluationData<dim, Number, true> &fe_eval)
+ {
+ const auto &shape_info = fe_eval.get_shape_info();
+ const auto &shape_data = shape_info.data.front();
+
+ const unsigned int dofs_per_comp_face =
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::fixed_power<dim - 1>(shape_data.fe_degree + 1);
+
+ Number *temp = fe_eval.get_scratch_data().begin();
+ Number *scratch_data = temp + 3 * n_components * dofs_per_comp_face;
+
+ FEFaceEvaluationImplIntegrateSelector<dim, Number>::
+ template integrate_in_face<fe_degree, n_q_points_1d>(
+ n_components, integration_flag, fe_eval, temp, scratch_data);
+
return false;
}
};
DeclException1(
ExcFEPointEvaluationAccessToUninitializedMappingField,
std::string,
- << "You are requesting information from an FEPointEvaluation "
+ << "You are requesting information from an FEPointEvaluationBase "
<< "object for which this kind of information has not been computed. "
<< "What information these objects compute is determined by the update_* "
<< "flags you pass to MappingInfo() in the Constructor. Here, "
/**
- * This class provides an interface to the evaluation of interpolated solution
- * values and gradients on cells on arbitrary reference point positions. These
- * points can change from cell to cell, both with respect to their quantity as
- * well to the location. The two typical use cases are evaluations on
- * non-matching grids and particle simulations.
- *
- * The use of this class is similar to FEValues or FEEvaluation: The class is
- * first initialized to a cell by calling `FEPointEvaluation::reinit(cell,
- * unit_points)`, with the main difference to the other concepts that the
- * underlying points in reference coordinates need to be passed along. Then,
- * upon call to evaluate() or integrate(), the user can compute information at
- * the give points. Eventually, the access functions get_value() or
- * get_gradient() allow to query this information at a specific point index.
- *
- * The functionality is similar to creating an FEValues object with a
- * Quadrature object on the `unit_points` on every cell separately and then
- * calling FEValues::get_function_values or FEValues::get_function_gradients,
- * and for some elements and mappings this is what actually happens
- * internally. For specific combinations of Mapping and FiniteElement
- * realizations, however, there is a much more efficient implementation that
- * avoids the memory allocation and other expensive start-up cost of
- * FEValues. Currently, the functionality is specialized for mappings derived
- * from MappingQ and MappingCartesian and for finite elements with tensor
- * product structure that work with the
- * @ref matrixfree
- * module. In those cases, the cost implied
- * by this class is similar (or sometimes even somewhat lower) than using
- * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`.
+ * Base class of FEPointEvaluation and FEFacePointEvaluation. This class needs
+ * usually not be called in user code and does not have any public
+ * constructor. The usage is through the class
+ * FEPointEvaluation/FEFacePointEvaluation instead.
*/
template <int n_components_,
int dim,
int spacedim = dim,
typename Number = double>
-class FEPointEvaluation
+class FEPointEvaluationBase
{
public:
static constexpr unsigned int dimension = dim;
using interface_vectorized_gradient_type =
typename ETT::interface_vectorized_gradient_type;
+protected:
/**
* Constructor.
*
* objects, this parameter allows to select a range of `n_components`
* components starting from this parameter.
*/
- FEPointEvaluation(const Mapping<dim> &mapping,
- const FiniteElement<dim> &fe,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component = 0);
+ FEPointEvaluationBase(const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component = 0);
/**
* Constructor to make the present class able to re-use the geometry
- * data also used by other `FEPointEvaluation` objects.
+ * data also used by other `FEPointEvaluationBase` objects.
*
* @param mapping_info The MappingInfo class describes the geometry-related
* data for evaluating finite-element solutions. This object enables to
* @param first_selected_component For multi-component FiniteElement
* objects, this parameter allows to select a range of `n_components`
* components starting from this parameter.
+ *
+ * @param is_interior Defines if interior or exterior. Only makes sense for
+ * faces.
*/
- FEPointEvaluation(
+ FEPointEvaluationBase(
NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
const FiniteElement<dim> &fe,
- const unsigned int first_selected_component = 0);
+ const unsigned int first_selected_component = 0,
+ const bool is_interior = true);
/**
* Copy constructor.
*/
- FEPointEvaluation(FEPointEvaluation &other) noexcept;
+ FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept;
/**
* Move constructor.
*/
- FEPointEvaluation(FEPointEvaluation &&other) noexcept;
+ FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept;
/**
* Destructor.
*/
- ~FEPointEvaluation();
-
- /**
- * Set up the mapping information for the given cell, e.g., by computing the
- * Jacobian of the mapping for the given points if gradients of the functions
- * are requested.
- *
- * @param[in] cell An iterator to the current cell
- *
- * @param[in] unit_points List of points in the reference locations of the
- * current cell where the FiniteElement object should be
- * evaluated/integrated in the evaluate() and integrate() functions.
- */
- void
- reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> &unit_points);
-
- /**
- * Reinitialize the evaluator to point to the correct precomputed mapping of
- * the single cell in the MappingInfo object.
- */
- void
- reinit();
-
- /**
- * Reinitialize the evaluator to point to the correct precomputed mapping of
- * the cell in the MappingInfo object.
- */
- void
- reinit(const unsigned int cell_index);
-
- /**
- * Reinitialize the evaluator to point to the correct precomputed mapping of
- * the face in the MappingInfo object.
- */
- void
- reinit(const unsigned int cell_index, const unsigned int face_number);
-
- /**
- * This function interpolates the finite element solution, represented by
- * `solution_values`, on the cell and `unit_points` passed to reinit().
- *
- * @param[in] solution_values This array is supposed to contain the unknown
- * values on the element read out by
- * `FEEvaluation::read_dof_values(global_vector)`.
- *
- * @param[in] evaluation_flags Flags specifying which quantities should be
- * evaluated at the points.
- */
- template <std::size_t stride_view>
- void
- evaluate(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
- /**
- * This function interpolates the finite element solution, represented by
- * `solution_values`, on the cell and `unit_points` passed to reinit().
- *
- * @param[in] solution_values This array is supposed to contain the unknown
- * values on the element as returned by `cell->get_dof_values(global_vector,
- * solution_values)`.
- *
- * @param[in] evaluation_flags Flags specifying which quantities should be
- * evaluated at the points.
- */
- void
- evaluate(const ArrayView<const ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
- /**
- * This function multiplies the quantities passed in by previous
- * submit_value() or submit_gradient() calls by the value or gradient of the
- * test functions, and performs summation over all given points multiplied be
- * the Jacobian determinant times the quadrature weight (JxW).
- *
- * @param[out] solution_values This array will contain the result of the
- * integral, which can be used during
- * `FEEvaluation::set_dof_values(global_vector)` or
- * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
- * that for multi-component systems where only some of the components are
- * selected by the present class, the entries in `solution_values` not touched
- * by this class will be set to zero.
- *
- * @param[in] integration_flags Flags specifying which quantities should be
- * integrated at the points.
- *
- * @param[in] sum_into_values Flag specifying if the integrated values
- * should be summed into the solution values. For the default value
- * `sum_into_values=false` every value of @p solution_values is zeroed out.
- *
- */
- template <std::size_t stride_view>
- void
- integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
-
- /**
- * This function multiplies the quantities passed in by previous
- * submit_value() or submit_gradient() calls by the value or gradient of the
- * test functions, and performs summation over all given points multiplied be
- * the Jacobian determinant times the quadrature weight (JxW).
- *
- * @param[out] solution_values This array will contain the result of the
- * integral, which can be used to during
- * `cell->set_dof_values(solution_values, global_vector)` or
- * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
- * that for multi-component systems where only some of the components are
- * selected by the present class, the entries in `solution_values` not touched
- * by this class will be set to zero.
- *
- * @param[in] integration_flags Flags specifying which quantities should be
- * integrated at the points.
- *
- * @param[in] sum_into_values Flag specifying if the integrated values
- * should be summed into the solution values. For the default value
- * `sum_into_values=false` every value of @p solution_values is zeroed out.
- *
- */
- void
- integrate(const ArrayView<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
-
- /**
- * This function multiplies the quantities passed in by previous
- * submit_value() or submit_gradient() calls by the value or gradient of the
- * test functions, and performs summation over all given points. This is
- * similar to the integration of a bilinear form in terms of the test
- * function, with the difference that this formula does not include a `JxW`
- * factor (in contrast to the integrate function of this class). This allows
- * the class to naturally embed point information (e.g. particles) into a
- * finite element formulation.
- *
- * @param[out] solution_values This array will contain the result of the
- * integral, which can be used during
- * `FEEvaluation::set_dof_values(global_vector)` or
- * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
- * that for multi-component systems where only some of the components are
- * selected by the present class, the entries in `solution_values` not touched
- * by this class will be set to zero.
- *
- * @param[in] integration_flags Flags specifying which quantities should be
- * integrated at the points.
- *
- * @param[in] sum_into_values Flag specifying if the integrated values
- * should be summed into the solution values. For the default value
- * `sum_into_values=false` every value of @p solution_values is zeroed out.
- *
- */
- template <std::size_t stride_view>
- void
- test_and_sum(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
-
- /**
- * This function multiplies the quantities passed in by previous
- * submit_value() or submit_gradient() calls by the value or gradient of the
- * test functions, and performs summation over all given points. This is
- * similar to the integration of a bilinear form in terms of the test
- * function, with the difference that this formula does not include a `JxW`
- * factor (in contrast to the integrate function of this class). This allows
- * the class to naturally embed point information (e.g. particles) into a
- * finite element formulation.
- *
- * @param[out] solution_values This array will contain the result of the
- * integral, which can be used during
- * `cell->set_dof_values(solution_values, global_vector)` or
- * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
- * that for multi-component systems where only some of the components are
- * selected by the present class, the entries in `solution_values` not touched
- * by this class will be set to zero.
- *
- * @param[in] integration_flags Flags specifying which quantities should be
- * integrated at the points.
- *
- * @param[in] sum_into_values Flag specifying if the integrated values
- * should be summed into the solution values. For the default value
- * `sum_into_values=false` every value of @p solution_values is zeroed out.
- *
- */
- void
- test_and_sum(const ArrayView<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values = false);
+ ~FEPointEvaluationBase();
+public:
/**
* Return the value at quadrature point number @p point_index after a call to
* FEPointEvaluation::evaluate() with EvaluationFlags::values set, or
* the value that has been stored there with a call to
- * FEPointEvaluation::submit_value(). If the object is vector-valued, a
+ * FEPointEvaluationBase::submit_value(). If the object is vector-valued, a
* vector-valued return argument is given.
*/
const value_type &
/**
* Write a value to the field containing the values on points
* with component point_index. Access to the same field as through
- * get_value(). If applied before the function FEPointEvaluation::integrate()
- * with EvaluationFlags::values set is called, this specifies the value
- * which is tested by all basis function on the current cell and
- * integrated over.
+ * get_value(). If applied before the function
+ * FEPointEvaluation::integrate() with EvaluationFlags::values set is
+ * called, this specifies the value which is tested by all basis function on
+ * the current cell and integrated over.
*/
void
submit_value(const value_type &value, const unsigned int point_index);
* Return the gradient in real coordinates at the point with index
* `point_index` after a call to FEPointEvaluation::evaluate() with
* EvaluationFlags::gradients set, or the gradient that has been stored there
- * with a call to FEPointEvaluation::submit_gradient(). The gradient in real
- * coordinates is obtained by taking the unit gradient (also accessible via
- * get_unit_gradient()) and applying the inverse Jacobian of the mapping. If
- * the object is vector-valued, a vector-valued return argument is given.
+ * with a call to FEPointEvaluationBase::submit_gradient(). The gradient in
+ * real coordinates is obtained by taking the unit gradient (also accessible
+ * via get_unit_gradient()) and applying the inverse Jacobian of the mapping.
+ * If the object is vector-valued, a vector-valued return argument is given.
*/
const gradient_type &
get_gradient(const unsigned int point_index) const;
* Write a contribution that is tested by the gradient to the field
* containing the values on points with the given `point_index`. Access to
* the same field as through get_gradient(). If applied before the function
- * FEPointEvaluation::integrate(EvaluationFlags::gradients) is called, this
- * specifies what is tested by all basis function gradients on the current
- * cell and integrated over.
+ * FEPointEvaluation::integrate(EvaluationFlags::gradients) is called,
+ * this specifies what is tested by all basis function gradients on the
+ * current cell and integrated over.
*/
void
submit_gradient(const gradient_type &, const unsigned int point_index);
Number
JxW(const unsigned int point_index) const;
- /**
- * Return the normal vector. This class or the MappingInfo object passed to
- * this function needs to be constructed with UpdateFlags containing
- * `update_normal_vectors`.
- */
- Tensor<1, spacedim, Number>
- normal_vector(const unsigned int point_index) const;
-
/**
* Return the position in real coordinates of the given point index among
* the points passed to reinit().
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
quadrature_point_indices() const;
-private:
+protected:
static constexpr std::size_t n_lanes_user_interface =
internal::VectorizedArrayTrait<Number>::width();
static constexpr std::size_t n_lanes_internal =
do_reinit();
/**
- * Resizes necessary data fields, reads in and renumbers solution values.
- * Interpolates onto face if face path is selected.
+ * Number of quadrature batches of the current cell/face.
*/
- template <bool is_face_path, bool is_linear, std::size_t stride_view>
- void
- prepare_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
+ const unsigned int n_q_batches;
/**
- * Evaluates the actual interpolation on the cell or face for a quadrature
- * batch.
+ * Number of quadrature points/batches of the current cell/face.
*/
- template <bool is_face_path, bool is_linear, std::size_t stride_view>
- void
- compute_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags,
- const unsigned int n_shapes,
- const unsigned int qb,
- vectorized_value_type &value,
- interface_vectorized_gradient_type &gradient);
+ const unsigned int n_q_points;
/**
- * Fast path of the evaluate function.
+ * Number of quadrature points of the current cell/face.
*/
- template <bool is_face_path, bool is_linear, std::size_t stride_view>
- void
- evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
+ const unsigned int n_q_points_scalar;
/**
- * Slow path of the evaluate function using FEValues.
- */
- template <std::size_t stride_view>
- void
- evaluate_slow(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
- /**
- * Integrates the product of the data passed in by submit_value() and
- * submit_gradient() with the values or gradients of test functions on the
- * cell or face for a given quadrature batch.
- */
- template <bool is_face_path, bool is_linear>
- void
- compute_integrate_fast(
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const unsigned int n_shapes,
- const unsigned int qb,
- const vectorized_value_type value,
- const interface_vectorized_gradient_type gradient,
- vectorized_value_type *solution_values_vectorized_linear);
-
- /**
- * Addition across the lanes of VectorizedArray as accumulated by the
- * compute_integrate_fast_function(), writing the sum into the result vector.
- * Applies face contributions to cell contributions for face path.
- */
- template <bool is_face_path, bool is_linear, std::size_t stride_view>
- void
- finish_integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- vectorized_value_type *solution_values_vectorized_linear,
- const bool sum_into_values);
-
- /**
- * Fast path of the integrate function.
- */
- template <bool do_JxW,
- bool is_face_path,
- bool is_linear,
- std::size_t stride_view>
- void
- integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values);
-
- /**
- * Slow path of the integrate function using FEValues.
- */
- template <bool do_JxW, std::size_t stride_view>
- void
- integrate_slow(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values);
-
- /**
- * Implementation of the integrate/test_and_sum function.
- */
- template <bool do_JxW, std::size_t stride_view>
- void
- do_integrate(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values);
-
- /**
- * Number of quadrature batches of the current cell/face.
- */
- const unsigned int n_q_batches;
-
- /**
- * Number of quadrature points/batches of the current cell/face.
- */
- const unsigned int n_q_points;
-
- /**
- * Number of quadrature points of the current cell/face.
- */
- const unsigned int n_q_points_scalar;
-
- /**
- * Pointer to the Mapping object passed to the constructor.
+ * Pointer to the Mapping object passed to the constructor.
*/
SmartPointer<const Mapping<dim, spacedim>> mapping;
AlignedVector<vectorized_value_type> solution_renumbered_vectorized;
/**
- * Temporary array for the use_face_path path (scalar).
+ * Temporary array for the face path (scalar).
*/
AlignedVector<ScalarNumber> scratch_data_scalar;
/**
* Pointer to first unit point batch of current face from MappingInfo,
- * set internally during do_reinit(). Needed for use_face_path path.
+ * set internally during do_reinit(). Needed for face path.
*/
const Point<dim - 1, VectorizedArrayType> *unit_point_faces_ptr;
unsigned int dofs_per_component_face;
/**
- * Bool indicating if use_face_path path should be chosen. Set during
- * do_reinit().
- */
- bool use_face_path;
-
- /**
- * Scalar ShapeInfo object needed for use_face_path path.
+ * Scalar ShapeInfo object needed for face path.
*/
internal::MatrixFreeFunctions::ShapeInfo<ScalarNumber> shape_info;
* reinit()) at the vectorized unit points on faces.
*/
AlignedVector<dealii::ndarray<VectorizedArrayType, 2, dim - 1>> shapes_faces;
+
+ const bool is_interior;
};
// ----------------------- template and inline function ----------------------
template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
- const Mapping<dim> &mapping,
- const FiniteElement<dim> &fe,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component)
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component)
: n_q_batches(numbers::invalid_unsigned_int)
, n_q_points(numbers::invalid_unsigned_int)
, n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping)
, fe(&fe)
, JxW_ptr(nullptr)
- , use_face_path(false)
, update_flags(update_flags)
, mapping_info_on_the_fly(
std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
, current_cell_index(numbers::invalid_unsigned_int)
, current_face_number(numbers::invalid_unsigned_int)
, is_reinitialized(false)
+ , is_interior(true)
{
setup(first_selected_component);
}
template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
- NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
- const FiniteElement<dim> &fe,
- const unsigned int first_selected_component)
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const unsigned int first_selected_component,
+ const bool is_interior)
: n_q_batches(numbers::invalid_unsigned_int)
, n_q_points(numbers::invalid_unsigned_int)
, n_q_points_scalar(numbers::invalid_unsigned_int)
, mapping(&mapping_info.get_mapping())
, fe(&fe)
, JxW_ptr(nullptr)
- , use_face_path(false)
, update_flags(mapping_info.get_update_flags())
, mapping_info(&mapping_info)
, current_cell_index(numbers::invalid_unsigned_int)
, current_face_number(numbers::invalid_unsigned_int)
, is_reinitialized(false)
+ , is_interior(is_interior)
{
setup(first_selected_component);
connection_is_reinitialized = mapping_info.connect_is_reinitialized(
template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
- FEPointEvaluation<n_components_, dim, spacedim, Number> &other) noexcept
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ FEPointEvaluationBase<n_components_, dim, spacedim, Number> &other) noexcept
: n_q_batches(other.n_q_batches)
, n_q_points(other.n_q_points)
, n_q_points_scalar(other.n_q_points_scalar)
, gradients(other.gradients)
, dofs_per_component(other.dofs_per_component)
, dofs_per_component_face(other.dofs_per_component_face)
- , use_face_path(false)
, component_in_base_element(other.component_in_base_element)
, nonzero_shape_function_component(other.nonzero_shape_function_component)
, update_flags(other.update_flags)
, is_reinitialized(false)
, shapes(other.shapes)
, shapes_faces(other.shapes_faces)
+ , is_interior(other.is_interior)
{
connection_is_reinitialized = mapping_info->connect_is_reinitialized(
[this]() { this->is_reinitialized = false; });
template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
- FEPointEvaluation<n_components_, dim, spacedim, Number> &&other) noexcept
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ FEPointEvaluationBase(
+ FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+ &&other) noexcept
: n_q_batches(other.n_q_batches)
, n_q_points(other.n_q_points)
, n_q_points_scalar(other.n_q_points_scalar)
, gradients(other.gradients)
, dofs_per_component(other.dofs_per_component)
, dofs_per_component_face(other.dofs_per_component_face)
- , use_face_path(false)
, component_in_base_element(other.component_in_base_element)
, nonzero_shape_function_component(other.nonzero_shape_function_component)
, update_flags(other.update_flags)
, is_reinitialized(false)
, shapes(other.shapes)
, shapes_faces(other.shapes_faces)
+ , is_interior(other.is_interior)
{
connection_is_reinitialized = mapping_info->connect_is_reinitialized(
[this]() { this->is_reinitialized = false; });
template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::~FEPointEvaluation()
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ ~FEPointEvaluationBase()
{
connection_is_reinitialized.disconnect();
}
template <int n_components_, int dim, int spacedim, typename Number>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::setup(
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::setup(
const unsigned int first_selected_component)
{
AssertIndexRange(first_selected_component + n_components,
poly[1].value(1.) == 1.) &&
(fe->n_components() == n_components);
- const unsigned int size_face = 2 * dofs_per_component_face;
- const unsigned int size_cell = dofs_per_component;
+ const unsigned int size_face = 3 * dofs_per_component_face * n_components;
+ const unsigned int size_cell = dofs_per_component * n_components;
scratch_data_scalar.resize(size_face + size_cell);
solution_renumbered.resize(dofs_per_component);
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const ArrayView<const Point<dim>> &unit_points)
-{
- // reinit is only allowed for mapping computation on the fly
- AssertThrow(mapping_info_on_the_fly.get() != nullptr, ExcNotImplemented());
-
- mapping_info->reinit(cell, unit_points);
-
- if (!fast_path)
- {
- fe_values = std::make_shared<FEValues<dim, spacedim>>(
- *mapping,
- *fe,
- Quadrature<dim>(
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
- update_flags);
- fe_values->reinit(cell);
- }
-
- if (use_linear_path)
- do_reinit<false, true>();
- else
- do_reinit<false, false>();
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
-{
- current_cell_index = numbers::invalid_unsigned_int;
- current_face_number = numbers::invalid_unsigned_int;
-
- if (use_linear_path)
- do_reinit<false, true>();
- else
- do_reinit<false, false>();
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
- const unsigned int cell_index)
-{
- current_cell_index = cell_index;
- current_face_number = numbers::invalid_unsigned_int;
-
- if (use_linear_path)
- do_reinit<false, true>();
- else
- do_reinit<false, false>();
-
- if (!fast_path)
- {
- std::vector<Point<dim>> unit_points(n_q_points_scalar);
-
- for (unsigned int v = 0; v < n_q_points_scalar; ++v)
- for (unsigned int d = 0; d < dim; ++d)
- unit_points[v][d] =
- unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
-
- fe_values = std::make_shared<FEValues<dim, spacedim>>(
- *mapping,
- *fe,
- Quadrature<dim>(
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
- update_flags);
-
- fe_values->reinit(mapping_info->get_cell_iterator(current_cell_index));
- }
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
- const unsigned int cell_index,
- const unsigned int face_number)
-{
- current_cell_index = cell_index;
- current_face_number = face_number;
-
- if (use_linear_path)
- do_reinit<true, true>();
- else
- do_reinit<true, false>();
-}
-
-
-
template <int n_components_, int dim, int spacedim, typename Number>
template <bool is_face, bool is_linear>
inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::do_reinit()
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::do_reinit()
{
const unsigned int geometry_index =
mapping_info->template compute_geometry_index_offset<is_face>(
return;
}
- // use face path if mapping_info in face state and number of quadrature points
- // is large enough
- use_face_path =
- is_face && (mapping_info->is_face_state() && n_q_points_scalar >= 6);
-
// set unit point pointer
const unsigned int unit_point_offset =
mapping_info->compute_unit_point_index_offset(geometry_index);
- if (use_face_path)
+ if (is_face)
unit_point_faces_ptr =
mapping_info->get_unit_point_faces(unit_point_offset);
else
if (update_flags_mapping & UpdateFlags::update_quadrature_points)
real_point_ptr = mapping_info->get_real_point(data_offset);
if (update_flags_mapping & UpdateFlags::update_jacobians)
- jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset);
+ jacobian_ptr =
+ mapping_info->get_jacobian(compressed_data_offset, is_interior);
if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
inverse_jacobian_ptr =
- mapping_info->get_inverse_jacobian(compressed_data_offset);
+ mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
if (update_flags_mapping & UpdateFlags::update_normal_vectors)
normal_ptr = mapping_info->get_normal_vector(data_offset);
if (update_flags_mapping & UpdateFlags::update_JxW_values)
JxW_ptr = mapping_info->get_JxW(data_offset);
#else
real_point_ptr = mapping_info->get_real_point(data_offset);
- jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset);
+ jacobian_ptr =
+ mapping_info->get_jacobian(compressed_data_offset, is_interior);
inverse_jacobian_ptr =
- mapping_info->get_inverse_jacobian(compressed_data_offset);
+ mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
normal_ptr = mapping_info->get_normal_vector(data_offset);
JxW_ptr = mapping_info->get_JxW(data_offset);
#endif
const std::size_t n_shapes = poly.size();
for (unsigned int qb = 0; qb < n_q_batches; ++qb)
- if (use_face_path)
+ if (is_face)
{
if (dim > 1)
{
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face_path, bool is_linear, std::size_t stride_view>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::prepare_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags)
+inline const typename FEPointEvaluationBase<n_components_,
+ dim,
+ spacedim,
+ Number>::value_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_value(
+ const unsigned int point_index) const
{
- const unsigned int dofs_per_comp =
- is_linear ? Utilities::pow(2, dim) : dofs_per_component;
+ AssertIndexRange(point_index, values.size());
+ return values[point_index];
+}
- for (unsigned int comp = 0; comp < n_components; ++comp)
- {
- const std::size_t offset =
- (component_in_base_element + comp) * dofs_per_comp;
- if (is_face_path)
- {
- const ScalarNumber *input;
- if (is_linear || renumber.empty())
- {
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- scratch_data_scalar[i] = solution_values[i + offset];
- input = scratch_data_scalar.data();
- }
- else
- {
- const unsigned int *renumber_ptr = renumber.data() + offset;
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- scratch_data_scalar[i] = solution_values[renumber_ptr[i]];
- input = scratch_data_scalar.data();
- }
- ScalarNumber *output = scratch_data_scalar.begin() + dofs_per_comp;
+template <int n_components_, int dim, int spacedim, typename Number>
+inline const typename FEPointEvaluationBase<n_components_,
+ dim,
+ spacedim,
+ Number>::gradient_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_gradient(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, gradients.size());
+ return gradients[point_index];
+}
- internal::FEFaceNormalEvaluationImpl<dim, -1, ScalarNumber>::
- template interpolate<true, false>(1,
- evaluation_flags,
- shape_info,
- input,
- output,
- current_face_number);
- const unsigned int dofs_per_comp_face =
- is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
- for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
- ETT::read_value(output[i], comp, solution_renumbered[i]);
- }
- else
- {
- if ((is_linear && n_components == 1) || renumber.empty())
- {
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- ETT::read_value(solution_values[i + offset],
- comp,
- solution_renumbered[i]);
- }
- else
- {
- const unsigned int *renumber_ptr = renumber.data() + offset;
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- ETT::read_value(solution_values[renumber_ptr[i]],
- comp,
- solution_renumbered[i]);
- }
- }
- }
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_value(
+ const value_type &value,
+ const unsigned int point_index)
+{
+ AssertIndexRange(point_index, n_q_points);
+ values[point_index] = value;
}
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face_path, bool is_linear, std::size_t stride_view>
inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags,
- const unsigned int n_shapes,
- const unsigned int qb,
- vectorized_value_type &value,
- interface_vectorized_gradient_type &gradient)
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_gradient(
+ const gradient_type &gradient,
+ const unsigned int point_index)
{
- if (is_face_path)
- {
- if (evaluation_flags & EvaluationFlags::gradients)
- {
- const std::array<vectorized_value_type, dim + 1> interpolated_value =
- is_linear ?
- internal::evaluate_tensor_product_value_and_gradient_linear<
- dim - 1,
- scalar_value_type,
- VectorizedArrayType,
- 2>(solution_renumbered.data(), unit_point_faces_ptr[qb]) :
- internal::evaluate_tensor_product_value_and_gradient_shapes<
- dim - 1,
- scalar_value_type,
- VectorizedArrayType,
- 2,
- false>(shapes_faces.data() + qb * n_shapes,
- n_shapes,
- solution_renumbered.data());
+ AssertIndexRange(point_index, n_q_points);
+ gradients[point_index] = gradient;
+}
- value = interpolated_value[dim - 1];
- // reorder derivative from tangential/normal derivatives into tensor
- // in physical coordinates
- if (current_face_number / 2 == 0)
- {
- gradient[0] = interpolated_value[dim];
- if (dim > 1)
- gradient[1] = interpolated_value[0];
- if (dim > 2)
- gradient[2] = interpolated_value[1];
- }
- else if (current_face_number / 2 == 1)
- {
- if (dim > 1)
- gradient[1] = interpolated_value[dim];
- if (dim == 3)
- {
- gradient[0] = interpolated_value[1];
- gradient[2] = interpolated_value[0];
- }
- else if (dim == 2)
- gradient[0] = interpolated_value[0];
- else
- Assert(false, ExcInternalError());
- }
- else if (current_face_number / 2 == 2)
- {
- if (dim > 2)
- {
- gradient[0] = interpolated_value[0];
- gradient[1] = interpolated_value[1];
- gradient[2] = interpolated_value[dim];
- }
- else
- Assert(false, ExcInternalError());
- }
- else
- Assert(false, ExcInternalError());
- }
- else
- {
- value = is_linear ? internal::evaluate_tensor_product_value_linear<
- dim - 1,
- scalar_value_type,
- VectorizedArrayType>(solution_renumbered.data(),
- unit_point_faces_ptr[qb]) :
- internal::evaluate_tensor_product_value_shapes<
- dim - 1,
- scalar_value_type,
- VectorizedArrayType,
- false>(shapes_faces.data() + qb * n_shapes,
- n_shapes,
- solution_renumbered.data());
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline DerivativeForm<1, dim, spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::jacobian(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, n_q_points);
+ Assert(jacobian_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_jacobians"));
+ return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions::
+ GeometryType::affine ?
+ 0 :
+ point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline DerivativeForm<1, spacedim, dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::inverse_jacobian(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, n_q_points);
+ Assert(inverse_jacobian_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_inverse_jacobians"));
+ return inverse_jacobian_ptr
+ [cell_type <=
+ ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ?
+ 0 :
+ point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Number
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::JxW(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, n_q_points);
+ Assert(JxW_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_JxW_values"));
+ return JxW_ptr[point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Point<spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::real_point(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, n_q_points);
+ Assert(real_point_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_quadrature_points"));
+ return real_point_ptr[point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Point<dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::unit_point(
+ const unsigned int point_index) const
+{
+ AssertIndexRange(point_index, n_q_points);
+ Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
+ Point<dim, Number> unit_point;
+ for (unsigned int d = 0; d < dim; ++d)
+ unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+ unit_point_ptr[point_index / stride][d], point_index % stride);
+ return unit_point;
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+ quadrature_point_indices() const
+{
+ return {0U, n_q_points};
+}
+
+
+
+/**
+ * This class provides an interface to the evaluation of interpolated solution
+ * values and gradients on cells on arbitrary reference point positions. These
+ * points can change from cell to cell, both with respect to their quantity as
+ * well to the location. The two typical use cases are evaluations on
+ * non-matching grids and particle simulations.
+ *
+ * The use of this class is similar to FEValues or FEEvaluation: The class is
+ * first initialized to a cell by calling `FEPointEvaluation::reinit(cell,
+ * unit_points)`, with the main difference to the other concepts that the
+ * underlying points in reference coordinates need to be passed along. Then,
+ * upon call to evaluate() or integrate(), the user can compute information at
+ * the give points. Eventually, the access functions get_value() or
+ * get_gradient() allow to query this information at a specific point index.
+ *
+ * The functionality is similar to creating an FEValues object with a
+ * Quadrature object on the `unit_points` on every cell separately and then
+ * calling FEValues::get_function_values or FEValues::get_function_gradients,
+ * and for some elements and mappings this is what actually happens
+ * internally. For specific combinations of Mapping and FiniteElement
+ * realizations, however, there is a much more efficient implementation that
+ * avoids the memory allocation and other expensive start-up cost of
+ * FEValues. Currently, the functionality is specialized for mappings derived
+ * from MappingQ and MappingCartesian and for finite elements with tensor
+ * product structure that work with the
+ * @ref matrixfree
+ * module. In those cases, the cost implied
+ * by this class is similar (or sometimes even somewhat lower) than using
+ * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`.
+ */
+template <int n_components_,
+ int dim,
+ int spacedim = dim,
+ typename Number = double>
+class FEPointEvaluation
+ : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+{
+public:
+ static constexpr unsigned int dimension = dim;
+ static constexpr unsigned int n_components = n_components_;
+
+ using number_type = Number;
+
+ using ScalarNumber =
+ typename internal::VectorizedArrayTrait<Number>::value_type;
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
+ using ETT = typename internal::FEPointEvaluation::
+ EvaluatorTypeTraits<dim, n_components, Number>;
+ using value_type = typename ETT::value_type;
+ using scalar_value_type = typename ETT::scalar_value_type;
+ using vectorized_value_type = typename ETT::vectorized_value_type;
+ using gradient_type = typename ETT::gradient_type;
+ using interface_vectorized_gradient_type =
+ typename ETT::interface_vectorized_gradient_type;
+
+ /**
+ * Constructor.
+ *
+ * @param mapping The Mapping class describing the actual geometry of a cell
+ * passed to the evaluate() function.
+ *
+ * @param fe The FiniteElement object that is used for the evaluation, which
+ * is typically the same on all cells to be evaluated.
+ *
+ * @param update_flags Specify the quantities to be computed by the mapping
+ * during the call of reinit(). During evaluate() or integrate(), this data
+ * is queried to produce the desired result (e.g., the gradient of a finite
+ * element solution).
+ *
+ * @param first_selected_component For multi-component FiniteElement
+ * objects, this parameter allows to select a range of `n_components`
+ * components starting from this parameter.
+ */
+ FEPointEvaluation(const Mapping<dim> &mapping,
+ const FiniteElement<dim> &fe,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component = 0)
+ : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+ mapping,
+ fe,
+ update_flags,
+ first_selected_component)
+ {}
+
+ /**
+ * Constructor to make the present class able to re-use the geometry
+ * data also used by other `FEPointEvaluation` objects.
+ *
+ * @param mapping_info The MappingInfo class describes the geometry-related
+ * data for evaluating finite-element solutions. This object enables to
+ * construct such an object on the outside, possibly re-using it between
+ * several objects or between several calls to the same cell and unit points.
+ *
+ * @param fe The FiniteElement object that is used for the evaluation, which
+ * is typically the same on all cells to be evaluated.
+ *
+ * @param first_selected_component For multi-component FiniteElement
+ * objects, this parameter allows to select a range of `n_components`
+ * components starting from this parameter.
+ */
+ FEPointEvaluation(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const unsigned int first_selected_component = 0)
+ : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+ mapping_info,
+ fe,
+ first_selected_component)
+ {}
+
+ /**
+ * Set up the mapping information for the given cell, e.g., by computing the
+ * Jacobian of the mapping for the given points if gradients of the functions
+ * are requested.
+ *
+ * @param[in] cell An iterator to the current cell
+ *
+ * @param[in] unit_points List of points in the reference locations of the
+ * current cell where the FiniteElement object should be
+ * evaluated/integrated in the evaluate() and integrate() functions.
+ */
+ void
+ reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const ArrayView<const Point<dim>> &unit_points);
+
+ /**
+ * Reinitialize the evaluator to point to the correct precomputed mapping of
+ * the single cell in the MappingInfo object.
+ */
+ void
+ reinit();
+
+ /**
+ * Reinitialize the evaluator to point to the correct precomputed mapping of
+ * the cell in the MappingInfo object.
+ */
+ void
+ reinit(const unsigned int cell_index);
+
+
+ /**
+ * This function interpolates the finite element solution, represented by
+ * `solution_values`, on the cell and `unit_points` passed to reinit().
+ *
+ * @param[in] solution_values This array is supposed to contain the unknown
+ * values on the element read out by
+ * `FEEvaluation::read_dof_values(global_vector)`.
+ *
+ * @param[in] evaluation_flags Flags specifying which quantities should be
+ * evaluated at the points.
+ */
+ template <std::size_t stride_view>
+ void
+ evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * This function interpolates the finite element solution, represented by
+ * `solution_values`, on the cell and `unit_points` passed to reinit().
+ *
+ * @param[in] solution_values This array is supposed to contain the unknown
+ * values on the element as returned by `cell->get_dof_values(global_vector,
+ * solution_values)`.
+ *
+ * @param[in] evaluation_flags Flags specifying which quantities should be
+ * evaluated at the points.
+ */
+ void
+ evaluate(const ArrayView<const ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used during
+ * `FEEvaluation::set_dof_values(global_vector)` or
+ * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ template <std::size_t stride_view>
+ void
+ integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used to during
+ * `cell->set_dof_values(solution_values, global_vector)` or
+ * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ void
+ integrate(const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points. This is
+ * similar to the integration of a bilinear form in terms of the test
+ * function, with the difference that this formula does not include a `JxW`
+ * factor (in contrast to the integrate function of this class). This allows
+ * the class to naturally embed point information (e.g. particles) into a
+ * finite element formulation.
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used during
+ * `FEEvaluation::set_dof_values(global_vector)` or
+ * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ template <std::size_t stride_view>
+ void
+ test_and_sum(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points. This is
+ * similar to the integration of a bilinear form in terms of the test
+ * function, with the difference that this formula does not include a `JxW`
+ * factor (in contrast to the integrate function of this class). This allows
+ * the class to naturally embed point information (e.g. particles) into a
+ * finite element formulation.
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used during
+ * `cell->set_dof_values(solution_values, global_vector)` or
+ * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ void
+ test_and_sum(const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * Return the normal vector. This class or the MappingInfo object passed to
+ * this function needs to be constructed with UpdateFlags containing
+ * `update_normal_vectors`.
+ */
+ Tensor<1, spacedim, Number>
+ normal_vector(const unsigned int point_index) const;
+
+private:
+ static constexpr std::size_t n_lanes_user_interface =
+ internal::VectorizedArrayTrait<Number>::width();
+ static constexpr std::size_t n_lanes_internal =
+ internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+ static constexpr std::size_t stride =
+ internal::VectorizedArrayTrait<Number>::stride();
+
+ /**
+ * Resizes necessary data fields, reads in and renumbers solution values.
+ * Interpolates onto face if face path is selected.
+ */
+ template <bool is_linear, std::size_t stride_view>
+ void
+ prepare_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values);
+
+ /**
+ * Evaluates the actual interpolation on the cell or face for a quadrature
+ * batch.
+ */
+ template <bool is_linear, std::size_t stride_view>
+ void
+ compute_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags,
+ const unsigned int n_shapes,
+ const unsigned int qb,
+ vectorized_value_type &value,
+ interface_vectorized_gradient_type &gradient);
+
+ /**
+ * Fast path of the evaluate function.
+ */
+ template <bool is_linear, std::size_t stride_view>
+ void
+ evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * Slow path of the evaluate function using FEValues.
+ */
+ template <std::size_t stride_view>
+ void
+ evaluate_slow(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * Integrates the product of the data passed in by submit_value() and
+ * submit_gradient() with the values or gradients of test functions on the
+ * cell or face for a given quadrature batch.
+ */
+ template <bool is_linear>
+ void
+ compute_integrate_fast(
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const unsigned int n_shapes,
+ const unsigned int qb,
+ const vectorized_value_type value,
+ const interface_vectorized_gradient_type gradient,
+ vectorized_value_type *solution_values_vectorized_linear);
+
+ /**
+ * Addition across the lanes of VectorizedArray as accumulated by the
+ * compute_integrate_fast_function(), writing the sum into the result vector.
+ * Applies face contributions to cell contributions for face path.
+ */
+ template <bool is_linear, std::size_t stride_view>
+ void
+ finish_integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ vectorized_value_type *solution_values_vectorized_linear,
+ const bool sum_into_values);
+
+ /**
+ * Fast path of the integrate function.
+ */
+ template <bool do_JxW, bool is_linear, std::size_t stride_view>
+ void
+ integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+
+ /**
+ * Slow path of the integrate function using FEValues.
+ */
+ template <bool do_JxW, std::size_t stride_view>
+ void
+ integrate_slow(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+
+ /**
+ * Implementation of the integrate/test_and_sum function.
+ */
+ template <bool do_JxW, std::size_t stride_view>
+ void
+ do_integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+};
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
+{
+ this->current_cell_index = numbers::invalid_unsigned_int;
+ this->current_face_number = numbers::invalid_unsigned_int;
+
+ if (this->use_linear_path)
+ this->template do_reinit<false, true>();
+ else
+ this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const ArrayView<const Point<dim>> &unit_points)
+{
+ // reinit is only allowed for mapping computation on the fly
+ AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
+ ExcNotImplemented());
+
+ this->mapping_info->reinit(cell, unit_points);
+
+ if (!this->fast_path)
+ {
+ this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
+ *this->mapping,
+ *this->fe,
+ Quadrature<dim>(
+ std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+ this->update_flags);
+ this->fe_values->reinit(cell);
+ }
+
+ if (this->use_linear_path)
+ this->template do_reinit<false, true>();
+ else
+ this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+ const unsigned int cell_index)
+{
+ this->current_cell_index = cell_index;
+ this->current_face_number = numbers::invalid_unsigned_int;
+
+ if (this->use_linear_path)
+ this->template do_reinit<false, true>();
+ else
+ this->template do_reinit<false, false>();
+
+ if (!this->fast_path)
+ {
+ std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
+
+ for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
+ for (unsigned int d = 0; d < dim; ++d)
+ unit_points[v][d] =
+ this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
+
+ this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
+ *this->mapping,
+ *this->fe,
+ Quadrature<dim>(
+ std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+ this->update_flags);
+
+ this->fe_values->reinit(
+ this->mapping_info->get_cell_iterator(this->current_cell_index));
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
+{
+ if (!this->is_reinitialized)
+ reinit();
+
+ if (this->n_q_points == 0)
+ return;
+
+ Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
+
+ if (!((evaluation_flags & EvaluationFlags::values) ||
+ (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
+ return;
+
+ AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
+ if (this->fast_path)
+ {
+ if (this->use_linear_path)
+ evaluate_fast<true>(solution_values, evaluation_flags);
+ else
+ evaluate_fast<false>(solution_values, evaluation_flags);
+ }
+ else
+ evaluate_slow(solution_values, evaluation_flags);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+ const ArrayView<const ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
+{
+ evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ evaluation_flags);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ do_integrate<true>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ do_integrate<false>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::prepare_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values)
+{
+ const unsigned int dofs_per_comp =
+ is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ {
+ const std::size_t offset =
+ (this->component_in_base_element + comp) * dofs_per_comp;
+
+ if ((is_linear && n_components == 1) || this->renumber.empty())
+ {
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ ETT::read_value(solution_values[i + offset],
+ comp,
+ this->solution_renumbered[i]);
+ }
+ else
+ {
+ const unsigned int *renumber_ptr = this->renumber.data() + offset;
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ ETT::read_value(solution_values[renumber_ptr[i]],
+ comp,
+ this->solution_renumbered[i]);
+ }
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags,
+ const unsigned int n_shapes,
+ const unsigned int qb,
+ vectorized_value_type &value,
+ interface_vectorized_gradient_type &gradient)
+{
+ if (evaluation_flags & EvaluationFlags::gradients)
+ {
+ std::array<vectorized_value_type, dim + 1> result;
+ if constexpr (is_linear)
+ {
+ if constexpr (n_components == 1)
+ result =
+ internal::evaluate_tensor_product_value_and_gradient_linear<
+ dim,
+ scalar_value_type,
+ VectorizedArrayType,
+ 1,
+ stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
+ else
+ result =
+ internal::evaluate_tensor_product_value_and_gradient_linear(
+ this->solution_renumbered.data(), this->unit_point_ptr[qb]);
+ }
+ else
+ result = internal::evaluate_tensor_product_value_and_gradient_shapes<
+ dim,
+ scalar_value_type,
+ VectorizedArrayType,
+ 1,
+ false>(this->shapes.data() + qb * n_shapes,
+ n_shapes,
+ this->solution_renumbered.data());
+ gradient[0] = result[0];
+ if (dim > 1)
+ gradient[1] = result[1];
+ if (dim > 2)
+ gradient[2] = result[2];
+ value = result[dim];
+ }
+ else
+ {
+ if constexpr (is_linear)
+ {
+ if constexpr (n_components == 1)
+ value = internal::evaluate_tensor_product_value_linear<
+ dim,
+ scalar_value_type,
+ VectorizedArrayType,
+ stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
+ else
+ value = internal::evaluate_tensor_product_value_linear(
+ this->solution_renumbered.data(), this->unit_point_ptr[qb]);
+ }
+ else
+ value =
+ internal::evaluate_tensor_product_value_shapes<dim,
+ scalar_value_type,
+ VectorizedArrayType,
+ false>(
+ this->shapes.data() + qb * n_shapes,
+ n_shapes,
+ this->solution_renumbered.data());
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
+{
+ if (!(is_linear && n_components == 1))
+ prepare_evaluate_fast<is_linear>(solution_values);
+
+ // loop over quadrature batches qb
+ const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
+
+ for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
+ {
+ vectorized_value_type value;
+ interface_vectorized_gradient_type gradient;
+
+ compute_evaluate_fast<is_linear>(
+ solution_values, evaluation_flags, n_shapes, qb, value, gradient);
+
+ if (evaluation_flags & EvaluationFlags::values)
+ {
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ ETT::set_value(value, v, this->values[offset]);
+ }
+ if (evaluation_flags & EvaluationFlags::gradients)
+ {
+ Assert(this->update_flags & update_gradients ||
+ this->update_flags & update_inverse_jacobians,
+ ExcNotInitialized());
+
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ {
+ gradient_type unit_gradient;
+ ETT::set_gradient(gradient, v, unit_gradient);
+ this->gradients[offset] =
+ this->cell_type <=
+ internal::MatrixFreeFunctions::GeometryType::cartesian ?
+ apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
+ unit_gradient) :
+ apply_transformation(
+ this
+ ->inverse_jacobian_ptr[this->cell_type <=
+ internal::MatrixFreeFunctions::
+ GeometryType::affine ?
+ 0 :
+ offset]
+ .transpose(),
+ unit_gradient);
+ }
+ }
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_slow(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
+{
+ // slow path with FEValues
+ Assert(this->fe_values.get() != nullptr,
+ ExcMessage(
+ "Not initialized. Please call FEPointEvaluation::reinit()!"));
+
+ const std::size_t n_points = this->fe_values->get_quadrature().size();
+
+ if (evaluation_flags & EvaluationFlags::values)
+ {
+ this->values.resize(this->n_q_points);
+ std::fill(this->values.begin(), this->values.end(), value_type());
+ for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
+ {
+ const ScalarNumber value = solution_values[i];
+ for (unsigned int d = 0; d < n_components; ++d)
+ if (this->nonzero_shape_function_component[i][d] &&
+ (this->fe->is_primitive(i) || this->fe->is_primitive()))
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ ETT::access(this->values[qb],
+ v,
+ d,
+ this->fe_values->shape_value(i, q + v) * value);
+ else if (this->nonzero_shape_function_component[i][d])
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ ETT::access(this->values[qb],
+ v,
+ d,
+ this->fe_values->shape_value_component(i,
+ q + v,
+ d) *
+ value);
+ }
+ }
+
+ if (evaluation_flags & EvaluationFlags::gradients)
+ {
+ this->gradients.resize(this->n_q_points);
+ std::fill(this->gradients.begin(),
+ this->gradients.end(),
+ gradient_type());
+ for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
+ {
+ const ScalarNumber value = solution_values[i];
+ for (unsigned int d = 0; d < n_components; ++d)
+ if (this->nonzero_shape_function_component[i][d] &&
+ (this->fe->is_primitive(i) || this->fe->is_primitive()))
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ ETT::access(this->gradients[qb],
+ v,
+ d,
+ this->fe_values->shape_grad(i, q + v) * value);
+ else if (this->nonzero_shape_function_component[i][d])
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ ETT::access(
+ this->gradients[qb],
+ v,
+ d,
+ this->fe_values->shape_grad_component(i, q + v, d) * value);
}
}
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_integrate_fast(
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const unsigned int n_shapes,
+ const unsigned int qb,
+ const vectorized_value_type value,
+ const interface_vectorized_gradient_type gradient,
+ vectorized_value_type *solution_values_vectorized_linear)
+{
+ if (integration_flags & EvaluationFlags::gradients)
+ internal::integrate_tensor_product_value_and_gradient<
+ is_linear,
+ dim,
+ VectorizedArrayType,
+ vectorized_value_type>(this->shapes.data() + qb * n_shapes,
+ n_shapes,
+ &value,
+ gradient,
+ is_linear ?
+ solution_values_vectorized_linear :
+ this->solution_renumbered_vectorized.data(),
+ this->unit_point_ptr[qb],
+ qb != 0);
else
+ internal::integrate_tensor_product_value<is_linear,
+ dim,
+ VectorizedArrayType,
+ vectorized_value_type>(
+ this->shapes.data() + qb * n_shapes,
+ n_shapes,
+ value,
+ is_linear ? solution_values_vectorized_linear :
+ this->solution_renumbered_vectorized.data(),
+ this->unit_point_ptr[qb],
+ qb != 0);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ vectorized_value_type *solution_values_vectorized_linear,
+ const bool sum_into_values)
+{
+ if (!sum_into_values && this->fe->n_components() > n_components)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
+
+ const unsigned int dofs_per_comp =
+ is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
{
- if (evaluation_flags & EvaluationFlags::gradients)
+ const std::size_t offset =
+ (this->component_in_base_element + comp) * dofs_per_comp;
+
+ if (is_linear || this->renumber.empty())
{
- std::array<vectorized_value_type, dim + 1> result;
- if constexpr (is_linear)
- {
- if constexpr (n_components == 1)
- result =
- internal::evaluate_tensor_product_value_and_gradient_linear<
- dim,
- scalar_value_type,
- VectorizedArrayType,
- 1,
- stride_view>(solution_values.data(), unit_point_ptr[qb]);
- else
- result =
- internal::evaluate_tensor_product_value_and_gradient_linear(
- solution_renumbered.data(), unit_point_ptr[qb]);
- }
- else
- result =
- internal::evaluate_tensor_product_value_and_gradient_shapes<
- dim,
- scalar_value_type,
- VectorizedArrayType,
- 1,
- false>(shapes.data() + qb * n_shapes,
- n_shapes,
- solution_renumbered.data());
- gradient[0] = result[0];
- if (dim > 1)
- gradient[1] = result[1];
- if (dim > 2)
- gradient[2] = result[2];
- value = result[dim];
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ if (sum_into_values)
+ solution_values[i + offset] +=
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear + i) :
+ this->solution_renumbered_vectorized[i]);
+ else
+ solution_values[i + offset] =
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear + i) :
+ this->solution_renumbered_vectorized[i]);
}
else
{
- if constexpr (is_linear)
- {
- if constexpr (n_components == 1)
- value = internal::evaluate_tensor_product_value_linear<
- dim,
- scalar_value_type,
- VectorizedArrayType,
- stride_view>(solution_values.data(), unit_point_ptr[qb]);
- else
- value = internal::evaluate_tensor_product_value_linear(
- solution_renumbered.data(), unit_point_ptr[qb]);
- }
- else
- value = internal::evaluate_tensor_product_value_shapes<
- dim,
- scalar_value_type,
- VectorizedArrayType,
- false>(shapes.data() + qb * n_shapes,
- n_shapes,
- solution_renumbered.data());
+ const unsigned int *renumber_ptr = this->renumber.data() + offset;
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ if (sum_into_values)
+ solution_values[renumber_ptr[i]] +=
+ ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
+ else
+ solution_values[renumber_ptr[i]] =
+ ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
+ }
+ }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW, bool is_linear, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ // zero out lanes of incomplete last quadrature point batch
+ if constexpr (stride == 1)
+ if (const unsigned int n_filled_lanes =
+ this->n_q_points_scalar & (n_lanes_internal - 1);
+ n_filled_lanes > 0)
+ {
+ if (integration_flags & EvaluationFlags::values)
+ for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+ ETT::set_zero_value(this->values.back(), v);
+ if (integration_flags & EvaluationFlags::gradients)
+ for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+ ETT::set_zero_gradient(this->gradients.back(), v);
+ }
+
+ std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
+ solution_values_vectorized_linear = {};
+
+ // loop over quadrature batches qb
+ const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
+
+ const bool cartesian_cell =
+ this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian;
+ const bool affine_cell =
+ this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine;
+ for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
+ {
+ vectorized_value_type value = {};
+ Tensor<1, dim, vectorized_value_type> gradient;
+
+ if (integration_flags & EvaluationFlags::values)
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ ETT::get_value(value,
+ v,
+ do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
+ this->values[offset]);
+
+ if (integration_flags & EvaluationFlags::gradients)
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ {
+ const auto grad_w =
+ do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
+ this->gradients[offset];
+ ETT::get_gradient(
+ gradient,
+ v,
+ cartesian_cell ?
+ apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
+ grad_w) :
+ apply_transformation(
+ this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
+ grad_w));
+ }
+
+ compute_integrate_fast<is_linear>(
+ integration_flags,
+ n_shapes,
+ qb,
+ value,
+ gradient,
+ solution_values_vectorized_linear.data());
+ }
+
+ // add between the lanes and write into the result
+ finish_integrate_fast<is_linear>(solution_values,
+ solution_values_vectorized_linear.data(),
+ sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW, std::size_t stride_view>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ // slow path with FEValues
+ Assert(this->fe_values.get() != nullptr,
+ ExcMessage(
+ "Not initialized. Please call FEPointEvaluation::reinit()!"));
+ if (!sum_into_values)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
+
+ const std::size_t n_points = this->fe_values->get_quadrature().size();
+
+ if (integration_flags & EvaluationFlags::values)
+ {
+ AssertIndexRange(this->n_q_points, this->values.size() + 1);
+ for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
+ {
+ for (unsigned int d = 0; d < n_components; ++d)
+ if (this->nonzero_shape_function_component[i][d] &&
+ (this->fe->is_primitive(i) || this->fe->is_primitive()))
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ solution_values[i] +=
+ this->fe_values->shape_value(i, q + v) *
+ ETT::access(this->values[qb], v, d) *
+ (do_JxW ? this->fe_values->JxW(q + v) : 1.);
+ else if (this->nonzero_shape_function_component[i][d])
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ solution_values[i] +=
+ this->fe_values->shape_value_component(i, q + v, d) *
+ ETT::access(this->values[qb], v, d) *
+ (do_JxW ? this->fe_values->JxW(q + v) : 1.);
+ }
+ }
+
+ if (integration_flags & EvaluationFlags::gradients)
+ {
+ AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
+ for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
+ {
+ for (unsigned int d = 0; d < n_components; ++d)
+ if (this->nonzero_shape_function_component[i][d] &&
+ (this->fe->is_primitive(i) || this->fe->is_primitive()))
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ solution_values[i] +=
+ this->fe_values->shape_grad(i, q + v) *
+ ETT::access(this->gradients[qb], v, d) *
+ (do_JxW ? this->fe_values->JxW(q + v) : 1.);
+ else if (this->nonzero_shape_function_component[i][d])
+ for (unsigned int qb = 0, q = 0; q < n_points;
+ ++qb, q += n_lanes_user_interface)
+ for (unsigned int v = 0;
+ v < n_lanes_user_interface && q + v < n_points;
+ ++v)
+ solution_values[i] +=
+ this->fe_values->shape_grad_component(i, q + v, d) *
+ ETT::access(this->gradients[qb], v, d) *
+ (do_JxW ? this->fe_values->JxW(q + v) : 1.);
}
}
}
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face_path, bool is_linear, std::size_t stride_view>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags)
+template <bool do_JxW, std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ if (!this->is_reinitialized)
+ reinit();
+
+ Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
+
+ if (this->n_q_points == 0 || // no evaluation points provided
+ !((integration_flags & EvaluationFlags::values) ||
+ (integration_flags &
+ EvaluationFlags::gradients))) // no integration flags
+ {
+ if (!sum_into_values)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
+ return;
+ }
+
+ Assert(
+ !do_JxW || this->JxW_ptr != nullptr,
+ ExcMessage(
+ "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
+
+ AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
+ if (this->fast_path)
+ {
+ if (this->use_linear_path)
+ integrate_fast<do_JxW, true>(solution_values,
+ integration_flags,
+ sum_into_values);
+ else
+ integrate_fast<do_JxW, false>(solution_values,
+ integration_flags,
+ sum_into_values);
+ }
+ else
+ integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Tensor<1, spacedim, Number>
+FEPointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
+ const unsigned int point_index) const
{
- if (!(is_linear && n_components == 1) || is_face_path)
- prepare_evaluate_fast<is_face_path, is_linear>(solution_values,
- evaluation_flags);
+ AssertIndexRange(point_index, this->n_q_points);
+ Assert(this->normal_ptr != nullptr,
+ internal::FEPointEvaluation::
+ ExcFEPointEvaluationAccessToUninitializedMappingField(
+ "update_normal_vectors"));
+ if (this->is_interior)
+ return this->normal_ptr[point_index];
+ else
+ return -this->normal_ptr[point_index];
+}
- // loop over quadrature batches qb
- const unsigned int n_shapes = is_linear ? 2 : poly.size();
- for (unsigned int qb = 0; qb < n_q_batches; ++qb)
- {
- vectorized_value_type value;
- interface_vectorized_gradient_type gradient;
- compute_evaluate_fast<is_face_path, is_linear>(
- solution_values, evaluation_flags, n_shapes, qb, value, gradient);
+/**
+ * This class provides an interface to the evaluation of interpolated solution
+ * values and gradients on faces on arbitrary reference point positions. These
+ * points can change from face to face, both with respect to their quantity as
+ * well to the location. A typical use case is evaluations on non-matching
+ * grids.
+ *
+ * The use of this class is similar to FEEvaluation: In the constructor, a
+ * reference to a NonMatching::MappingInfo object is passed, where the
+ * quadrature points in reference position is stored together with the mapping
+ * information. The class is then reinitialized to a cell by calling
+ * `FEFacePointEvaluation::reinit(face_index)` or
+ * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to
+ * evaluate() or integrate(), the user can compute information at the given
+ * points. Eventually, the access functions get_value() or get_gradient() allow
+ * to query this information at a specific point index.
+ */
+template <int n_components_,
+ int dim,
+ int spacedim = dim,
+ typename Number = double>
+class FEFacePointEvaluation
+ : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+{
+public:
+ static constexpr unsigned int dimension = dim;
+ static constexpr unsigned int n_components = n_components_;
- if (evaluation_flags & EvaluationFlags::values)
- {
- for (unsigned int v = 0, offset = qb * stride;
- v < stride && (stride == 1 || offset < n_q_points_scalar);
- ++v, ++offset)
- ETT::set_value(value, v, values[offset]);
- }
- if (evaluation_flags & EvaluationFlags::gradients)
- {
- Assert(update_flags & update_gradients ||
- update_flags & update_inverse_jacobians,
- ExcNotInitialized());
+ using number_type = Number;
- for (unsigned int v = 0, offset = qb * stride;
- v < stride && (stride == 1 || offset < n_q_points_scalar);
- ++v, ++offset)
- {
- gradient_type unit_gradient;
- ETT::set_gradient(gradient, v, unit_gradient);
- gradients[offset] =
- cell_type <=
- internal::MatrixFreeFunctions::GeometryType::cartesian ?
- apply_diagonal_transformation(inverse_jacobian_ptr[0],
- unit_gradient) :
- apply_transformation(
- inverse_jacobian_ptr
- [cell_type <=
- internal::MatrixFreeFunctions::GeometryType::affine ?
- 0 :
- offset]
- .transpose(),
- unit_gradient);
- }
- }
- }
-}
+ using ScalarNumber =
+ typename internal::VectorizedArrayTrait<Number>::value_type;
+ using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+ Number>::vectorized_value_type;
+ using ETT = typename internal::FEPointEvaluation::
+ EvaluatorTypeTraits<dim, n_components, Number>;
+ using value_type = typename ETT::value_type;
+ using scalar_value_type = typename ETT::scalar_value_type;
+ using vectorized_value_type = typename ETT::vectorized_value_type;
+ using gradient_type = typename ETT::gradient_type;
+ using interface_vectorized_gradient_type =
+ typename ETT::interface_vectorized_gradient_type;
+
+ /**
+ * Constructor. Allows to select if interior or exterior face is selected.
+ */
+ FEFacePointEvaluation(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const bool is_interior = true,
+ const unsigned int first_selected_component = 0);
+
+ /**
+ * Reinitialize the evaluator to point to the correct precomputed mapping of
+ * the face in the MappingInfo object. Used in element-centric loops (ECL).
+ */
+ void
+ reinit(const unsigned int cell_index, const unsigned int face_number);
+
+ /**
+ * Reinitialize the evaluator to point to the correct precomputed mapping of
+ * the face in the MappingInfo object. Used in face-centric loops (FCL).
+ */
+ void
+ reinit(const unsigned int face_index);
+
+ /**
+ * This function interpolates the finite element solution, represented by
+ * `solution_values`, on the cell and `unit_points` passed to reinit().
+ *
+ * @param[in] solution_values This array is supposed to contain the unknown
+ * values on the element read out by
+ * `FEEvaluation::read_dof_values(global_vector)`.
+ *
+ * @param[in] evaluation_flags Flags specifying which quantities should be
+ * evaluated at the points.
+ */
+ template <std::size_t stride_view>
+ void
+ evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * This function interpolates the finite element solution, represented by
+ * `solution_values`, on the cell and `unit_points` passed to reinit().
+ *
+ * @param[in] solution_values This array is supposed to contain the unknown
+ * values on the element as returned by `cell->get_dof_values(global_vector,
+ * solution_values)`.
+ *
+ * @param[in] evaluation_flags Flags specifying which quantities should be
+ * evaluated at the points.
+ */
+ void
+ evaluate(const ArrayView<const ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used during
+ * `FEEvaluation::set_dof_values(global_vector)` or
+ * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ template <std::size_t stride_view>
+ void
+ integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used to during
+ * `cell->set_dof_values(solution_values, global_vector)` or
+ * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ void
+ integrate(const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used during
+ * `FEEvaluation::set_dof_values(global_vector)` or
+ * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ template <std::size_t stride_view>
+ void
+ test_and_sum(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * This function multiplies the quantities passed in by previous
+ * submit_value() or submit_gradient() calls by the value or gradient of the
+ * test functions, and performs summation over all given points multiplied be
+ * the Jacobian determinant times the quadrature weight (JxW).
+ *
+ * @param[out] solution_values This array will contain the result of the
+ * integral, which can be used to during
+ * `cell->set_dof_values(solution_values, global_vector)` or
+ * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+ * that for multi-component systems where only some of the components are
+ * selected by the present class, the entries in `solution_values` not touched
+ * by this class will be set to zero.
+ *
+ * @param[in] integration_flags Flags specifying which quantities should be
+ * integrated at the points.
+ *
+ * @param[in] sum_into_values Flag specifying if the integrated values
+ * should be summed into the solution values. Defaults to false.
+ *
+ */
+ void
+ test_and_sum(const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * Evaluate values and gradients in face for the selected face (lane) of the
+ * batch. Default stride into the face dofs is width of
+ * VectorizedArray<selected_floating_point_type> which is the default
+ * vectorization over faces for FEFaceEvaluation.
+ */
+ template <int stride_face_dof = VectorizedArrayType::size()>
+ void
+ evaluate_in_face(const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * Integrate values and gradients in face for the selected face (lane) of the
+ * batch. Default stride into the face dofs is width of
+ * VectorizedArray<selected_floating_point_type> which is the default
+ * vectorization over faces for FEFaceEvaluation.
+ */
+ template <int stride_face_dof = VectorizedArrayType::size()>
+ void
+ integrate_in_face(ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values = false);
+
+ /**
+ * Return the normal vector. This class or the MappingInfo object passed to
+ * this function needs to be constructed with UpdateFlags containing
+ * `update_normal_vectors`.
+ */
+ Tensor<1, spacedim, Number>
+ normal_vector(const unsigned int point_index) const;
+
+private:
+ static constexpr std::size_t n_lanes_user_interface =
+ internal::VectorizedArrayTrait<Number>::width();
+ static constexpr std::size_t n_lanes_internal =
+ internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+ static constexpr std::size_t stride =
+ internal::VectorizedArrayTrait<Number>::stride();
+
+ template <bool is_linear, std::size_t stride_view>
+ void
+ do_evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ template <bool do_JxW, bool is_linear, std::size_t stride_view>
+ void
+ do_integrate(
+ const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+
+ /**
+ * Actually does the evaluation templated on the chosen code path (linear or
+ * higher order).
+ */
+ template <bool is_linear, int stride_face_dof>
+ void
+ do_evaluate_in_face(const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags);
+
+ /**
+ * Actually does the integration templated on the chosen code path (linear or
+ * higher order).
+ */
+ template <bool do_JxW, bool is_linear, int stride_face_dof>
+ void
+ do_integrate_in_face(
+ ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values);
+};
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
+ FEFacePointEvaluation(
+ NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+ const FiniteElement<dim> &fe,
+ const bool is_interior,
+ const unsigned int first_selected_component)
+ : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+ mapping_info,
+ fe,
+ first_selected_component,
+ is_interior)
+{}
template <int n_components_, int dim, int spacedim, typename Number>
-template <std::size_t stride_view>
inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_slow(
- const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &evaluation_flags)
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+ const unsigned int cell_index,
+ const unsigned int face_number)
{
- // slow path with FEValues
- Assert(fe_values.get() != nullptr,
- ExcMessage(
- "Not initialized. Please call FEPointEvaluation::reinit()!"));
+ this->current_cell_index = cell_index;
+ this->current_face_number = face_number;
- const std::size_t n_points = fe_values->get_quadrature().size();
+ if (this->use_linear_path)
+ this->template do_reinit<true, true>();
+ else
+ this->template do_reinit<true, false>();
+}
- if (evaluation_flags & EvaluationFlags::values)
- {
- values.resize(n_q_points);
- std::fill(values.begin(), values.end(), value_type());
- for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
- {
- const ScalarNumber value = solution_values[i];
- for (unsigned int d = 0; d < n_components; ++d)
- if (nonzero_shape_function_component[i][d] &&
- (fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- ETT::access(values[qb],
- v,
- d,
- fe_values->shape_value(i, q + v) * value);
- else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- ETT::access(values[qb],
- v,
- d,
- fe_values->shape_value_component(i, q + v, d) *
- value);
- }
- }
- if (evaluation_flags & EvaluationFlags::gradients)
- {
- gradients.resize(n_q_points);
- std::fill(gradients.begin(), gradients.end(), gradient_type());
- for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
- {
- const ScalarNumber value = solution_values[i];
- for (unsigned int d = 0; d < n_components; ++d)
- if (nonzero_shape_function_component[i][d] &&
- (fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- ETT::access(gradients[qb],
- v,
- d,
- fe_values->shape_grad(i, q + v) * value);
- else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- ETT::access(gradients[qb],
- v,
- d,
- fe_values->shape_grad_component(i, q + v, d) *
- value);
- }
- }
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+ const unsigned int face_index)
+{
+ this->current_cell_index = face_index;
+ this->current_face_number =
+ this->mapping_info->get_face_number(face_index, this->is_interior);
+
+ if (this->use_linear_path)
+ this->template do_reinit<true, true>();
+ else
+ this->template do_reinit<true, false>();
}
template <int n_components_, int dim, int spacedim, typename Number>
template <std::size_t stride_view>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &evaluation_flags)
{
- if (!is_reinitialized)
- reinit();
+ Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!"));
- if (n_q_points == 0)
+ if (this->n_q_points == 0)
return;
Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
(evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
return;
- AssertDimension(solution_values.size(), fe->dofs_per_cell);
- if (fast_path)
- {
- if (use_face_path)
- {
- if (use_linear_path)
- evaluate_fast<true, true>(solution_values, evaluation_flags);
- else
- evaluate_fast<true, false>(solution_values, evaluation_flags);
- }
- else
- {
- if (use_linear_path)
- evaluate_fast<false, true>(solution_values, evaluation_flags);
- else
- evaluate_fast<false, false>(solution_values, evaluation_flags);
- }
- }
+ AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
+
+ if (this->use_linear_path)
+ do_evaluate<true>(solution_values, evaluation_flags);
else
- evaluate_slow(solution_values, evaluation_flags);
+ do_evaluate<false>(solution_values, evaluation_flags);
}
template <int n_components_, int dim, int spacedim, typename Number>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
const ArrayView<const ScalarNumber> &solution_values,
const EvaluationFlags::EvaluationFlags &evaluation_flags)
{
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face_path, bool is_linear>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_integrate_fast(
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const unsigned int n_shapes,
- const unsigned int qb,
- const vectorized_value_type value,
- const interface_vectorized_gradient_type gradient,
- vectorized_value_type *solution_values_vectorized_linear)
-{
- if (is_face_path)
- {
- if (integration_flags & EvaluationFlags::gradients)
- {
- std::array<vectorized_value_type, 2> value_face = {};
- Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
-
- value_face[0] = value;
- // fill derivative in physical coordinates into tangential/normal
- // derivatives
- if (current_face_number / 2 == 0)
- {
- value_face[1] = gradient[0];
- if (dim > 1)
- gradient_in_face[0] = gradient[1];
- if (dim > 2)
- gradient_in_face[1] = gradient[2];
- }
- else if (current_face_number / 2 == 1)
- {
- if (dim > 1)
- value_face[1] = gradient[1];
- if (dim == 3)
- {
- gradient_in_face[0] = gradient[2];
- gradient_in_face[1] = gradient[0];
- }
- else if (dim == 2)
- gradient_in_face[0] = gradient[0];
- else
- Assert(false, ExcInternalError());
- }
- else if (current_face_number / 2 == 2)
- {
- if (dim > 2)
- {
- value_face[1] = gradient[2];
- gradient_in_face[0] = gradient[0];
- gradient_in_face[1] = gradient[1];
- }
- else
- Assert(false, ExcInternalError());
- }
- else
- Assert(false, ExcInternalError());
-
- internal::integrate_tensor_product_value_and_gradient<
- is_linear,
- dim - 1,
- VectorizedArrayType,
- vectorized_value_type,
- 2>(shapes_faces.data() + qb * n_shapes,
- n_shapes,
- value_face.data(),
- gradient_in_face,
- is_linear ? solution_values_vectorized_linear :
- solution_renumbered_vectorized.data(),
- unit_point_faces_ptr[qb],
- qb != 0);
- }
- else
- internal::integrate_tensor_product_value<is_linear,
- dim - 1,
- VectorizedArrayType,
- vectorized_value_type>(
- shapes_faces.data() + qb * n_shapes,
- n_shapes,
- value,
- is_linear ? solution_values_vectorized_linear :
- solution_renumbered_vectorized.data(),
- unit_point_faces_ptr[qb],
- qb != 0);
- }
- else
- {
- if (integration_flags & EvaluationFlags::gradients)
- internal::integrate_tensor_product_value_and_gradient<
- is_linear,
- dim,
- VectorizedArrayType,
- vectorized_value_type>(shapes.data() + qb * n_shapes,
- n_shapes,
- &value,
- gradient,
- is_linear ?
- solution_values_vectorized_linear :
- solution_renumbered_vectorized.data(),
- unit_point_ptr[qb],
- qb != 0);
- else
- internal::integrate_tensor_product_value<is_linear,
- dim,
- VectorizedArrayType,
- vectorized_value_type>(
- shapes.data() + qb * n_shapes,
- n_shapes,
- value,
- is_linear ? solution_values_vectorized_linear :
- solution_renumbered_vectorized.data(),
- unit_point_ptr[qb],
- qb != 0);
- }
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face_path, bool is_linear, std::size_t stride_view>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- vectorized_value_type *solution_values_vectorized_linear,
- const bool sum_into_values)
+template <bool is_linear, std::size_t stride_view>
+void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::do_evaluate(
+ const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
{
- if (!sum_into_values && fe->n_components() > n_components)
- for (unsigned int i = 0; i < solution_values.size(); ++i)
- solution_values[i] = 0;
-
const unsigned int dofs_per_comp =
- is_linear ? Utilities::pow(2, dim) : dofs_per_component;
-
- for (unsigned int comp = 0; comp < n_components; ++comp)
- {
- const std::size_t offset =
- (component_in_base_element + comp) * dofs_per_comp;
-
- if (is_face_path)
- {
- const unsigned int dofs_per_comp_face =
- is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
-
- const unsigned int size_input = 2 * dofs_per_comp_face;
- ScalarNumber *input = scratch_data_scalar.begin();
- ScalarNumber *output = input + size_input;
-
- for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
- input[i] =
- ETT::sum_value(comp,
- is_linear ?
- *(solution_values_vectorized_linear + i) :
- solution_renumbered_vectorized[i]);
-
- internal::FEFaceNormalEvaluationImpl<dim, -1, ScalarNumber>::
- template interpolate<false, false>(1,
- integration_flags,
- shape_info,
- input,
- output,
- current_face_number);
-
- if (is_linear || renumber.empty())
- {
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- if (sum_into_values)
- solution_values[i + offset] += output[i];
- else
- solution_values[i + offset] = output[i];
- }
- else
- {
- const unsigned int *renumber_ptr = renumber.data() + offset;
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- if (sum_into_values)
- solution_values[renumber_ptr[i]] += output[i];
- else
- solution_values[renumber_ptr[i]] = output[i];
- }
- }
- else
- {
- if (is_linear || renumber.empty())
- {
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- if (sum_into_values)
- solution_values[i + offset] +=
- ETT::sum_value(comp,
- is_linear ?
- *(solution_values_vectorized_linear + i) :
- solution_renumbered_vectorized[i]);
- else
- solution_values[i + offset] =
- ETT::sum_value(comp,
- is_linear ?
- *(solution_values_vectorized_linear + i) :
- solution_renumbered_vectorized[i]);
- }
- else
- {
- const unsigned int *renumber_ptr = renumber.data() + offset;
- for (unsigned int i = 0; i < dofs_per_comp; ++i)
- if (sum_into_values)
- solution_values[renumber_ptr[i]] +=
- ETT::sum_value(comp, solution_renumbered_vectorized[i]);
- else
- solution_values[renumber_ptr[i]] =
- ETT::sum_value(comp, solution_renumbered_vectorized[i]);
- }
- }
- }
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW,
- bool is_face_path,
- bool is_linear,
- std::size_t stride_view>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values)
-{
- // zero out lanes of incomplete last quadrature point batch
- if constexpr (stride == 1)
- if (const unsigned int n_filled_lanes =
- n_q_points_scalar & (n_lanes_internal - 1);
- n_filled_lanes > 0)
- {
- if (integration_flags & EvaluationFlags::values)
- for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
- ETT::set_zero_value(values.back(), v);
- if (integration_flags & EvaluationFlags::gradients)
- for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
- ETT::set_zero_gradient(gradients.back(), v);
- }
-
- std::array<vectorized_value_type,
- is_linear ? (is_face_path ? 2 * Utilities::pow(2, dim - 1) :
- Utilities::pow(2, dim)) :
- 0>
- solution_values_vectorized_linear = {};
-
- // loop over quadrature batches qb
- const unsigned int n_shapes = is_linear ? 2 : poly.size();
+ is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
- const bool cartesian_cell =
- cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian;
- const bool affine_cell =
- cell_type <= internal::MatrixFreeFunctions::GeometryType::affine;
- for (unsigned int qb = 0; qb < n_q_batches; ++qb)
+ const ScalarNumber *input;
+ if (stride_view == 1 && this->component_in_base_element == 0 &&
+ (is_linear || this->renumber.empty()))
+ input = solution_values.data();
+ else
{
- vectorized_value_type value = {};
- Tensor<1, dim, vectorized_value_type> gradient;
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ {
+ const std::size_t offset =
+ (this->component_in_base_element + comp) * dofs_per_comp;
- if (integration_flags & EvaluationFlags::values)
- for (unsigned int v = 0, offset = qb * stride;
- v < stride && (stride == 1 || offset < n_q_points_scalar);
- ++v, ++offset)
- ETT::get_value(value,
- v,
- do_JxW ? values[offset] * JxW_ptr[offset] :
- values[offset]);
+ if (is_linear || this->renumber.empty())
+ {
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ this->scratch_data_scalar[i + comp * dofs_per_comp] =
+ solution_values[i + offset];
+ }
+ else
+ {
+ const unsigned int *renumber_ptr = this->renumber.data() + offset;
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ this->scratch_data_scalar[i + comp * dofs_per_comp] =
+ solution_values[renumber_ptr[i]];
+ }
+ }
+ input = this->scratch_data_scalar.data();
+ }
- if (integration_flags & EvaluationFlags::gradients)
- for (unsigned int v = 0, offset = qb * stride;
- v < stride && (stride == 1 || offset < n_q_points_scalar);
- ++v, ++offset)
- {
- const auto grad_w =
- do_JxW ? gradients[offset] * JxW_ptr[offset] : gradients[offset];
- ETT::get_gradient(
- gradient,
- v,
- cartesian_cell ?
- apply_diagonal_transformation(inverse_jacobian_ptr[0], grad_w) :
- apply_transformation(
- inverse_jacobian_ptr[affine_cell ? 0 : offset], grad_w));
- }
+ ScalarNumber *output =
+ this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
- compute_integrate_fast<is_face_path, is_linear>(
- integration_flags,
- n_shapes,
- qb,
- value,
- gradient,
- solution_values_vectorized_linear.data());
- }
+ internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
+ template interpolate<true, false>(n_components,
+ evaluation_flags,
+ this->shape_info,
+ input,
+ output,
+ this->current_face_number);
- // add between the lanes and write into the result
- finish_integrate_fast<is_face_path, is_linear>(
- solution_values,
- integration_flags,
- solution_values_vectorized_linear.data(),
- sum_into_values);
+ do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
}
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW, std::size_t stride_view>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
+template <std::size_t stride_view>
+void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::integrate(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
const bool sum_into_values)
{
- // slow path with FEValues
- Assert(fe_values.get() != nullptr,
- ExcMessage(
- "Not initialized. Please call FEPointEvaluation::reinit()!"));
- if (!sum_into_values)
- for (unsigned int i = 0; i < solution_values.size(); ++i)
- solution_values[i] = 0;
+ Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!"));
- const std::size_t n_points = fe_values->get_quadrature().size();
+ Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
- if (integration_flags & EvaluationFlags::values)
+ if (this->n_q_points == 0 || // no evaluation points provided
+ !((integration_flags & EvaluationFlags::values) ||
+ (integration_flags &
+ EvaluationFlags::gradients))) // no integration flags
{
- AssertIndexRange(n_q_points, values.size() + 1);
- for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
- {
- for (unsigned int d = 0; d < n_components; ++d)
- if (nonzero_shape_function_component[i][d] &&
- (fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- solution_values[i] += fe_values->shape_value(i, q + v) *
- ETT::access(values[qb], v, d) *
- (do_JxW ? fe_values->JxW(q + v) : 1.);
- else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- solution_values[i] +=
- fe_values->shape_value_component(i, q + v, d) *
- ETT::access(values[qb], v, d) *
- (do_JxW ? fe_values->JxW(q + v) : 1.);
- }
+ if (!sum_into_values)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
+ return;
}
- if (integration_flags & EvaluationFlags::gradients)
- {
- AssertIndexRange(n_q_points, gradients.size() + 1);
- for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
- {
- for (unsigned int d = 0; d < n_components; ++d)
- if (nonzero_shape_function_component[i][d] &&
- (fe->is_primitive(i) || fe->is_primitive()))
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- solution_values[i] += fe_values->shape_grad(i, q + v) *
- ETT::access(gradients[qb], v, d) *
- (do_JxW ? fe_values->JxW(q + v) : 1.);
- else if (nonzero_shape_function_component[i][d])
- for (unsigned int qb = 0, q = 0; q < n_points;
- ++qb, q += n_lanes_user_interface)
- for (unsigned int v = 0;
- v < n_lanes_user_interface && q + v < n_points;
- ++v)
- solution_values[i] +=
- fe_values->shape_grad_component(i, q + v, d) *
- ETT::access(gradients[qb], v, d) *
- (do_JxW ? fe_values->JxW(q + v) : 1.);
- }
- }
+ AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
+
+ if (this->use_linear_path)
+ do_integrate<true, true>(solution_values,
+ integration_flags,
+ sum_into_values);
+ else
+ do_integrate<true, false>(solution_values,
+ integration_flags,
+ sum_into_values);
}
template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW, std::size_t stride_view>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
const bool sum_into_values)
{
- if (!is_reinitialized)
- reinit();
+ Assert(this->is_reinitialized, ExcMessage("Is not reinitalized!"));
Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
- if (n_q_points == 0 || // no evaluation points provided
+ if (this->n_q_points == 0 || // no evaluation points provided
!((integration_flags & EvaluationFlags::values) ||
(integration_flags &
EvaluationFlags::gradients))) // no integration flags
return;
}
- Assert(
- !do_JxW || JxW_ptr != nullptr,
- ExcMessage(
- "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
+ AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
- AssertDimension(solution_values.size(), fe->dofs_per_cell);
- if (fast_path)
- {
- if (use_face_path)
- {
- if (use_linear_path)
- integrate_fast<do_JxW, true, true>(solution_values,
- integration_flags,
- sum_into_values);
- else
- integrate_fast<do_JxW, true, false>(solution_values,
- integration_flags,
- sum_into_values);
- }
- else
- {
- if (use_linear_path)
- integrate_fast<do_JxW, false, true>(solution_values,
- integration_flags,
- sum_into_values);
- else
- integrate_fast<do_JxW, false, false>(solution_values,
- integration_flags,
- sum_into_values);
- }
- }
+ if (this->use_linear_path)
+ do_integrate<false, true>(solution_values,
+ integration_flags,
+ sum_into_values);
else
- integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
+ do_integrate<false, false>(solution_values,
+ integration_flags,
+ sum_into_values);
}
template <int n_components_, int dim, int spacedim, typename Number>
-template <std::size_t stride_view>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
+ const ArrayView<ScalarNumber> &solution_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
+{
+ test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+ solution_values.size()),
+ integration_flags,
+ sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW, bool is_linear, std::size_t stride_view>
+void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
const StridedArrayView<ScalarNumber, stride_view> &solution_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
const bool sum_into_values)
{
- do_integrate<true>(solution_values, integration_flags, sum_into_values);
+ if (!sum_into_values && this->fe->n_components() > n_components)
+ for (unsigned int i = 0; i < solution_values.size(); ++i)
+ solution_values[i] = 0;
+
+ do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
+ integration_flags,
+ false);
+
+ ScalarNumber *input = this->scratch_data_scalar.begin();
+
+ if (stride_view == 1 && this->component_in_base_element == 0 &&
+ (is_linear || this->renumber.empty()))
+ {
+ if (sum_into_values)
+ internal::
+ FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
+ template interpolate<false, true>(n_components,
+ integration_flags,
+ this->shape_info,
+ input,
+ solution_values.data(),
+ this->current_face_number);
+ else
+ internal::
+ FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
+ template interpolate<false, false>(n_components,
+ integration_flags,
+ this->shape_info,
+ input,
+ solution_values.data(),
+ this->current_face_number);
+ }
+ else
+ {
+ const unsigned int dofs_per_comp_face =
+ is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
+
+ const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
+ ScalarNumber *output = input + size_input;
+
+ internal::
+ FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
+ template interpolate<false, false>(n_components,
+ integration_flags,
+ this->shape_info,
+ input,
+ output,
+ this->current_face_number);
+
+ const unsigned int dofs_per_comp =
+ is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
+
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ {
+ const std::size_t offset =
+ (this->component_in_base_element + comp) * dofs_per_comp;
+
+ if (is_linear || this->renumber.empty())
+ {
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ if (sum_into_values)
+ solution_values[i + offset] +=
+ output[i + comp * dofs_per_comp];
+ else
+ solution_values[i + offset] =
+ output[i + comp * dofs_per_comp];
+ }
+ else
+ {
+ const unsigned int *renumber_ptr = this->renumber.data() + offset;
+ for (unsigned int i = 0; i < dofs_per_comp; ++i)
+ if (sum_into_values)
+ solution_values[renumber_ptr[i]] +=
+ output[i + comp * dofs_per_comp];
+ else
+ solution_values[renumber_ptr[i]] =
+ output[i + comp * dofs_per_comp];
+ }
+ }
+ }
}
template <int n_components_, int dim, int spacedim, typename Number>
+template <int stride_face_dof>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
- const ArrayView<ScalarNumber> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values)
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::evaluate_in_face(
+ const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
{
- integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
- solution_values.size()),
- integration_flags,
- sum_into_values);
+ if (this->use_linear_path)
+ do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
+ evaluation_flags);
+ else
+ do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
+ evaluation_flags);
}
-template <int n_components_, int dim, int spacedim, typename Number>
-template <std::size_t stride_view>
-void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
- const StridedArrayView<ScalarNumber, stride_view> &solution_values,
- const EvaluationFlags::EvaluationFlags &integration_flags,
- const bool sum_into_values)
-{
- do_integrate<false>(solution_values, integration_flags, sum_into_values);
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_linear, int stride_face_dof>
+inline void
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
+ do_evaluate_in_face(const ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &evaluation_flags)
+{
+ const scalar_value_type *face_dof_values_ptr;
+ if constexpr (n_components == 1)
+ face_dof_values_ptr = face_dof_values;
+ else
+ {
+ const unsigned int dofs_per_comp_face =
+ is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
+ ETT::read_value(face_dof_values[i + comp * 3 * dofs_per_comp_face],
+ comp,
+ this->solution_renumbered[i]);
+
+ face_dof_values_ptr = this->solution_renumbered.data();
+ }
+
+ // loop over quadrature batches qb
+ const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
+
+ for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
+ {
+ vectorized_value_type value;
+ interface_vectorized_gradient_type gradient;
+
+ if (evaluation_flags & EvaluationFlags::gradients)
+ {
+ const std::array<vectorized_value_type, dim + 1> interpolated_value =
+ is_linear ?
+ internal::evaluate_tensor_product_value_and_gradient_linear<
+ dim - 1,
+ scalar_value_type,
+ VectorizedArrayType,
+ 2,
+ stride_face_dof>(face_dof_values_ptr,
+ this->unit_point_faces_ptr[qb]) :
+ internal::evaluate_tensor_product_value_and_gradient_shapes<
+ dim - 1,
+ scalar_value_type,
+ VectorizedArrayType,
+ 2,
+ false,
+ stride_face_dof>(this->shapes_faces.data() + qb * n_shapes,
+ n_shapes,
+ face_dof_values_ptr);
+
+ value = interpolated_value[dim - 1];
+ // reorder derivative from tangential/normal derivatives into tensor
+ // in physical coordinates
+ if (this->current_face_number / 2 == 0)
+ {
+ gradient[0] = interpolated_value[dim];
+ if (dim > 1)
+ gradient[1] = interpolated_value[0];
+ if (dim > 2)
+ gradient[2] = interpolated_value[1];
+ }
+ else if (this->current_face_number / 2 == 1)
+ {
+ if (dim > 1)
+ gradient[1] = interpolated_value[dim];
+ if (dim == 3)
+ {
+ gradient[0] = interpolated_value[1];
+ gradient[2] = interpolated_value[0];
+ }
+ else if (dim == 2)
+ gradient[0] = interpolated_value[0];
+ else
+ Assert(false, ExcInternalError());
+ }
+ else if (this->current_face_number / 2 == 2)
+ {
+ if (dim > 2)
+ {
+ gradient[0] = interpolated_value[0];
+ gradient[1] = interpolated_value[1];
+ gradient[2] = interpolated_value[dim];
+ }
+ else
+ Assert(false, ExcInternalError());
+ }
+ else
+ Assert(false, ExcInternalError());
+ }
+ else
+ {
+ value =
+ is_linear ?
+ internal::evaluate_tensor_product_value_linear<
+ dim - 1,
+ scalar_value_type,
+ VectorizedArrayType,
+ stride_face_dof>(face_dof_values_ptr,
+ this->unit_point_faces_ptr[qb]) :
+ internal::evaluate_tensor_product_value_shapes<
+ dim - 1,
+ scalar_value_type,
+ VectorizedArrayType,
+ false,
+ stride_face_dof>(this->shapes_faces.data() + qb * n_shapes,
+ n_shapes,
+ face_dof_values_ptr);
+ }
+
+ if (evaluation_flags & EvaluationFlags::values)
+ {
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ ETT::set_value(value, v, this->values[offset]);
+ }
+ if (evaluation_flags & EvaluationFlags::gradients)
+ {
+ Assert(this->update_flags & update_gradients ||
+ this->update_flags & update_inverse_jacobians,
+ ExcNotInitialized());
+
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ {
+ gradient_type unit_gradient;
+ ETT::set_gradient(gradient, v, unit_gradient);
+ this->gradients[offset] =
+ this->cell_type <=
+ internal::MatrixFreeFunctions::GeometryType::cartesian ?
+ apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
+ unit_gradient) :
+ apply_transformation(
+ this
+ ->inverse_jacobian_ptr[this->cell_type <=
+ internal::MatrixFreeFunctions::
+ GeometryType::affine ?
+ 0 :
+ offset]
+ .transpose(),
+ unit_gradient);
+ }
+ }
+ }
}
template <int n_components_, int dim, int spacedim, typename Number>
+template <int stride_face_dof>
void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
- const ArrayView<ScalarNumber> &solution_values,
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::integrate_in_face(
+ ScalarNumber *face_dof_values,
const EvaluationFlags::EvaluationFlags &integration_flags,
const bool sum_into_values)
{
- test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
- solution_values.size()),
- integration_flags,
- sum_into_values);
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluation<n_components_, dim, spacedim, Number>::
- value_type &
- FEPointEvaluation<n_components_, dim, spacedim, Number>::get_value(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, values.size());
- return values[point_index];
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluation<n_components_, dim, spacedim, Number>::
- gradient_type &
- FEPointEvaluation<n_components_, dim, spacedim, Number>::get_gradient(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, gradients.size());
- return gradients[point_index];
+ if (this->use_linear_path)
+ do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
+ integration_flags,
+ sum_into_values);
+ else
+ do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
+ integration_flags,
+ sum_into_values);
}
template <int n_components_, int dim, int spacedim, typename Number>
+template <bool do_JxW, bool is_linear, int stride_face_dof>
inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::submit_value(
- const value_type &value,
- const unsigned int point_index)
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
+ do_integrate_in_face(
+ ScalarNumber *face_dof_values,
+ const EvaluationFlags::EvaluationFlags &integration_flags,
+ const bool sum_into_values)
{
- AssertIndexRange(point_index, n_q_points);
- values[point_index] = value;
-}
-
-
+ // zero out lanes of incomplete last quadrature point batch
+ if constexpr (stride == 1)
+ if (const unsigned int n_filled_lanes =
+ this->n_q_points_scalar & (n_lanes_internal - 1);
+ n_filled_lanes > 0)
+ {
+ if (integration_flags & EvaluationFlags::values)
+ for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+ ETT::set_zero_value(this->values.back(), v);
+ if (integration_flags & EvaluationFlags::gradients)
+ for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+ ETT::set_zero_gradient(this->gradients.back(), v);
+ }
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::submit_gradient(
- const gradient_type &gradient,
- const unsigned int point_index)
-{
- AssertIndexRange(point_index, n_q_points);
- gradients[point_index] = gradient;
-}
+ std::array<vectorized_value_type,
+ is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
+ solution_values_vectorized_linear = {};
+ // loop over quadrature batches qb
+ const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
+ const bool cartesian_cell =
+ this->cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian;
+ const bool affine_cell =
+ this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine;
+ for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
+ {
+ vectorized_value_type value = {};
+ Tensor<1, dim, vectorized_value_type> gradient;
-template <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, dim, spacedim, Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::jacobian(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, n_q_points);
- Assert(jacobian_ptr != nullptr,
- internal::FEPointEvaluation::
- ExcFEPointEvaluationAccessToUninitializedMappingField(
- "update_jacobians"));
- return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions::
- GeometryType::affine ?
- 0 :
- point_index];
-}
+ if (integration_flags & EvaluationFlags::values)
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ ETT::get_value(value,
+ v,
+ do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
+ this->values[offset]);
+ if (integration_flags & EvaluationFlags::gradients)
+ for (unsigned int v = 0, offset = qb * stride;
+ v < stride && (stride == 1 || offset < this->n_q_points_scalar);
+ ++v, ++offset)
+ {
+ const auto grad_w =
+ do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
+ this->gradients[offset];
+ ETT::get_gradient(
+ gradient,
+ v,
+ cartesian_cell ?
+ apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
+ grad_w) :
+ apply_transformation(
+ this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
+ grad_w));
+ }
+ if (integration_flags & EvaluationFlags::gradients)
+ {
+ std::array<vectorized_value_type, 2> value_face = {};
+ Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
-template <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, spacedim, dim, Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::inverse_jacobian(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, n_q_points);
- Assert(inverse_jacobian_ptr != nullptr,
- internal::FEPointEvaluation::
- ExcFEPointEvaluationAccessToUninitializedMappingField(
- "update_inverse_jacobians"));
- return inverse_jacobian_ptr
- [cell_type <=
- ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ?
- 0 :
- point_index];
-}
+ value_face[0] = value;
+ // fill derivative in physical coordinates into tangential/normal
+ // derivatives
+ if (this->current_face_number / 2 == 0)
+ {
+ value_face[1] = gradient[0];
+ if (dim > 1)
+ gradient_in_face[0] = gradient[1];
+ if (dim > 2)
+ gradient_in_face[1] = gradient[2];
+ }
+ else if (this->current_face_number / 2 == 1)
+ {
+ if (dim > 1)
+ value_face[1] = gradient[1];
+ if (dim == 3)
+ {
+ gradient_in_face[0] = gradient[2];
+ gradient_in_face[1] = gradient[0];
+ }
+ else if (dim == 2)
+ gradient_in_face[0] = gradient[0];
+ else
+ Assert(false, ExcInternalError());
+ }
+ else if (this->current_face_number / 2 == 2)
+ {
+ if (dim > 2)
+ {
+ value_face[1] = gradient[2];
+ gradient_in_face[0] = gradient[0];
+ gradient_in_face[1] = gradient[1];
+ }
+ else
+ Assert(false, ExcInternalError());
+ }
+ else
+ Assert(false, ExcInternalError());
+ internal::integrate_tensor_product_value_and_gradient<
+ is_linear,
+ dim - 1,
+ VectorizedArrayType,
+ vectorized_value_type,
+ 2>(this->shapes_faces.data() + qb * n_shapes,
+ n_shapes,
+ value_face.data(),
+ gradient_in_face,
+ is_linear ? solution_values_vectorized_linear.data() :
+ this->solution_renumbered_vectorized.data(),
+ this->unit_point_faces_ptr[qb],
+ qb != 0);
+ }
+ else
+ internal::integrate_tensor_product_value<is_linear,
+ dim - 1,
+ VectorizedArrayType,
+ vectorized_value_type>(
+ this->shapes_faces.data() + qb * n_shapes,
+ n_shapes,
+ value,
+ is_linear ? solution_values_vectorized_linear.data() :
+ this->solution_renumbered_vectorized.data(),
+ this->unit_point_faces_ptr[qb],
+ qb != 0);
+ }
+ const unsigned int dofs_per_comp_face =
+ is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
-template <int n_components_, int dim, int spacedim, typename Number>
-inline Number
-FEPointEvaluation<n_components_, dim, spacedim, Number>::JxW(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, n_q_points);
- Assert(JxW_ptr != nullptr,
- internal::FEPointEvaluation::
- ExcFEPointEvaluationAccessToUninitializedMappingField(
- "update_JxW_values"));
- return JxW_ptr[point_index];
+ for (unsigned int comp = 0; comp < n_components; ++comp)
+ for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
+ if (sum_into_values)
+ face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
+ stride_face_dof] +=
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear.data() + i) :
+ this->solution_renumbered_vectorized[i]);
+ else
+ face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
+ ETT::sum_value(comp,
+ is_linear ?
+ *(solution_values_vectorized_linear.data() + i) :
+ this->solution_renumbered_vectorized[i]);
}
template <int n_components_, int dim, int spacedim, typename Number>
inline Tensor<1, spacedim, Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
+FEFacePointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
const unsigned int point_index) const
{
- AssertIndexRange(point_index, n_q_points);
- Assert(normal_ptr != nullptr,
+ AssertIndexRange(point_index, this->n_q_points);
+ Assert(this->normal_ptr != nullptr,
internal::FEPointEvaluation::
ExcFEPointEvaluationAccessToUninitializedMappingField(
"update_normal_vectors"));
- return normal_ptr[point_index];
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline Point<spacedim, Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::real_point(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, n_q_points);
- Assert(real_point_ptr != nullptr,
- internal::FEPointEvaluation::
- ExcFEPointEvaluationAccessToUninitializedMappingField(
- "update_quadrature_points"));
- return real_point_ptr[point_index];
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline Point<dim, Number>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::unit_point(
- const unsigned int point_index) const
-{
- AssertIndexRange(point_index, n_q_points);
- Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
- Point<dim, Number> unit_point;
- for (unsigned int d = 0; d < dim; ++d)
- unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
- unit_point_ptr[point_index / stride][d], point_index % stride);
- return unit_point;
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
-FEPointEvaluation<n_components_, dim, spacedim, Number>::
- quadrature_point_indices() const
-{
- return {0U, n_q_points};
+ if (this->cell_type <= dealii::internal::MatrixFreeFunctions::affine)
+ {
+ Tensor<1, spacedim, Number> normal;
+ for (unsigned int d = 0; d < dim; ++d)
+ normal[d] =
+ internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
+ 0);
+ if (this->is_interior)
+ return normal;
+ else
+ return -normal;
+ }
+ else
+ {
+ if (this->is_interior)
+ return this->normal_ptr[point_index];
+ else
+ return -(this->normal_ptr[point_index]);
+ }
}
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff