17#ifndef dealii_matrix_free_fe_evaluation_h
18#define dealii_matrix_free_fe_evaluation_h
55 <<
"You are requesting information from an FEEvaluation/FEFaceEvaluation "
56 <<
"object for which this kind of information has not been computed. What "
57 <<
"information these objects compute is determined by the update_* flags you "
58 <<
"pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, "
59 <<
"the operation you are attempting requires the <" << arg1
60 <<
"> flag to be set, but it was apparently not specified "
61 <<
"upon initialization.");
66 int n_q_points_1d = fe_degree + 1,
67 int n_components_ = 1,
68 typename Number = double,
100template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
104 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
105 "Type of Number and of VectorizedArrayType do not match.");
152 JxW(
const unsigned int q_point)
const;
200 const VectorizedArrayType & value)
const;
206 template <
typename T>
207 std::array<
T, VectorizedArrayType::size()>
215 template <
typename T>
219 const std::array<
T, VectorizedArrayType::size()> & value)
const;
225 std::array<
unsigned int, VectorizedArrayType::size()>
232 std::array<
unsigned int, VectorizedArrayType::size()>
241 const std::vector<unsigned int> &
294 const unsigned int dof_no,
295 const unsigned int first_selected_component,
297 const unsigned int fe_degree,
298 const unsigned int n_q_points,
302 const unsigned int face_type);
313 const unsigned int first_selected_component,
370 (is_face ? dim - 1 : dim),
392 (is_face ? dim - 1 : dim),
490 std::shared_ptr<internal::MatrixFreeFunctions::
491 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
496 template <
int,
int,
int,
int,
typename,
typename>
542 bool is_face =
false,
591 template <
typename VectorType>
623 template <
typename VectorType>
626 const unsigned int first_index = 0);
659 template <
typename VectorType>
663 const unsigned int first_index = 0,
664 const std::bitset<VectorizedArrayType::size()> &mask =
665 std::bitset<VectorizedArrayType::size()>().flip())
const;
705 template <
typename VectorType>
708 const unsigned int first_index = 0,
709 const std::bitset<VectorizedArrayType::size()> &mask =
710 std::bitset<VectorizedArrayType::size()>().flip())
const;
715 template <
typename VectorType>
719 const unsigned int first_index = 0,
720 const std::bitset<VectorizedArrayType::size()> &mask =
721 std::bitset<VectorizedArrayType::size()>().flip())
const;
856 const unsigned int q_point);
929 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
949 const unsigned int q_point);
970 const unsigned int q_point);
986 const unsigned int q_point);
1023 const VectorizedArrayType *
1034 VectorizedArrayType *
1047 const VectorizedArrayType *
1060 VectorizedArrayType *
1074 const VectorizedArrayType *
1088 VectorizedArrayType *
1103 const VectorizedArrayType *
1118 VectorizedArrayType *
1140 const unsigned int dof_no,
1143 const unsigned int fe_degree,
1144 const unsigned int n_q_points,
1148 const unsigned int face_type);
1217 template <
typename VectorType,
typename VectorOperation>
1221 const std::array<VectorType *, n_components_> &vectors,
1224 n_components_> & vectors_sm,
1225 const std::bitset<VectorizedArrayType::size()> &mask,
1226 const bool apply_constraints =
true)
const;
1235 template <
typename VectorType,
typename VectorOperation>
1239 const std::array<VectorType *, n_components_> &vectors,
1242 n_components_> & vectors_sm,
1243 const std::bitset<VectorizedArrayType::size()> &mask)
const;
1252 template <
typename VectorType,
typename VectorOperation>
1256 const std::array<VectorType *, n_components_> &vectors)
const;
1400 VectorizedArrayType>
1403 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1404 "Type of Number and of VectorizedArrayType do not match.");
1426 const unsigned int dof_no,
1429 const unsigned int fe_degree,
1430 const unsigned int n_q_points,
1471template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
1476 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1477 "Type of Number and of VectorizedArrayType do not match.");
1516 const unsigned int q_point);
1541 const unsigned int q_point);
1577 const unsigned int dof_no,
1580 const unsigned int fe_degree,
1581 const unsigned int n_q_points,
1623template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
1628 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1629 "Type of Number and of VectorizedArrayType do not match.");
1666 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1698 const unsigned int q_point);
1710 const unsigned int q_point);
1723 const unsigned int q_point);
1731 const unsigned int q_point);
1743 const unsigned int dof_no,
1746 const unsigned int dofs_per_cell,
1747 const unsigned int n_q_points,
1787template <
typename Number,
bool is_face,
typename VectorizedArrayType>
1792 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1793 "Type of Number and of VectorizedArrayType do not match.");
1868 const unsigned int q_point);
1875 const unsigned int q_point);
1882 const unsigned int q_point);
1918 const unsigned int dof_no,
1921 const unsigned int fe_degree,
1922 const unsigned int n_q_points,
2513 typename VectorizedArrayType>
2518 VectorizedArrayType>
2521 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2522 "Type of Number and of VectorizedArrayType do not match.");
2636 const unsigned int dof_no = 0,
2637 const unsigned int quad_no = 0,
2650 const std::pair<unsigned int, unsigned int> & range,
2651 const unsigned int dof_no = 0,
2652 const unsigned int quad_no = 0,
2752 template <
bool level_dof_access>
2774 const unsigned int give_n_q_points_1d);
2794 const bool evaluate_gradients,
2795 const bool evaluate_hessians =
false);
2819 const bool evaluate_values,
2820 const bool evaluate_gradients,
2821 const bool evaluate_hessians =
false);
2836 template <
typename VectorType>
2844 template <
typename VectorType>
2847 const bool evaluate_values,
2848 const bool evaluate_gradients,
2849 const bool evaluate_hessians =
false);
2869 integrate(
const bool integrate_values,
const bool integrate_gradients);
2884 VectorizedArrayType * values_array);
2891 const bool integrate_gradients,
2892 VectorizedArrayType *values_array);
2907 template <
typename VectorType>
2910 VectorType & output_vector);
2915 template <
typename VectorType>
2918 const bool integrate_gradients,
2919 VectorType &output_vector);
3002 int n_q_points_1d = fe_degree + 1,
3003 int n_components_ = 1,
3004 typename Number = double,
3010 VectorizedArrayType>
3013 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
3014 "Type of Number and of VectorizedArrayType do not match.");
3143 const unsigned int dof_no = 0,
3144 const unsigned int quad_no = 0,
3159 const std::pair<unsigned int, unsigned int> & range,
3161 const unsigned int dof_no = 0,
3162 const unsigned int quad_no = 0,
3176 reinit(
const unsigned int face_batch_number);
3186 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
3193 const unsigned int give_n_q_points_1d);
3212 evaluate(
const bool evaluate_values,
const bool evaluate_gradients);
3235 const bool evaluate_values,
3236 const bool evaluate_gradients);
3249 template <
typename VectorType>
3257 template <
typename VectorType>
3260 const bool evaluate_values,
3261 const bool evaluate_gradients);
3279 integrate(
const bool integrate_values,
const bool integrate_gradients);
3291 VectorizedArrayType * values_array);
3298 const bool integrate_gradients,
3299 VectorizedArrayType *values_array);
3312 template <
typename VectorType>
3315 VectorType & output_vector);
3320 template <
typename VectorType>
3323 const bool integrate_gradients,
3324 VectorType &output_vector);
3363 std::array<
unsigned int, VectorizedArrayType::size()>
3369 std::array<
unsigned int, VectorizedArrayType::size()>
3377 namespace MatrixFreeFunctions
3381 template <
int dim,
int degree>
3390 template <
int degree>
3393 static constexpr unsigned int value = degree + 1;
3406template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3410 const unsigned int dof_no,
3411 const unsigned int first_selected_component,
3412 const unsigned int quad_no_in,
3413 const unsigned int fe_degree,
3414 const unsigned int n_q_points,
3415 const bool is_interior_face,
3416 const unsigned int active_fe_index_in,
3417 const unsigned int active_quad_index_in,
3418 const unsigned int face_type)
3419 : scratch_data_array(data_in.acquire_scratch_data())
3420 , quad_no(quad_no_in)
3421 , matrix_info(&data_in)
3422 , dof_info(&data_in.get_dof_info(dof_no))
3425 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3426 data_in.get_mapping_info(),
3429 data_in.get_dof_info(dof_no).fe_index_from_degree(
3430 first_selected_component,
3433 active_fe_index_in :
3435 , active_quad_index(
3437 (mapping_data->quad_index_from_n_q_points(n_q_points)) :
3439 active_quad_index_in :
3441 mapping_data->descriptor.size() - 1)))
3443 &mapping_data->descriptor
3445 (active_quad_index *
std::
max<unsigned
int>(1, dim - 1) +
3448 , n_quadrature_points(descriptor->n_q_points)
3449 , data(&data_in.get_shape_info(
3452 dof_info->component_to_base_index[first_selected_component],
3457 , normal_vectors(nullptr)
3458 , normal_x_jacobian(nullptr)
3459 , quadrature_weights(descriptor->quadrature_weights.
begin())
3461 , is_interior_face(is_interior_face)
3465 internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3466 internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3467 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3472 VectorizedArrayType::size());
3478template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3485 const unsigned int first_selected_component,
3488 : scratch_data_array(new
AlignedVector<VectorizedArrayType>())
3492 , descriptor(nullptr)
3493 , n_quadrature_points(
3495 , matrix_info(nullptr)
3497 , mapping_data(nullptr)
3500 data(new
internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
3503 fe.component_to_base_index(first_selected_component).
first))
3506 , normal_vectors(nullptr)
3507 , normal_x_jacobian(nullptr)
3508 , quadrature_weights(nullptr)
3511 , is_interior_face(true)
3512 , dof_access_index(
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3516 if (other !=
nullptr &&
3521 std::make_shared<internal::MatrixFreeFunctions::
3522 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3523 mapping, quadrature, update_flags);
3526 mapping_data = &mapped_geometry->get_data_storage();
3527 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3528 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3533template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3538 : scratch_data_array(other.matrix_info == nullptr ?
3540 other.matrix_info->acquire_scratch_data())
3541 , quad_no(other.quad_no)
3542 , active_fe_index(other.active_fe_index)
3543 , active_quad_index(other.active_quad_index)
3544 , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
3545 , n_quadrature_points(other.n_quadrature_points)
3546 , matrix_info(other.matrix_info)
3547 , dof_info(other.dof_info)
3548 , mapping_data(other.mapping_data)
3549 , data(other.matrix_info == nullptr ?
3550 new
internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
3555 , normal_vectors(nullptr)
3556 , normal_x_jacobian(nullptr)
3557 , quadrature_weights(other.descriptor == nullptr ?
3559 descriptor->quadrature_weights.
begin())
3562 , is_interior_face(other.is_interior_face)
3563 , dof_access_index(other.dof_access_index)
3568 mapped_geometry = std::make_shared<
3569 internal::MatrixFreeFunctions::
3570 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3574 mapping_data = &mapped_geometry->get_data_storage();
3577 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3578 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3584template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3594 if (matrix_info ==
nullptr)
3597 delete scratch_data_array;
3601 matrix_info->release_scratch_data(scratch_data_array);
3617 scratch_data_array = matrix_info->acquire_scratch_data();
3620 quadrature_weights =
3621 (descriptor !=
nullptr ? descriptor->quadrature_weights.begin() :
nullptr);
3630 mapped_geometry = std::make_shared<
3631 internal::MatrixFreeFunctions::
3632 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3637 mapping_data = &mapped_geometry->get_data_storage();
3638 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3639 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3647template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3651 if (matrix_info !=
nullptr)
3655 matrix_info->release_scratch_data(scratch_data_array);
3662 delete scratch_data_array;
3666 scratch_data_array =
nullptr;
3671template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3676 if (matrix_info ==
nullptr)
3681 return this->mapping_data->data_index_offsets[cell];
3687template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3698template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3709template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3714 Assert(dof_info !=
nullptr,
3716 "FEEvaluation was not initialized with a MatrixFree object!"));
3722template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3728 Assert(normal_vectors !=
nullptr,
3730 "update_normal_vectors"));
3732 return normal_vectors[0];
3734 return normal_vectors[q_point];
3739template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3742 const unsigned int q_point)
const
3745 Assert(J_value !=
nullptr,
3747 "update_values|update_gradients"));
3751 return J_value[0] * this->quadrature_weights[q_point];
3754 return J_value[q_point];
3759template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3765 Assert(this->jacobian !=
nullptr,
3767 "update_gradients"));
3771 return jacobian[q_point];
3776template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3777inline std::array<
unsigned int, VectorizedArrayType::size()>
3783 const unsigned int n_lanes = VectorizedArrayType::size();
3784 std::array<unsigned int, n_lanes> cells;
3787 for (
unsigned int i = 0; i < n_lanes; ++i)
3790 if ((is_face ==
false) ||
3792 this->dof_access_index ==
3794 this->is_interior_face))
3797 for (
unsigned int i = 0; i < n_lanes; ++i)
3798 cells[i] = cell * n_lanes + i;
3801 this->dof_access_index ==
3803 this->is_interior_face ==
false)
3809 for (
unsigned int i = 0; i < n_lanes; i++)
3812 const unsigned int cell_this = this->cell * n_lanes + i;
3814 unsigned int face_index =
3815 this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3823 auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3824 .cells_interior[face_index % n_lanes];
3825 auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3826 .cells_exterior[face_index % n_lanes];
3829 if (cell_m == cell_this)
3831 else if (cell_p == cell_this)
3838 const unsigned int *cells_ =
3840 &this->matrix_info->get_face_info(cell).cells_interior[0] :
3841 &this->matrix_info->get_face_info(cell).cells_exterior[0];
3842 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3844 cells[i] = cells_[i];
3856 typename VectorizedArrayType,
3857 typename VectorizedArrayType2,
3858 typename GlobalVectorType,
3864 GlobalVectorType & array,
3865 VectorizedArrayType2 & out,
3877 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3880 array[cells[i] / VectorizedArrayType::size()]
3881 [cells[i] % VectorizedArrayType::size()]);
3887template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3888std::array<
unsigned int, VectorizedArrayType::size()>
3892 const unsigned int v_len = VectorizedArrayType::size();
3893 std::array<
unsigned int, VectorizedArrayType::size()> cells;
3896 for (
unsigned int i = 0; i < v_len; ++i)
3900 this->dof_access_index ==
3902 this->is_interior_face ==
false)
3905 for (
unsigned int i = 0; i < v_len; i++)
3908 const unsigned int cell_this = this->cell * v_len + i;
3920 .cells_interior[fn % v_len];
3922 .cells_exterior[fn % v_len];
3925 if (cell_m == cell_this)
3927 else if (cell_p == cell_this)
3933 for (
unsigned int i = 0; i < v_len; ++i)
3934 cells[i] = cell * v_len + i;
3942template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3943inline VectorizedArrayType
3947 VectorizedArrayType out = Number(1.);
3948 internal::process_cell_data(
3949 *
this, this->matrix_info, array, out, [](
auto &local,
const auto &global) {
3957template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3961 const VectorizedArrayType & in)
const
3963 internal::process_cell_data(
3964 *
this, this->matrix_info, array, in, [](
const auto &local,
auto &global) {
3971template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3972template <
typename T>
3973inline std::array<
T, VectorizedArrayType::size()>
3975 const AlignedVector<std::array<
T, VectorizedArrayType::size()>> &array)
const
3977 std::array<
T, VectorizedArrayType::size()> out;
3978 internal::process_cell_data(
3979 *
this, this->matrix_info, array, out, [](
auto &local,
const auto &global) {
3987template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3988template <
typename T>
3991 AlignedVector<std::array<
T, VectorizedArrayType::size()>> &array,
3992 const std::array<
T, VectorizedArrayType::size()> & in)
const
3994 internal::process_cell_data(
3995 *
this, this->matrix_info, array, in, [](
const auto &local,
auto &global) {
4002template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4003inline const std::vector<unsigned int> &
4007 return data->lexicographic_numbering;
4012template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4018 const_cast<VectorizedArrayType *
>(scratch_data),
4019 scratch_data_array->end() - scratch_data);
4024template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4029 return this->quad_no;
4034template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4039 if (is_face && this->dof_access_index ==
4048template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4053 return active_fe_index;
4058template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4063 return active_quad_index;
4068template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4073 Assert(matrix_info !=
nullptr,
4075 "FEEvaluation was not initialized with a MatrixFree object!"));
4076 return *matrix_info;
4086 typename VectorizedArrayType>
4091 VectorizedArrayType>::
4093 const unsigned int dof_no,
4094 const unsigned int first_selected_component,
4095 const unsigned int quad_no_in,
4096 const unsigned int fe_degree,
4097 const unsigned int n_q_points,
4098 const bool is_interior_face,
4099 const unsigned int active_fe_index,
4100 const unsigned int active_quad_index,
4101 const unsigned int face_type)
4105 first_selected_component,
4113 , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
4114 , dof_values_initialized(false)
4115 , values_quad_initialized(false)
4116 , gradients_quad_initialized(false)
4117 , hessians_quad_initialized(false)
4118 , values_quad_submitted(false)
4119 , gradients_quad_submitted(false)
4120 , first_selected_component(first_selected_component)
4122 set_data_pointers();
4124 this->dof_info->start_components.back() == 1 ||
4125 static_cast<int>(n_components_) <=
4127 this->dof_info->start_components
4128 [this->dof_info->component_to_base_index[first_selected_component] +
4130 first_selected_component,
4132 "You tried to construct a vector-valued evaluator with " +
4134 " components. However, "
4135 "the current base element has only " +
4137 this->dof_info->start_components
4138 [this->dof_info->component_to_base_index[first_selected_component] +
4140 first_selected_component) +
4141 " components left when starting from local element index " +
4143 first_selected_component -
4144 this->dof_info->start_components
4145 [this->dof_info->component_to_base_index[first_selected_component]]) +
4146 " (global index " +
std::to_string(first_selected_component) +
")"));
4158 typename VectorizedArrayType>
4163 VectorizedArrayType>::
4169 const unsigned int first_selected_component,
4177 first_selected_component,
4179 , n_fe_components(n_components_)
4180 , dof_values_initialized(false)
4181 , values_quad_initialized(false)
4182 , gradients_quad_initialized(false)
4183 , hessians_quad_initialized(false)
4184 , values_quad_submitted(false)
4185 , gradients_quad_submitted(false)
4188 , first_selected_component(first_selected_component)
4190 set_data_pointers();
4192 const unsigned int base_element_number =
4196 first_selected_component >=
4198 ExcMessage(
"The underlying element must at least contain as many "
4199 "components as requested by this class"));
4200 (void)base_element_number;
4209 typename VectorizedArrayType>
4214 VectorizedArrayType>::
4219 VectorizedArrayType> &other)
4221 , n_fe_components(other.n_fe_components)
4222 , dof_values_initialized(false)
4223 , values_quad_initialized(false)
4224 , gradients_quad_initialized(false)
4225 , hessians_quad_initialized(false)
4226 , values_quad_submitted(false)
4227 , gradients_quad_submitted(false)
4228 , first_selected_component(other.first_selected_component)
4230 set_data_pointers();
4239 typename VectorizedArrayType>
4244 VectorizedArrayType> &
4250 VectorizedArrayType> &other)
4255 AssertDimension(first_selected_component, other.first_selected_component);
4266 typename VectorizedArrayType>
4273 const unsigned int tensor_dofs_per_component =
4274 Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4275 const unsigned int dofs_per_component =
4276 this->data->dofs_per_component_on_cell;
4277 const unsigned int n_quadrature_points = this->n_quadrature_points;
4279 const unsigned int shift =
4280 std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4282 2 * n_quadrature_points;
4283 const unsigned int allocated_size =
4284 shift + n_components_ * dofs_per_component +
4285 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4286 this->scratch_data_array->resize_fast(allocated_size);
4289 for (
unsigned int c = 0; c < n_components_; ++c)
4292 this->scratch_data_array->begin() + c * dofs_per_component;
4295 this->scratch_data_array->begin() + n_components * dofs_per_component;
4296 gradients_quad = this->scratch_data_array->begin() +
4297 n_components * (dofs_per_component + n_quadrature_points);
4299 this->scratch_data_array->begin() +
4300 n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4301 this->scratch_data =
4302 this->scratch_data_array->begin() + n_components_ * dofs_per_component +
4303 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4313 template <
typename VectorType,
bool>
4314 struct BlockVectorSelector
4317 template <
typename VectorType>
4318 struct BlockVectorSelector<VectorType, true>
4320 using BaseVectorType =
typename VectorType::BlockType;
4322 static BaseVectorType *
4323 get_vector_component(VectorType &vec,
const unsigned int component)
4326 return &vec.block(component);
4330 template <
typename VectorType>
4331 struct BlockVectorSelector<VectorType, false>
4333 using BaseVectorType = VectorType;
4335 static BaseVectorType *
4336 get_vector_component(VectorType &vec,
const unsigned int component)
4355 template <
typename VectorType>
4356 struct BlockVectorSelector<
std::vector<VectorType>, false>
4358 using BaseVectorType = VectorType;
4360 static BaseVectorType *
4361 get_vector_component(std::vector<VectorType> &vec,
4362 const unsigned int component)
4365 return &vec[component];
4369 template <
typename VectorType>
4370 struct BlockVectorSelector<
std::vector<VectorType *>, false>
4372 using BaseVectorType = VectorType;
4374 static BaseVectorType *
4375 get_vector_component(std::vector<VectorType *> &vec,
4376 const unsigned int component)
4379 return vec[component];
4390 typename VectorizedArrayType>
4391template <
typename VectorType,
typename VectorOperation>
4396 const std::array<VectorType *, n_components_> &src,
4399 n_components_> & src_sm,
4400 const std::bitset<VectorizedArrayType::size()> &mask,
4401 const bool apply_constraints)
const
4406 if (this->matrix_info ==
nullptr)
4408 read_write_operation_global(operation, src);
4414 if (n_fe_components == 1)
4415 for (
unsigned int comp = 0; comp < n_components; ++comp)
4417 Assert(src[comp] !=
nullptr,
4418 ExcMessage(
"The finite element underlying this FEEvaluation "
4419 "object is scalar, but you requested " +
4421 " components via the template argument in "
4422 "FEEvaluation. In that case, you must pass an "
4423 "std::vector<VectorType> or a BlockVector to " +
4424 "read_dof_values and distribute_local_to_global."));
4436 this->dof_info->index_storage_variants[this->dof_access_index].size());
4437 if (this->dof_info->index_storage_variants
4438 [is_face ? this->dof_access_index :
4443 read_write_operation_contiguous(operation, src, src_sm, mask);
4449 constexpr unsigned int n_lanes = VectorizedArrayType::size();
4450 Assert(mask.count() == n_lanes,
4452 "non-contiguous DoF storage"));
4454 std::integral_constant<
bool,
4455 internal::is_vectorizable<VectorType, Number>::value>
4458 const unsigned int dofs_per_component =
4459 this->data->dofs_per_component_on_cell;
4460 if (this->dof_info->index_storage_variants
4461 [is_face ? this->dof_access_index :
4466 const unsigned int *dof_indices =
4467 this->dof_info->dof_indices_interleaved.data() +
4468 this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4471 ->component_dof_indices_offset[this->active_fe_index]
4472 [this->first_selected_component] *
4474 if (n_components == 1 || n_fe_components == 1)
4475 for (
unsigned int i = 0; i < dofs_per_component;
4476 ++i, dof_indices += n_lanes)
4477 for (
unsigned int comp = 0; comp < n_components; ++comp)
4478 operation.process_dof_gather(dof_indices,
4481 values_dofs[comp][i],
4484 for (
unsigned int comp = 0; comp < n_components; ++comp)
4485 for (
unsigned int i = 0; i < dofs_per_component;
4486 ++i, dof_indices += n_lanes)
4487 operation.process_dof_gather(
4488 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4492 const unsigned int * dof_indices[n_lanes];
4493 VectorizedArrayType **values_dofs =
4494 const_cast<VectorizedArrayType **
>(&this->values_dofs[0]);
4498 unsigned int cells_copied[n_lanes];
4499 const unsigned int *cells;
4500 unsigned int n_vectorization_actual =
4502 ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4503 bool has_constraints =
false;
4504 const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4507 if (this->dof_access_index ==
4509 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4510 cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4512 this->dof_access_index ==
4515 (this->is_interior_face ?
4516 &this->matrix_info->
get_face_info(this->cell).cells_interior[0] :
4517 &this->matrix_info->
get_face_info(this->cell).cells_exterior[0]);
4518 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4520 Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4522 const std::pair<unsigned int, unsigned int> *my_index_start =
4523 &this->dof_info->row_starts[cells[v] * n_fe_components +
4524 first_selected_component];
4529 if (my_index_start[n_components_read].
second !=
4530 my_index_start[0].
second)
4531 has_constraints =
true;
4534 this->dof_info->dof_indices.data() + my_index_start[0].first;
4536 for (
unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4537 dof_indices[v] =
nullptr;
4542 this->dof_info->row_starts.size());
4543 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4545 const std::pair<unsigned int, unsigned int> *my_index_start =
4547 ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4548 first_selected_component];
4549 if (my_index_start[n_components_read].
second !=
4550 my_index_start[0].
second)
4551 has_constraints =
true;
4553 my_index_start[0].
first ||
4554 my_index_start[0].first < this->dof_info->dof_indices.size(),
4556 my_index_start[0].
first,
4557 this->dof_info->dof_indices.size()));
4559 this->dof_info->dof_indices.data() + my_index_start[0].first;
4561 for (
unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4562 dof_indices[v] =
nullptr;
4567 if (!has_constraints)
4569 if (n_vectorization_actual < n_lanes)
4570 for (
unsigned int comp = 0; comp < n_components; ++comp)
4571 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4572 operation.process_empty(values_dofs[comp][i]);
4573 if (n_components == 1 || n_fe_components == 1)
4575 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4576 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4577 for (
unsigned int comp = 0; comp < n_components; ++comp)
4578 operation.process_dof(dof_indices[v][i],
4580 values_dofs[comp][i][v]);
4584 for (
unsigned int comp = 0; comp < n_components; ++comp)
4585 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4586 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4587 operation.process_dof(
4588 dof_indices[v][comp * dofs_per_component + i],
4590 values_dofs[comp][i][v]);
4600 if (n_vectorization_actual < n_lanes)
4601 for (
unsigned int comp = 0; comp < n_components; ++comp)
4602 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4603 operation.process_empty(values_dofs[comp][i]);
4604 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4607 is_face ? cells[v] : this->cell * n_lanes + v;
4608 const unsigned int cell_dof_index =
4609 cell_index * n_fe_components + first_selected_component;
4610 const unsigned int n_components_read =
4611 n_fe_components > 1 ? n_components : 1;
4612 unsigned int index_indicators =
4613 this->dof_info->row_starts[cell_dof_index].second;
4614 unsigned int next_index_indicators =
4615 this->dof_info->row_starts[cell_dof_index + 1].second;
4619 if (apply_constraints ==
false &&
4620 this->dof_info->row_starts[cell_dof_index].second !=
4621 this->dof_info->row_starts[cell_dof_index + n_components_read]
4628 this->dof_info->plain_dof_indices.data() +
4630 ->component_dof_indices_offset[this->active_fe_index]
4631 [this->first_selected_component] +
4632 this->dof_info->row_starts_plain_indices[
cell_index];
4633 next_index_indicators = index_indicators;
4636 if (n_components == 1 || n_fe_components == 1)
4638 unsigned int ind_local = 0;
4639 for (; index_indicators != next_index_indicators; ++index_indicators)
4641 const std::pair<unsigned short, unsigned short> indicator =
4642 this->dof_info->constraint_indicator[index_indicators];
4644 for (
unsigned int j = 0; j < indicator.first; ++j)
4645 for (
unsigned int comp = 0; comp < n_components; ++comp)
4646 operation.process_dof(dof_indices[v][j],
4648 values_dofs[comp][ind_local + j][v]);
4650 ind_local += indicator.first;
4651 dof_indices[v] += indicator.first;
4655 Number
value[n_components];
4656 for (
unsigned int comp = 0; comp < n_components; ++comp)
4657 operation.pre_constraints(values_dofs[comp][ind_local][v],
4660 const Number *data_val =
4662 const Number *end_pool =
4664 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4665 for (
unsigned int comp = 0; comp < n_components; ++comp)
4666 operation.process_constraint(*dof_indices[v],
4671 for (
unsigned int comp = 0; comp < n_components; ++comp)
4672 operation.post_constraints(value[comp],
4673 values_dofs[comp][ind_local][v]);
4679 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4680 for (
unsigned int comp = 0; comp < n_components; ++comp)
4681 operation.process_dof(*dof_indices[v],
4683 values_dofs[comp][ind_local][v]);
4692 for (
unsigned int comp = 0; comp < n_components; ++comp)
4694 unsigned int ind_local = 0;
4697 for (; index_indicators != next_index_indicators;
4700 const std::pair<unsigned short, unsigned short> indicator =
4701 this->dof_info->constraint_indicator[index_indicators];
4704 for (
unsigned int j = 0; j < indicator.first; ++j)
4705 operation.process_dof(dof_indices[v][j],
4707 values_dofs[comp][ind_local + j][v]);
4708 ind_local += indicator.first;
4709 dof_indices[v] += indicator.first;
4714 operation.pre_constraints(values_dofs[comp][ind_local][v],
4717 const Number *data_val =
4719 const Number *end_pool =
4722 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4723 operation.process_constraint(*dof_indices[v],
4728 operation.post_constraints(value,
4729 values_dofs[comp][ind_local][v]);
4736 for (; ind_local < dofs_per_component;
4737 ++dof_indices[v], ++ind_local)
4740 operation.process_dof(*dof_indices[v],
4742 values_dofs[comp][ind_local][v]);
4745 if (apply_constraints ==
true && comp + 1 < n_components)
4746 next_index_indicators =
4747 this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4759 typename VectorizedArrayType>
4760template <
typename VectorType,
typename VectorOperation>
4765 const std::array<VectorType *, n_components_> &src)
const
4769 unsigned int index =
4770 first_selected_component * this->data->dofs_per_component_on_cell;
4771 for (
unsigned int comp = 0; comp < n_components; ++comp)
4773 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4776 operation.process_empty(values_dofs[comp][i]);
4777 operation.process_dof_global(
4778 local_dof_indices[this->data->lexicographic_numbering[index]],
4780 values_dofs[comp][i][0]);
4791 typename VectorizedArrayType>
4792template <
typename VectorType,
typename VectorOperation>
4797 const std::array<VectorType *, n_components_> &src,
4800 n_components_> & vectors_sm,
4801 const std::bitset<VectorizedArrayType::size()> &mask)
const
4810 std::integral_constant<
bool,
4811 internal::is_vectorizable<VectorType, Number>::value>
4814 is_face ? this->dof_access_index :
4816 const unsigned int n_lanes = mask.count();
4818 const std::vector<unsigned int> &dof_indices_cont =
4819 this->dof_info->dof_indices_contiguous[ind];
4823 if ((this->dof_info->index_storage_variants[ind][this->cell] ==
4825 interleaved_contiguous &&
4826 n_lanes == VectorizedArrayType::size()) &&
4828 this->dof_access_index ==
4830 this->is_interior_face ==
false))
4832 const unsigned int dof_index =
4833 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4834 this->dof_info->component_dof_indices_offset[this->active_fe_index]
4835 [first_selected_component] *
4836 VectorizedArrayType::size();
4837 if (n_components == 1 || n_fe_components == 1)
4838 for (
unsigned int comp = 0; comp < n_components; ++comp)
4839 operation.process_dofs_vectorized(
4840 this->data->dofs_per_component_on_cell,
4846 operation.process_dofs_vectorized(
4847 this->data->dofs_per_component_on_cell * n_components,
4855 std::array<
unsigned int, VectorizedArrayType::size()> cells =
4856 this->get_cell_or_face_ids();
4860 const unsigned int n_filled_lanes =
4861 this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4864 this->dof_access_index ==
4866 this->is_interior_face ==
false;
4868 if (vectors_sm[0] !=
nullptr)
4870 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
4871 std::array<
typename VectorType::value_type *,
4872 VectorizedArrayType::size()>
4875 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4879 Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
4881 ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
4883 this->dof_info->dof_indices_contiguous_sm[ind].size(),
4887 this->dof_info->dof_indices_contiguous_sm[ind].size()));
4890 this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
4893 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
4894 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
4895 this->dof_info->component_dof_indices_offset
4896 [this->active_fe_index][this->first_selected_component]);
4898 vector_ptrs[v] =
nullptr;
4900 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
4902 vector_ptrs[v] =
nullptr;
4907 if (n_filled_lanes == VectorizedArrayType::size() &&
4908 n_lanes == VectorizedArrayType::size() && !is_ecl)
4910 if (n_components == 1 || n_fe_components == 1)
4912 for (
unsigned int comp = 0; comp < n_components; ++comp)
4914 auto vector_ptrs = compute_vector_ptrs(comp);
4915 operation.process_dofs_vectorized_transpose(
4916 this->data->dofs_per_component_on_cell,
4924 auto vector_ptrs = compute_vector_ptrs(0);
4925 operation.process_dofs_vectorized_transpose(
4926 this->data->dofs_per_component_on_cell * n_components,
4933 for (
unsigned int comp = 0; comp < n_components; ++comp)
4935 auto vector_ptrs = compute_vector_ptrs(
4936 (n_components == 1 || n_fe_components == 1) ? comp : 0);
4938 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4940 operation.process_empty(values_dofs[comp][i]);
4942 if (n_components == 1 || n_fe_components == 1)
4944 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4945 if (mask[v] ==
true)
4946 for (
unsigned int i = 0;
4947 i < this->data->dofs_per_component_on_cell;
4949 operation.process_dof(vector_ptrs[v][i],
4950 values_dofs[comp][i][v]);
4954 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4955 if (mask[v] ==
true)
4956 for (
unsigned int i = 0;
4957 i < this->data->dofs_per_component_on_cell;
4959 operation.process_dof(
4961 [i + comp * this->data
4962 ->dofs_per_component_on_cell],
4963 values_dofs[comp][i][v]);
4969 unsigned int dof_indices[VectorizedArrayType::size()];
4971 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4975 dof_indices_cont[cells[v]] +
4977 ->component_dof_indices_offset[this->active_fe_index]
4978 [this->first_selected_component] *
4979 this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
4982 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4987 if (n_filled_lanes == VectorizedArrayType::size() &&
4988 n_lanes == VectorizedArrayType::size() && !is_ecl)
4990 if (this->dof_info->index_storage_variants[ind][this->cell] ==
4994 if (n_components == 1 || n_fe_components == 1)
4995 for (
unsigned int comp = 0; comp < n_components; ++comp)
4996 operation.process_dofs_vectorized_transpose(
4997 this->data->dofs_per_component_on_cell,
5003 operation.process_dofs_vectorized_transpose(
5004 this->data->dofs_per_component_on_cell * n_components,
5010 else if (this->dof_info->index_storage_variants[ind][this->cell] ==
5012 interleaved_contiguous_strided)
5014 if (n_components == 1 || n_fe_components == 1)
5015 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5018 for (
unsigned int comp = 0; comp < n_components; ++comp)
5019 operation.process_dof_gather(dof_indices,
5021 i * VectorizedArrayType::size(),
5022 values_dofs[comp][i],
5026 for (
unsigned int comp = 0; comp < n_components; ++comp)
5027 for (
unsigned int i = 0;
5028 i < this->data->dofs_per_component_on_cell;
5031 operation.process_dof_gather(
5034 (comp * this->data->dofs_per_component_on_cell + i) *
5035 VectorizedArrayType::size(),
5036 values_dofs[comp][i],
5042 Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
5044 IndexStorageVariants::interleaved_contiguous_mixed_strides,
5046 const unsigned int *offsets =
5047 &this->dof_info->dof_indices_interleave_strides
5048 [ind][VectorizedArrayType::size() * this->cell];
5049 if (n_components == 1 || n_fe_components == 1)
5050 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5053 for (
unsigned int comp = 0; comp < n_components; ++comp)
5054 operation.process_dof_gather(dof_indices,
5057 values_dofs[comp][i],
5060 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5061 dof_indices[v] += offsets[v];
5064 for (
unsigned int comp = 0; comp < n_components; ++comp)
5065 for (
unsigned int i = 0;
5066 i < this->data->dofs_per_component_on_cell;
5069 operation.process_dof_gather(dof_indices,
5072 values_dofs[comp][i],
5075 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5076 dof_indices[v] += offsets[v];
5081 for (
unsigned int comp = 0; comp < n_components; ++comp)
5083 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5085 operation.process_empty(values_dofs[comp][i]);
5086 if (this->dof_info->index_storage_variants[ind][this->cell] ==
5090 if (n_components == 1 || n_fe_components == 1)
5092 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5093 if (mask[v] ==
true)
5094 for (
unsigned int i = 0;
5095 i < this->data->dofs_per_component_on_cell;
5097 operation.process_dof(dof_indices[v] + i,
5099 values_dofs[comp][i][v]);
5103 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5104 if (mask[v] ==
true)
5105 for (
unsigned int i = 0;
5106 i < this->data->dofs_per_component_on_cell;
5108 operation.process_dof(
5109 dof_indices[v] + i +
5110 comp * this->data->dofs_per_component_on_cell,
5112 values_dofs[comp][i][v]);
5117 const unsigned int *offsets =
5118 &this->dof_info->dof_indices_interleave_strides
5119 [ind][VectorizedArrayType::size() * this->cell];
5120 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5122 if (n_components == 1 || n_fe_components == 1)
5123 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5125 if (mask[v] ==
true)
5126 for (
unsigned int i = 0;
5127 i < this->data->dofs_per_component_on_cell;
5129 operation.process_dof(dof_indices[v] + i * offsets[v],
5131 values_dofs[comp][i][v]);
5135 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5136 if (mask[v] ==
true)
5137 for (
unsigned int i = 0;
5138 i < this->data->dofs_per_component_on_cell;
5140 operation.process_dof(
5142 (i + comp * this->data->dofs_per_component_on_cell) *
5145 values_dofs[comp][i][v]);
5153 template <
typename Number,
5154 typename VectorType,
5155 typename std::enable_if<!IsBlockVector<VectorType>::value,
5156 VectorType>::type * =
nullptr>
5157 decltype(std::declval<VectorType>().begin())
5158 get_beginning(VectorType &vec)
5163 template <
typename Number,
5164 typename VectorType,
5165 typename std::enable_if<IsBlockVector<VectorType>::value,
5166 VectorType>::type * =
nullptr>
5167 typename VectorType::value_type *
5168 get_beginning(VectorType &)
5173 template <
typename VectorType,
5174 typename std::enable_if<has_shared_vector_data<VectorType>::value,
5175 VectorType>::type * =
nullptr>
5176 const std::vector<ArrayView<const typename VectorType::value_type>> *
5177 get_shared_vector_data(VectorType & vec,
5178 const bool is_valid_mode_for_sm,
5179 const unsigned int active_fe_index,
5183 if (is_valid_mode_for_sm &&
5186 active_fe_index == 0)
5187 return &vec.shared_vector_data();
5192 template <
typename VectorType,
5193 typename std::enable_if<!has_shared_vector_data<VectorType>::value,
5194 VectorType>::type * =
nullptr>
5195 const std::vector<ArrayView<const typename VectorType::value_type>> *
5196 get_shared_vector_data(VectorType &,
5204 template <
int n_components,
typename VectorType>
5206 std::array<
typename internal::BlockVectorSelector<
5207 typename std::remove_const<VectorType>::type,
5209 value>::BaseVectorType *,
5212 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
5213 typename std::remove_const<VectorType>::type,
5215 BaseVectorType::value_type>> *,
5217 get_vector_data(VectorType & src,
5218 const unsigned int first_index,
5219 const bool is_valid_mode_for_sm,
5220 const unsigned int active_fe_index,
5226 std::array<
typename internal::BlockVectorSelector<
5227 typename std::remove_const<VectorType>::type,
5229 value>::BaseVectorType *,
5233 ArrayView<
const typename internal::BlockVectorSelector<
5234 typename std::remove_const<VectorType>::type,
5236 value>::BaseVectorType::value_type>> *,
5240 for (
unsigned int d = 0;
d < n_components; ++
d)
5241 src_data.first[
d] = internal::BlockVectorSelector<
5242 typename std::remove_const<VectorType>::type,
5243 IsBlockVector<
typename std::remove_const<VectorType>::type>::value>::
5244 get_vector_component(
5245 const_cast<typename std::remove_const<VectorType>::type &
>(src),
5248 for (
unsigned int d = 0;
d < n_components; ++
d)
5249 src_data.second[
d] = get_shared_vector_data(*src_data.first[
d],
5250 is_valid_mode_for_sm,
5264 typename VectorizedArrayType>
5265template <
typename VectorType>
5268 read_dof_values(
const VectorType &src,
const unsigned int first_index)
5270 const auto src_data = internal::get_vector_data<n_components_>(
5273 this->dof_access_index ==
5275 this->active_fe_index,
5279 read_write_operation(reader,
5282 std::bitset<VectorizedArrayType::size()>().flip(),
5286 dof_values_initialized =
true;
5296 typename VectorizedArrayType>
5297template <
typename VectorType>
5302 const auto src_data = internal::get_vector_data<n_components_>(
5305 this->dof_access_index ==
5307 this->active_fe_index,
5311 read_write_operation(reader,
5314 std::bitset<VectorizedArrayType::size()>().flip(),
5318 dof_values_initialized =
true;
5328 typename VectorizedArrayType>
5329template <
typename VectorType>
5334 const unsigned int first_index,
5335 const std::bitset<VectorizedArrayType::size()> &mask)
const
5338 Assert(dof_values_initialized ==
true,
5342 const auto dst_data = internal::get_vector_data<n_components_>(
5345 this->dof_access_index ==
5347 this->active_fe_index,
5352 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
5361 typename VectorizedArrayType>
5362template <
typename VectorType>
5366 const unsigned int first_index,
5367 const std::bitset<VectorizedArrayType::size()> &mask)
const
5370 Assert(dof_values_initialized ==
true,
5374 const auto dst_data = internal::get_vector_data<n_components_>(
5377 this->dof_access_index ==
5379 this->active_fe_index,
5383 read_write_operation(setter, dst_data.first, dst_data.second, mask);
5392 typename VectorizedArrayType>
5393template <
typename VectorType>
5398 const unsigned int first_index,
5399 const std::bitset<VectorizedArrayType::size()> &mask)
const
5402 Assert(dof_values_initialized ==
true,
5406 const auto dst_data = internal::get_vector_data<n_components_>(
5409 this->dof_access_index ==
5411 this->active_fe_index,
5415 read_write_operation(setter, dst_data.first, dst_data.second, mask,
false);
5428 typename VectorizedArrayType>
5429inline const VectorizedArrayType *
5433 return &values_dofs[0][0];
5442 typename VectorizedArrayType>
5443inline VectorizedArrayType *
5448 dof_values_initialized =
true;
5450 return &values_dofs[0][0];
5459 typename VectorizedArrayType>
5460inline const VectorizedArrayType *
5476 typename VectorizedArrayType>
5477inline VectorizedArrayType *
5482 values_quad_initialized =
true;
5483 values_quad_submitted =
true;
5494 typename VectorizedArrayType>
5495inline const VectorizedArrayType *
5500 Assert(gradients_quad_initialized || gradients_quad_submitted,
5503 return gradients_quad;
5512 typename VectorizedArrayType>
5513inline VectorizedArrayType *
5518 gradients_quad_submitted =
true;
5519 gradients_quad_initialized =
true;
5521 return gradients_quad;
5530 typename VectorizedArrayType>
5531inline const VectorizedArrayType *
5538 return hessians_quad;
5547 typename VectorizedArrayType>
5548inline VectorizedArrayType *
5553 hessians_quad_initialized =
true;
5555 return hessians_quad;
5564 typename VectorizedArrayType>
5569 return first_selected_component;
5578 typename VectorizedArrayType>
5585 for (
unsigned int comp = 0; comp < n_components; comp++)
5586 return_value[comp] = this->values_dofs[comp][dof];
5587 return return_value;
5596 typename VectorizedArrayType>
5599 get_value(
const unsigned int q_point)
const
5602 Assert(this->values_quad_initialized ==
true,
5607 const std::size_t nqp = this->n_quadrature_points;
5609 for (
unsigned int comp = 0; comp < n_components; comp++)
5610 return_value[comp] = values_quad[comp * nqp + q_point];
5611 return return_value;
5620 typename VectorizedArrayType>
5627 Assert(this->gradients_quad_initialized ==
true,
5632 Assert(this->jacobian !=
nullptr,
5634 "update_gradients"));
5635 const std::size_t nqp = this->n_quadrature_points;
5641 for (
unsigned int d = 0;
d < dim; ++
d)
5642 for (
unsigned int comp = 0; comp < n_components; comp++)
5643 grad_out[comp][
d] = gradients_quad[(comp * dim +
d) * nqp + q_point] *
5644 this->jacobian[0][
d][
d];
5653 for (
unsigned int comp = 0; comp < n_components; comp++)
5654 for (
unsigned int d = 0;
d < dim; ++
d)
5657 jac[
d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5658 for (
unsigned int e = 1;
e < dim; ++
e)
5659 grad_out[comp][
d] +=
5660 jac[
d][
e] * gradients_quad[(comp * dim +
e) * nqp + q_point];
5672 typename VectorizedArrayType>
5679 Assert(this->gradients_quad_initialized ==
true,
5683 Assert(this->normal_x_jacobian !=
nullptr,
5685 "update_gradients"));
5687 const std::size_t nqp = this->n_quadrature_points;
5691 for (
unsigned int comp = 0; comp < n_components; comp++)
5692 grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5693 (this->normal_x_jacobian[0][dim - 1]);
5696 const std::size_t index =
5698 for (
unsigned int comp = 0; comp < n_components; comp++)
5700 grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5701 this->normal_x_jacobian[index][0];
5702 for (
unsigned int d = 1;
d < dim; ++
d)
5703 grad_out[comp] += gradients_quad[(comp * dim +
d) * nqp + q_point] *
5704 this->normal_x_jacobian[index][
d];
5716 template <
typename VectorizedArrayType>
5719 const VectorizedArrayType *
const hessians,
5721 VectorizedArrayType (&tmp)[1][1])
5723 tmp[0][0] = jac[0][0] *
hessians[0];
5726 template <
typename VectorizedArrayType>
5729 const VectorizedArrayType *
const hessians,
5730 const unsigned int nqp,
5731 VectorizedArrayType (&tmp)[2][2])
5733 for (
unsigned int d = 0;
d < 2; ++
d)
5741 template <
typename VectorizedArrayType>
5744 const VectorizedArrayType *
const hessians,
5745 const unsigned int nqp,
5746 VectorizedArrayType (&tmp)[3][3])
5748 for (
unsigned int d = 0;
d < 3; ++
d)
5769 typename VectorizedArrayType>
5776 Assert(this->hessians_quad_initialized ==
true,
5781 Assert(this->jacobian !=
nullptr,
5791 const std::size_t nqp = this->n_quadrature_points;
5792 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5797 for (
unsigned int comp = 0; comp < n_components; comp++)
5799 for (
unsigned int d = 0;
d < dim; ++
d)
5800 hessian_out[comp][
d][
d] =
5801 hessians_quad[(comp * hdim +
d) * nqp + q_point] *
5802 (jac[
d][
d] * jac[
d][
d]);
5808 hessian_out[comp][0][1] =
5809 hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5810 (jac[0][0] * jac[1][1]);
5813 hessian_out[comp][0][1] =
5814 hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5815 (jac[0][0] * jac[1][1]);
5816 hessian_out[comp][0][2] =
5817 hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5818 (jac[0][0] * jac[2][2]);
5819 hessian_out[comp][1][2] =
5820 hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5821 (jac[1][1] * jac[2][2]);
5826 for (
unsigned int d = 0;
d < dim; ++
d)
5827 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5828 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5834 for (
unsigned int comp = 0; comp < n_components; comp++)
5836 VectorizedArrayType tmp[dim][dim];
5837 internal::hessian_unit_times_jac(
5838 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5841 for (
unsigned int d = 0;
d < dim; ++
d)
5842 for (
unsigned int e =
d;
e < dim; ++
e)
5844 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
5845 for (
unsigned int f = 1; f < dim; ++f)
5846 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
5853 for (
unsigned int d = 0;
d < dim; ++
d)
5854 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5855 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5861 const auto &jac_grad =
5862 this->mapping_data->jacobian_gradients
5863 [1 - this->is_interior_face]
5864 [this->mapping_data->data_index_offsets[this->cell] + q_point];
5865 for (
unsigned int comp = 0; comp < n_components; comp++)
5869 VectorizedArrayType tmp[dim][dim];
5870 internal::hessian_unit_times_jac(
5871 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5874 for (
unsigned int d = 0;
d < dim; ++
d)
5875 for (
unsigned int e =
d;
e < dim; ++
e)
5877 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
5878 for (
unsigned int f = 1; f < dim; ++f)
5879 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
5883 for (
unsigned int d = 0;
d < dim; ++
d)
5884 for (
unsigned int e = 0;
e < dim; ++
e)
5885 hessian_out[comp][
d][
d] +=
5887 gradients_quad[(comp * dim +
e) * nqp + q_point];
5890 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
5891 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++count)
5892 for (
unsigned int f = 0; f < dim; ++f)
5893 hessian_out[comp][
d][
e] +=
5894 jac_grad[count][f] *
5895 gradients_quad[(comp * dim + f) * nqp + q_point];
5898 for (
unsigned int d = 0;
d < dim; ++
d)
5899 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5900 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5912 typename VectorizedArrayType>
5919 Assert(this->hessians_quad_initialized ==
true,
5930 const std::size_t nqp = this->n_quadrature_points;
5931 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5937 for (
unsigned int comp = 0; comp < n_components; comp++)
5938 for (
unsigned int d = 0;
d < dim; ++
d)
5939 hessian_out[comp][
d] =
5940 hessians_quad[(comp * hdim +
d) * nqp + q_point] *
5941 (jac[
d][
d] * jac[
d][
d]);
5946 for (
unsigned int comp = 0; comp < n_components; comp++)
5950 VectorizedArrayType tmp[dim][dim];
5951 internal::hessian_unit_times_jac(
5952 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5956 for (
unsigned int d = 0;
d < dim; ++
d)
5958 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
5959 for (
unsigned int f = 1; f < dim; ++f)
5960 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
5969 this->mapping_data->jacobian_gradients
5970 [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5971 for (
unsigned int comp = 0; comp < n_components; comp++)
5975 VectorizedArrayType tmp[dim][dim];
5976 internal::hessian_unit_times_jac(
5977 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5981 for (
unsigned int d = 0;
d < dim; ++
d)
5983 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
5984 for (
unsigned int f = 1; f < dim; ++f)
5985 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
5988 for (
unsigned int d = 0;
d < dim; ++
d)
5989 for (
unsigned int e = 0;
e < dim; ++
e)
5990 hessian_out[comp][
d] +=
5992 gradients_quad[(comp * dim +
e) * nqp + q_point];
6004 typename VectorizedArrayType>
6011 Assert(this->hessians_quad_initialized ==
true,
6017 const auto hess_diag = get_hessian_diagonal(q_point);
6018 for (
unsigned int comp = 0; comp < n_components; ++comp)
6020 laplacian_out[comp] = hess_diag[comp][0];
6021 for (
unsigned int d = 1;
d < dim; ++
d)
6022 laplacian_out[comp] += hess_diag[comp][
d];
6024 return laplacian_out;
6033 typename VectorizedArrayType>
6037 const unsigned int dof)
6040 this->dof_values_initialized =
true;
6043 for (
unsigned int comp = 0; comp < n_components; comp++)
6044 this->values_dofs[comp][dof] = val_in[comp];
6053 typename VectorizedArrayType>
6057 const unsigned int q_point)
6061 Assert(this->J_value !=
nullptr,
6065 this->values_quad_submitted =
true;
6068 const std::size_t nqp = this->n_quadrature_points;
6071 const VectorizedArrayType JxW =
6072 this->J_value[0] * this->quadrature_weights[q_point];
6073 for (
unsigned int comp = 0; comp < n_components; ++comp)
6074 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6078 const VectorizedArrayType JxW = this->J_value[q_point];
6079 for (
unsigned int comp = 0; comp < n_components; ++comp)
6080 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6090 typename VectorizedArrayType>
6095 const unsigned int q_point)
6099 Assert(this->J_value !=
nullptr,
6101 "update_gradients"));
6102 Assert(this->jacobian !=
nullptr,
6104 "update_gradients"));
6106 this->gradients_quad_submitted =
true;
6109 const std::size_t nqp = this->n_quadrature_points;
6112 const VectorizedArrayType JxW =
6113 this->J_value[0] * this->quadrature_weights[q_point];
6114 for (
unsigned int d = 0;
d < dim; ++
d)
6116 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
6117 for (
unsigned int comp = 0; comp < n_components; comp++)
6118 gradients_quad[(comp * dim +
d) * nqp + q_point] =
6119 grad_in[comp][
d] * factor;
6126 this->jacobian[q_point] :
6128 const VectorizedArrayType JxW =
6130 this->J_value[q_point] :
6131 this->J_value[0] * this->quadrature_weights[q_point];
6132 for (
unsigned int comp = 0; comp < n_components; ++comp)
6133 for (
unsigned int d = 0;
d < dim; ++
d)
6135 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
6136 for (
unsigned int e = 1;
e < dim; ++
e)
6137 new_val += (jac[
e][
d] * grad_in[comp][
e]);
6138 gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val * JxW;
6149 typename VectorizedArrayType>
6154 const unsigned int q_point)
6157 Assert(this->normal_x_jacobian !=
nullptr,
6159 "update_gradients"));
6161 this->gradients_quad_submitted =
true;
6164 const std::size_t nqp = this->n_quadrature_points;
6166 for (
unsigned int comp = 0; comp < n_components; comp++)
6168 for (
unsigned int d = 0;
d < dim - 1; ++
d)
6169 gradients_quad[(comp * dim +
d) * nqp + q_point] =
6170 VectorizedArrayType();
6171 gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
6173 (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
6174 this->quadrature_weights[q_point]);
6178 const unsigned int index =
6181 this->normal_x_jacobian[index];
6182 for (
unsigned int comp = 0; comp < n_components; comp++)
6184 VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
6186 factor = factor * this->quadrature_weights[q_point];
6187 for (
unsigned int d = 0;
d < dim; ++
d)
6188 gradients_quad[(comp * dim +
d) * nqp + q_point] = factor * jac[
d];
6199 typename VectorizedArrayType>
6206 Assert(this->values_quad_submitted ==
true,
6211 const std::size_t nqp = this->n_quadrature_points;
6212 for (
unsigned int q = 0; q < nqp; ++q)
6213 for (
unsigned int comp = 0; comp < n_components; ++comp)
6214 return_value[comp] += this->values_quad[comp * nqp + q];
6215 return (return_value);
6227 typename VectorizedArrayType>
6232 VectorizedArrayType>::
6235 const unsigned int dof_no,
6236 const unsigned int first_selected_component,
6237 const unsigned int quad_no_in,
6238 const unsigned int fe_degree,
6239 const unsigned int n_q_points,
6240 const bool is_interior_face,
6241 const unsigned int active_fe_index,
6242 const unsigned int active_quad_index,
6243 const unsigned int face_type)
6244 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6247 first_selected_component,
6263 typename VectorizedArrayType>
6268 VectorizedArrayType>::
6274 const unsigned int first_selected_component,
6277 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6282 first_selected_component,
6292 typename VectorizedArrayType>
6297 VectorizedArrayType>::
6302 VectorizedArrayType> &other)
6303 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6313 typename VectorizedArrayType>
6318 VectorizedArrayType> &
6324 VectorizedArrayType> &other)
6330 VectorizedArrayType>::operator=(other);
6339template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6343 const unsigned int dof_no,
6344 const unsigned int first_selected_component,
6345 const unsigned int quad_no_in,
6346 const unsigned int fe_degree,
6347 const unsigned int n_q_points,
6348 const bool is_interior_face,
6349 const unsigned int active_fe_index,
6350 const unsigned int active_quad_index,
6351 const unsigned int face_type)
6355 first_selected_component,
6367template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6374 const unsigned int first_selected_component,
6382 first_selected_component,
6388template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6398template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6410template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6413 const unsigned int dof)
const
6416 return this->values_dofs[0][dof];
6421template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6424 const unsigned int q_point)
const
6427 Assert(this->values_quad_initialized ==
true,
6431 return this->values_quad[q_point];
6436template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6441 return BaseClass::get_normal_derivative(q_point)[0];
6446template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6449 const unsigned int q_point)
const
6455 Assert(this->gradients_quad_initialized ==
true,
6460 Assert(this->jacobian !=
nullptr,
6462 "update_gradients"));
6466 const std::size_t nqp = this->n_quadrature_points;
6469 for (
unsigned int d = 0;
d < dim; ++
d)
6471 this->gradients_quad[
d * nqp + q_point] * this->jacobian[0][
d][
d];
6480 for (
unsigned int d = 0;
d < dim; ++
d)
6482 grad_out[
d] = jac[
d][0] * this->gradients_quad[q_point];
6483 for (
unsigned int e = 1;
e < dim; ++
e)
6484 grad_out[
d] += jac[
d][
e] * this->gradients_quad[
e * nqp + q_point];
6492template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6495 const unsigned int q_point)
const
6497 return BaseClass::get_hessian(q_point)[0];
6502template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6507 return BaseClass::get_hessian_diagonal(q_point)[0];
6512template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6513inline VectorizedArrayType
6515 const unsigned int q_point)
const
6517 return BaseClass::get_laplacian(q_point)[0];
6522template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6528 this->dof_values_initialized =
true;
6531 this->values_dofs[0][dof] = val_in;
6536template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6539 const VectorizedArrayType val_in,
6540 const unsigned int q_point)
6544 Assert(this->J_value !=
nullptr,
6548 this->values_quad_submitted =
true;
6553 const VectorizedArrayType JxW =
6554 this->J_value[0] * this->quadrature_weights[q_point];
6555 this->values_quad[q_point] = val_in * JxW;
6559 this->values_quad[q_point] = val_in * this->J_value[q_point];
6565template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6569 const unsigned int q_point)
6571 submit_value(val_in[0], q_point);
6576template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6580 const unsigned int q_point)
6584 BaseClass::submit_normal_derivative(grad, q_point);
6589template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6593 const unsigned int q_point)
6597 Assert(this->J_value !=
nullptr,
6599 "update_gradients"));
6600 Assert(this->jacobian !=
nullptr,
6602 "update_gradients"));
6604 this->gradients_quad_submitted =
true;
6607 const std::size_t nqp = this->n_quadrature_points;
6610 const VectorizedArrayType JxW =
6611 this->J_value[0] * this->quadrature_weights[q_point];
6612 for (
unsigned int d = 0;
d < dim; ++
d)
6613 this->gradients_quad[
d * nqp + q_point] =
6614 (grad_in[
d] * this->jacobian[0][
d][
d] * JxW);
6621 this->jacobian[q_point] :
6623 const VectorizedArrayType JxW =
6625 this->J_value[q_point] :
6626 this->J_value[0] * this->quadrature_weights[q_point];
6627 for (
unsigned int d = 0;
d < dim; ++
d)
6629 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
6630 for (
unsigned int e = 1;
e < dim; ++
e)
6631 new_val += jac[
e][
d] * grad_in[
e];
6632 this->gradients_quad[
d * nqp + q_point] = new_val * JxW;
6639template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6640inline VectorizedArrayType
6644 return BaseClass::integrate_value()[0];
6652template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6656 const unsigned int dof_no,
6657 const unsigned int first_selected_component,
6658 const unsigned int quad_no_in,
6659 const unsigned int fe_degree,
6660 const unsigned int n_q_points,
6661 const bool is_interior_face,
6662 const unsigned int active_fe_index,
6663 const unsigned int active_quad_index,
6664 const unsigned int face_type)
6668 first_selected_component,
6680template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6687 const unsigned int first_selected_component,
6695 first_selected_component,
6701template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6711template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6724template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6729 return BaseClass::get_gradient(q_point);
6734template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6740 Assert(this->gradients_quad_initialized ==
true,
6744 Assert(this->jacobian !=
nullptr,
6746 "update_gradients"));
6748 VectorizedArrayType divergence;
6749 const std::size_t nqp = this->n_quadrature_points;
6754 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6755 for (
unsigned int d = 1;
d < dim; ++
d)
6756 divergence += this->gradients_quad[(dim *
d +
d) * nqp + q_point] *
6757 this->jacobian[0][
d][
d];
6764 this->jacobian[q_point] :
6766 divergence = jac[0][0] * this->gradients_quad[q_point];
6767 for (
unsigned int e = 1;
e < dim; ++
e)
6768 divergence += jac[0][
e] * this->gradients_quad[
e * nqp + q_point];
6769 for (
unsigned int d = 1;
d < dim; ++
d)
6770 for (
unsigned int e = 0;
e < dim; ++
e)
6772 jac[
d][
e] * this->gradients_quad[(
d * dim +
e) * nqp + q_point];
6779template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6785 const auto grad = get_gradient(q_point);
6786 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6787 VectorizedArrayType half = Number(0.5);
6788 for (
unsigned int d = 0;
d < dim; ++
d)
6789 symmetrized[
d] = grad[
d][
d];
6795 symmetrized[2] = grad[0][1] + grad[1][0];
6796 symmetrized[2] *= half;
6799 symmetrized[3] = grad[0][1] + grad[1][0];
6800 symmetrized[3] *= half;
6801 symmetrized[4] = grad[0][2] + grad[2][0];
6802 symmetrized[4] *= half;
6803 symmetrized[5] = grad[1][2] + grad[2][1];
6804 symmetrized[5] *= half;
6814template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6816 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6818 const unsigned int q_point)
const
6822 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6828 "Computing the curl in 1d is not a useful operation"));
6831 curl[0] = grad[1][0] - grad[0][1];
6834 curl[0] = grad[2][1] - grad[1][2];
6835 curl[1] = grad[0][2] - grad[2][0];
6836 curl[2] = grad[1][0] - grad[0][1];
6846template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6851 return BaseClass::get_hessian_diagonal(q_point);
6856template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6859 const unsigned int q_point)
const
6862 Assert(this->hessians_quad_initialized ==
true,
6866 return BaseClass::get_hessian(q_point);
6871template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6875 const unsigned int q_point)
6877 BaseClass::submit_gradient(grad_in, q_point);
6882template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6887 const unsigned int q_point)
6889 BaseClass::submit_gradient(grad_in, q_point);
6894template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6898 const unsigned int q_point)
6902 Assert(this->J_value !=
nullptr,
6904 "update_gradients"));
6905 Assert(this->jacobian !=
nullptr,
6907 "update_gradients"));
6909 this->gradients_quad_submitted =
true;
6912 const std::size_t nqp = this->n_quadrature_points;
6915 const VectorizedArrayType fac =
6916 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6917 for (
unsigned int d = 0;
d < dim; ++
d)
6919 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6920 (fac * this->jacobian[0][
d][
d]);
6921 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6923 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6924 VectorizedArrayType();
6925 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6926 VectorizedArrayType();
6934 this->jacobian[q_point] :
6936 const VectorizedArrayType fac =
6938 this->J_value[q_point] :
6939 this->J_value[0] * this->quadrature_weights[q_point]) *
6941 for (
unsigned int d = 0;
d < dim; ++
d)
6943 for (
unsigned int e = 0;
e < dim; ++
e)
6944 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6952template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6957 const unsigned int q_point)
6964 Assert(this->J_value !=
nullptr,
6966 "update_gradients"));
6967 Assert(this->jacobian !=
nullptr,
6969 "update_gradients"));
6971 this->gradients_quad_submitted =
true;
6974 const std::size_t nqp = this->n_quadrature_points;
6977 const VectorizedArrayType JxW =
6978 this->J_value[0] * this->quadrature_weights[q_point];
6979 for (
unsigned int d = 0;
d < dim; ++
d)
6980 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6982 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6983 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6985 const VectorizedArrayType
value =
6987 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6988 value * this->jacobian[0][
d][
d];
6989 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6990 value * this->jacobian[0][
e][
e];
6996 const VectorizedArrayType JxW =
6998 this->J_value[q_point] :
6999 this->J_value[0] * this->quadrature_weights[q_point];
7002 this->jacobian[q_point] :
7004 VectorizedArrayType weighted[dim][dim];
7005 for (
unsigned int i = 0; i < dim; ++i)
7007 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
7008 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
7010 const VectorizedArrayType
value =
7012 weighted[i][j] =
value;
7013 weighted[j][i] =
value;
7015 for (
unsigned int comp = 0; comp < dim; ++comp)
7016 for (
unsigned int d = 0;
d < dim; ++
d)
7018 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
7019 for (
unsigned int e = 1;
e < dim; ++
e)
7020 new_val += jac[
e][
d] * weighted[comp][
e];
7021 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
7028template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
7032 const unsigned int q_point)
7040 "Testing by the curl in 1d is not a useful operation"));
7043 grad[1][0] = curl[0];
7044 grad[0][1] = -curl[0];
7047 grad[2][1] = curl[0];
7048 grad[1][2] = -curl[0];
7049 grad[0][2] = curl[1];
7050 grad[2][0] = -curl[1];
7051 grad[1][0] = curl[2];
7052 grad[0][1] = -curl[2];
7057 submit_gradient(grad, q_point);
7064template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7067 const unsigned int dof_no,
7068 const unsigned int first_selected_component,
7069 const unsigned int quad_no_in,
7070 const unsigned int fe_degree,
7071 const unsigned int n_q_points,
7072 const bool is_interior_face,
7073 const unsigned int active_fe_index,
7074 const unsigned int active_quad_index,
7075 const unsigned int face_type)
7079 first_selected_component,
7091template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7098 const unsigned int first_selected_component,
7105 first_selected_component,
7111template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7120template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7132template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7135 const unsigned int dof)
const
7138 return this->values_dofs[0][dof];
7143template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7146 const unsigned int q_point)
const
7149 Assert(this->values_quad_initialized ==
true,
7153 return this->values_quad[q_point];
7158template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7161 const unsigned int q_point)
const
7167 Assert(this->gradients_quad_initialized ==
true,
7174 this->jacobian[q_point] :
7178 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
7185template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7188 const unsigned int q_point)
const
7190 return get_gradient(q_point)[0];
7195template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7200 return BaseClass::get_normal_derivative(q_point)[0];
7205template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7208 const unsigned int q_point)
const
7210 return BaseClass::get_hessian(q_point)[0];
7215template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7220 return BaseClass::get_hessian_diagonal(q_point)[0];
7225template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7228 const unsigned int q_point)
const
7230 return BaseClass::get_laplacian(q_point)[0];
7235template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7241 this->dof_values_initialized =
true;
7244 this->values_dofs[0][dof] = val_in;
7249template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7252 const VectorizedArrayType val_in,
7253 const unsigned int q_point)
7258 this->values_quad_submitted =
true;
7263 const VectorizedArrayType JxW = this->J_value[q_point];
7264 this->values_quad[q_point] = val_in * JxW;
7268 const VectorizedArrayType JxW =
7269 this->J_value[0] * this->quadrature_weights[q_point];
7270 this->values_quad[q_point] = val_in * JxW;
7276template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7280 const unsigned int q_point)
7282 submit_value(val_in[0], q_point);
7287template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7291 const unsigned int q_point)
7293 submit_gradient(grad_in[0], q_point);
7298template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7301 const VectorizedArrayType grad_in,
7302 const unsigned int q_point)
7307 this->gradients_quad_submitted =
true;
7312 this->jacobian[q_point] :
7314 const VectorizedArrayType JxW =
7316 this->J_value[q_point] :
7317 this->J_value[0] * this->quadrature_weights[q_point];
7319 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7324template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7328 const unsigned int q_point)
7330 submit_gradient(grad_in[0][0], q_point);
7335template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7339 const unsigned int q_point)
7343 BaseClass::submit_normal_derivative(grad, q_point);
7348template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7352 const unsigned int q_point)
7354 BaseClass::submit_normal_derivative(grad_in, q_point);
7359template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7360inline VectorizedArrayType
7364 return BaseClass::integrate_value()[0];
7377 typename VectorizedArrayType>
7383 VectorizedArrayType>::
7385 const unsigned int fe_no,
7386 const unsigned int quad_no,
7387 const unsigned int first_selected_component,
7388 const unsigned int active_fe_index,
7389 const unsigned int active_quad_index)
7390 : BaseClass(data_in,
7392 first_selected_component,
7399 , dofs_per_component(this->data->dofs_per_component_on_cell)
7400 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7401 , n_q_points(this->data->n_q_points)
7403 check_template_arguments(fe_no, 0);
7413 typename VectorizedArrayType>
7419 VectorizedArrayType>::
7421 const std::pair<unsigned int, unsigned int> & range,
7422 const unsigned int dof_no,
7423 const unsigned int quad_no,
7424 const unsigned int first_selected_component)
7428 first_selected_component,
7429 matrix_free.get_cell_active_fe_index(range))
7439 typename VectorizedArrayType>
7445 VectorizedArrayType>::
7450 const unsigned int first_selected_component)
7451 : BaseClass(mapping,
7455 first_selected_component,
7457 , dofs_per_component(this->data->dofs_per_component_on_cell)
7458 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7459 , n_q_points(this->data->n_q_points)
7471 typename VectorizedArrayType>
7477 VectorizedArrayType>::
7481 const unsigned int first_selected_component)
7486 first_selected_component,
7488 , dofs_per_component(this->data->dofs_per_component_on_cell)
7489 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7490 , n_q_points(this->data->n_q_points)
7502 typename VectorizedArrayType>
7508 VectorizedArrayType>::
7512 const unsigned int first_selected_component)
7513 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7515 other.mapped_geometry->get_quadrature(),
7516 other.mapped_geometry->get_fe_values().get_update_flags(),
7517 first_selected_component,
7519 , dofs_per_component(this->data->dofs_per_component_on_cell)
7520 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7521 , n_q_points(this->data->n_q_points)
7533 typename VectorizedArrayType>
7542 , dofs_per_component(this->data->dofs_per_component_on_cell)
7543 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7544 , n_q_points(this->data->n_q_points)
7556 typename VectorizedArrayType>
7562 VectorizedArrayType> &
7568 VectorizedArrayType>::operator=(
const FEEvaluation &other)
7570 BaseClass::operator=(other);
7582 typename VectorizedArrayType>
7589 VectorizedArrayType>::
7590 check_template_arguments(
const unsigned int dof_no,
7591 const unsigned int first_selected_component)
7594 (void)first_selected_component;
7600 static_cast<unsigned int>(fe_degree) !=
7601 this->data->data.front().fe_degree) ||
7602 n_q_points != this->n_quadrature_points)
7604 std::string message =
7605 "-------------------------------------------------------\n";
7606 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7607 message +=
" Called --> FEEvaluation<dim,";
7611 message +=
",Number>(data";
7627 if (
static_cast<unsigned int>(fe_degree) ==
7628 this->data->data.front().fe_degree)
7630 proposed_dof_comp = dof_no;
7631 proposed_fe_comp = first_selected_component;
7634 for (
unsigned int no = 0; no < this->matrix_info->
n_components();
7636 for (
unsigned int nf = 0;
7639 if (this->matrix_info
7640 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7642 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7644 proposed_dof_comp = no;
7645 proposed_fe_comp = nf;
7649 this->mapping_data->descriptor[this->active_quad_index]
7651 proposed_quad_comp = this->quad_no;
7653 for (
unsigned int no = 0;
7658 .descriptor[this->active_quad_index]
7659 .n_q_points == n_q_points)
7661 proposed_quad_comp = no;
7668 if (proposed_dof_comp != first_selected_component)
7669 message +=
"Wrong vector component selection:\n";
7671 message +=
"Wrong quadrature formula selection:\n";
7672 message +=
" Did you mean FEEvaluation<dim,";
7676 message +=
",Number>(data";
7685 std::string correct_pos;
7686 if (proposed_dof_comp != dof_no)
7687 correct_pos =
" ^ ";
7690 if (proposed_quad_comp != this->quad_no)
7691 correct_pos +=
" ^ ";
7694 if (proposed_fe_comp != first_selected_component)
7695 correct_pos +=
" ^\n";
7697 correct_pos +=
" \n";
7703 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7704 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7705 message +=
"Wrong template arguments:\n";
7706 message +=
" Did you mean FEEvaluation<dim,";
7711 message +=
",Number>(data";
7719 std::string correct_pos;
7720 if (this->data->data.front().fe_degree !=
7721 static_cast<unsigned int>(fe_degree))
7725 if (proposed_n_q_points_1d != n_q_points_1d)
7726 correct_pos +=
" ^\n";
7728 correct_pos +=
" \n";
7729 message +=
" " + correct_pos;
7731 Assert(
static_cast<unsigned int>(fe_degree) ==
7732 this->data->data.front().fe_degree &&
7733 n_q_points == this->n_quadrature_points,
7739 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7750 typename VectorizedArrayType>
7759 Assert(this->mapped_geometry ==
nullptr,
7760 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7761 " Integer indexing is not possible"));
7762 if (this->mapped_geometry !=
nullptr)
7771 const unsigned int offsets =
7772 this->mapping_data->data_index_offsets[
cell_index];
7773 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7774 this->J_value = &this->mapping_data->JxW_values[offsets];
7777 this->dof_values_initialized =
false;
7778 this->values_quad_initialized =
false;
7779 this->gradients_quad_initialized =
false;
7780 this->hessians_quad_initialized =
false;
7791 typename VectorizedArrayType>
7792template <
bool level_dof_access>
7799 VectorizedArrayType>
::
7802 Assert(this->matrix_info ==
nullptr,
7803 ExcMessage(
"Cannot use initialization from cell iterator if "
7804 "initialized from MatrixFree object. Use variant for "
7805 "on the fly computation with arguments as for FEValues "
7808 this->mapped_geometry->reinit(
7810 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7811 if (level_dof_access)
7812 cell->get_mg_dof_indices(this->local_dof_indices);
7814 cell->get_dof_indices(this->local_dof_indices);
7824 typename VectorizedArrayType>
7831 VectorizedArrayType>
::
7834 Assert(this->matrix_info == 0,
7835 ExcMessage(
"Cannot use initialization from cell iterator if "
7836 "initialized from MatrixFree object. Use variant for "
7837 "on the fly computation with arguments as for FEValues "
7840 this->mapped_geometry->reinit(cell);
7850 typename VectorizedArrayType>
7857 VectorizedArrayType>::quadrature_point(
const unsigned int q)
const
7859 if (this->matrix_info ==
nullptr)
7861 Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7864 "update_quadrature_points"));
7868 Assert(this->mapping_data->quadrature_point_offsets.empty() ==
false,
7870 "update_quadrature_points"));
7876 &this->mapping_data->quadrature_points
7877 [this->mapping_data->quadrature_point_offsets[this->cell]];
7889 for (
unsigned int d = 0;
d < dim; ++
d)
7890 point[
d] += jac[
d][
d] *
static_cast<Number
>(
7891 this->descriptor->quadrature.point(q)[
d]);
7893 for (
unsigned int d = 0;
d < dim; ++
d)
7894 for (
unsigned int e = 0;
e < dim; ++
e)
7895 point[
d] += jac[
d][
e] *
static_cast<Number
>(
7896 this->descriptor->quadrature.point(q)[
e]);
7910 typename VectorizedArrayType>
7917 VectorizedArrayType>::evaluate(
const bool evaluate_values,
7918 const bool evaluate_gradients,
7919 const bool evaluate_hessians)
7922 Assert(this->dof_values_initialized ==
true,
7925 evaluate(this->values_dofs[0],
7937 typename VectorizedArrayType>
7944 VectorizedArrayType>::
7948 Assert(this->dof_values_initialized ==
true,
7951 evaluate(this->values_dofs[0], evaluation_flags);
7961 typename VectorizedArrayType>
7968 VectorizedArrayType>::evaluate(
const VectorizedArrayType
7970 const bool evaluate_values,
7971 const bool evaluate_gradients,
7972 const bool evaluate_hessians)
7981 evaluate(values_array, flag);
7991 typename VectorizedArrayType>
7998 VectorizedArrayType>::
7999 evaluate(
const VectorizedArrayType * values_array,
8007 const_cast<VectorizedArrayType *
>(values_array),
8009 this->gradients_quad,
8010 this->hessians_quad,
8011 this->scratch_data);
8017 const_cast<VectorizedArrayType *
>(values_array),
8019 this->gradients_quad,
8020 this->hessians_quad,
8021 this->scratch_data);
8025 this->values_quad_initialized =
true;
8027 this->gradients_quad_initialized =
true;
8029 this->hessians_quad_initialized =
true;
8040 typename VectorizedArrayType>
8041template <
typename VectorType>
8049 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
8050 const bool evaluate_values,
8051 const bool evaluate_gradients,
8052 const bool evaluate_hessians)
8061 gather_evaluate(input_vector, flag);
8070 template <
typename Number,
8071 typename VectorizedArrayType,
8072 typename VectorType,
8074 typename std::enable_if<
8075 internal::has_begin<VectorType>::value &&
8076 std::is_same<decltype(std::declval<VectorType>().begin()),
8078 VectorType>::type * =
nullptr>
8080 try_gather_evaluate_inplace(
8082 const VectorType & input_vector,
8083 const unsigned int cell,
8084 const unsigned int active_fe_index,
8085 const unsigned int first_selected_component,
8093 if (std::is_same<typename VectorType::value_type, Number>::value &&
8097 interleaved_contiguous &&
8098 reinterpret_cast<std::size_t
>(
8099 input_vector.begin() +
8102 [cell * VectorizedArrayType::size()]) %
8103 sizeof(VectorizedArrayType) ==
8106 const VectorizedArrayType *vec_values =
8107 reinterpret_cast<const VectorizedArrayType *
>(
8108 input_vector.begin() +
8111 [cell * VectorizedArrayType::size()] +
8113 [first_selected_component] *
8114 VectorizedArrayType::size());
8116 phi->evaluate(vec_values, evaluation_flag);
8127 template <
typename Number,
8128 typename VectorizedArrayType,
8129 typename VectorType,
8131 typename std::enable_if<
8132 !internal::has_begin<VectorType>::value ||
8133 !std::is_same<decltype(std::declval<VectorType>().begin()),
8135 VectorType>::type * =
nullptr>
8137 try_gather_evaluate_inplace(
T,
8155 typename VectorizedArrayType,
8156 typename VectorType,
8157 typename std::enable_if<
8158 internal::has_begin<VectorType>::value &&
8159 std::is_same<decltype(std::declval<VectorType>().begin()),
8161 VectorType>::type * =
nullptr>
8163 try_integrate_scatter_inplace(
8164 VectorType & destination,
8165 const unsigned int cell,
8166 const unsigned int n_components,
8167 const unsigned int active_fe_index,
8168 const unsigned int first_selected_component,
8170 VectorizedArrayType * values_quad,
8171 VectorizedArrayType * gradients_quad,
8172 VectorizedArrayType * scratch_data,
8181 if (std::is_same<typename VectorType::value_type, Number>::value &&
8185 interleaved_contiguous &&
8186 reinterpret_cast<std::size_t
>(
8187 destination.begin() +
8190 [cell * VectorizedArrayType::size()]) %
8191 sizeof(VectorizedArrayType) ==
8194 VectorizedArrayType *vec_values =
8195 reinterpret_cast<VectorizedArrayType *
>(
8196 destination.begin() +
8199 [cell * VectorizedArrayType::size()] +
8201 [first_selected_component] *
8202 VectorizedArrayType::size());
8237 typename VectorizedArrayType,
8238 typename VectorType,
8239 typename std::enable_if<
8240 !internal::has_begin<VectorType>::value ||
8241 !std::is_same<decltype(std::declval<VectorType>().begin()),
8243 VectorType>::type * =
nullptr>
8245 try_integrate_scatter_inplace(
8252 const VectorizedArrayType *,
8253 const VectorizedArrayType *,
8254 const VectorizedArrayType *,
8269 typename VectorizedArrayType>
8270template <
typename VectorType>
8277 VectorizedArrayType>::
8278 gather_evaluate(
const VectorType & input_vector,
8281 if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
8285 this->active_fe_index,
8286 this->first_selected_component,
8288 evaluation_flag) ==
false)
8290 this->read_dof_values(input_vector);
8291 evaluate(this->begin_dof_values(), evaluation_flag);
8302 typename VectorizedArrayType>
8309 VectorizedArrayType>::integrate(
const bool integrate_values,
8310 const bool integrate_gradients)
8312 integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8315 this->dof_values_initialized =
true;
8326 typename VectorizedArrayType>
8333 VectorizedArrayType>::
8336 integrate(integration_flag, this->values_dofs[0]);
8339 this->dof_values_initialized =
true;
8350 typename VectorizedArrayType>
8357 VectorizedArrayType>::integrate(
const bool integrate_values,
8358 const bool integrate_gradients,
8359 VectorizedArrayType *values_array)
8365 integrate(flag, values_array);
8375 typename VectorizedArrayType>
8382 VectorizedArrayType>::
8384 VectorizedArrayType * values_array)
8388 Assert(this->values_quad_submitted ==
true,
8391 Assert(this->gradients_quad_submitted ==
true,
8394 Assert(this->matrix_info !=
nullptr ||
8395 this->mapped_geometry->is_initialized(),
8402 "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8411 this->gradients_quad,
8421 this->gradients_quad,
8426 this->dof_values_initialized =
true;
8437 typename VectorizedArrayType>
8438template <
typename VectorType>
8446 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8447 const bool integrate_gradients,
8448 VectorType &destination)
8455 integrate_scatter(flag, destination);
8465 typename VectorizedArrayType>
8466template <
typename VectorType>
8473 VectorizedArrayType>::
8475 VectorType & destination)
8477 if (internal::try_integrate_scatter_inplace<dim,
8481 VectorizedArrayType>(
8485 this->active_fe_index,
8486 this->first_selected_component,
8489 this->gradients_quad,
8492 integration_flag) ==
false)
8494 integrate(integration_flag, this->begin_dof_values());
8495 this->distribute_local_to_global(destination);
8510 typename VectorizedArrayType>
8516 VectorizedArrayType>::
8519 const bool is_interior_face,
8520 const unsigned int dof_no,
8521 const unsigned int quad_no,
8522 const unsigned int first_selected_component,
8523 const unsigned int active_fe_index,
8524 const unsigned int active_quad_index,
8525 const unsigned int face_type)
8526 : BaseClass(matrix_free,
8528 first_selected_component,
8536 , dofs_per_component(this->data->dofs_per_component_on_cell)
8537 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
8538 , n_q_points(this->n_quadrature_points)
8548 typename VectorizedArrayType>
8554 VectorizedArrayType>::
8557 const std::pair<unsigned int, unsigned int> & range,
8558 const bool is_interior_face,
8559 const unsigned int dof_no,
8560 const unsigned int quad_no,
8561 const unsigned int first_selected_component)
8566 first_selected_component,
8567 matrix_free.get_face_active_fe_index(range,
8570 matrix_free.get_face_info(range.
first).face_type)
8580 typename VectorizedArrayType>
8587 VectorizedArrayType>
::reinit(
const unsigned int face_index)
8589 Assert(this->mapped_geometry ==
nullptr,
8590 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8591 " Integer indexing is not possible"));
8592 if (this->mapped_geometry !=
nullptr)
8595 this->cell = face_index;
8596 this->dof_access_index =
8597 this->is_interior_face ?
8602 VectorizedArrayType::size()> &faces =
8607 this->matrix_info->get_task_info().boundary_partition_data.back())
8608 Assert(this->is_interior_face,
8610 "Boundary faces do not have a neighbor. When looping over "
8611 "boundary faces use FEFaceEvaluation with the parameter "
8612 "is_interior_face set to true. "));
8616 this->subface_index = this->is_interior_face ==
true ?
8625 this->face_orientation =
8630 this->cell_type = this->matrix_info->
get_mapping_info().face_type[face_index];
8631 const unsigned int offsets =
8632 this->mapping_data->data_index_offsets[face_index];
8633 this->J_value = &this->mapping_data->JxW_values[offsets];
8634 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8636 &this->mapping_data->jacobians[!this->is_interior_face][offsets];
8637 this->normal_x_jacobian =
8639 ->normals_times_jacobians[!this->is_interior_face][offsets];
8642 this->dof_values_initialized =
false;
8643 this->values_quad_initialized =
false;
8644 this->gradients_quad_initialized =
false;
8645 this->hessians_quad_initialized =
false;
8656 typename VectorizedArrayType>
8664 const unsigned int face_number)
8670 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8674 Assert(this->mapped_geometry ==
nullptr,
8675 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8676 " Integer indexing is not possible"));
8677 if (this->mapped_geometry !=
nullptr)
8683 this->face_orientation = 0;
8685 this->face_no = face_number;
8686 this->dof_access_index =
8689 const unsigned int offsets =
8691 .face_data_by_cells[this->quad_no]
8696 .face_data_by_cells[this->quad_no]
8697 .JxW_values.size());
8699 .face_data_by_cells[this->quad_no]
8700 .JxW_values[offsets];
8702 .face_data_by_cells[this->quad_no]
8703 .normal_vectors[offsets];
8705 .face_data_by_cells[this->quad_no]
8706 .jacobians[!this->is_interior_face][offsets];
8707 this->normal_x_jacobian =
8709 .face_data_by_cells[this->quad_no]
8710 .normals_times_jacobians[!this->is_interior_face][offsets];
8713 this->dof_values_initialized =
false;
8714 this->values_quad_initialized =
false;
8715 this->gradients_quad_initialized =
false;
8716 this->hessians_quad_initialized =
false;
8727 typename VectorizedArrayType>
8734 VectorizedArrayType>::evaluate(
const bool evaluate_values,
8735 const bool evaluate_gradients)
8741 evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients);
8751 typename VectorizedArrayType>
8758 VectorizedArrayType>::
8765 evaluate(this->values_dofs[0], evaluation_flag);
8775 typename VectorizedArrayType>
8782 VectorizedArrayType>::evaluate(
const VectorizedArrayType
8784 const bool evaluate_values,
8785 const bool evaluate_gradients)
8792 evaluate(values_array, flag);
8802 typename VectorizedArrayType>
8809 VectorizedArrayType>::
8810 evaluate(
const VectorizedArrayType * values_array,
8817 "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8823 if (this->dof_access_index ==
8825 this->is_interior_face ==
false)
8827 const auto face_nos = this->compute_face_no_data();
8828 const auto face_orientations = this->compute_face_orientations();
8833 for (
unsigned int v = 1; v < VectorizedArrayType::size(); ++v)
8843 template run<fe_degree, n_q_points_1d>(
8847 this->begin_values(),
8848 this->begin_gradients(),
8853 this->subface_index,
8854 face_orientations[0],
8855 this->descriptor->face_orientations);
8861 VectorizedArrayType>::
8862 template run<fe_degree, n_q_points_1d>(
8866 this->begin_values(),
8867 this->begin_gradients(),
8872 this->subface_index,
8873 this->face_orientation,
8874 this->descriptor->face_orientations);
8880 this->begin_values(),
8881 this->begin_gradients(),
8886 this->subface_index,
8887 this->face_orientation,
8888 this->descriptor->face_orientations);
8893 this->values_quad_initialized =
true;
8895 this->gradients_quad_initialized =
true;
8906 typename VectorizedArrayType>
8913 VectorizedArrayType>::
8916 integrate(evaluation_flag, this->values_dofs[0]);
8919 this->dof_values_initialized =
true;
8930 typename VectorizedArrayType>
8937 VectorizedArrayType>::integrate(
const bool integrate_values,
8938 const bool integrate_gradients)
8940 integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8943 this->dof_values_initialized =
true;
8954 typename VectorizedArrayType>
8961 VectorizedArrayType>::integrate(
const bool integrate_values,
8962 const bool integrate_gradients,
8971 integrate(flag, values_array);
8981 typename VectorizedArrayType>
8988 VectorizedArrayType>::
8990 VectorizedArrayType * values_array)
8996 "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
9004 template run<fe_degree, n_q_points_1d>(
9008 this->begin_values(),
9009 this->begin_gradients(),
9014 this->subface_index,
9015 this->face_orientation,
9016 this->descriptor->face_orientations);
9022 this->begin_values(),
9023 this->begin_gradients(),
9028 this->subface_index,
9029 this->face_orientation,
9030 this->descriptor->face_orientations);
9040 typename VectorizedArrayType>
9041template <
typename VectorType>
9049 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
9050 const bool evaluate_values,
9051 const bool evaluate_gradients)
9058 gather_evaluate(input_vector, flag);
9068 typename VectorizedArrayType>
9069template <
typename VectorType>
9076 VectorizedArrayType>::
9077 gather_evaluate(
const VectorType & input_vector,
9084 "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
9086 const auto fu = [&]() {
9087 const auto shared_vector_data = internal::get_shared_vector_data(
9089 this->dof_access_index ==
9091 this->active_fe_index,
9094 if (this->dof_access_index ==
9096 this->is_interior_face ==
false)
9098 const auto cells = this->get_cell_or_face_ids();
9099 const auto face_nos = this->compute_face_no_data();
9100 const auto face_orientations = this->compute_face_orientations();
9105 VectorizedArrayType>::template
run<fe_degree,
9108 VectorizedArrayType::size(),
9109 internal::get_beginning<Number>(input_vector),
9113 this->begin_values(),
9114 this->begin_gradients(),
9118 this->active_fe_index,
9119 this->first_selected_component,
9122 this->subface_index,
9123 this->dof_access_index,
9125 this->descriptor->face_orientations);
9131 std::array<
unsigned int, VectorizedArrayType::size()> cells_ = {};
9132 std::array<
unsigned int, VectorizedArrayType::size()> face_no_ = {};
9133 std::array<
unsigned int, VectorizedArrayType::size()>
9134 face_orientation_ = {};
9136 cells_[0] = this->cell;
9137 face_no_[0] = this->face_no;
9138 face_orientation_[0] = this->face_orientation;
9145 VectorizedArrayType>::template
run<fe_degree,
9149 internal::get_beginning<Number>(input_vector),
9153 this->begin_values(),
9154 this->begin_gradients(),
9158 this->active_fe_index,
9159 this->first_selected_component,
9162 this->subface_index,
9163 this->dof_access_index,
9165 this->descriptor->face_orientations);
9172 internal::get_beginning<Number>(input_vector),
9176 this->begin_values(),
9177 this->begin_gradients(),
9181 this->active_fe_index,
9182 this->first_selected_component,
9185 this->subface_index,
9186 this->dof_access_index,
9188 this->descriptor->face_orientations);
9194 this->read_dof_values(input_vector);
9195 this->evaluate(evaluation_flag);
9200 this->values_quad_initialized =
true;
9202 this->gradients_quad_initialized =
true;
9213 typename VectorizedArrayType>
9214template <
typename VectorType>
9222 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
9223 const bool integrate_gradients,
9224 VectorType &destination)
9231 integrate_scatter(flag, destination);
9241 typename VectorizedArrayType>
9242template <
typename VectorType>
9249 VectorizedArrayType>::
9251 VectorType & destination)
9253 Assert((this->dof_access_index ==
9255 this->is_interior_face ==
false) ==
false,
9258 const auto shared_vector_data = internal::get_shared_vector_data(
9260 this->dof_access_index ==
9262 this->active_fe_index,
9267 std::array<
unsigned int, VectorizedArrayType::size()> cells_ = {};
9268 std::array<
unsigned int, VectorizedArrayType::size()> face_no_ = {};
9269 std::array<
unsigned int, VectorizedArrayType::size()> face_orientation_ = {};
9271 cells_[0] = this->cell;
9272 face_no_[0] = this->face_no;
9273 face_orientation_[0] = this->face_orientation;
9280 VectorizedArrayType>::template
run<fe_degree,
9284 internal::get_beginning<Number>(destination),
9288 this->begin_dof_values(),
9289 this->begin_values(),
9290 this->begin_gradients(),
9294 this->active_fe_index,
9295 this->first_selected_component,
9298 this->subface_index,
9299 this->dof_access_index,
9301 this->descriptor->face_orientations))
9307 this->distribute_local_to_global(destination);
9313 integrate_scatter(n_components,
9315 internal::get_beginning<Number>(destination),
9319 this->begin_dof_values(),
9320 this->begin_values(),
9321 this->begin_gradients(),
9325 this->active_fe_index,
9326 this->first_selected_component,
9329 this->subface_index,
9330 this->dof_access_index,
9332 this->descriptor->face_orientations))
9334 this->distribute_local_to_global(destination);
9346 typename VectorizedArrayType>
9353 VectorizedArrayType>::quadrature_point(
const unsigned int q)
9357 if (this->dof_access_index < 2)
9359 Assert(this->mapping_data->quadrature_point_offsets.empty() ==
false,
9361 "update_quadrature_points"));
9363 this->mapping_data->quadrature_point_offsets.size());
9364 return this->mapping_data->quadrature_points
9365 [this->mapping_data->quadrature_point_offsets[this->cell] + q];
9370 .face_data_by_cells[this->quad_no]
9371 .quadrature_point_offsets.empty() ==
false,
9373 "update_quadrature_points"));
9374 const unsigned int index =
9378 .face_data_by_cells[this->quad_no]
9379 .quadrature_point_offsets.size());
9381 .face_data_by_cells[this->quad_no]
9383 .face_data_by_cells[this->quad_no]
9384 .quadrature_point_offsets[index] +
9396 typename VectorizedArrayType>
9397std::array<
unsigned int, VectorizedArrayType::size()>
9403 VectorizedArrayType>::compute_face_no_data()
9405 std::array<
unsigned int, VectorizedArrayType::size()> face_no_data;
9407 if (this->dof_access_index !=
9409 this->is_interior_face ==
true)
9411 std::fill(face_no_data.begin(),
9412 face_no_data.begin() +
9413 this->dof_info->n_vectorization_lanes_filled
9420 std::fill(face_no_data.begin(),
9424 for (
unsigned int i = 0;
9431 const unsigned int cell_this =
9432 this->cell * VectorizedArrayType::size() + i;
9434 const unsigned int face_index =
9446 .cells_interior[face_index % VectorizedArrayType::size()];
9450 (cell_m == cell_this) ?
9452 ->get_face_info(face_index / VectorizedArrayType::size())
9455 ->get_face_info(face_index / VectorizedArrayType::size())
9460 return face_no_data;
9470 typename VectorizedArrayType>
9471std::array<
unsigned int, VectorizedArrayType::size()>
9477 VectorizedArrayType>::compute_face_orientations()
9479 std::array<
unsigned int, VectorizedArrayType::size()> face_no_data;
9481 if (this->dof_access_index !=
9483 this->is_interior_face ==
true)
9489 std::fill(face_no_data.begin(),
9495 for (
unsigned int i = 0;
9502 const unsigned int cell_this =
9503 this->cell * VectorizedArrayType::size() + i;
9505 const unsigned int face_index =
9507 this->cell, this->face_no, i);
9512 const unsigned int macro =
9513 face_index / VectorizedArrayType::size();
9514 const unsigned int lane =
9515 face_index % VectorizedArrayType::size();
9517 const auto &faces = this->matrix_info->
get_face_info(macro);
9522 const bool is_interior_face = cell_m != cell_this;
9527 if (is_interior_face != fo_interior_face)
9531 static const std::array<unsigned int, 8> table{
9532 {0, 1, 0, 3, 6, 5, 4, 7}};
9534 face_orientation = table[face_orientation];
9538 face_no_data[i] = face_orientation;
9544 face_no_data.begin(),
9545 face_no_data.begin() +
9546 this->dof_info->n_vectorization_lanes_filled
9553 return face_no_data;
9563 typename VectorizedArrayType>
9570 VectorizedArrayType>::
9571 fast_evaluation_supported(
const unsigned int given_degree,
9572 const unsigned int give_n_q_points_1d)
9574 return fe_degree == -1 ?
9587 typename VectorizedArrayType>
9594 VectorizedArrayType>::
9595 fast_evaluation_supported(
const unsigned int given_degree,
9596 const unsigned int give_n_q_points_1d)
9598 return fe_degree == -1 ?
FEEvaluationAccess(const FEEvaluationAccess &other)
value_type get_value(const unsigned int q_point) const
value_type integrate_value() const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
Tensor< 2, 1, VectorizedArrayType > get_hessian(unsigned int q_point) const
value_type get_laplacian(const unsigned int q_point) const
void submit_value(const gradient_type val_in, const unsigned int q_point)
value_type get_normal_derivative(const unsigned int q_point) const
void submit_normal_derivative(const gradient_type grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
FEEvaluationAccess(const Mapping< 1 > &mapping, const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< 1, Number, is_face, VectorizedArrayType > *other)
void submit_gradient(const value_type grad_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
VectorizedArrayType value_type
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > grad_in, const unsigned int q_point)
value_type get_divergence(const unsigned int q_point) const
value_type get_dof_value(const unsigned int dof) const
FEEvaluationAccess(const MatrixFree< 1, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< dim, Number, is_face, VectorizedArrayType > *other)
value_type get_value(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
Tensor< 2, dim, VectorizedArrayType > get_hessian(unsigned int q_point) const
value_type get_dof_value(const unsigned int dof) const
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
FEEvaluationAccess(const FEEvaluationAccess &other)
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
value_type get_laplacian(const unsigned int q_point) const
VectorizedArrayType value_type
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
value_type get_normal_derivative(const unsigned int q_point) const
value_type integrate_value() const
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const FEEvaluationAccess &other)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
Tensor< 3, dim, VectorizedArrayType > get_hessian(const unsigned int q_point) const
VectorizedArrayType get_divergence(const unsigned int q_point) const
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 1, dim, Tensor< 1, dim, VectorizedArrayType > > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< dim, Number, is_face, VectorizedArrayType > *other)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
Tensor< 1, n_components_, VectorizedArrayType > value_type
static constexpr unsigned int dimension
Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > gradient_type
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< dim, Number, is_face, VectorizedArrayType > *other)
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const FEEvaluationAccess &other)
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
static constexpr unsigned int n_components
std::array< T, VectorizedArrayType::size()> read_cell_data(const AlignedVector< std::array< T, VectorizedArrayType::size()> > &array) const
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
Tensor< 1, dim, VectorizedArrayType > get_normal_vector(const unsigned int q_point) const
FEEvaluationBaseData & operator=(const FEEvaluationBaseData &other)
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
const internal::MatrixFreeFunctions::MappingInfoStorage<(is_face?dim-1:dim), dim, Number, VectorizedArrayType >::QuadratureDescriptor * descriptor
void set_cell_data(AlignedVector< VectorizedArrayType > &array, const VectorizedArrayType &value) const
const unsigned int active_quad_index
unsigned int subface_index
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info() const
FEEvaluationBaseData(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
const Tensor< 1, dim, VectorizedArrayType > * normal_x_jacobian
internal::MatrixFreeFunctions::GeometryType get_cell_type() const
const std::vector< unsigned int > & get_internal_dof_numbering() const
std::array< unsigned int, VectorizedArrayType::size()> get_cell_ids() const
Tensor< 2, dim, VectorizedArrayType > inverse_jacobian(const unsigned int q_point) const
unsigned int get_mapping_data_index_offset() const
const Tensor< 2, dim, VectorizedArrayType > * jacobian
unsigned int face_orientation
internal::MatrixFreeFunctions::GeometryType cell_type
const unsigned int quad_no
VectorizedArrayType JxW(const unsigned int q_point) const
const Number * quadrature_weights
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > * data
unsigned int get_current_cell_index() const
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number, VectorizedArrayType > > mapped_geometry
unsigned int get_active_quadrature_index() const
const unsigned int n_quadrature_points
static constexpr unsigned int dimension
void set_cell_data(AlignedVector< std::array< T, VectorizedArrayType::size()> > &array, const std::array< T, VectorizedArrayType::size()> &value) const
VectorizedArrayType read_cell_data(const AlignedVector< VectorizedArrayType > &array) const
const Tensor< 1, dim, VectorizedArrayType > * normal_vectors
const internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 :dim), dim, Number, VectorizedArrayType > * mapping_data
FEEvaluationBaseData(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< dim, Number, is_face, VectorizedArrayType > *other)
unsigned int get_quadrature_index() const
FEEvaluationBaseData(const FEEvaluationBaseData &other)
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
const VectorizedArrayType * J_value
AlignedVector< VectorizedArrayType > * scratch_data_array
const internal::MatrixFreeFunctions::DoFInfo * dof_info
ArrayView< VectorizedArrayType > get_scratch_data() const
std::array< unsigned int, VectorizedArrayType::size()> get_cell_or_face_ids() const
unsigned int get_active_fe_index() const
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_info
VectorizedArrayType * scratch_data
const unsigned int active_fe_index
VectorizedArrayType * gradients_quad
value_type get_dof_value(const unsigned int dof) const
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > get_hessian(const unsigned int q_point) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
gradient_type get_gradient(const unsigned int q_point) const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
const VectorizedArrayType * begin_hessians() const
static constexpr unsigned int dimension
bool values_quad_initialized
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
VectorizedArrayType * begin_values()
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
const unsigned int first_selected_component
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
bool dof_values_initialized
bool hessians_quad_initialized
const VectorizedArrayType * begin_dof_values() const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0)
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
VectorizedArrayType * begin_hessians()
unsigned int get_first_selected_component() const
const VectorizedArrayType * begin_values() const
bool gradients_quad_submitted
VectorizedArrayType * begin_dof_values()
VectorizedArrayType get_divergence(const unsigned int q_point) const
VectorizedArrayType * begin_gradients()
FEEvaluationBase & operator=(const FEEvaluationBase &other)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
VectorizedArrayType * values_quad
static constexpr unsigned int n_components
gradient_type get_hessian_diagonal(const unsigned int q_point) const
VectorizedArrayType * values_dofs[n_components]
const VectorizedArrayType * begin_gradients() const
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask, const bool apply_constraints=true) const
VectorizedArrayType * hessians_quad
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationBaseData< dim, Number, is_face, VectorizedArrayType > *other)
const unsigned int n_fe_components
bool values_quad_submitted
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
bool gradients_quad_initialized
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
Point< dim, VectorizedArrayType > quadrature_point(const unsigned int q_point) const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
void integrate(const bool integrate_values, const bool integrate_gradients)
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
void integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationBaseData< dim, Number, false, VectorizedArrayType > &other, const unsigned int first_selected_component=0)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
FEEvaluation(const FEEvaluation &other)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
static constexpr unsigned int static_n_q_points_cell
static constexpr unsigned int tensor_dofs_per_cell
std::array< unsigned int, VectorizedArrayType::size()> compute_face_no_data()
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
static constexpr unsigned int static_dofs_per_component
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
static constexpr unsigned int static_n_q_points
std::array< unsigned int, VectorizedArrayType::size()> compute_face_orientations()
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags evaluation_flag)
Point< dim, VectorizedArrayType > quadrature_point(const unsigned int q_point) const
void integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag, VectorType &output_vector)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate(const EvaluationFlags::EvaluationFlags evaluation_flag, VectorizedArrayType *values_array)
void integrate(const bool integrate_values, const bool integrate_gradients)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
@ update_quadrature_points
Transformed quadrature points.
#define DeclException0(Exception0)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
#define Assert(cond, exc)
std::string to_string(const T &t)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
EvaluationFlags
The EvaluationFlags enum.
@ general
No special properties.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
constexpr T pow(const T base, const int iexp)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
void check_vector_compatibility(const VectorType &vec, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
static const unsigned int invalid_unsigned_int
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *hessians_quad, Number *scratch_data)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const internal::MatrixFreeFunctions::ShapeInfo< Number > &shape_info, Number *values_dofs_actual, Number *values_quad, Number *gradients_quad, Number *scratch_data, const bool sum_into_values_array=false)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, VectorizedArrayType *values_dofs_actual, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *hessians_quad, VectorizedArrayType *scratch_data)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, VectorizedArrayType *values_dofs_actual, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool sum_into_values_array)
static bool gather_evaluate(const unsigned int n_components, const std::size_t n_face_orientations, const Number *src_ptr, const std::vector< ArrayView< const Number > > *sm_ptr, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const MatrixFreeFunctions::DoFInfo &dof_info, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int active_fe_index, const unsigned int first_selected_component, const std::array< unsigned int, VectorizedArrayType::size()> cells, const std::array< unsigned int, VectorizedArrayType::size()> face_nos, const unsigned int subface_index, const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index, const std::array< unsigned int, VectorizedArrayType::size()> face_orientations, const Table< 2, unsigned int > &orientation_map)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, const VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool evaluate_values, const bool evaluate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static void integrate(const unsigned int n_components, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &data, VectorizedArrayType *values_array, VectorizedArrayType *values_quad, VectorizedArrayType *gradients_quad, VectorizedArrayType *scratch_data, const bool integrate_values, const bool integrate_gradients, const unsigned int face_no, const unsigned int subface_index, const unsigned int face_orientation, const Table< 2, unsigned int > &orientation_map)
static constexpr unsigned int value
@ dof_access_face_exterior
@ dof_access_face_interior
std::vector< std::vector< unsigned int > > component_dof_indices_offset
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
unsigned char subface_index
unsigned char interior_face_no
std::array< unsigned int, vectorization_width > cells_interior
types::boundary_id exterior_face_no
unsigned char face_orientation
std::vector< unsigned int > cell_partition_data
std::vector< unsigned int > face_partition_data