17 #ifndef dealii_matrix_free_fe_evaluation_h
18 #define dealii_matrix_free_fe_evaluation_h
45 #include <type_traits>
93 typename VectorizedArrayType>
106 static constexpr
unsigned int n_lanes = VectorizedArrayType::size();
144 template <
typename VectorType>
147 const VectorType &src,
148 const unsigned int first_index = 0,
149 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
179 template <
typename VectorType>
182 const VectorType &src,
183 const unsigned int first_index = 0,
184 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
217 template <
typename VectorType>
221 const unsigned int first_index = 0,
222 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
262 template <
typename VectorType>
266 const unsigned int first_index = 0,
267 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
272 template <
typename VectorType>
276 const unsigned int first_index = 0,
277 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
412 const unsigned int q_point);
500 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
520 const unsigned int q_point);
541 const unsigned int q_point);
557 const unsigned int q_point);
599 const unsigned int dof_no,
602 const unsigned int fe_degree,
603 const unsigned int n_q_points,
607 const unsigned int face_type);
681 template <
typename VectorType,
typename VectorOperation>
685 const std::array<VectorType *, n_components_> &vectors,
688 n_components_> &vectors_sm,
689 const std::bitset<n_lanes> &mask,
690 const bool apply_constraints =
true)
const;
699 template <
typename VectorType,
typename VectorOperation>
703 const std::array<VectorType *, n_components_> &vectors,
706 n_components_> &vectors_sm,
707 const std::bitset<n_lanes> &mask)
const;
716 template <
typename VectorType,
typename VectorOperation>
720 const std::array<VectorType *, n_components_> &vectors)
const;
766 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
767 "Type of Number and of VectorizedArrayType do not match.");
789 const unsigned int dof_no,
792 const unsigned int fe_degree,
793 const unsigned int n_q_points,
833 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
838 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
839 "Type of Number and of VectorizedArrayType do not match.");
879 const unsigned int q_point);
904 const unsigned int q_point);
946 const unsigned int dof_no,
949 const unsigned int fe_degree,
950 const unsigned int n_q_points,
991 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
996 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
997 "Type of Number and of VectorizedArrayType do not match.");
1040 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1060 const unsigned int q_point);
1079 const unsigned int q_point);
1091 const unsigned int q_point);
1104 const unsigned int q_point);
1112 const unsigned int q_point);
1124 const unsigned int dof_no,
1127 const unsigned int dofs_per_cell,
1128 const unsigned int n_q_points,
1167 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
1172 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1173 "Type of Number and of VectorizedArrayType do not match.");
1249 const unsigned int q_point);
1256 const unsigned int q_point);
1263 const unsigned int q_point);
1305 const unsigned int dof_no,
1308 const unsigned int fe_degree,
1309 const unsigned int n_q_points,
1909 typename VectorizedArrayType>
1914 VectorizedArrayType>
1917 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1918 "Type of Number and of VectorizedArrayType do not match.");
1960 static constexpr
unsigned int n_lanes = VectorizedArrayType::size();
2045 const unsigned int dof_no = 0,
2046 const unsigned int quad_no = 0,
2059 const std::pair<unsigned int, unsigned int> &range,
2060 const unsigned int dof_no = 0,
2061 const unsigned int quad_no = 0,
2171 template <
bool level_dof_access>
2193 const unsigned int given_n_q_points_1d);
2236 template <
typename VectorType>
2267 VectorizedArrayType *values_array,
2268 const bool sum_into_values =
false);
2283 template <
typename VectorType>
2286 VectorType &output_vector);
2370 int n_q_points_1d = fe_degree + 1,
2371 int n_components_ = 1,
2372 typename Number =
double,
2378 VectorizedArrayType>
2381 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
2382 "Type of Number and of VectorizedArrayType do not match.");
2424 static constexpr
unsigned int n_lanes = VectorizedArrayType::size();
2526 const unsigned int dof_no = 0,
2527 const unsigned int quad_no = 0,
2542 const std::pair<unsigned int, unsigned int> &range,
2544 const unsigned int dof_no = 0,
2545 const unsigned int quad_no = 0,
2559 reinit(
const unsigned int face_batch_number);
2569 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2576 const unsigned int given_n_q_points_1d);
2618 template <
typename VectorType>
2645 VectorizedArrayType *values_array);
2658 template <
typename VectorType>
2661 VectorType &output_vector);
2666 template <
typename VectorType>
2669 const bool integrate_gradients,
2670 VectorType &output_vector);
2732 namespace MatrixFreeFunctions
2736 template <
int dim,
int degree>
2745 template <
int degree>
2748 static constexpr
unsigned int value = degree + 1;
2764 template <
bool is_face,
2767 typename VectorizedArrayType>
2770 extract_initialization_data(
2772 const unsigned int dof_no,
2773 const unsigned int first_selected_component,
2774 const unsigned int quad_no,
2775 const unsigned int fe_degree,
2776 const unsigned int n_q_points,
2777 const unsigned int active_fe_index_given,
2778 const unsigned int active_quad_index_given,
2779 const unsigned int face_type)
2782 InitializationData init_data;
2786 &internal::MatrixFreeFunctions::
2787 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2795 active_fe_index_given :
2800 active_quad_index_given :
2831 typename VectorizedArrayType>
2836 VectorizedArrayType>::
2839 const unsigned int dof_no,
2840 const unsigned int first_selected_component,
2841 const unsigned int quad_no,
2842 const unsigned int fe_degree,
2843 const unsigned int n_q_points,
2844 const bool is_interior_face,
2845 const unsigned int active_fe_index,
2846 const unsigned int active_quad_index,
2847 const unsigned int face_type)
2849 internal::extract_initialization_data<is_face>(matrix_free,
2851 first_selected_component,
2860 first_selected_component)
2861 , scratch_data_array(matrix_free.acquire_scratch_data())
2862 , matrix_free(&matrix_free)
2864 this->set_data_pointers(scratch_data_array, n_components_);
2866 this->dof_info->start_components.back() == 1 ||
2867 static_cast<int>(n_components_) <=
2869 this->dof_info->start_components
2870 [this->dof_info->component_to_base_index[first_selected_component] +
2872 first_selected_component,
2874 "You tried to construct a vector-valued evaluator with " +
2876 " components. However, "
2877 "the current base element has only " +
2879 this->dof_info->start_components
2880 [this->dof_info->component_to_base_index[first_selected_component] +
2882 first_selected_component) +
2883 " components left when starting from local element index " +
2885 first_selected_component -
2886 this->dof_info->start_components
2887 [this->dof_info->component_to_base_index[first_selected_component]]) +
2888 " (global index " +
std::to_string(first_selected_component) +
")"));
2900 typename VectorizedArrayType>
2905 VectorizedArrayType>::
2911 const unsigned int first_selected_component,
2915 other->mapped_geometry->get_quadrature() == quadrature ?
2916 other->mapped_geometry :
2917 std::make_shared<
internal::MatrixFreeFunctions::
2918 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2923 first_selected_component)
2924 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2925 , matrix_free(nullptr)
2927 const unsigned int base_element_number =
2931 first_selected_component >=
2933 ExcMessage(
"The underlying element must at least contain as many "
2934 "components as requested by this class"));
2935 (void)base_element_number;
2940 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2944 this->set_data_pointers(scratch_data_array, n_components_);
2953 typename VectorizedArrayType>
2958 VectorizedArrayType>::
2963 VectorizedArrayType> &other)
2965 , scratch_data_array(other.matrix_free == nullptr ?
2967 other.matrix_free->acquire_scratch_data())
2968 , matrix_free(other.matrix_free)
2970 if (other.matrix_free ==
nullptr)
2978 this->mapped_geometry =
2979 std::make_shared<internal::MatrixFreeFunctions::
2980 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2984 this->mapping_data = &this->mapped_geometry->get_data_storage();
2988 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2990 this->mapped_geometry->get_data_storage().JxW_values.begin();
2991 this->jacobian_gradients =
2992 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2993 this->jacobian_gradients_non_inverse =
2994 this->mapped_geometry->get_data_storage()
2995 .jacobian_gradients_non_inverse[0]
2998 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3001 this->set_data_pointers(scratch_data_array, n_components_);
3010 typename VectorizedArrayType>
3015 VectorizedArrayType> &
3021 VectorizedArrayType> &other)
3024 if (matrix_free ==
nullptr)
3027 delete scratch_data_array;
3036 matrix_free = other.matrix_free;
3038 if (other.matrix_free ==
nullptr)
3047 this->mapped_geometry =
3048 std::make_shared<internal::MatrixFreeFunctions::
3049 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3054 this->mapping_data = &this->mapped_geometry->get_data_storage();
3056 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3058 this->mapped_geometry->get_data_storage().JxW_values.begin();
3059 this->jacobian_gradients =
3060 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3061 this->jacobian_gradients_non_inverse =
3062 this->mapped_geometry->get_data_storage()
3063 .jacobian_gradients_non_inverse[0]
3066 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3073 this->set_data_pointers(scratch_data_array, n_components_);
3084 typename VectorizedArrayType>
3089 VectorizedArrayType>::~FEEvaluationBase()
3091 if (matrix_free !=
nullptr)
3102 delete scratch_data_array;
3113 typename VectorizedArrayType>
3118 Assert(matrix_free !=
nullptr,
3120 "FEEvaluation was not initialized with a MatrixFree object!"));
3121 return *matrix_free;
3130 template <
typename VectorType,
bool>
3131 struct ConstBlockVectorSelector;
3133 template <
typename VectorType>
3134 struct ConstBlockVectorSelector<VectorType, true>
3136 using BaseVectorType =
const typename VectorType::BlockType;
3139 template <
typename VectorType>
3140 struct ConstBlockVectorSelector<VectorType, false>
3142 using BaseVectorType =
typename VectorType::BlockType;
3148 template <
typename VectorType,
bool>
3149 struct BlockVectorSelector;
3151 template <
typename VectorType>
3152 struct BlockVectorSelector<VectorType, true>
3154 using BaseVectorType =
typename ConstBlockVectorSelector<
3156 std::is_const_v<VectorType>>::BaseVectorType;
3158 static BaseVectorType *
3159 get_vector_component(VectorType &vec,
const unsigned int component)
3162 return &vec.block(component);
3166 template <
typename VectorType>
3167 struct BlockVectorSelector<VectorType, false>
3169 using BaseVectorType = VectorType;
3171 static BaseVectorType *
3172 get_vector_component(VectorType &vec,
const unsigned int component)
3191 template <
typename VectorType>
3192 struct BlockVectorSelector<std::vector<VectorType>, false>
3194 using BaseVectorType = VectorType;
3196 static BaseVectorType *
3197 get_vector_component(std::vector<VectorType> &vec,
3198 const unsigned int component)
3201 return &vec[component];
3205 template <
typename VectorType>
3206 struct BlockVectorSelector<const std::vector<VectorType>, false>
3208 using BaseVectorType =
const VectorType;
3210 static const BaseVectorType *
3211 get_vector_component(
const std::vector<VectorType> &vec,
3212 const unsigned int component)
3215 return &vec[component];
3219 template <
typename VectorType>
3220 struct BlockVectorSelector<std::vector<VectorType *>, false>
3222 using BaseVectorType = VectorType;
3224 static BaseVectorType *
3225 get_vector_component(std::vector<VectorType *> &vec,
3226 const unsigned int component)
3229 return vec[component];
3233 template <
typename VectorType>
3234 struct BlockVectorSelector<const std::vector<VectorType *>, false>
3236 using BaseVectorType =
const VectorType;
3238 static const BaseVectorType *
3239 get_vector_component(
const std::vector<VectorType *> &vec,
3240 const unsigned int component)
3243 return vec[component];
3254 typename VectorizedArrayType>
3255 template <
typename VectorType,
typename VectorOperation>
3260 const std::array<VectorType *, n_components_> &src,
3263 n_components_> &src_sm,
3264 const std::bitset<n_lanes> &mask,
3265 const bool apply_constraints)
const
3270 if (this->matrix_free ==
nullptr)
3272 read_write_operation_global(operation, src);
3279 if (this->n_fe_components == 1)
3280 for (
unsigned int comp = 0; comp < n_components; ++comp)
3282 Assert(src[comp] !=
nullptr,
3283 ExcMessage(
"The finite element underlying this FEEvaluation "
3284 "object is scalar, but you requested " +
3286 " components via the template argument in "
3287 "FEEvaluation. In that case, you must pass an "
3288 "std::vector<VectorType> or a BlockVector to " +
3289 "read_dof_values and distribute_local_to_global."));
3301 const bool accesses_exterior_dofs =
3302 this->dof_access_index ==
3304 this->is_interior_face() ==
false;
3314 bool is_contiguous =
true;
3316 if (accesses_exterior_dofs)
3318 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3319 const unsigned int n_filled_lanes =
3324 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3325 if (mask[v] ==
true &&
3328 [cells[v] / n_lanes] <
3331 is_contiguous =
false;
3335 this->dof_access_index :
3338 IndexStorageVariants::contiguous)
3340 is_contiguous =
false;
3345 read_write_operation_contiguous(operation, src, src_sm, mask);
3352 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
3354 const bool masking_is_active =
mask.count() < n_lanes;
3355 if (masking_is_active)
3356 for (
unsigned int v = 0; v < n_lanes; ++v)
3357 if (mask[v] ==
false)
3360 bool has_hn_constraints =
false;
3362 if (is_face ==
false)
3368 [this->first_selected_component])
3369 for (
unsigned int v = 0; v < n_lanes; ++v)
3374 has_hn_constraints =
true;
3377 std::integral_constant<
bool,
3378 internal::is_vectorizable<VectorType, Number>::value>
3381 const bool use_vectorized_path =
3382 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
3384 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3385 std::array<VectorizedArrayType *, n_components> values_dofs;
3386 for (
unsigned int c = 0; c < n_components; ++c)
3387 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3388 c * dofs_per_component;
3392 [is_face ? this->dof_access_index :
3395 IndexStorageVariants::interleaved &&
3396 use_vectorized_path)
3398 const unsigned int *dof_indices =
3400 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
3404 [this->first_selected_component] *
3407 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3408 if (n_components == 1 || this->n_fe_components == 1)
3409 for (
unsigned int comp = 0; comp < n_components; ++comp)
3411 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
3414 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3416 if (n_components == 1 || this->n_fe_components == 1)
3417 for (
unsigned int i = 0; i < dofs_per_component;
3418 ++i, dof_indices += n_lanes)
3419 for (
unsigned int comp = 0; comp < n_components; ++comp)
3420 operation.process_dof_gather(dof_indices,
3424 values_dofs[comp][i],
3427 for (
unsigned int comp = 0; comp < n_components; ++comp)
3428 for (
unsigned int i = 0; i < dofs_per_component;
3429 ++i, dof_indices += n_lanes)
3430 operation.process_dof_gather(dof_indices,
3434 values_dofs[comp][i],
3441 std::array<const unsigned int *, n_lanes> dof_indices;
3442 dof_indices.fill(
nullptr);
3447 bool has_constraints =
false;
3448 const unsigned int n_components_read =
3449 this->n_fe_components > 1 ? n_components : 1;
3453 for (
unsigned int v = 0; v < n_lanes; ++v)
3459 const std::pair<unsigned int, unsigned int> *my_index_start =
3460 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3461 this->first_selected_component];
3466 if (my_index_start[n_components_read].
second !=
3467 my_index_start[0].
second)
3468 has_constraints =
true;
3471 dof_info.
dof_indices.data() + my_index_start[0].first;
3476 for (
unsigned int v = 0; v < n_lanes; ++v)
3481 const std::pair<unsigned int, unsigned int> *my_index_start =
3482 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3483 this->first_selected_component];
3484 if (my_index_start[n_components_read].
second !=
3485 my_index_start[0].
second)
3486 has_constraints =
true;
3494 [this->active_fe_index][this->first_selected_component])
3495 has_hn_constraints =
true;
3498 my_index_start[0].
first ||
3501 my_index_start[0].
first,
3504 dof_info.
dof_indices.data() + my_index_start[0].first;
3508 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3509 return i != numbers::invalid_unsigned_int;
3511 for (
unsigned int comp = 0; comp < n_components; ++comp)
3512 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3513 operation.process_empty(values_dofs[comp][i]);
3517 if (!has_constraints && apply_constraints)
3519 if (n_components == 1 || this->n_fe_components == 1)
3521 for (
unsigned int v = 0; v < n_lanes; ++v)
3526 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3527 for (
unsigned int comp = 0; comp < n_components; ++comp)
3528 operation.process_dof(dof_indices[v][i],
3530 values_dofs[comp][i][v]);
3535 for (
unsigned int comp = 0; comp < n_components; ++comp)
3536 for (
unsigned int v = 0; v < n_lanes; ++v)
3541 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3542 operation.process_dof(
3543 dof_indices[v][comp * dofs_per_component + i],
3545 values_dofs[comp][i][v]);
3557 for (
unsigned int v = 0; v < n_lanes; ++v)
3563 const unsigned int cell_dof_index =
3564 cell_index * this->n_fe_components + this->first_selected_component;
3565 const unsigned int n_components_read =
3566 this->n_fe_components > 1 ? n_components : 1;
3567 unsigned int index_indicators =
3569 unsigned int next_index_indicators =
3570 dof_info.
row_starts[cell_dof_index + 1].second;
3574 if (apply_constraints ==
false &&
3575 (dof_info.
row_starts[cell_dof_index].second !=
3576 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3583 [this->active_fe_index][this->first_selected_component])))
3592 [this->first_selected_component] +
3594 next_index_indicators = index_indicators;
3597 if (n_components == 1 || this->n_fe_components == 1)
3599 unsigned int ind_local = 0;
3600 for (; index_indicators != next_index_indicators; ++index_indicators)
3602 const std::pair<unsigned short, unsigned short> indicator =
3605 for (
unsigned int j = 0; j < indicator.first; ++j)
3606 for (
unsigned int comp = 0; comp < n_components; ++comp)
3607 operation.process_dof(dof_indices[v][j],
3609 values_dofs[comp][ind_local + j][v]);
3611 ind_local += indicator.first;
3612 dof_indices[v] += indicator.first;
3616 Number
value[n_components];
3617 for (
unsigned int comp = 0; comp < n_components; ++comp)
3618 operation.pre_constraints(values_dofs[comp][ind_local][v],
3621 const Number *data_val =
3623 const Number *end_pool =
3625 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3626 for (
unsigned int comp = 0; comp < n_components; ++comp)
3627 operation.process_constraint(*dof_indices[v],
3632 for (
unsigned int comp = 0; comp < n_components; ++comp)
3633 operation.post_constraints(value[comp],
3634 values_dofs[comp][ind_local][v]);
3640 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3641 for (
unsigned int comp = 0; comp < n_components; ++comp)
3642 operation.process_dof(*dof_indices[v],
3644 values_dofs[comp][ind_local][v]);
3653 for (
unsigned int comp = 0; comp < n_components; ++comp)
3655 unsigned int ind_local = 0;
3658 for (; index_indicators != next_index_indicators;
3661 const std::pair<unsigned short, unsigned short> indicator =
3665 for (
unsigned int j = 0; j < indicator.first; ++j)
3666 operation.process_dof(dof_indices[v][j],
3668 values_dofs[comp][ind_local + j][v]);
3669 ind_local += indicator.first;
3670 dof_indices[v] += indicator.first;
3675 operation.pre_constraints(values_dofs[comp][ind_local][v],
3678 const Number *data_val =
3680 const Number *end_pool =
3683 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3684 operation.process_constraint(*dof_indices[v],
3689 operation.post_constraints(value,
3690 values_dofs[comp][ind_local][v]);
3697 for (; ind_local < dofs_per_component;
3698 ++dof_indices[v], ++ind_local)
3701 operation.process_dof(*dof_indices[v],
3703 values_dofs[comp][ind_local][v]);
3706 if (apply_constraints ==
true && comp + 1 < n_components)
3707 next_index_indicators =
3708 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3720 typename VectorizedArrayType>
3721 template <
typename VectorType,
typename VectorOperation>
3726 const std::array<VectorType *, n_components_> &src)
const
3730 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3731 unsigned int index = this->first_selected_component * dofs_per_component;
3732 for (
unsigned int comp = 0; comp < n_components; ++comp)
3734 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3736 operation.process_empty(
3737 this->values_dofs[comp * dofs_per_component + i]);
3738 operation.process_dof_global(
3739 local_dof_indices[this->data->lexicographic_numbering[index]],
3741 this->values_dofs[comp * dofs_per_component + i][0]);
3752 typename VectorizedArrayType>
3753 template <
typename VectorType,
typename VectorOperation>
3758 const std::array<VectorType *, n_components_> &src,
3761 n_components_> &vectors_sm,
3762 const std::bitset<n_lanes> &mask)
const
3771 std::integral_constant<
bool,
3772 internal::is_vectorizable<VectorType, Number>::value>
3775 is_face ? this->dof_access_index :
3777 const unsigned int n_active_lanes =
mask.count();
3780 const std::vector<unsigned int> &dof_indices_cont =
3783 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3784 std::array<VectorizedArrayType *, n_components> values_dofs;
3785 for (
unsigned int c = 0; c < n_components; ++c)
3786 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3787 c * dofs_per_component;
3791 const bool accesses_exterior_dofs =
3792 this->dof_access_index ==
3794 this->is_interior_face() ==
false;
3800 interleaved_contiguous &&
3801 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3803 const unsigned int dof_index =
3804 dof_indices_cont[this->cell * n_lanes] +
3807 [this->first_selected_component] *
3809 if (n_components == 1 || this->n_fe_components == 1)
3810 for (
unsigned int comp = 0; comp < n_components; ++comp)
3811 operation.process_dofs_vectorized(dofs_per_component,
3817 operation.process_dofs_vectorized(dofs_per_component * n_components,
3825 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3829 const unsigned int n_filled_lanes =
3832 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3833 n_active_lanes == n_lanes &&
3834 !accesses_exterior_dofs;
3836 if (vectors_sm[0] !=
nullptr)
3838 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3839 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3841 const auto upper_bound =
3842 std::min<unsigned int>(n_filled_lanes, n_lanes);
3843 for (
unsigned int v = 0; v < upper_bound; ++v)
3845 if (mask[v] ==
false)
3847 vector_ptrs[v] =
nullptr;
3867 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3868 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3870 [this->active_fe_index][this->first_selected_component]);
3872 vector_ptrs[v] =
nullptr;
3874 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3875 vector_ptrs[v] =
nullptr;
3880 if (use_vectorized_path)
3882 if (n_components == 1 || this->n_fe_components == 1)
3884 for (
unsigned int comp = 0; comp < n_components; ++comp)
3886 auto vector_ptrs = compute_vector_ptrs(comp);
3887 operation.process_dofs_vectorized_transpose(
3896 auto vector_ptrs = compute_vector_ptrs(0);
3897 operation.process_dofs_vectorized_transpose(dofs_per_component *
3905 for (
unsigned int comp = 0; comp < n_components; ++comp)
3907 auto vector_ptrs = compute_vector_ptrs(
3908 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3910 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3911 operation.process_empty(values_dofs[comp][i]);
3913 if (n_components == 1 || this->n_fe_components == 1)
3915 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3916 if (mask[v] ==
true)
3917 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3918 operation.process_dof(vector_ptrs[v][i],
3919 values_dofs[comp][i][v]);
3923 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3924 if (mask[v] ==
true)
3925 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3926 operation.process_dof(
3927 vector_ptrs[v][i + comp * dofs_per_component],
3928 values_dofs[comp][i][v]);
3934 std::array<unsigned int, n_lanes> dof_indices;
3935 std::fill(dof_indices.begin(),
3940 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3944 if (mask[v] ==
true)
3946 dof_indices_cont[cells[v]] +
3949 [this->first_selected_component] *
3955 if (use_vectorized_path)
3961 if (n_components == 1 || this->n_fe_components == 1)
3962 for (
unsigned int comp = 0; comp < n_components; ++comp)
3963 operation.process_dofs_vectorized_transpose(dofs_per_component,
3969 operation.process_dofs_vectorized_transpose(dofs_per_component *
3978 interleaved_contiguous_strided)
3980 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3981 if (n_components == 1 || this->n_fe_components == 1)
3982 for (
unsigned int comp = 0; comp < n_components; ++comp)
3983 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3984 src[comp]->
begin());
3987 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3989 if (n_components == 1 || this->n_fe_components == 1)
3990 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3992 for (
unsigned int comp = 0; comp < n_components; ++comp)
3993 operation.process_dof_gather(dof_indices.data(),
3996 src_ptrs[comp] + i * n_lanes,
3997 values_dofs[comp][i],
4001 for (
unsigned int comp = 0; comp < n_components; ++comp)
4002 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4004 operation.process_dof_gather(
4007 (comp * dofs_per_component + i) * n_lanes,
4008 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
4009 values_dofs[comp][i],
4017 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4019 std::array<typename VectorType::value_type *, n_components> src_ptrs;
4020 if (n_components == 1 || this->n_fe_components == 1)
4021 for (
unsigned int comp = 0; comp < n_components; ++comp)
4022 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
4023 src[comp]->
begin());
4026 const_cast<typename VectorType::value_type *
>(src[0]->begin());
4028 const unsigned int *offsets =
4030 if (n_components == 1 || this->n_fe_components == 1)
4031 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4033 for (
unsigned int comp = 0; comp < n_components; ++comp)
4034 operation.process_dof_gather(dof_indices.data(),
4038 values_dofs[comp][i],
4041 for (
unsigned int v = 0; v < n_lanes; ++v)
4042 dof_indices[v] += offsets[v];
4045 for (
unsigned int comp = 0; comp < n_components; ++comp)
4046 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4048 operation.process_dof_gather(dof_indices.data(),
4052 values_dofs[comp][i],
4055 for (
unsigned int v = 0; v < n_lanes; ++v)
4056 dof_indices[v] += offsets[v];
4061 for (
unsigned int comp = 0; comp < n_components; ++comp)
4063 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4064 operation.process_empty(values_dofs[comp][i]);
4065 if (accesses_exterior_dofs)
4067 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4068 if (mask[v] ==
true)
4071 [ind][cells[v] / VectorizedArrayType::size()] ==
4075 if (n_components == 1 || this->n_fe_components == 1)
4077 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4078 operation.process_dof(dof_indices[v] + i,
4080 values_dofs[comp][i][v]);
4084 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4085 operation.process_dof(dof_indices[v] + i +
4086 comp * dofs_per_component,
4088 values_dofs[comp][i][v]);
4093 const unsigned int offset =
4096 if (n_components == 1 || this->n_fe_components == 1)
4098 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4099 operation.process_dof(dof_indices[v] + i * offset,
4101 values_dofs[comp][i][v]);
4105 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4106 operation.process_dof(
4108 (i + comp * dofs_per_component) * offset,
4110 values_dofs[comp][i][v]);
4121 if (n_components == 1 || this->n_fe_components == 1)
4123 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4124 if (mask[v] ==
true)
4125 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4126 operation.process_dof(dof_indices[v] + i,
4128 values_dofs[comp][i][v]);
4132 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4133 if (mask[v] ==
true)
4134 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4135 operation.process_dof(dof_indices[v] + i +
4136 comp * dofs_per_component,
4138 values_dofs[comp][i][v]);
4143 const unsigned int *offsets =
4145 [ind][VectorizedArrayType::size() * this->cell];
4146 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4148 if (n_components == 1 || this->n_fe_components == 1)
4149 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4151 if (mask[v] ==
true)
4152 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4153 operation.process_dof(dof_indices[v] + i * offsets[v],
4155 values_dofs[comp][i][v]);
4159 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4160 if (mask[v] ==
true)
4161 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4162 operation.process_dof(
4164 (i + comp * dofs_per_component) * offsets[v],
4166 values_dofs[comp][i][v]);
4177 typename VectorType,
4178 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4179 decltype(std::declval<VectorType>().
begin())
4180 get_beginning(VectorType &vec)
4187 typename VectorType,
4188 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4189 typename VectorType::value_type *
4190 get_beginning(VectorType &)
4195 template <
typename VectorType,
4196 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4198 const std::vector<ArrayView<const typename VectorType::value_type>> *
4199 get_shared_vector_data(VectorType *vec,
4200 const bool is_valid_mode_for_sm,
4201 const unsigned int active_fe_index,
4205 if (is_valid_mode_for_sm &&
4208 active_fe_index == 0)
4209 return &vec->shared_vector_data();
4214 template <
typename VectorType,
4215 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4217 const std::vector<ArrayView<const typename VectorType::value_type>> *
4218 get_shared_vector_data(VectorType *,
4226 template <
int n_components,
typename VectorType>
4228 std::array<
typename internal::BlockVectorSelector<
4233 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
4237 get_vector_data(VectorType &src,
4238 const unsigned int first_index,
4239 const bool is_valid_mode_for_sm,
4240 const unsigned int active_fe_index,
4246 std::array<
typename internal::BlockVectorSelector<
4252 ArrayView<
const typename internal::BlockVectorSelector<
4258 for (
unsigned int d = 0;
d < n_components; ++
d)
4259 src_data.first[
d] = internal::BlockVectorSelector<
4265 for (
unsigned int d = 0;
d < n_components; ++
d)
4266 src_data.second[
d] = get_shared_vector_data(
4267 const_cast<typename internal::BlockVectorSelector<
4268 std::remove_const_t<VectorType>,
4270 *>(src_data.first[
d]),
4271 is_valid_mode_for_sm,
4285 typename VectorizedArrayType>
4290 if (this->dof_info ==
nullptr ||
4292 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
4293 this->dof_info->hanging_node_constraint_masks_comp
4294 [this->active_fe_index][this->first_selected_component] ==
false)
4297 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4300 bool hn_available =
false;
4302 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
4304 for (
unsigned int v = 0; v < n_lanes; ++v)
4316 constraint_mask[v] =
mask;
4322 if (hn_available ==
false)
4327 this->data->data.front().fe_degree,
4328 this->get_shape_info(),
4340 typename VectorizedArrayType>
4341 template <
typename VectorType>
4345 const unsigned int first_index,
4346 const std::bitset<n_lanes> &mask)
4348 const auto src_data = internal::get_vector_data<n_components_>(
4351 this->dof_access_index ==
4353 this->active_fe_index,
4357 read_write_operation(reader, src_data.first, src_data.second, mask,
true);
4359 apply_hanging_node_constraints(
false);
4362 this->dof_values_initialized =
true;
4372 typename VectorizedArrayType>
4373 template <
typename VectorType>
4377 const unsigned int first_index,
4378 const std::bitset<n_lanes> &mask)
4380 const auto src_data = internal::get_vector_data<n_components_>(
4383 this->dof_access_index ==
4385 this->active_fe_index,
4389 read_write_operation(reader, src_data.first, src_data.second, mask,
false);
4392 this->dof_values_initialized =
true;
4402 typename VectorizedArrayType>
4403 template <
typename VectorType>
4407 const unsigned int first_index,
4408 const std::bitset<n_lanes> &mask)
const
4411 Assert(this->dof_values_initialized ==
true,
4415 apply_hanging_node_constraints(
true);
4417 const auto dst_data = internal::get_vector_data<n_components_>(
4420 this->dof_access_index ==
4422 this->active_fe_index,
4427 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4436 typename VectorizedArrayType>
4437 template <
typename VectorType>
4441 const unsigned int first_index,
4442 const std::bitset<n_lanes> &mask)
const
4445 Assert(this->dof_values_initialized ==
true,
4449 const auto dst_data = internal::get_vector_data<n_components_>(
4452 this->dof_access_index ==
4454 this->active_fe_index,
4458 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4467 typename VectorizedArrayType>
4468 template <
typename VectorType>
4472 const unsigned int first_index,
4473 const std::bitset<n_lanes> &mask)
const
4476 Assert(this->dof_values_initialized ==
true,
4480 const auto dst_data = internal::get_vector_data<n_components_>(
4483 this->dof_access_index ==
4485 this->active_fe_index,
4489 read_write_operation(setter, dst_data.first, dst_data.second, mask,
false);
4502 typename VectorizedArrayType>
4508 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4510 for (
unsigned int comp = 0; comp < n_components; ++comp)
4511 return_value[comp] = this->values_dofs[comp * dofs + dof];
4512 return return_value;
4521 typename VectorizedArrayType>
4524 get_value(
const unsigned int q_point)
const
4527 Assert(this->values_quad_initialized ==
true,
4532 const std::size_t nqp = this->n_quadrature_points;
4534 for (
unsigned int comp = 0; comp < n_components; ++comp)
4535 return_value[comp] = this->values_quad[comp * nqp + q_point];
4536 return return_value;
4545 typename VectorizedArrayType>
4552 Assert(this->gradients_quad_initialized ==
true,
4557 Assert(this->jacobian !=
nullptr,
4559 "update_gradients"));
4560 const std::size_t nqp = this->n_quadrature_points;
4566 for (
unsigned int comp = 0; comp < n_components; ++comp)
4567 for (
unsigned int d = 0;
d < dim; ++
d)
4569 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4570 this->jacobian[0][
d][
d];
4579 for (
unsigned int comp = 0; comp < n_components; ++comp)
4580 for (
unsigned int d = 0;
d < dim; ++
d)
4583 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4584 for (
unsigned int e = 1;
e < dim; ++
e)
4585 grad_out[comp][
d] +=
4587 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4599 typename VectorizedArrayType>
4606 Assert(this->gradients_quad_initialized ==
true,
4610 Assert(this->normal_x_jacobian !=
nullptr,
4612 "update_gradients"));
4614 const std::size_t nqp = this->n_quadrature_points;
4618 for (
unsigned int comp = 0; comp < n_components; ++comp)
4620 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4621 (this->normal_x_jacobian[0][dim - 1]);
4624 const std::size_t
index =
4626 for (
unsigned int comp = 0; comp < n_components; ++comp)
4628 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4629 this->normal_x_jacobian[index][0];
4630 for (
unsigned int d = 1;
d < dim; ++
d)
4632 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4633 this->normal_x_jacobian[
index][
d];
4645 template <
typename VectorizedArrayType>
4648 const VectorizedArrayType *
const hessians,
4650 VectorizedArrayType (&tmp)[1][1])
4652 tmp[0][0] = jac[0][0] *
hessians[0];
4655 template <
typename VectorizedArrayType>
4658 const VectorizedArrayType *
const hessians,
4659 const unsigned int nqp,
4660 VectorizedArrayType (&tmp)[2][2])
4662 for (
unsigned int d = 0;
d < 2; ++
d)
4670 template <
typename VectorizedArrayType>
4673 const VectorizedArrayType *
const hessians,
4674 const unsigned int nqp,
4675 VectorizedArrayType (&tmp)[3][3])
4677 for (
unsigned int d = 0;
d < 3; ++
d)
4698 typename VectorizedArrayType>
4704 Assert(this->hessians_quad_initialized ==
true,
4709 Assert(this->jacobian !=
nullptr,
4719 const std::size_t nqp = this->n_quadrature_points;
4720 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4725 for (
unsigned int comp = 0; comp < n_components; ++comp)
4727 for (
unsigned int d = 0;
d < dim; ++
d)
4728 hessian_out[comp][
d][
d] =
4729 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4730 (jac[
d][
d] * jac[
d][
d]);
4736 hessian_out[comp][0][1] =
4737 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4738 (jac[0][0] * jac[1][1]);
4741 hessian_out[comp][0][1] =
4742 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4743 (jac[0][0] * jac[1][1]);
4744 hessian_out[comp][0][2] =
4745 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4746 (jac[0][0] * jac[2][2]);
4747 hessian_out[comp][1][2] =
4748 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4749 (jac[1][1] * jac[2][2]);
4754 for (
unsigned int d = 0;
d < dim; ++
d)
4755 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4756 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4762 for (
unsigned int comp = 0; comp < n_components; ++comp)
4764 VectorizedArrayType tmp[dim][dim];
4765 internal::hessian_unit_times_jac(
4766 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4769 for (
unsigned int d = 0;
d < dim; ++
d)
4770 for (
unsigned int e =
d;
e < dim; ++
e)
4772 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4773 for (
unsigned int f = 1; f < dim; ++f)
4774 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4781 for (
unsigned int d = 0;
d < dim; ++
d)
4782 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4783 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4789 const auto &jac_grad = this->jacobian_gradients[q_point];
4790 for (
unsigned int comp = 0; comp < n_components; ++comp)
4792 VectorizedArrayType tmp[dim][dim];
4793 internal::hessian_unit_times_jac(
4794 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4797 for (
unsigned int d = 0;
d < dim; ++
d)
4798 for (
unsigned int e =
d;
e < dim; ++
e)
4800 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4801 for (
unsigned int f = 1; f < dim; ++f)
4802 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4806 for (
unsigned int d = 0;
d < dim; ++
d)
4807 for (
unsigned int e = 0;
e < dim; ++
e)
4808 hessian_out[comp][
d][
d] +=
4810 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4813 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4814 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++count)
4815 for (
unsigned int f = 0; f < dim; ++f)
4816 hessian_out[comp][
d][
e] +=
4817 jac_grad[count][f] *
4818 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4821 for (
unsigned int d = 0;
d < dim; ++
d)
4822 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4823 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4835 typename VectorizedArrayType>
4842 Assert(this->hessians_quad_initialized ==
true,
4853 const std::size_t nqp = this->n_quadrature_points;
4854 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4860 for (
unsigned int comp = 0; comp < n_components; ++comp)
4861 for (
unsigned int d = 0;
d < dim; ++
d)
4862 hessian_out[comp][
d] =
4863 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4864 (jac[
d][
d] * jac[
d][
d]);
4869 for (
unsigned int comp = 0; comp < n_components; ++comp)
4873 VectorizedArrayType tmp[dim][dim];
4874 internal::hessian_unit_times_jac(
4875 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4879 for (
unsigned int d = 0;
d < dim; ++
d)
4881 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4882 for (
unsigned int f = 1; f < dim; ++f)
4883 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4890 const auto &jac_grad = this->jacobian_gradients[q_point];
4891 for (
unsigned int comp = 0; comp < n_components; ++comp)
4895 VectorizedArrayType tmp[dim][dim];
4896 internal::hessian_unit_times_jac(
4897 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4901 for (
unsigned int d = 0;
d < dim; ++
d)
4903 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4904 for (
unsigned int f = 1; f < dim; ++f)
4905 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4908 for (
unsigned int d = 0;
d < dim; ++
d)
4909 for (
unsigned int e = 0;
e < dim; ++
e)
4910 hessian_out[comp][
d] +=
4912 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4924 typename VectorizedArrayType>
4931 Assert(this->hessians_quad_initialized ==
true,
4937 const auto hess_diag = get_hessian_diagonal(q_point);
4938 for (
unsigned int comp = 0; comp < n_components; ++comp)
4940 laplacian_out[comp] = hess_diag[comp][0];
4941 for (
unsigned int d = 1;
d < dim; ++
d)
4942 laplacian_out[comp] += hess_diag[comp][
d];
4944 return laplacian_out;
4953 typename VectorizedArrayType>
4957 const unsigned int dof)
4960 this->dof_values_initialized =
true;
4962 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4964 for (
unsigned int comp = 0; comp < n_components; ++comp)
4965 this->values_dofs[comp * dofs + dof] = val_in[comp];
4974 typename VectorizedArrayType>
4978 const unsigned int q_point)
4984 Assert(this->J_value !=
nullptr,
4988 this->values_quad_submitted =
true;
4991 const std::size_t nqp = this->n_quadrature_points;
4994 const VectorizedArrayType JxW =
4995 this->J_value[0] * this->quadrature_weights[q_point];
4996 for (
unsigned int comp = 0; comp < n_components; ++comp)
4997 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5001 const VectorizedArrayType JxW = this->J_value[q_point];
5002 for (
unsigned int comp = 0; comp < n_components; ++comp)
5003 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5013 typename VectorizedArrayType>
5018 const unsigned int q_point)
5024 Assert(this->J_value !=
nullptr,
5026 "update_gradients"));
5027 Assert(this->jacobian !=
nullptr,
5029 "update_gradients"));
5031 this->gradients_quad_submitted =
true;
5034 const std::size_t nqp_d = this->n_quadrature_points * dim;
5035 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5039 const VectorizedArrayType JxW =
5040 this->J_value[0] * this->quadrature_weights[q_point];
5041 std::array<VectorizedArrayType, dim> jac;
5042 for (
unsigned int d = 0;
d < dim; ++
d)
5043 jac[
d] = this->jacobian[0][
d][
d];
5044 for (
unsigned int d = 0;
d < dim; ++
d)
5046 const VectorizedArrayType factor = jac[
d] * JxW;
5047 for (
unsigned int comp = 0; comp < n_components; ++comp)
5048 gradients[comp * nqp_d +
d] = grad_in[comp][
d] * factor;
5055 this->jacobian[q_point] :
5057 const VectorizedArrayType JxW =
5059 this->J_value[q_point] :
5060 this->J_value[0] * this->quadrature_weights[q_point];
5061 for (
unsigned int comp = 0; comp < n_components; ++comp)
5062 for (
unsigned int d = 0;
d < dim; ++
d)
5064 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5065 for (
unsigned int e = 1;
e < dim; ++
e)
5066 new_val += (jac[
e][
d] * grad_in[comp][
e]);
5078 typename VectorizedArrayType>
5083 const unsigned int q_point)
5086 Assert(this->normal_x_jacobian !=
nullptr,
5088 "update_gradients"));
5090 this->gradients_quad_submitted =
true;
5093 const std::size_t nqp_d = this->n_quadrature_points * dim;
5094 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5098 const VectorizedArrayType JxW_jac = this->J_value[0] *
5099 this->quadrature_weights[q_point] *
5100 this->normal_x_jacobian[0][dim - 1];
5101 for (
unsigned int comp = 0; comp < n_components; ++comp)
5103 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5104 gradients[comp * nqp_d +
d] = VectorizedArrayType();
5105 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5110 const unsigned int index =
5113 this->normal_x_jacobian[
index];
5114 const VectorizedArrayType JxW =
5116 this->J_value[index] * this->quadrature_weights[q_point] :
5117 this->J_value[index];
5118 for (
unsigned int comp = 0; comp < n_components; ++comp)
5120 for (
unsigned int d = 0;
d < dim; ++
d)
5121 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[
d];
5132 typename VectorizedArrayType>
5138 const unsigned int q_point)
5144 Assert(this->J_value !=
nullptr,
5146 "update_hessians"));
5147 Assert(this->jacobian !=
nullptr,
5149 "update_hessians"));
5151 this->hessians_quad_submitted =
true;
5155 const std::size_t nqp = this->n_quadrature_points;
5156 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
5159 const VectorizedArrayType JxW =
5160 this->J_value[0] * this->quadrature_weights[q_point];
5163 for (
unsigned int d = 0;
d < dim; ++
d)
5165 const auto jac_d = this->jacobian[0][
d][
d];
5166 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5167 for (
unsigned int comp = 0; comp < n_components; ++comp)
5168 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5169 hessian_in[comp][
d][
d] * factor;
5173 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5174 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5176 const auto jac_d = this->jacobian[0][
d][
d];
5177 const auto jac_e = this->jacobian[0][
e][
e];
5178 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5179 for (
unsigned int comp = 0; comp < n_components; ++comp)
5180 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5181 (hessian_in[comp][
d][
e] + hessian_in[comp][
e][
d]) * factor;
5188 const VectorizedArrayType JxW =
5189 this->J_value[0] * this->quadrature_weights[q_point];
5190 for (
unsigned int comp = 0; comp < n_components; ++comp)
5193 VectorizedArrayType tmp[dim][dim];
5194 for (
unsigned int i = 0; i < dim; ++i)
5195 for (
unsigned int j = 0; j < dim; ++j)
5197 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5198 for (
unsigned int k = 1; k < dim; ++k)
5199 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5203 VectorizedArrayType tmp2[dim][dim];
5204 for (
unsigned int i = 0; i < dim; ++i)
5205 for (
unsigned int j = 0; j < dim; ++j)
5207 tmp2[i][j] = jac[0][i] * tmp[0][j];
5208 for (
unsigned int k = 1; k < dim; ++k)
5209 tmp2[i][j] += jac[k][i] * tmp[k][j];
5213 for (
unsigned int d = 0;
d < dim; ++
d)
5214 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5218 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5219 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5220 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5221 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5227 const VectorizedArrayType JxW = this->J_value[q_point];
5228 const auto &jac_grad = this->jacobian_gradients[q_point];
5229 for (
unsigned int comp = 0; comp < n_components; ++comp)
5232 VectorizedArrayType tmp[dim][dim];
5233 for (
unsigned int i = 0; i < dim; ++i)
5234 for (
unsigned int j = 0; j < dim; ++j)
5236 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5237 for (
unsigned int k = 1; k < dim; ++k)
5238 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5242 VectorizedArrayType tmp2[dim][dim];
5243 for (
unsigned int i = 0; i < dim; ++i)
5244 for (
unsigned int j = 0; j < dim; ++j)
5246 tmp2[i][j] = jac[0][i] * tmp[0][j];
5247 for (
unsigned int k = 1; k < dim; ++k)
5248 tmp2[i][j] += jac[k][i] * tmp[k][j];
5252 for (
unsigned int d = 0;
d < dim; ++
d)
5253 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5257 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5258 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5259 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5260 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5263 for (
unsigned int d = 0;
d < dim; ++
d)
5265 VectorizedArrayType
sum = 0;
5266 for (
unsigned int e = 0;
e < dim; ++
e)
5267 sum += hessian_in[comp][
e][
e] * jac_grad[
e][
d];
5268 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5269 for (
unsigned int f =
e + 1; f < dim; ++f, ++count)
5270 sum += (hessian_in[comp][
e][f] + hessian_in[comp][f][
e]) *
5272 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5285 typename VectorizedArrayType>
5292 Assert(this->values_quad_submitted ==
true,
5297 const std::size_t nqp = this->n_quadrature_points;
5298 for (
unsigned int q = 0; q < nqp; ++q)
5299 for (
unsigned int comp = 0; comp < n_components; ++comp)
5300 return_value[comp] += this->values_quad[comp * nqp + q];
5301 return (return_value);
5313 typename VectorizedArrayType>
5318 VectorizedArrayType>::
5321 const unsigned int dof_no,
5322 const unsigned int first_selected_component,
5323 const unsigned int quad_no,
5324 const unsigned int fe_degree,
5325 const unsigned int n_q_points,
5326 const bool is_interior_face,
5327 const unsigned int active_fe_index,
5328 const unsigned int active_quad_index,
5329 const unsigned int face_type)
5330 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5333 first_selected_component,
5349 typename VectorizedArrayType>
5354 VectorizedArrayType>::
5360 const unsigned int first_selected_component,
5362 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5367 first_selected_component,
5377 typename VectorizedArrayType>
5382 VectorizedArrayType>::
5387 VectorizedArrayType> &other)
5388 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5398 typename VectorizedArrayType>
5403 VectorizedArrayType> &
5409 VectorizedArrayType> &other)
5415 VectorizedArrayType>::operator=(other);
5424 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5428 const unsigned int dof_no,
5429 const unsigned int first_selected_component,
5430 const unsigned int quad_no,
5431 const unsigned int fe_degree,
5432 const unsigned int n_q_points,
5433 const bool is_interior_face,
5434 const unsigned int active_fe_index,
5435 const unsigned int active_quad_index,
5436 const unsigned int face_type)
5440 first_selected_component,
5452 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5459 const unsigned int first_selected_component,
5466 first_selected_component,
5472 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5482 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5488 ->
FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5495 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5498 const unsigned int dof)
const
5501 return this->values_dofs[dof];
5506 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5509 const unsigned int q_point)
const
5512 Assert(this->values_quad_initialized ==
true,
5516 return this->values_quad[q_point];
5521 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5526 return BaseClass::get_normal_derivative(q_point)[0];
5531 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5534 const unsigned int q_point_in)
const
5540 Assert(this->gradients_quad_initialized ==
true,
5545 Assert(this->jacobian !=
nullptr,
5547 "update_gradients"));
5551 const std::size_t q_point = q_point_in;
5554 for (
unsigned int d = 0;
d < dim; ++
d)
5556 this->gradients_quad[dim * q_point +
d] * this->jacobian[0][
d][
d];
5565 for (
unsigned int d = 0;
d < dim; ++
d)
5567 grad_out[
d] = jac[
d][0] * this->gradients_quad[dim * q_point];
5568 for (
unsigned int e = 1;
e < dim; ++
e)
5569 grad_out[
d] += jac[
d][
e] * this->gradients_quad[dim * q_point +
e];
5577 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5580 const unsigned int q_point)
const
5582 return BaseClass::get_hessian(q_point)[0];
5587 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5592 return BaseClass::get_hessian_diagonal(q_point)[0];
5597 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5598 inline VectorizedArrayType
5600 const unsigned int q_point)
const
5602 return BaseClass::get_laplacian(q_point)[0];
5607 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5613 this->dof_values_initialized =
true;
5616 this->values_dofs[dof] = val_in;
5621 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5624 const VectorizedArrayType val_in,
5625 const unsigned int q_point)
5631 Assert(this->J_value !=
nullptr,
5635 this->values_quad_submitted =
true;
5640 const VectorizedArrayType JxW =
5641 this->J_value[0] * this->quadrature_weights[q_point];
5642 this->values_quad[q_point] = val_in * JxW;
5646 this->values_quad[q_point] = val_in * this->J_value[q_point];
5652 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5656 const unsigned int q_point)
5658 submit_value(val_in[0], q_point);
5663 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5667 const unsigned int q_point)
5671 BaseClass::submit_normal_derivative(grad, q_point);
5676 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5680 const unsigned int q_point_in)
5686 Assert(this->J_value !=
nullptr,
5688 "update_gradients"));
5689 Assert(this->jacobian !=
nullptr,
5691 "update_gradients"));
5693 this->gradients_quad_submitted =
true;
5696 const std::size_t q_point = q_point_in;
5697 VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point;
5700 const VectorizedArrayType JxW =
5701 this->J_value[0] * this->quadrature_weights[q_point];
5705 std::array<VectorizedArrayType, dim> jac;
5706 for (
unsigned int d = 0;
d < dim; ++
d)
5707 jac[
d] = this->jacobian[0][
d][
d];
5709 for (
unsigned int d = 0;
d < dim; ++
d)
5710 grad_ptr[
d] = grad_in[
d] * jac[
d] * JxW;
5717 this->jacobian[q_point] :
5719 const VectorizedArrayType JxW =
5721 this->J_value[q_point] :
5722 this->J_value[0] * this->quadrature_weights[q_point];
5723 for (
unsigned int d = 0;
d < dim; ++
d)
5725 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5726 for (
unsigned int e = 1;
e < dim; ++
e)
5727 new_val += jac[
e][
d] * grad_in[
e];
5728 grad_ptr[
d] = new_val * JxW;
5735 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5739 const unsigned int q_point)
5743 BaseClass::submit_hessian(hessian, q_point);
5748 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5749 inline VectorizedArrayType
5753 return BaseClass::integrate_value()[0];
5761 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5765 const unsigned int dof_no,
5766 const unsigned int first_selected_component,
5767 const unsigned int quad_no,
5768 const unsigned int fe_degree,
5769 const unsigned int n_q_points,
5770 const bool is_interior_face,
5771 const unsigned int active_fe_index,
5772 const unsigned int active_quad_index,
5773 const unsigned int face_type)
5777 first_selected_component,
5789 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5796 const unsigned int first_selected_component,
5803 first_selected_component,
5809 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5819 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5831 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5834 const unsigned int q_point)
const
5836 if (this->data->element_type ==
5841 Assert(this->values_quad_initialized ==
true,
5846 Assert(this->J_value !=
nullptr,
5849 const std::size_t nqp = this->n_quadrature_points;
5857 const VectorizedArrayType inv_det =
5858 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5859 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5860 this->jacobian[0][2][2];
5863 for (
unsigned int comp = 0; comp < n_components; ++comp)
5864 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5865 jac[comp][comp] * inv_det;
5872 this->jacobian[q_point] :
5880 const VectorizedArrayType inv_det =
5881 (is_face && dim == 2 && this->get_face_no() < 2) ?
5885 for (
unsigned int comp = 0; comp < n_components; ++comp)
5887 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
5888 for (
unsigned int e = 1;
e < dim; ++
e)
5890 this->values_quad[
e * nqp + q_point] * jac[comp][
e];
5891 value_out[comp] *= inv_det;
5899 return BaseClass::get_value(q_point);
5905 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5910 if (this->data->element_type ==
5915 Assert(this->gradients_quad_initialized ==
true,
5920 Assert(this->jacobian !=
nullptr,
5922 "update_gradients"));
5923 const std::size_t nqp = this->n_quadrature_points;
5924 const std::size_t nqp_d = nqp * dim;
5927 this->gradients_quad + q_point * dim;
5937 const VectorizedArrayType inv_det =
5938 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5939 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5940 this->jacobian[0][2][2];
5943 for (
unsigned int d = 0;
d < dim; ++
d)
5944 for (
unsigned int comp = 0; comp < n_components; ++comp)
5946 inv_t_jac[
d][
d] * (jac[comp][comp] * inv_det);
5956 const VectorizedArrayType inv_det =
5957 (is_face && dim == 2 && this->get_face_no() < 2) ?
5961 VectorizedArrayType tmp[dim][dim];
5963 for (
unsigned int d = 0;
d < dim; ++
d)
5964 for (
unsigned int e = 0;
e < dim; ++
e)
5967 for (
unsigned int f = 1; f < dim; ++f)
5970 for (
unsigned int comp = 0; comp < n_components; ++comp)
5971 for (
unsigned int d = 0;
d < dim; ++
d)
5973 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
5974 for (
unsigned int f = 1; f < dim; ++f)
5975 res += jac[comp][f] * tmp[
d][f];
5977 grad_out[comp][
d] = res * inv_det;
5988 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
5990 "update_hessians"));
5992 const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
5994 this->jacobian[q_point];
5997 const VectorizedArrayType inv_det =
5998 (is_face && dim == 2 && this->get_face_no() < 2) ?
6004 VectorizedArrayType tmp[dim][dim];
6005 for (
unsigned int d = 0;
d < dim; ++
d)
6006 for (
unsigned int e = 0;
e < dim; ++
e)
6009 for (
unsigned int f = 1; f < dim; ++f)
6015 for (
unsigned int d = 0;
d < dim; ++
d)
6017 for (
unsigned int e = 0;
e < dim; ++
e)
6019 jac_grad[
e][
d] * this->values_quad[
e * nqp + q_point];
6020 for (
unsigned int f = 0, r = dim; f < dim; ++f)
6021 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
6024 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
6026 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
6031 for (
unsigned int d = 0;
d < dim; ++
d)
6032 for (
unsigned int e = 0;
e < dim; ++
e)
6034 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
6035 for (
unsigned int f = 1; f < dim; ++f)
6036 res += tmp[f][
d] * inv_t_jac[
e][f];
6037 grad_out[
d][
e] = res;
6043 VectorizedArrayType tmp3[dim], tmp4[dim];
6044 for (
unsigned int d = 0;
d < dim; ++
d)
6046 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
6047 for (
unsigned int e = 1;
e < dim; ++
e)
6048 tmp3[
d] += inv_t_jac[
e][
d] * jac_grad[
d][
e];
6050 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
6051 for (
unsigned int f =
e + 1; f < dim; ++k, ++f)
6052 for (
unsigned int d = 0;
d < dim; ++
d)
6054 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
6055 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
6057 for (
unsigned int d = 0;
d < dim; ++
d)
6059 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
6060 for (
unsigned int e = 1;
e < dim; ++
e)
6061 tmp4[
d] += tmp3[
e] * inv_t_jac[
d][
e];
6064 VectorizedArrayType tmp2[dim];
6065 for (
unsigned int d = 0;
d < dim; ++
d)
6067 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
6068 for (
unsigned e = 1;
e < dim; ++
e)
6069 tmp2[
d] += t_jac[
e][
d] * this->values_quad[
e * nqp + q_point];
6072 for (
unsigned int d = 0;
d < dim; ++
d)
6073 for (
unsigned int e = 0;
e < dim; ++
e)
6075 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
6079 grad_out[
d][
e] *= inv_det;
6086 return BaseClass::get_gradient(q_point);
6092 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6098 Assert(this->gradients_quad_initialized ==
true,
6102 Assert(this->jacobian !=
nullptr,
6104 "update_gradients"));
6106 VectorizedArrayType divergence;
6107 const std::size_t nqp = this->n_quadrature_points;
6109 if (this->data->element_type ==
6112 VectorizedArrayType inv_det =
6115 this->jacobian[0][0][0] *
6116 ((dim == 2) ? this->jacobian[0][1][1] :
6117 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
6125 if (is_face && dim == 2 && this->get_face_no() < 2)
6129 divergence = this->gradients_quad[q_point * dim];
6130 for (
unsigned int d = 1;
d < dim; ++
d)
6131 divergence += this->gradients_quad[(
d * nqp + q_point) * dim +
d];
6132 divergence *= inv_det;
6141 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
6142 for (
unsigned int d = 1;
d < dim; ++
d)
6143 divergence += this->gradients_quad[(
d * nqp + q_point) * dim +
d] *
6144 this->jacobian[0][
d][
d];
6151 this->jacobian[q_point] :
6153 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
6154 for (
unsigned int e = 1;
e < dim; ++
e)
6155 divergence += jac[0][
e] * this->gradients_quad[q_point * dim +
e];
6156 for (
unsigned int d = 1;
d < dim; ++
d)
6157 for (
unsigned int e = 0;
e < dim; ++
e)
6159 jac[
d][
e] * this->gradients_quad[(
d * nqp + q_point) * dim +
e];
6167 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6173 const auto grad = get_gradient(q_point);
6174 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6175 VectorizedArrayType half = Number(0.5);
6176 for (
unsigned int d = 0;
d < dim; ++
d)
6177 symmetrized[
d] = grad[
d][
d];
6183 symmetrized[2] = grad[0][1] + grad[1][0];
6184 symmetrized[2] *= half;
6187 symmetrized[3] = grad[0][1] + grad[1][0];
6188 symmetrized[3] *= half;
6189 symmetrized[4] = grad[0][2] + grad[2][0];
6190 symmetrized[4] *= half;
6191 symmetrized[5] = grad[1][2] + grad[2][1];
6192 symmetrized[5] *= half;
6202 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6204 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6206 const unsigned int q_point)
const
6210 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6216 "Computing the curl in 1d is not a useful operation"));
6219 curl[0] = grad[1][0] - grad[0][1];
6222 curl[0] = grad[2][1] - grad[1][2];
6223 curl[1] = grad[0][2] - grad[2][0];
6224 curl[2] = grad[1][0] - grad[0][1];
6234 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6239 return BaseClass::get_hessian_diagonal(q_point);
6244 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6247 const unsigned int q_point)
const
6250 Assert(this->hessians_quad_initialized ==
true,
6254 return BaseClass::get_hessian(q_point);
6258 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6262 const unsigned int q_point)
6264 if (this->data->element_type ==
6269 Assert(this->J_value !=
nullptr,
6274 this->values_quad_submitted =
true;
6277 VectorizedArrayType *
values = this->values_quad + q_point;
6278 const std::size_t nqp = this->n_quadrature_points;
6284 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6286 for (
unsigned int comp = 0; comp < n_components; ++comp)
6287 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
6294 this->jacobian[q_point] :
6299 const VectorizedArrayType fac =
6301 this->quadrature_weights[q_point] :
6303 this->J_value[q_point] :
6304 this->J_value[0] * this->quadrature_weights[q_point]) *
6305 ((dim == 2 && this->get_face_no() < 2) ?
6314 for (
unsigned int comp = 0; comp < n_components; ++comp)
6316 values[comp * nqp] = val_in[0] * jac[0][comp];
6317 for (
unsigned int e = 1;
e < dim; ++
e)
6318 values[comp * nqp] += val_in[
e] * jac[
e][comp];
6319 values[comp * nqp] *= fac;
6326 BaseClass::submit_value(val_in, q_point);
6332 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6336 const unsigned int q_point)
6338 if (this->data->element_type ==
6347 Assert(this->J_value !=
nullptr,
6349 "update_gradients"));
6350 Assert(this->jacobian !=
nullptr,
6352 "update_gradients"));
6354 this->gradients_quad_submitted =
true;
6357 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
6358 VectorizedArrayType *
values = this->values_from_gradients_quad + q_point;
6359 const std::size_t nqp = this->n_quadrature_points;
6360 const std::size_t nqp_d = nqp * dim;
6369 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6370 for (
unsigned int d = 0;
d < dim; ++
d)
6371 for (
unsigned int comp = 0; comp < n_components; ++comp)
6373 grad_in[comp][
d] * inv_t_jac[
d][
d] * (jac[comp][comp] * weight);
6384 const VectorizedArrayType fac =
6385 (!is_face) ? this->quadrature_weights[q_point] :
6386 this->J_value[0] * this->quadrature_weights[q_point] *
6387 ((dim == 2 && this->get_face_no() < 2) ?
6392 VectorizedArrayType tmp[dim][dim];
6393 for (
unsigned int d = 0;
d < dim; ++
d)
6394 for (
unsigned int e = 0;
e < dim; ++
e)
6396 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
6397 for (
unsigned int f = 1; f < dim; ++f)
6398 tmp[
d][
e] += inv_t_jac[f][
d] * grad_in[
e][f];
6400 for (
unsigned int comp = 0; comp < n_components; ++comp)
6401 for (
unsigned int d = 0;
d < dim; ++
d)
6403 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
6404 for (
unsigned int f = 1; f < dim; ++f)
6405 res += jac[f][comp] * tmp[
d][f];
6414 const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
6416 this->jacobian[q_point];
6420 const VectorizedArrayType fac =
6422 this->quadrature_weights[q_point] :
6423 this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6431 VectorizedArrayType tmp3[dim], tmp4[dim];
6432 for (
unsigned int d = 0;
d < dim; ++
d)
6434 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
6435 for (
unsigned int e = 1;
e < dim; ++
e)
6436 tmp3[
d] += inv_t_jac[
e][
d] * jac_grad[
d][
e];
6438 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
6439 for (
unsigned int f =
e + 1; f < dim; ++k, ++f)
6440 for (
unsigned int d = 0;
d < dim; ++
d)
6442 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
6443 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
6445 for (
unsigned int d = 0;
d < dim; ++
d)
6447 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
6448 for (
unsigned int e = 1;
e < dim; ++
e)
6449 tmp4[
d] += tmp3[
e] * inv_t_jac[
d][
e];
6455 VectorizedArrayType tmp[dim][dim];
6458 for (
unsigned int d = 0;
d < dim; ++
d)
6459 for (
unsigned int e = 0;
e < dim; ++
e)
6461 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
6462 for (
unsigned int f = 1; f < dim; ++f)
6463 tmp[
d][
e] += inv_t_jac[f][
d] * grad_in_scaled[
e][f];
6466 for (
unsigned int d = 0;
d < dim; ++
d)
6467 for (
unsigned int e = 0;
e < dim; ++
e)
6469 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
6470 for (
unsigned int f = 1; f < dim; ++f)
6471 res += t_jac[
d][f] * tmp[
e][f];
6478 VectorizedArrayType
value[dim];
6479 for (
unsigned int d = 0;
d < dim; ++
d)
6481 value[
d] = tmp[
d][0] * jac_grad[
d][0];
6482 for (
unsigned int e = 1;
e < dim; ++
e)
6483 value[
d] += tmp[
d][
e] * jac_grad[
d][
e];
6485 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
6486 for (
unsigned int f =
e + 1; f < dim; ++k, ++f)
6487 for (
unsigned int d = 0;
d < dim; ++
d)
6489 value[
e] += tmp[f][
d] * jac_grad[k][
d];
6490 value[f] += tmp[
e][
d] * jac_grad[k][
d];
6495 for (
unsigned int d = 0;
d < dim; ++
d)
6497 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
6498 for (
unsigned int e = 1;
e < dim; ++
e)
6499 tmp2 += grad_in_scaled[
d][
e] * tmp4[
e];
6500 for (
unsigned int e = 0;
e < dim; ++
e)
6501 value[
e] -= t_jac[
e][
d] * tmp2;
6504 for (
unsigned int d = 0;
d < dim; ++
d)
6510 BaseClass::submit_gradient(grad_in, q_point);
6516 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6521 const unsigned int q_point)
6523 if (this->data->element_type ==
6533 BaseClass::submit_gradient(grad_in, q_point);
6539 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6543 const unsigned int q_point)
6549 Assert(this->J_value !=
nullptr,
6551 "update_gradients"));
6552 Assert(this->jacobian !=
nullptr,
6554 "update_gradients"));
6556 this->gradients_quad_submitted =
true;
6559 const std::size_t nqp_d = this->n_quadrature_points * dim;
6560 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
6562 if (this->data->element_type ==
6569 const VectorizedArrayType fac =
6571 this->quadrature_weights[q_point] * div_in :
6573 this->J_value[q_point] :
6574 this->J_value[0] * this->quadrature_weights[q_point]) *
6577 this->jacobian[this->cell_type >
6581 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6583 for (
unsigned int d = 0;
d < dim; ++
d)
6585 for (
unsigned int e = 0;
e < dim; ++
e)
6588 this->divergence_is_requested =
true;
6595 const VectorizedArrayType fac =
6596 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6597 for (
unsigned int d = 0;
d < dim; ++
d)
6599 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
6600 for (
unsigned int e = 0;
e < dim; ++
e)
6608 this->jacobian[q_point] :
6610 const VectorizedArrayType fac =
6612 this->J_value[q_point] :
6613 this->J_value[0] * this->quadrature_weights[q_point]) *
6615 for (
unsigned int d = 0;
d < dim; ++
d)
6617 for (
unsigned int e = 0;
e < dim; ++
e)
6626 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6631 const unsigned int q_point)
6634 this->data->element_type !=
6645 Assert(this->J_value !=
nullptr,
6647 "update_gradients"));
6648 Assert(this->jacobian !=
nullptr,
6650 "update_gradients"));
6652 this->gradients_quad_submitted =
true;
6655 const std::size_t nqp_d = this->n_quadrature_points * dim;
6656 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
6659 const VectorizedArrayType JxW =
6660 this->J_value[0] * this->quadrature_weights[q_point];
6662 for (
unsigned int d = 0;
d < dim; ++
d)
6665 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6666 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6668 const VectorizedArrayType
value =
6677 const VectorizedArrayType JxW =
6679 this->J_value[q_point] :
6680 this->J_value[0] * this->quadrature_weights[q_point];
6683 this->jacobian[q_point] :
6685 VectorizedArrayType weighted[dim][dim];
6686 for (
unsigned int i = 0; i < dim; ++i)
6688 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6689 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6691 const VectorizedArrayType
value =
6693 weighted[i][j] =
value;
6694 weighted[j][i] =
value;
6696 for (
unsigned int comp = 0; comp < dim; ++comp)
6697 for (
unsigned int d = 0;
d < dim; ++
d)
6699 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6700 for (
unsigned int e = 1;
e < dim; ++
e)
6701 new_val += jac[
e][
d] * weighted[comp][
e];
6709 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6713 const unsigned int q_point)
6721 "Testing by the curl in 1d is not a useful operation"));
6724 grad[1][0] = curl[0];
6725 grad[0][1] = -curl[0];
6728 grad[2][1] = curl[0];
6729 grad[1][2] = -curl[0];
6730 grad[0][2] = curl[1];
6731 grad[2][0] = -curl[1];
6732 grad[1][0] = curl[2];
6733 grad[0][1] = -curl[2];
6738 submit_gradient(grad, q_point);
6745 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6749 const unsigned int dof_no,
6750 const unsigned int first_selected_component,
6751 const unsigned int quad_no,
6752 const unsigned int fe_degree,
6753 const unsigned int n_q_points,
6754 const bool is_interior_face,
6755 const unsigned int active_fe_index,
6756 const unsigned int active_quad_index,
6757 const unsigned int face_type)
6761 first_selected_component,
6773 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6780 const unsigned int first_selected_component,
6787 first_selected_component,
6793 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6802 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6814 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6817 const unsigned int dof)
const
6820 return this->values_dofs[dof];
6825 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6828 const unsigned int q_point)
const
6831 Assert(this->values_quad_initialized ==
true,
6835 return this->values_quad[q_point];
6840 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6843 const unsigned int q_point)
const
6849 Assert(this->gradients_quad_initialized ==
true,
6856 this->jacobian[q_point] :
6860 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6867 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6870 const unsigned int q_point)
const
6872 return get_gradient(q_point)[0];
6877 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6882 return BaseClass::get_normal_derivative(q_point)[0];
6887 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6890 const unsigned int q_point)
const
6892 return BaseClass::get_hessian(q_point)[0];
6897 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6902 return BaseClass::get_hessian_diagonal(q_point)[0];
6907 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6910 const unsigned int q_point)
const
6912 return BaseClass::get_laplacian(q_point)[0];
6917 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6923 this->dof_values_initialized =
true;
6926 this->values_dofs[dof] = val_in;
6931 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6934 const VectorizedArrayType val_in,
6935 const unsigned int q_point)
6942 this->values_quad_submitted =
true;
6947 const VectorizedArrayType JxW = this->J_value[q_point];
6948 this->values_quad[q_point] = val_in * JxW;
6952 const VectorizedArrayType JxW =
6953 this->J_value[0] * this->quadrature_weights[q_point];
6954 this->values_quad[q_point] = val_in * JxW;
6960 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6964 const unsigned int q_point)
6966 submit_value(val_in[0], q_point);
6971 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6975 const unsigned int q_point)
6977 submit_gradient(grad_in[0], q_point);
6982 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6985 const VectorizedArrayType grad_in,
6986 const unsigned int q_point)
6993 this->gradients_quad_submitted =
true;
6998 this->jacobian[q_point] :
7000 const VectorizedArrayType JxW =
7002 this->J_value[q_point] :
7003 this->J_value[0] * this->quadrature_weights[q_point];
7005 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7010 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7014 const unsigned int q_point)
7016 submit_gradient(grad_in[0][0], q_point);
7021 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7025 const unsigned int q_point)
7029 BaseClass::submit_normal_derivative(grad, q_point);
7034 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7038 const unsigned int q_point)
7040 BaseClass::submit_normal_derivative(grad_in, q_point);
7044 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7048 const unsigned int q_point)
7052 BaseClass::submit_hessian(hessian, q_point);
7056 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7057 inline VectorizedArrayType
7061 return BaseClass::integrate_value()[0];
7074 typename VectorizedArrayType>
7080 VectorizedArrayType>::
7082 const unsigned int fe_no,
7083 const unsigned int quad_no,
7084 const unsigned int first_selected_component,
7085 const unsigned int active_fe_index,
7086 const unsigned int active_quad_index)
7087 : BaseClass(matrix_free,
7089 first_selected_component,
7096 , dofs_per_component(this->data->dofs_per_component_on_cell)
7097 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7098 , n_q_points(this->data->n_q_points)
7100 check_template_arguments(fe_no, 0);
7110 typename VectorizedArrayType>
7116 VectorizedArrayType>::
7118 const std::pair<unsigned int, unsigned int> &range,
7119 const unsigned int dof_no,
7120 const unsigned int quad_no,
7121 const unsigned int first_selected_component)
7125 first_selected_component,
7126 matrix_free.get_cell_active_fe_index(range))
7136 typename VectorizedArrayType>
7142 VectorizedArrayType>::
7147 const unsigned int first_selected_component)
7148 : BaseClass(mapping,
7152 first_selected_component,
7154 , dofs_per_component(this->data->dofs_per_component_on_cell)
7155 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7156 , n_q_points(this->data->n_q_points)
7168 typename VectorizedArrayType>
7174 VectorizedArrayType>::
7178 const unsigned int first_selected_component)
7183 first_selected_component,
7185 , dofs_per_component(this->data->dofs_per_component_on_cell)
7186 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7187 , n_q_points(this->data->n_q_points)
7199 typename VectorizedArrayType>
7205 VectorizedArrayType>::
7208 const unsigned int first_selected_component)
7209 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7211 other.mapped_geometry->get_quadrature(),
7212 other.mapped_geometry->get_fe_values().get_update_flags(),
7213 first_selected_component,
7215 , dofs_per_component(this->data->dofs_per_component_on_cell)
7216 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7217 , n_q_points(this->data->n_q_points)
7229 typename VectorizedArrayType>
7238 , dofs_per_component(this->data->dofs_per_component_on_cell)
7239 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7240 , n_q_points(this->data->n_q_points)
7252 typename VectorizedArrayType>
7258 VectorizedArrayType> &
7264 VectorizedArrayType>::operator=(
const FEEvaluation &other)
7266 BaseClass::operator=(other);
7278 typename VectorizedArrayType>
7285 VectorizedArrayType>::
7286 check_template_arguments(
const unsigned int dof_no,
7287 const unsigned int first_selected_component)
7290 (void)first_selected_component;
7293 this->data->dofs_per_component_on_cell > 0,
7295 "There is nothing useful you can do with an FEEvaluation object with "
7296 "FE_Nothing, i.e., without DoFs! If you have passed to "
7297 "MatrixFree::reinit() a collection of finite elements also containing "
7298 "FE_Nothing, please check - before creating FEEvaluation - the category "
7299 "of the current range by calling either "
7300 "MatrixFree::get_cell_range_category(range) or "
7301 "MatrixFree::get_face_range_category(range). The returned category "
7302 "is the index of the active FE, which you can use to exclude "
7309 static_cast<unsigned int>(fe_degree) !=
7310 this->data->data.front().fe_degree) ||
7311 n_q_points != this->n_quadrature_points)
7313 std::string message =
7314 "-------------------------------------------------------\n";
7315 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7316 message +=
" Called --> FEEvaluation<dim,";
7320 message +=
",Number>(data";
7336 if (
static_cast<unsigned int>(fe_degree) ==
7337 this->data->data.front().fe_degree)
7339 proposed_dof_comp = dof_no;
7340 proposed_fe_comp = first_selected_component;
7343 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
7345 for (
unsigned int nf = 0;
7348 if (this->matrix_free
7349 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7351 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7353 proposed_dof_comp = no;
7354 proposed_fe_comp = nf;
7358 this->mapping_data->descriptor[this->active_quad_index]
7360 proposed_quad_comp = this->quad_no;
7362 for (
unsigned int no = 0;
7367 .descriptor[this->active_quad_index]
7368 .n_q_points == n_q_points)
7370 proposed_quad_comp = no;
7377 if (proposed_dof_comp != first_selected_component)
7378 message +=
"Wrong vector component selection:\n";
7380 message +=
"Wrong quadrature formula selection:\n";
7381 message +=
" Did you mean FEEvaluation<dim,";
7385 message +=
",Number>(data";
7394 std::string correct_pos;
7395 if (proposed_dof_comp != dof_no)
7396 correct_pos =
" ^ ";
7399 if (proposed_quad_comp != this->quad_no)
7400 correct_pos +=
" ^ ";
7403 if (proposed_fe_comp != first_selected_component)
7404 correct_pos +=
" ^\n";