16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
92 typename VectorizedArrayType>
99 std::conditional_t<n_components_ == 1,
106 n_components_ == dim,
113 n_components_ == dim,
118 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
156 template <
typename VectorType>
159 const VectorType &src,
160 const unsigned int first_index = 0,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
191 template <
typename VectorType>
194 const VectorType &src,
195 const unsigned int first_index = 0,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
229 template <
typename VectorType>
233 const unsigned int first_index = 0,
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
274 template <
typename VectorType>
278 const unsigned int first_index = 0,
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
284 template <
typename VectorType>
288 const unsigned int first_index = 0,
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
365 typename = std::enable_if_t<n_components == n_components_local>>
368 const unsigned int q_point);
419 template <
int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
423 const unsigned int q_point);
442 const unsigned int q_point);
518 const unsigned int q_point);
527 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
546 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
549 const unsigned int q_point);
559 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
578 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
582 const unsigned int q_point);
592 template <
int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
610 template <
int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
614 const unsigned int q_point);
655 const unsigned int dof_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
663 const unsigned int face_type);
737 template <
typename VectorType,
typename VectorOperation>
741 const std::array<VectorType *, n_components_> &vectors,
744 n_components_> &vectors_sm,
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints =
true)
const;
755 template <
typename VectorType,
typename VectorOperation>
759 const std::array<VectorType *, n_components_> &vectors,
762 n_components_> &vectors_sm,
763 const std::bitset<n_lanes> &mask)
const;
772 template <
typename VectorType,
typename VectorOperation>
776 const std::array<VectorType *, n_components_> &vectors)
const;
1380 typename VectorizedArrayType>
1385 VectorizedArrayType>
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1530 const std::pair<unsigned int, unsigned int> &range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1642 template <
bool level_dof_access>
1664 const unsigned int given_n_q_points_1d);
1707 template <
typename VectorType>
1737 VectorizedArrayType *values_array,
1738 const bool sum_into_values =
false);
1753 template <
typename VectorType>
1756 VectorType &output_vector);
1840 int n_q_points_1d = fe_degree + 1,
1841 int n_components_ = 1,
1842 typename Number = double,
1848 VectorizedArrayType>
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
2012 const std::pair<unsigned int, unsigned int> &range,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2029 reinit(
const unsigned int face_batch_number);
2039 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2046 const unsigned int given_n_q_points_1d);
2110 template <
typename VectorType>
2126 const bool sum_into_values =
false);
2139 VectorizedArrayType *values_array,
2140 const bool sum_into_values =
false);
2157 const bool sum_into_values =
false);
2165 VectorizedArrayType *values_array,
2166 const bool sum_into_values =
false);
2179 template <
typename VectorType>
2182 VectorType &output_vector);
2187 template <
typename VectorType>
2190 const bool integrate_gradients,
2191 VectorType &output_vector);
2267 namespace MatrixFreeFunctions
2271 template <
int dim,
int degree>
2280 template <
int degree>
2283 static constexpr unsigned int value = degree + 1;
2299 template <
bool is_face,
2302 typename VectorizedArrayType>
2305 extract_initialization_data(
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2312 const unsigned int active_fe_index_given,
2313 const unsigned int active_quad_index_given,
2314 const unsigned int face_type)
2317 InitializationData init_data;
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2330 active_fe_index_given :
2336 active_quad_index_given :
2337 std::min<unsigned int>(
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2369 typename VectorizedArrayType>
2374 VectorizedArrayType>::
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2387 internal::extract_initialization_data<is_face>(matrix_free,
2389 first_selected_component,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2402 this->set_data_pointers(scratch_data_array, n_components_);
2404 this->dof_info->start_components.back() == 1 ||
2405 static_cast<int>(n_components_) <=
2407 this->dof_info->start_components
2408 [this->dof_info->component_to_base_index[first_selected_component] +
2410 first_selected_component,
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2417 this->dof_info->start_components
2418 [this->dof_info->component_to_base_index[first_selected_component] +
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) +
")"));
2438 typename VectorizedArrayType>
2443 VectorizedArrayType>::
2449 const unsigned int first_selected_component,
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2461 first_selected_component)
2462 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2463 , matrix_free(nullptr)
2465 const unsigned int base_element_number =
2469 first_selected_component >=
2471 ExcMessage(
"The underlying element must at least contain as many "
2472 "components as requested by this class"));
2473 (void)base_element_number;
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2479 fe.component_to_base_index(first_selected_component).
first);
2481 this->set_data_pointers(scratch_data_array, n_components_);
2490 typename VectorizedArrayType>
2495 VectorizedArrayType>::
2500 VectorizedArrayType> &other)
2502 , scratch_data_array(other.matrix_free == nullptr ?
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(other.matrix_free)
2507 if (other.matrix_free ==
nullptr)
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2521 if constexpr (is_face ==
false)
2522 this->mapping_data = &this->mapped_geometry->get_data_storage();
2526 "face evaluators is not currently "
2532 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2534 this->mapped_geometry->get_data_storage().JxW_values.begin();
2535 this->jacobian_gradients =
2536 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2537 this->jacobian_gradients_non_inverse =
2538 this->mapped_geometry->get_data_storage()
2539 .jacobian_gradients_non_inverse[0]
2542 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2545 this->set_data_pointers(scratch_data_array, n_components_);
2554 typename VectorizedArrayType>
2559 VectorizedArrayType> &
2565 VectorizedArrayType> &other)
2568 if (matrix_free ==
nullptr)
2571 delete scratch_data_array;
2580 matrix_free = other.matrix_free;
2582 if (other.matrix_free ==
nullptr)
2590 this->mapped_geometry =
2591 std::make_shared<internal::MatrixFreeFunctions::
2592 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2597 if constexpr (is_face ==
false)
2598 this->mapping_data = &this->mapped_geometry->get_data_storage();
2602 "face evaluators is not currently "
2607 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2609 this->mapped_geometry->get_data_storage().JxW_values.begin();
2610 this->jacobian_gradients =
2611 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2612 this->jacobian_gradients_non_inverse =
2613 this->mapped_geometry->get_data_storage()
2614 .jacobian_gradients_non_inverse[0]
2617 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2624 this->set_data_pointers(scratch_data_array, n_components_);
2635 typename VectorizedArrayType>
2640 VectorizedArrayType>::~FEEvaluationBase()
2642 if (matrix_free !=
nullptr)
2653 delete scratch_data_array;
2664 typename VectorizedArrayType>
2669 Assert(matrix_free !=
nullptr,
2671 "FEEvaluation was not initialized with a MatrixFree object!"));
2672 return *matrix_free;
2681 template <
typename VectorType,
bool>
2682 struct ConstBlockVectorSelector;
2684 template <
typename VectorType>
2685 struct ConstBlockVectorSelector<
VectorType, true>
2687 using BaseVectorType =
const typename VectorType::BlockType;
2690 template <
typename VectorType>
2691 struct ConstBlockVectorSelector<
VectorType, false>
2693 using BaseVectorType =
typename VectorType::BlockType;
2699 template <
typename VectorType,
bool>
2700 struct BlockVectorSelector;
2702 template <
typename VectorType>
2705 using BaseVectorType =
typename ConstBlockVectorSelector<
2707 std::is_const_v<VectorType>>::BaseVectorType;
2709 static BaseVectorType *
2710 get_vector_component(VectorType &vec,
const unsigned int component)
2713 return &vec.block(component);
2717 template <
typename VectorType>
2718 struct BlockVectorSelector<
VectorType, false>
2722 static BaseVectorType *
2723 get_vector_component(VectorType &vec,
const unsigned int component)
2742 template <
typename VectorType>
2743 struct BlockVectorSelector<
std::vector<VectorType>, false>
2747 static BaseVectorType *
2748 get_vector_component(std::vector<VectorType> &vec,
2749 const unsigned int component)
2752 return &vec[component];
2756 template <
typename VectorType>
2757 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2761 static const BaseVectorType *
2762 get_vector_component(
const std::vector<VectorType> &vec,
2763 const unsigned int component)
2766 return &vec[component];
2770 template <
typename VectorType>
2771 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2775 static BaseVectorType *
2776 get_vector_component(std::vector<VectorType *> &vec,
2777 const unsigned int component)
2780 return vec[component];
2784 template <
typename VectorType>
2785 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2789 static const BaseVectorType *
2790 get_vector_component(
const std::vector<VectorType *> &vec,
2791 const unsigned int component)
2794 return vec[component];
2798 template <
typename VectorType, std::
size_t N>
2799 struct BlockVectorSelector<
std::array<VectorType *, N>, false>
2803 static BaseVectorType *
2804 get_vector_component(std::array<VectorType *, N> &vec,
2805 const unsigned int component)
2808 return vec[component];
2819 typename VectorizedArrayType>
2820template <
typename VectorType,
typename VectorOperation>
2825 const std::array<VectorType *, n_components_> &src,
2828 n_components_> &src_sm,
2829 const std::bitset<n_lanes> &
mask,
2830 const bool apply_constraints)
const
2835 if (this->matrix_free ==
nullptr)
2837 read_write_operation_global(operation, src);
2844 if (this->n_fe_components == 1)
2845 for (
unsigned int comp = 0; comp < n_components; ++comp)
2847 Assert(src[comp] !=
nullptr,
2848 ExcMessage(
"The finite element underlying this FEEvaluation "
2849 "object is scalar, but you requested " +
2850 std::to_string(n_components) +
2851 " components via the template argument in "
2852 "FEEvaluation. In that case, you must pass an "
2853 "std::vector<VectorType> or a BlockVector to " +
2854 "read_dof_values and distribute_local_to_global."));
2866 const bool accesses_exterior_dofs =
2867 this->dof_access_index ==
2869 this->is_interior_face() ==
false;
2879 bool is_contiguous =
true;
2881 if (accesses_exterior_dofs)
2883 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2884 const unsigned int n_filled_lanes =
2889 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2890 if (
mask[v] ==
true &&
2893 [cells[v] / n_lanes] <
2896 is_contiguous = false;
2900 this->dof_access_index :
2901 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2902 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2905 is_contiguous =
false;
2910 read_write_operation_contiguous(operation, src, src_sm,
mask);
2917 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2919 const bool masking_is_active =
mask.count() < n_lanes;
2920 if (masking_is_active)
2921 for (
unsigned int v = 0; v < n_lanes; ++v)
2922 if (
mask[v] ==
false)
2925 bool has_hn_constraints =
false;
2927 if (is_face ==
false)
2933 [this->first_selected_component])
2934 for (
unsigned int v = 0; v < n_lanes; ++v)
2939 has_hn_constraints = true;
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2945 const bool use_vectorized_path =
2946 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2948 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
2949 std::array<VectorizedArrayType *, n_components> values_dofs;
2950 for (
unsigned int c = 0; c < n_components; ++c)
2951 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2952 c * dofs_per_component;
2956 [is_face ? this->dof_access_index :
2957 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2958 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2959 IndexStorageVariants::interleaved &&
2960 use_vectorized_path)
2962 const unsigned int *dof_indices =
2964 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2968 [this->first_selected_component] *
2971 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2972 if (n_components == 1 || this->n_fe_components == 1)
2973 for (
unsigned int comp = 0; comp < n_components; ++comp)
2975 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2978 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2980 if (n_components == 1 || this->n_fe_components == 1)
2981 for (
unsigned int i = 0; i < dofs_per_component;
2982 ++i, dof_indices += n_lanes)
2983 for (
unsigned int comp = 0; comp < n_components; ++comp)
2984 operation.process_dof_gather(dof_indices,
2988 values_dofs[comp][i],
2991 for (
unsigned int comp = 0; comp < n_components; ++comp)
2992 for (
unsigned int i = 0; i < dofs_per_component;
2993 ++i, dof_indices += n_lanes)
2994 operation.process_dof_gather(dof_indices,
2998 values_dofs[comp][i],
3005 std::array<const unsigned int *, n_lanes> dof_indices;
3006 dof_indices.fill(
nullptr);
3011 bool has_constraints =
false;
3012 const unsigned int n_components_read =
3013 this->n_fe_components > 1 ? n_components : 1;
3017 for (
unsigned int v = 0; v < n_lanes; ++v)
3023 const std::pair<unsigned int, unsigned int> *my_index_start =
3024 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3025 this->first_selected_component];
3030 if (my_index_start[n_components_read].
second !=
3031 my_index_start[0].
second)
3032 has_constraints =
true;
3035 dof_info.
dof_indices.data() + my_index_start[0].first;
3040 for (
unsigned int v = 0; v < n_lanes; ++v)
3045 const std::pair<unsigned int, unsigned int> *my_index_start =
3046 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3047 this->first_selected_component];
3048 if (my_index_start[n_components_read].
second !=
3049 my_index_start[0].
second)
3050 has_constraints =
true;
3057 dof_info.hanging_node_constraint_masks_comp
3058 [this->active_fe_index][this->first_selected_component])
3059 has_hn_constraints = true;
3062 my_index_start[0].
first ||
3065 my_index_start[0].
first,
3068 dof_info.
dof_indices.data() + my_index_start[0].first;
3072 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3073 return i != numbers::invalid_unsigned_int;
3075 for (
unsigned int comp = 0; comp < n_components; ++comp)
3076 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3077 operation.process_empty(values_dofs[comp][i]);
3081 if (!has_constraints && apply_constraints)
3083 if (n_components == 1 || this->n_fe_components == 1)
3085 for (
unsigned int v = 0; v < n_lanes; ++v)
3090 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3091 for (
unsigned int comp = 0; comp < n_components; ++comp)
3092 operation.process_dof(dof_indices[v][i],
3094 values_dofs[comp][i][v]);
3099 for (
unsigned int comp = 0; comp < n_components; ++comp)
3100 for (
unsigned int v = 0; v < n_lanes; ++v)
3105 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3106 operation.process_dof(
3107 dof_indices[v][comp * dofs_per_component + i],
3109 values_dofs[comp][i][v]);
3121 for (
unsigned int v = 0; v < n_lanes; ++v)
3127 const unsigned int cell_dof_index =
3128 cell_index * this->n_fe_components + this->first_selected_component;
3129 const unsigned int n_components_read =
3130 this->n_fe_components > 1 ? n_components : 1;
3131 unsigned int index_indicators =
3133 unsigned int next_index_indicators =
3134 dof_info.
row_starts[cell_dof_index + 1].second;
3138 if (apply_constraints ==
false &&
3139 (dof_info.
row_starts[cell_dof_index].second !=
3140 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3146 dof_info.hanging_node_constraint_masks_comp
3147 [this->active_fe_index][this->first_selected_component])))
3156 [this->first_selected_component] +
3158 next_index_indicators = index_indicators;
3161 if (n_components == 1 || this->n_fe_components == 1)
3163 unsigned int ind_local = 0;
3164 for (; index_indicators != next_index_indicators; ++index_indicators)
3166 const std::pair<unsigned short, unsigned short> indicator =
3169 for (
unsigned int j = 0; j < indicator.first; ++j)
3170 for (
unsigned int comp = 0; comp < n_components; ++comp)
3171 operation.process_dof(dof_indices[v][j],
3173 values_dofs[comp][ind_local + j][v]);
3175 ind_local += indicator.first;
3176 dof_indices[v] += indicator.first;
3180 Number
value[n_components];
3181 for (
unsigned int comp = 0; comp < n_components; ++comp)
3182 operation.pre_constraints(values_dofs[comp][ind_local][v],
3185 const Number *data_val =
3187 const Number *end_pool =
3189 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3190 for (
unsigned int comp = 0; comp < n_components; ++comp)
3191 operation.process_constraint(*dof_indices[v],
3196 for (
unsigned int comp = 0; comp < n_components; ++comp)
3197 operation.post_constraints(
value[comp],
3198 values_dofs[comp][ind_local][v]);
3204 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3205 for (
unsigned int comp = 0; comp < n_components; ++comp)
3206 operation.process_dof(*dof_indices[v],
3208 values_dofs[comp][ind_local][v]);
3217 for (
unsigned int comp = 0; comp < n_components; ++comp)
3219 unsigned int ind_local = 0;
3222 for (; index_indicators != next_index_indicators;
3225 const std::pair<unsigned short, unsigned short> indicator =
3229 for (
unsigned int j = 0; j < indicator.first; ++j)
3230 operation.process_dof(dof_indices[v][j],
3232 values_dofs[comp][ind_local + j][v]);
3233 ind_local += indicator.first;
3234 dof_indices[v] += indicator.first;
3239 operation.pre_constraints(values_dofs[comp][ind_local][v],
3242 const Number *data_val =
3244 const Number *end_pool =
3247 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3248 operation.process_constraint(*dof_indices[v],
3253 operation.post_constraints(
value,
3254 values_dofs[comp][ind_local][v]);
3261 for (; ind_local < dofs_per_component;
3262 ++dof_indices[v], ++ind_local)
3265 operation.process_dof(*dof_indices[v],
3267 values_dofs[comp][ind_local][v]);
3270 if (apply_constraints ==
true && comp + 1 < n_components)
3271 next_index_indicators =
3272 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3284 typename VectorizedArrayType>
3285template <
typename VectorType,
typename VectorOperation>
3290 const std::array<VectorType *, n_components_> &src)
const
3294 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3295 unsigned int index = this->first_selected_component * dofs_per_component;
3296 for (
unsigned int comp = 0; comp < n_components; ++comp)
3298 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3300 operation.process_empty(
3301 this->values_dofs[comp * dofs_per_component + i]);
3302 operation.process_dof_global(
3303 local_dof_indices[this->
data->lexicographic_numbering[index]],
3305 this->values_dofs[comp * dofs_per_component + i][0]);
3316 typename VectorizedArrayType>
3317template <
typename VectorType,
typename VectorOperation>
3322 const std::array<VectorType *, n_components_> &src,
3325 n_components_> &vectors_sm,
3326 const std::bitset<n_lanes> &
mask)
const
3335 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3338 is_face ? this->dof_access_index :
3340 const unsigned int n_active_lanes =
mask.count();
3343 const std::vector<unsigned int> &dof_indices_cont =
3346 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3347 std::array<VectorizedArrayType *, n_components> values_dofs{{
nullptr}};
3348 for (
unsigned int c = 0; c < n_components; ++c)
3349 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3350 c * dofs_per_component;
3354 const bool accesses_exterior_dofs =
3355 this->dof_access_index ==
3357 this->is_interior_face() ==
false;
3363 interleaved_contiguous &&
3364 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3366 const unsigned int dof_index =
3367 dof_indices_cont[this->cell * n_lanes] +
3370 [this->first_selected_component] *
3372 if (n_components == 1 || this->n_fe_components == 1)
3373 for (
unsigned int comp = 0; comp < n_components; ++comp)
3374 operation.process_dofs_vectorized(dofs_per_component,
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3388 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3392 const unsigned int n_filled_lanes =
3395 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3396 n_active_lanes == n_lanes &&
3397 !accesses_exterior_dofs;
3399 if (vectors_sm[0] !=
nullptr)
3401 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3402 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3405 const auto upper_bound =
3406 std::min<unsigned int>(n_filled_lanes, n_lanes);
3407 for (
unsigned int v = 0; v < upper_bound; ++v)
3409 if (
mask[v] ==
false)
3411 vector_ptrs[v] =
nullptr;
3431 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3432 vectors_sm[comp]->operator[](temp.first).
data() + temp.second +
3434 [this->active_fe_index][this->first_selected_component]);
3436 vector_ptrs[v] =
nullptr;
3438 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3439 vector_ptrs[v] =
nullptr;
3444 if (use_vectorized_path)
3446 if (n_components == 1 || this->n_fe_components == 1)
3448 for (
unsigned int comp = 0; comp < n_components; ++comp)
3450 auto vector_ptrs = compute_vector_ptrs(comp);
3451 operation.process_dofs_vectorized_transpose(
3460 auto vector_ptrs = compute_vector_ptrs(0);
3461 operation.process_dofs_vectorized_transpose(dofs_per_component *
3469 for (
unsigned int comp = 0; comp < n_components; ++comp)
3471 auto vector_ptrs = compute_vector_ptrs(
3472 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3474 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3475 operation.process_empty(values_dofs[comp][i]);
3477 if (n_components == 1 || this->n_fe_components == 1)
3479 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3480 if (
mask[v] ==
true)
3481 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3482 operation.process_dof(vector_ptrs[v][i],
3483 values_dofs[comp][i][v]);
3487 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3488 if (
mask[v] ==
true)
3489 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3490 operation.process_dof(
3491 vector_ptrs[v][i + comp * dofs_per_component],
3492 values_dofs[comp][i][v]);
3498 std::array<unsigned int, n_lanes> dof_indices{
3501 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3505 if (
mask[v] ==
true)
3507 dof_indices_cont[cells[v]] +
3510 [this->first_selected_component] *
3516 if (use_vectorized_path)
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (
unsigned int comp = 0; comp < n_components; ++comp)
3524 operation.process_dofs_vectorized_transpose(dofs_per_component,
3530 operation.process_dofs_vectorized_transpose(dofs_per_component *
3539 interleaved_contiguous_strided)
3541 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (
unsigned int comp = 0; comp < n_components; ++comp)
3545 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3546 src[comp]->
begin());
3549 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3551 if (n_components == 1 || this->n_fe_components == 1)
3552 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3554 for (
unsigned int comp = 0; comp < n_components; ++comp)
3555 operation.process_dof_gather(dof_indices.data(),
3558 src_ptrs[comp] + i * n_lanes,
3559 values_dofs[comp][i],
3563 for (
unsigned int comp = 0; comp < n_components; ++comp)
3564 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3566 operation.process_dof_gather(
3569 (comp * dofs_per_component + i) * n_lanes,
3570 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3571 values_dofs[comp][i],
3579 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3581 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3583 if (n_components == 1 || this->n_fe_components == 1)
3584 for (
unsigned int comp = 0; comp < n_components; ++comp)
3585 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3586 src[comp]->
begin());
3589 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3591 const unsigned int *offsets =
3593 if (n_components == 1 || this->n_fe_components == 1)
3594 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3596 for (
unsigned int comp = 0; comp < n_components; ++comp)
3597 operation.process_dof_gather(dof_indices.data(),
3601 values_dofs[comp][i],
3604 for (
unsigned int v = 0; v < n_lanes; ++v)
3605 dof_indices[v] += offsets[v];
3608 for (
unsigned int comp = 0; comp < n_components; ++comp)
3609 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3611 operation.process_dof_gather(dof_indices.data(),
3615 values_dofs[comp][i],
3618 for (
unsigned int v = 0; v < n_lanes; ++v)
3619 dof_indices[v] += offsets[v];
3624 for (
unsigned int comp = 0; comp < n_components; ++comp)
3626 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3627 operation.process_empty(values_dofs[comp][i]);
3628 if (accesses_exterior_dofs)
3630 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3631 if (
mask[v] ==
true)
3634 [ind][cells[v] / VectorizedArrayType::size()] ==
3638 if (n_components == 1 || this->n_fe_components == 1)
3640 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3641 operation.process_dof(dof_indices[v] + i,
3643 values_dofs[comp][i][v]);
3647 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3651 values_dofs[comp][i][v]);
3656 const unsigned int offset =
3659 if (n_components == 1 || this->n_fe_components == 1)
3661 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3662 operation.process_dof(dof_indices[v] + i * offset,
3664 values_dofs[comp][i][v]);
3668 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3669 operation.process_dof(
3671 (i + comp * dofs_per_component) * offset,
3673 values_dofs[comp][i][v]);
3684 if (n_components == 1 || this->n_fe_components == 1)
3686 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3687 if (
mask[v] ==
true)
3688 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3689 operation.process_dof(dof_indices[v] + i,
3691 values_dofs[comp][i][v]);
3695 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3696 if (
mask[v] ==
true)
3697 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3698 operation.process_dof(dof_indices[v] + i +
3699 comp * dofs_per_component,
3701 values_dofs[comp][i][v]);
3706 const unsigned int *offsets =
3708 [ind][VectorizedArrayType::size() * this->cell];
3709 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3711 if (n_components == 1 || this->n_fe_components == 1)
3712 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3714 if (
mask[v] ==
true)
3715 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3716 operation.process_dof(dof_indices[v] + i * offsets[v],
3718 values_dofs[comp][i][v]);
3722 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3723 if (
mask[v] ==
true)
3724 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3725 operation.process_dof(
3727 (i + comp * dofs_per_component) * offsets[v],
3729 values_dofs[comp][i][v]);
3741 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3742 decltype(std::declval<VectorType>().begin())
3743 get_beginning(VectorType &vec)
3751 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3752 typename VectorType::value_type *
3753 get_beginning(VectorType &)
3759 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3761 const std::vector<ArrayView<const typename VectorType::value_type>> *
3762 get_shared_vector_data(VectorType *vec,
3763 const bool is_valid_mode_for_sm,
3764 const unsigned int active_fe_index,
3768 if (is_valid_mode_for_sm &&
3771 active_fe_index == 0)
3772 return &vec->shared_vector_data();
3778 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3780 const std::vector<ArrayView<const typename VectorType::value_type>> *
3781 get_shared_vector_data(VectorType *,
3789 template <
int n_components,
typename VectorType>
3791 std::array<
typename internal::BlockVectorSelector<
3796 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3800 get_vector_data(VectorType &src,
3801 const unsigned int first_index,
3802 const bool is_valid_mode_for_sm,
3803 const unsigned int active_fe_index,
3809 std::array<
typename internal::BlockVectorSelector<
3815 ArrayView<
const typename internal::BlockVectorSelector<
3821 for (
unsigned int d = 0;
d < n_components; ++
d)
3822 src_data.first[d] = internal::BlockVectorSelector<
3828 for (
unsigned int d = 0;
d < n_components; ++
d)
3829 src_data.second[d] = get_shared_vector_data(
3830 const_cast<typename internal::BlockVectorSelector<
3831 std::remove_const_t<VectorType>,
3833 *>(src_data.first[d]),
3834 is_valid_mode_for_sm,
3848 typename VectorizedArrayType>
3853 if (this->dof_info ==
nullptr ||
3855 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3856 this->dof_info->hanging_node_constraint_masks_comp
3857 [this->active_fe_index][this->first_selected_component] ==
false)
3860 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3861 constraint_mask{{internal::MatrixFreeFunctions::
3862 unconstrained_compressed_constraint_kind}};
3864 bool hn_available =
false;
3866 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3868 for (
unsigned int v = 0; v < n_lanes; ++v)
3880 constraint_mask[v] =
mask;
3886 if (hn_available ==
false)
3891 this->
data->data.front().fe_degree,
3892 this->get_shape_info(),
3904 typename VectorizedArrayType>
3905template <
typename VectorType>
3909 const unsigned int first_index,
3910 const std::bitset<n_lanes> &
mask)
3912 const auto src_data = internal::get_vector_data<n_components_>(
3915 this->dof_access_index ==
3917 this->active_fe_index,
3921 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3923 apply_hanging_node_constraints(
false);
3926 this->dof_values_initialized =
true;
3936 typename VectorizedArrayType>
3937template <
typename VectorType>
3941 const unsigned int first_index,
3942 const std::bitset<n_lanes> &
mask)
3944 const auto src_data = internal::get_vector_data<n_components_>(
3947 this->dof_access_index ==
3949 this->active_fe_index,
3953 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3956 this->dof_values_initialized =
true;
3966 typename VectorizedArrayType>
3967template <
typename VectorType>
3971 const unsigned int first_index,
3972 const std::bitset<n_lanes> &
mask)
const
3975 Assert(this->dof_values_initialized ==
true,
3979 apply_hanging_node_constraints(
true);
3981 const auto dst_data = internal::get_vector_data<n_components_>(
3984 this->dof_access_index ==
3986 this->active_fe_index,
3991 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
4000 typename VectorizedArrayType>
4001template <
typename VectorType>
4005 const unsigned int first_index,
4006 const std::bitset<n_lanes> &
mask)
const
4009 Assert(this->dof_values_initialized ==
true,
4013 const auto dst_data = internal::get_vector_data<n_components_>(
4016 this->dof_access_index ==
4018 this->active_fe_index,
4022 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4031 typename VectorizedArrayType>
4032template <
typename VectorType>
4036 const unsigned int first_index,
4037 const std::bitset<n_lanes> &
mask)
const
4040 Assert(this->dof_values_initialized ==
true,
4044 const auto dst_data = internal::get_vector_data<n_components_>(
4047 this->dof_access_index ==
4049 this->active_fe_index,
4053 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4066 typename VectorizedArrayType>
4072 VectorizedArrayType>::value_type
4077 if constexpr (n_components == 1)
4078 return this->values_dofs[dof];
4081 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4082 Tensor<1, n_components_, VectorizedArrayType> return_value;
4083 for (
unsigned int comp = 0; comp < n_components; ++comp)
4084 return_value[comp] = this->values_dofs[comp * dofs + dof];
4085 return return_value;
4095 typename VectorizedArrayType>
4101 VectorizedArrayType>::value_type
4103 get_value(
const unsigned int q_point)
const
4106 Assert(this->values_quad_initialized ==
true,
4111 if constexpr (n_components == 1)
4112 return this->values_quad[q_point];
4115 if (n_components == dim &&
4116 this->
data->element_type ==
4121 Assert(this->values_quad_initialized ==
true,
4126 Assert(this->J_value !=
nullptr,
4129 const std::size_t nqp = this->n_quadrature_points;
4137 const VectorizedArrayType inv_det =
4138 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4139 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4140 this->jacobian[0][2][2];
4143 for (
unsigned int comp = 0; comp < n_components; ++comp)
4144 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4145 jac[comp][comp] * inv_det;
4152 this->jacobian[q_point] :
4161 const VectorizedArrayType inv_det =
4162 (is_face && dim == 2 && this->get_face_no() < 2) ?
4166 for (
unsigned int comp = 0; comp < n_components; ++comp)
4168 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4169 for (
unsigned int e = 1;
e < dim; ++
e)
4171 this->values_quad[e * nqp + q_point] * jac[comp][e];
4172 value_out[comp] *= inv_det;
4179 const std::size_t nqp = this->n_quadrature_points;
4181 for (
unsigned int comp = 0; comp < n_components; ++comp)
4182 return_value[comp] = this->values_quad[comp * nqp + q_point];
4183 return return_value;
4194 typename VectorizedArrayType>
4200 VectorizedArrayType>::gradient_type
4205 Assert(this->gradients_quad_initialized ==
true,
4210 Assert(this->jacobian !=
nullptr,
4212 "update_gradients"));
4213 const std::size_t nqp = this->n_quadrature_points;
4215 if constexpr (n_components == dim && dim > 1)
4217 if (this->
data->element_type ==
4222 Assert(this->gradients_quad_initialized ==
true,
4227 Assert(this->jacobian !=
nullptr,
4229 "update_gradients"));
4230 const std::size_t nqp = this->n_quadrature_points;
4231 const std::size_t nqp_d = nqp * dim;
4234 this->gradients_quad + q_point * dim;
4245 const VectorizedArrayType inv_det =
4246 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4247 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4248 this->jacobian[0][2][2];
4251 for (
unsigned int d = 0;
d < dim; ++
d)
4252 for (
unsigned int comp = 0; comp < n_components; ++comp)
4253 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4255 (jac[comp][comp] * inv_det);
4267 const VectorizedArrayType inv_det =
4268 (is_face && dim == 2 && this->get_face_no() < 2) ?
4272 VectorizedArrayType tmp[dim][dim];
4274 for (
unsigned int d = 0;
d < dim; ++
d)
4275 for (
unsigned int e = 0;
e < dim; ++
e)
4278 for (
unsigned int f = 1; f < dim; ++f)
4279 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4281 for (
unsigned int comp = 0; comp < n_components; ++comp)
4282 for (
unsigned int d = 0;
d < dim; ++
d)
4284 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4285 for (
unsigned int f = 1; f < dim; ++f)
4286 res += jac[comp][f] * tmp[d][f];
4288 grad_out[comp][
d] = res * inv_det;
4299 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4301 "update_hessians"));
4303 const auto jac_grad =
4304 this->jacobian_gradients_non_inverse[q_point];
4306 this->jacobian[q_point];
4310 const VectorizedArrayType inv_det =
4311 (is_face && dim == 2 && this->get_face_no() < 2) ?
4318 VectorizedArrayType tmp[dim][dim];
4319 for (
unsigned int d = 0;
d < dim; ++
d)
4320 for (
unsigned int e = 0;
e < dim; ++
e)
4323 for (
unsigned int f = 1; f < dim; ++f)
4324 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4329 for (
unsigned int d = 0;
d < dim; ++
d)
4331 for (
unsigned int e = 0;
e < dim; ++
e)
4333 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4334 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4335 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4338 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4340 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4345 for (
unsigned int d = 0;
d < dim; ++
d)
4346 for (
unsigned int e = 0;
e < dim; ++
e)
4348 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4349 for (
unsigned int f = 1; f < dim; ++f)
4350 res += tmp[f][d] * inv_t_jac[e][f];
4351 grad_out[
d][
e] = res;
4357 VectorizedArrayType tmp3[dim], tmp4[dim];
4358 for (
unsigned int d = 0;
d < dim; ++
d)
4360 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4361 for (
unsigned int e = 1;
e < dim; ++
e)
4362 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4364 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4365 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4366 for (
unsigned int d = 0;
d < dim; ++
d)
4368 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4369 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4371 for (
unsigned int d = 0;
d < dim; ++
d)
4373 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4374 for (
unsigned int e = 1;
e < dim; ++
e)
4375 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4378 VectorizedArrayType tmp2[dim];
4379 for (
unsigned int d = 0;
d < dim; ++
d)
4381 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4382 for (
unsigned e = 1;
e < dim; ++
e)
4384 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4387 for (
unsigned int d = 0;
d < dim; ++
d)
4388 for (
unsigned int e = 0;
e < dim; ++
e)
4390 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4394 grad_out[
d][
e] *= inv_det;
4405 for (
unsigned int comp = 0; comp < n_components; ++comp)
4406 for (
unsigned int d = 0;
d < dim; ++
d)
4408 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4409 this->jacobian[0][
d][
d];
4418 for (
unsigned int comp = 0; comp < n_components; ++comp)
4419 for (
unsigned int d = 0;
d < dim; ++
d)
4422 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4423 for (
unsigned int e = 1;
e < dim; ++
e)
4424 grad_out[comp][d] +=
4426 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4429 if constexpr (n_components == 1)
4441 typename VectorizedArrayType>
4447 VectorizedArrayType>::value_type
4453 Assert(this->gradients_quad_initialized ==
true,
4457 Assert(this->normal_x_jacobian !=
nullptr,
4459 "update_gradients"));
4461 const std::size_t nqp = this->n_quadrature_points;
4465 for (
unsigned int comp = 0; comp < n_components; ++comp)
4467 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4468 (this->normal_x_jacobian[0][dim - 1]);
4471 const std::size_t
index =
4473 for (
unsigned int comp = 0; comp < n_components; ++comp)
4475 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4476 this->normal_x_jacobian[
index][0];
4477 for (
unsigned int d = 1;
d < dim; ++
d)
4479 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4480 this->normal_x_jacobian[
index][
d];
4483 if constexpr (n_components == 1)
4495 template <
typename VectorizedArrayType>
4498 const VectorizedArrayType *
const hessians,
4500 VectorizedArrayType (&tmp)[1][1])
4502 tmp[0][0] = jac[0][0] *
hessians[0];
4505 template <
typename VectorizedArrayType>
4508 const VectorizedArrayType *
const hessians,
4509 const unsigned int nqp,
4510 VectorizedArrayType (&tmp)[2][2])
4512 for (
unsigned int d = 0;
d < 2; ++
d)
4520 template <
typename VectorizedArrayType>
4523 const VectorizedArrayType *
const hessians,
4524 const unsigned int nqp,
4525 VectorizedArrayType (&tmp)[3][3])
4527 for (
unsigned int d = 0;
d < 3; ++
d)
4548 typename VectorizedArrayType>
4553 VectorizedArrayType>::hessian_type
4558 Assert(this->hessians_quad_initialized ==
true,
4563 Assert(this->jacobian !=
nullptr,
4573 const std::size_t nqp = this->n_quadrature_points;
4574 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4579 for (
unsigned int comp = 0; comp < n_components; ++comp)
4581 for (
unsigned int d = 0;
d < dim; ++
d)
4582 hessian_out[comp][d][d] =
4583 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4584 (jac[
d][
d] * jac[
d][
d]);
4590 hessian_out[comp][0][1] =
4591 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4592 (jac[0][0] * jac[1][1]);
4595 hessian_out[comp][0][1] =
4596 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4597 (jac[0][0] * jac[1][1]);
4598 hessian_out[comp][0][2] =
4599 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4600 (jac[0][0] * jac[2][2]);
4601 hessian_out[comp][1][2] =
4602 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4603 (jac[1][1] * jac[2][2]);
4608 for (
unsigned int d = 0;
d < dim; ++
d)
4609 for (
unsigned int e = d + 1;
e < dim; ++
e)
4610 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4616 for (
unsigned int comp = 0; comp < n_components; ++comp)
4618 VectorizedArrayType tmp[dim][dim];
4619 internal::hessian_unit_times_jac(
4620 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4623 for (
unsigned int d = 0;
d < dim; ++
d)
4624 for (
unsigned int e = d;
e < dim; ++
e)
4626 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4627 for (
unsigned int f = 1; f < dim; ++f)
4628 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4635 for (
unsigned int d = 0;
d < dim; ++
d)
4636 for (
unsigned int e = d + 1;
e < dim; ++
e)
4637 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4643 const auto &jac_grad = this->jacobian_gradients[q_point];
4644 for (
unsigned int comp = 0; comp < n_components; ++comp)
4646 VectorizedArrayType tmp[dim][dim];
4647 internal::hessian_unit_times_jac(
4648 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4651 for (
unsigned int d = 0;
d < dim; ++
d)
4652 for (
unsigned int e = d;
e < dim; ++
e)
4654 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4655 for (
unsigned int f = 1; f < dim; ++f)
4656 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4660 for (
unsigned int d = 0;
d < dim; ++
d)
4661 for (
unsigned int e = 0;
e < dim; ++
e)
4662 hessian_out[comp][d][d] +=
4664 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4667 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4668 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4669 for (
unsigned int f = 0; f < dim; ++f)
4670 hessian_out[comp][d][e] +=
4671 jac_grad[count][f] *
4672 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4675 for (
unsigned int d = 0;
d < dim; ++
d)
4676 for (
unsigned int e = d + 1;
e < dim; ++
e)
4677 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4680 if constexpr (n_components == 1)
4681 return hessian_out[0];
4692 typename VectorizedArrayType>
4697 VectorizedArrayType>::gradient_type
4703 Assert(this->hessians_quad_initialized ==
true,
4714 const std::size_t nqp = this->n_quadrature_points;
4715 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4721 for (
unsigned int comp = 0; comp < n_components; ++comp)
4722 for (
unsigned int d = 0;
d < dim; ++
d)
4723 hessian_out[comp][d] =
4724 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4725 (jac[
d][
d] * jac[
d][
d]);
4730 for (
unsigned int comp = 0; comp < n_components; ++comp)
4734 VectorizedArrayType tmp[dim][dim];
4735 internal::hessian_unit_times_jac(
4736 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4740 for (
unsigned int d = 0;
d < dim; ++
d)
4742 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4743 for (
unsigned int f = 1; f < dim; ++f)
4744 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4751 const auto &jac_grad = this->jacobian_gradients[q_point];
4752 for (
unsigned int comp = 0; comp < n_components; ++comp)
4756 VectorizedArrayType tmp[dim][dim];
4757 internal::hessian_unit_times_jac(
4758 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4762 for (
unsigned int d = 0;
d < dim; ++
d)
4764 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4765 for (
unsigned int f = 1; f < dim; ++f)
4766 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4769 for (
unsigned int d = 0;
d < dim; ++
d)
4770 for (
unsigned int e = 0;
e < dim; ++
e)
4771 hessian_out[comp][d] +=
4773 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4777 if constexpr (n_components == 1)
4778 return hessian_out[0];
4789 typename VectorizedArrayType>
4794 VectorizedArrayType>::value_type
4800 Assert(this->hessians_quad_initialized ==
true,
4805 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4806 if constexpr (n_components == 1)
4808 VectorizedArrayType
sum = hess_diag[0];
4809 for (
unsigned int d = 1;
d < dim; ++
d)
4810 sum += hess_diag[d];
4816 for (
unsigned int comp = 0; comp < n_components; ++comp)
4818 laplacian_out[comp] = hess_diag[comp][0];
4819 for (
unsigned int d = 1;
d < dim; ++
d)
4820 laplacian_out[comp] += hess_diag[comp][d];
4822 return laplacian_out;
4832 typename VectorizedArrayType>
4837 VectorizedArrayType>::value_type
4842 Assert(this->hessians_quad_initialized ==
true,
4847 Assert(this->normal_x_jacobian !=
nullptr,
4849 "update_hessians"));
4853 const std::size_t nqp = this->n_quadrature_points;
4854 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4858 const auto nxj = this->normal_x_jacobian[0];
4860 for (
unsigned int comp = 0; comp < n_components; ++comp)
4862 for (
unsigned int d = 0;
d < dim; ++
d)
4863 hessian_out[comp] +=
4864 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4865 (nxj[
d]) * (nxj[d]);
4872 hessian_out[comp] +=
4873 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4877 hessian_out[comp] +=
4878 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4880 hessian_out[comp] +=
4881 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4883 hessian_out[comp] +=
4884 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4895 const auto normal = this->normal_vector(q_point);
4896 const auto hessian = get_hessian(q_point);
4898 if constexpr (n_components == 1)
4899 hessian_out[0] =
hessian * normal * normal;
4901 for (
unsigned int comp = 0; comp < n_components; ++comp)
4902 hessian_out[comp] =
hessian[comp] * normal * normal;
4904 if constexpr (n_components == 1)
4905 return hessian_out[0];
4916 typename VectorizedArrayType>
4922 this->dof_values_initialized =
true;
4924 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4926 for (
unsigned int comp = 0; comp < n_components; ++comp)
4927 if constexpr (n_components == 1)
4928 this->values_dofs[comp * dofs + dof] = val_in;
4930 this->values_dofs[comp * dofs + dof] = val_in[comp];
4939 typename VectorizedArrayType>
4942 submit_value(
const value_type val_in,
const unsigned int q_point)
4948 Assert(this->J_value !=
nullptr,
4952 this->values_quad_submitted =
true;
4955 const std::size_t nqp = this->n_quadrature_points;
4956 VectorizedArrayType *
values = this->values_quad + q_point;
4958 const VectorizedArrayType JxW =
4960 this->J_value[0] * this->quadrature_weights[q_point] :
4961 this->J_value[q_point];
4962 if constexpr (n_components == 1)
4963 values[0] = val_in * JxW;
4966 if (n_components == dim &&
4967 this->
data->element_type ==
4972 Assert(this->J_value !=
nullptr,
4977 this->values_quad_submitted =
true;
4980 VectorizedArrayType *
values = this->values_quad + q_point;
4981 const std::size_t nqp = this->n_quadrature_points;
4987 const VectorizedArrayType weight =
4988 this->quadrature_weights[q_point];
4990 for (
unsigned int comp = 0; comp < n_components; ++comp)
4991 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4998 this->jacobian[q_point] :
5003 const VectorizedArrayType fac =
5005 this->quadrature_weights[q_point] :
5007 this->J_value[q_point] :
5008 this->J_value[0] * this->quadrature_weights[q_point]) *
5009 ((dim == 2 && this->get_face_no() < 2) ?
5018 for (
unsigned int comp = 0; comp < n_components; ++comp)
5020 values[comp * nqp] = val_in[0] * jac[0][comp];
5021 for (
unsigned int e = 1;
e < dim; ++
e)
5022 values[comp * nqp] += val_in[e] * jac[e][comp];
5023 values[comp * nqp] *= fac;
5028 for (
unsigned int comp = 0; comp < n_components; ++comp)
5029 values[comp * nqp] = val_in[comp] * JxW;
5039 typename VectorizedArrayType>
5040template <
int,
typename>
5044 const unsigned int q_point)
5046 static_assert(n_components == 1,
5047 "Do not try to modify the default template parameters used for"
5048 " selectively enabling this function via std::enable_if!");
5049 submit_value(val_in[0], q_point);
5058 typename VectorizedArrayType>
5061 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
5067 Assert(this->J_value !=
nullptr,
5069 "update_gradients"));
5070 Assert(this->jacobian !=
nullptr,
5072 "update_gradients"));
5074 this->gradients_quad_submitted =
true;
5077 if constexpr (dim > 1 && n_components == dim)
5079 if (this->
data->element_type ==
5088 Assert(this->J_value !=
nullptr,
5090 "update_gradients"));
5091 Assert(this->jacobian !=
nullptr,
5093 "update_gradients"));
5095 this->gradients_quad_submitted =
true;
5098 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5099 VectorizedArrayType *
values =
5100 this->values_from_gradients_quad + q_point;
5101 const std::size_t nqp = this->n_quadrature_points;
5102 const std::size_t nqp_d = nqp * dim;
5112 const VectorizedArrayType weight =
5113 this->quadrature_weights[q_point];
5114 for (
unsigned int d = 0;
d < dim; ++
d)
5115 for (
unsigned int comp = 0; comp < n_components; ++comp)
5116 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5118 (jac[comp][comp] * weight);
5130 const VectorizedArrayType fac =
5132 this->quadrature_weights[q_point] :
5133 this->J_value[0] * this->quadrature_weights[q_point] *
5134 ((dim == 2 && this->get_face_no() < 2) ?
5139 VectorizedArrayType tmp[dim][dim];
5140 for (
unsigned int d = 0;
d < dim; ++
d)
5141 for (
unsigned int e = 0;
e < dim; ++
e)
5143 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5144 for (
unsigned int f = 1; f < dim; ++f)
5145 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5147 for (
unsigned int comp = 0; comp < n_components; ++comp)
5148 for (
unsigned int d = 0;
d < dim; ++
d)
5150 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5151 for (
unsigned int f = 1; f < dim; ++f)
5152 res += jac[f][comp] * tmp[d][f];
5161 const auto jac_grad =
5162 this->jacobian_gradients_non_inverse[q_point];
5164 this->jacobian[q_point];
5168 const VectorizedArrayType fac =
5169 (!is_face) ? this->quadrature_weights[q_point] :
5170 this->J_value[q_point] *
5171 ((dim == 2 && this->get_face_no() < 2) ?
5180 VectorizedArrayType tmp3[dim], tmp4[dim];
5181 for (
unsigned int d = 0;
d < dim; ++
d)
5183 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5184 for (
unsigned int e = 1;
e < dim; ++
e)
5185 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5187 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5188 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5189 for (
unsigned int d = 0;
d < dim; ++
d)
5191 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5192 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5194 for (
unsigned int d = 0;
d < dim; ++
d)
5196 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5197 for (
unsigned int e = 1;
e < dim; ++
e)
5198 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5204 VectorizedArrayType tmp[dim][dim];
5207 for (
unsigned int d = 0;
d < dim; ++
d)
5208 for (
unsigned int e = 0;
e < dim; ++
e)
5210 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5211 for (
unsigned int f = 1; f < dim; ++f)
5212 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5215 for (
unsigned int d = 0;
d < dim; ++
d)
5216 for (
unsigned int e = 0;
e < dim; ++
e)
5218 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5219 for (
unsigned int f = 1; f < dim; ++f)
5220 res += t_jac[d][f] * tmp[e][f];
5227 VectorizedArrayType
value[dim];
5228 for (
unsigned int d = 0;
d < dim; ++
d)
5230 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5231 for (
unsigned int e = 1;
e < dim; ++
e)
5232 value[d] += tmp[d][e] * jac_grad[d][e];
5234 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5235 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5236 for (
unsigned int d = 0;
d < dim; ++
d)
5238 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5239 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5244 for (
unsigned int d = 0;
d < dim; ++
d)
5246 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5247 for (
unsigned int e = 1;
e < dim; ++
e)
5248 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5249 for (
unsigned int e = 0;
e < dim; ++
e)
5250 value[e] -= t_jac[e][d] * tmp2;
5253 for (
unsigned int d = 0;
d < dim; ++
d)
5254 values[d * nqp] =
value[d];
5260 const std::size_t nqp_d = this->n_quadrature_points * dim;
5261 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5265 const VectorizedArrayType JxW =
5266 this->J_value[0] * this->quadrature_weights[q_point];
5272 std::array<VectorizedArrayType, dim> jac;
5273 for (
unsigned int d = 0;
d < dim; ++
d)
5274 jac[d] = this->jacobian[0][d][d];
5276 for (
unsigned int d = 0;
d < dim; ++
d)
5278 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5279 if constexpr (n_components == 1)
5280 gradients[d] = grad_in[d] * factor;
5282 for (
unsigned int comp = 0; comp < n_components; ++comp)
5283 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5290 this->jacobian[q_point] :
5292 const VectorizedArrayType JxW =
5294 this->J_value[q_point] :
5295 this->J_value[0] * this->quadrature_weights[q_point];
5296 if constexpr (n_components == 1)
5297 for (
unsigned int d = 0;
d < dim; ++
d)
5299 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5300 for (
unsigned int e = 1;
e < dim; ++
e)
5301 new_val += (jac[e][d] * grad_in[e]);
5305 for (
unsigned int comp = 0; comp < n_components; ++comp)
5306 for (
unsigned int d = 0;
d < dim; ++
d)
5308 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5309 for (
unsigned int e = 1;
e < dim; ++
e)
5310 new_val += (jac[e][d] * grad_in[comp][e]);
5322 typename VectorizedArrayType>
5323template <
int,
typename>
5327 const unsigned int q_point)
5329 static_assert(n_components == 1 && dim == 1,
5330 "Do not try to modify the default template parameters used for"
5331 " selectively enabling this function via std::enable_if!");
5332 submit_gradient(grad_in[0], q_point);
5341 typename VectorizedArrayType>
5347 Assert(this->normal_x_jacobian !=
nullptr,
5349 "update_gradients"));
5351 this->gradients_quad_submitted =
true;
5354 const std::size_t nqp_d = this->n_quadrature_points * dim;
5355 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5359 const VectorizedArrayType JxW_jac = this->J_value[0] *
5360 this->quadrature_weights[q_point] *
5361 this->normal_x_jacobian[0][dim - 1];
5362 for (
unsigned int comp = 0; comp < n_components; ++comp)
5364 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5365 gradients[comp * nqp_d + d] = VectorizedArrayType();
5366 if constexpr (n_components == 1)
5367 gradients[dim - 1] = grad_in * JxW_jac;
5369 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5374 const unsigned int index =
5377 this->normal_x_jacobian[
index];
5378 const VectorizedArrayType JxW =
5380 this->J_value[
index] * this->quadrature_weights[q_point] :
5381 this->J_value[
index];
5382 for (
unsigned int comp = 0; comp < n_components; ++comp)
5383 for (
unsigned int d = 0;
d < dim; ++
d)
5384 if constexpr (n_components == 1)
5387 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5397 typename VectorizedArrayType>
5400 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5406 Assert(this->J_value !=
nullptr,
5408 "update_hessians"));
5409 Assert(this->jacobian !=
nullptr,
5411 "update_hessians"));
5413 this->hessians_quad_submitted =
true;
5417 const std::size_t nqp = this->n_quadrature_points;
5418 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5421 const VectorizedArrayType JxW =
5422 this->J_value[0] * this->quadrature_weights[q_point];
5425 for (
unsigned int d = 0;
d < dim; ++
d)
5427 const auto jac_d = this->jacobian[0][
d][
d];
5428 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5429 for (
unsigned int comp = 0; comp < n_components; ++comp)
5430 if constexpr (n_components == 1)
5431 this->hessians_quad[
d * nqp + q_point] =
5432 hessian_in[
d][
d] * factor;
5434 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5435 hessian_in[comp][d][d] * factor;
5439 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5440 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5442 const auto jac_d = this->jacobian[0][
d][
d];
5443 const auto jac_e = this->jacobian[0][
e][
e];
5444 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5445 for (
unsigned int comp = 0; comp < n_components; ++comp)
5446 if constexpr (n_components == 1)
5447 this->hessians_quad[off_dia * nqp + q_point] =
5448 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5450 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5451 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5458 const VectorizedArrayType JxW =
5459 this->J_value[0] * this->quadrature_weights[q_point];
5460 for (
unsigned int comp = 0; comp < n_components; ++comp)
5463 if constexpr (n_components == 1)
5464 hessian_c = hessian_in;
5466 hessian_c = hessian_in[comp];
5469 VectorizedArrayType tmp[dim][dim];
5470 for (
unsigned int i = 0; i < dim; ++i)
5471 for (
unsigned int j = 0; j < dim; ++j)
5473 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5474 for (
unsigned int k = 1; k < dim; ++k)
5475 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5479 VectorizedArrayType tmp2[dim][dim];
5480 for (
unsigned int i = 0; i < dim; ++i)
5481 for (
unsigned int j = 0; j < dim; ++j)
5483 tmp2[i][j] = jac[0][i] * tmp[0][j];
5484 for (
unsigned int k = 1; k < dim; ++k)
5485 tmp2[i][j] += jac[k][i] * tmp[k][j];
5489 for (
unsigned int d = 0;
d < dim; ++
d)
5490 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5494 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5495 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5496 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5497 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5503 const VectorizedArrayType JxW = this->J_value[q_point];
5504 const auto &jac_grad = this->jacobian_gradients[q_point];
5505 for (
unsigned int comp = 0; comp < n_components; ++comp)
5508 if constexpr (n_components == 1)
5509 hessian_c = hessian_in;
5511 hessian_c = hessian_in[comp];
5514 VectorizedArrayType tmp[dim][dim];
5515 for (
unsigned int i = 0; i < dim; ++i)
5516 for (
unsigned int j = 0; j < dim; ++j)
5518 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5519 for (
unsigned int k = 1; k < dim; ++k)
5520 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5524 VectorizedArrayType tmp2[dim][dim];
5525 for (
unsigned int i = 0; i < dim; ++i)
5526 for (
unsigned int j = 0; j < dim; ++j)
5528 tmp2[i][j] = jac[0][i] * tmp[0][j];
5529 for (
unsigned int k = 1; k < dim; ++k)
5530 tmp2[i][j] += jac[k][i] * tmp[k][j];
5534 for (
unsigned int d = 0;
d < dim; ++
d)
5535 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5539 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5540 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5541 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5542 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5545 for (
unsigned int d = 0;
d < dim; ++
d)
5547 VectorizedArrayType
sum = 0;
5548 for (
unsigned int e = 0;
e < dim; ++
e)
5549 sum += hessian_c[e][e] * jac_grad[e][d];
5550 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5551 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5553 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5554 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5567 typename VectorizedArrayType>
5571 const unsigned int q_point)
5577 Assert(this->J_value !=
nullptr,
5579 "update_hessians"));
5580 Assert(this->jacobian !=
nullptr,
5582 "update_hessians"));
5584 this->hessians_quad_submitted =
true;
5588 const std::size_t nqp = this->n_quadrature_points;
5589 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5592 const VectorizedArrayType JxW =
5593 this->J_value[0] * this->quadrature_weights[q_point];
5595 const auto nxj = this->normal_x_jacobian[0];
5598 for (
unsigned int d = 0;
d < dim; ++
d)
5600 const auto nxj_d = nxj[
d];
5601 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5602 for (
unsigned int comp = 0; comp < n_components; ++comp)
5603 if constexpr (n_components == 1)
5604 this->hessians_quad[
d * nqp + q_point] =
5605 normal_hessian_in * factor;
5607 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5608 normal_hessian_in[comp] * factor;
5612 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5613 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5615 const auto jac_d = nxj[
d];
5616 const auto jac_e = nxj[
e];
5617 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5618 for (
unsigned int comp = 0; comp < n_components; ++comp)
5619 if constexpr (n_components == 1)
5620 this->hessians_quad[off_dia * nqp + q_point] =
5621 2. * normal_hessian_in * factor;
5623 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5624 2. * normal_hessian_in[comp] * factor;
5629 const auto normal = this->normal_vector(q_point);
5630 const auto normal_projector =
outer_product(normal, normal);
5631 if constexpr (n_components == 1)
5632 submit_hessian(normal_hessian_in * normal_projector, q_point);
5636 for (
unsigned int comp = 0; comp < n_components; ++comp)
5637 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5638 submit_hessian(tmp, q_point);
5649 typename VectorizedArrayType>
5654 VectorizedArrayType>::value_type
5660 Assert(this->values_quad_submitted ==
true,
5665 const std::size_t nqp = this->n_quadrature_points;
5666 for (
unsigned int q = 0; q < nqp; ++q)
5667 for (
unsigned int comp = 0; comp < n_components; ++comp)
5668 return_value[comp] += this->values_quad[comp * nqp + q];
5669 if constexpr (n_components == 1)
5670 return return_value[0];
5672 return return_value;
5681 typename VectorizedArrayType>
5682template <
int,
typename>
5687 static_assert(n_components == dim,
5688 "Do not try to modify the default template parameters used for"
5689 " selectively enabling this function via std::enable_if!");
5692 Assert(this->gradients_quad_initialized ==
true,
5696 Assert(this->jacobian !=
nullptr,
5698 "update_gradients"));
5700 VectorizedArrayType divergence;
5701 const std::size_t nqp = this->n_quadrature_points;
5704 this->
data->element_type ==
5707 VectorizedArrayType inv_det =
5710 this->jacobian[0][0][0] *
5711 ((dim == 2) ? this->jacobian[0][1][1] :
5712 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5720 if (is_face && dim == 2 && this->get_face_no() < 2)
5724 divergence = this->gradients_quad[q_point * dim];
5725 for (
unsigned int d = 1;
d < dim; ++
d)
5726 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5727 divergence *= inv_det;
5736 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5737 for (
unsigned int d = 1;
d < dim; ++
d)
5738 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5739 this->jacobian[0][
d][
d];
5746 this->jacobian[q_point] :
5748 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5749 for (
unsigned int e = 1;
e < dim; ++
e)
5750 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5751 for (
unsigned int d = 1;
d < dim; ++
d)
5752 for (
unsigned int e = 0;
e < dim; ++
e)
5754 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5766 typename VectorizedArrayType>
5767template <
int,
typename>
5772 static_assert(n_components == dim,
5773 "Do not try to modify the default template parameters used for"
5774 " selectively enabling this function via std::enable_if!");
5777 const auto grad = get_gradient(q_point);
5778 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5779 VectorizedArrayType half = Number(0.5);
5780 for (
unsigned int d = 0;
d < dim; ++
d)
5781 symmetrized[d] = grad[d][d];
5787 symmetrized[2] = grad[0][1] + grad[1][0];
5788 symmetrized[2] *= half;
5791 symmetrized[3] = grad[0][1] + grad[1][0];
5792 symmetrized[3] *= half;
5793 symmetrized[4] = grad[0][2] + grad[2][0];
5794 symmetrized[4] *= half;
5795 symmetrized[5] = grad[1][2] + grad[2][1];
5796 symmetrized[5] *= half;
5810 typename VectorizedArrayType>
5811template <
int,
typename>
5813 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5815 get_curl(const unsigned
int q_point) const
5817 static_assert(dim > 1 && n_components == dim,
5818 "Do not try to modify the default template parameters used for"
5819 " selectively enabling this function via std::enable_if!");
5823 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5827 curl[0] = grad[1][0] - grad[0][1];
5830 curl[0] = grad[2][1] - grad[1][2];
5831 curl[1] = grad[0][2] - grad[2][0];
5832 curl[2] = grad[1][0] - grad[0][1];
5846 typename VectorizedArrayType>
5847template <
int,
typename>
5851 const unsigned int q_point)
5853 static_assert(n_components == dim,
5854 "Do not try to modify the default template parameters used for"
5855 " selectively enabling this function via std::enable_if!");
5861 Assert(this->J_value !=
nullptr,
5863 "update_gradients"));
5864 Assert(this->jacobian !=
nullptr,
5866 "update_gradients"));
5868 this->gradients_quad_submitted =
true;
5871 const std::size_t nqp_d = this->n_quadrature_points * dim;
5872 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5874 if (this->
data->element_type ==
5881 const VectorizedArrayType fac =
5883 this->quadrature_weights[q_point] * div_in :
5885 this->J_value[q_point] :
5886 this->J_value[0] * this->quadrature_weights[q_point]) *
5889 this->jacobian[this->cell_type >
5893 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5895 for (
unsigned int d = 0;
d < dim; ++
d)
5897 for (
unsigned int e = 0;
e < dim; ++
e)
5898 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5900 this->divergence_is_requested =
true;
5907 const VectorizedArrayType fac =
5908 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5909 for (
unsigned int d = 0;
d < dim; ++
d)
5911 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5912 for (
unsigned int e = 0;
e < dim; ++
e)
5913 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5920 this->jacobian[q_point] :
5922 const VectorizedArrayType fac =
5924 this->J_value[q_point] :
5925 this->J_value[0] * this->quadrature_weights[q_point]) *
5927 for (
unsigned int d = 0;
d < dim; ++
d)
5929 for (
unsigned int e = 0;
e < dim; ++
e)
5930 gradients[d * nqp_d + e] = jac[d][e] * fac;
5942 typename VectorizedArrayType>
5943template <
int,
typename>
5948 const unsigned int q_point)
5950 static_assert(n_components == dim,
5951 "Do not try to modify the default template parameters used for"
5952 " selectively enabling this function via std::enable_if!");
5955 this->
data->element_type !=
5966 Assert(this->J_value !=
nullptr,
5968 "update_gradients"));
5969 Assert(this->jacobian !=
nullptr,
5971 "update_gradients"));
5973 this->gradients_quad_submitted =
true;
5976 const std::size_t nqp_d = this->n_quadrature_points * dim;
5977 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5980 const VectorizedArrayType JxW =
5981 this->J_value[0] * this->quadrature_weights[q_point];
5983 for (
unsigned int d = 0;
d < dim; ++
d)
5984 gradients[d * nqp_d + d] =
5986 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5987 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5989 const VectorizedArrayType
value =
5998 const VectorizedArrayType JxW =
6000 this->J_value[q_point] :
6001 this->J_value[0] * this->quadrature_weights[q_point];
6004 this->jacobian[q_point] :
6006 VectorizedArrayType weighted[dim][dim];
6007 for (
unsigned int i = 0; i < dim; ++i)
6009 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6010 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6012 const VectorizedArrayType
value =
6014 weighted[i][j] =
value;
6015 weighted[j][i] =
value;
6017 for (
unsigned int comp = 0; comp < dim; ++comp)
6018 for (
unsigned int d = 0;
d < dim; ++
d)
6020 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6021 for (
unsigned int e = 1;
e < dim; ++
e)
6022 new_val += jac[e][d] * weighted[comp][e];
6034 typename VectorizedArrayType>
6035template <
int,
typename>
6039 const unsigned int q_point)
6041 static_assert(n_components == dim,
6042 "Do not try to modify the default template parameters used for"
6043 " selectively enabling this function via std::enable_if!");
6049 grad[1][0] = curl[0];
6050 grad[0][1] = -curl[0];
6053 grad[2][1] = curl[0];
6054 grad[1][2] = -curl[0];
6055 grad[0][2] = curl[1];
6056 grad[2][0] = -curl[1];
6057 grad[1][0] = curl[2];
6058 grad[0][1] = -curl[2];
6063 submit_gradient(grad, q_point);
6076 typename VectorizedArrayType>
6082 VectorizedArrayType>::
6084 const unsigned int fe_no,
6085 const unsigned int quad_no,
6086 const unsigned int first_selected_component,
6087 const unsigned int active_fe_index,
6088 const unsigned int active_quad_index)
6089 : BaseClass(matrix_free,
6091 first_selected_component,
6098 numbers::invalid_unsigned_int )
6099 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6100 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6101 , n_q_points(this->
data->n_q_points)
6103 check_template_arguments(fe_no, 0);
6113 typename VectorizedArrayType>
6119 VectorizedArrayType>::
6121 const std::pair<unsigned int, unsigned int> &range,
6122 const unsigned int dof_no,
6123 const unsigned int quad_no,
6124 const unsigned int first_selected_component)
6128 first_selected_component,
6129 matrix_free.get_cell_active_fe_index(range))
6139 typename VectorizedArrayType>
6145 VectorizedArrayType>::
6150 const unsigned int first_selected_component)
6151 : BaseClass(mapping,
6155 first_selected_component,
6157 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6158 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6159 , n_q_points(this->
data->n_q_points)
6171 typename VectorizedArrayType>
6177 VectorizedArrayType>::
6181 const unsigned int first_selected_component)
6186 first_selected_component,
6188 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6189 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6190 , n_q_points(this->
data->n_q_points)
6202 typename VectorizedArrayType>
6208 VectorizedArrayType>::
6211 const unsigned int first_selected_component)
6212 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6214 other.mapped_geometry->get_quadrature(),
6215 other.mapped_geometry->get_fe_values().get_update_flags(),
6216 first_selected_component,
6218 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6219 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6220 , n_q_points(this->
data->n_q_points)
6232 typename VectorizedArrayType>
6241 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6242 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6243 , n_q_points(this->
data->n_q_points)
6255 typename VectorizedArrayType>
6261 VectorizedArrayType> &
6267 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6269 BaseClass::operator=(other);
6281 typename VectorizedArrayType>
6288 VectorizedArrayType>::
6289 check_template_arguments(
const unsigned int dof_no,
6290 const unsigned int first_selected_component)
6293 (void)first_selected_component;
6296 this->
data->dofs_per_component_on_cell > 0,
6298 "There is nothing useful you can do with an FEEvaluation object with "
6299 "FE_Nothing, i.e., without DoFs! If you have passed to "
6300 "MatrixFree::reinit() a collection of finite elements also containing "
6301 "FE_Nothing, please check - before creating FEEvaluation - the category "
6302 "of the current range by calling either "
6303 "MatrixFree::get_cell_range_category(range) or "
6304 "MatrixFree::get_face_range_category(range). The returned category "
6305 "is the index of the active FE, which you can use to exclude "
6312 static_cast<unsigned int>(fe_degree) !=
6313 this->
data->data.front().fe_degree) ||
6314 n_q_points != this->n_quadrature_points)
6316 std::string message =
6317 "-------------------------------------------------------\n";
6318 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6319 message +=
" Called --> FEEvaluation<dim,";
6323 message +=
",Number>(data";
6339 if (
static_cast<unsigned int>(fe_degree) ==
6340 this->
data->data.front().fe_degree)
6342 proposed_dof_comp = dof_no;
6343 proposed_fe_comp = first_selected_component;
6346 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6348 for (
unsigned int nf = 0;
6351 if (this->matrix_free
6352 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6354 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6356 proposed_dof_comp = no;
6357 proposed_fe_comp = nf;
6361 this->mapping_data->descriptor[this->active_quad_index]
6363 proposed_quad_comp = this->quad_no;
6365 for (
unsigned int no = 0;
6370 .descriptor[this->active_quad_index]
6371 .n_q_points == n_q_points)
6373 proposed_quad_comp = no;
6380 if (proposed_dof_comp != first_selected_component)
6381 message +=
"Wrong vector component selection:\n";
6383 message +=
"Wrong quadrature formula selection:\n";
6384 message +=
" Did you mean FEEvaluation<dim,";
6388 message +=
",Number>(data";
6397 std::string correct_pos;
6398 if (proposed_dof_comp != dof_no)
6399 correct_pos =
" ^ ";
6402 if (proposed_quad_comp != this->quad_no)
6403 correct_pos +=
" ^ ";
6406 if (proposed_fe_comp != first_selected_component)
6407 correct_pos +=
" ^\n";
6409 correct_pos +=
" \n";
6415 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6416 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6417 message +=
"Wrong template arguments:\n";
6418 message +=
" Did you mean FEEvaluation<dim,";
6423 message +=
",Number>(data";
6431 std::string correct_pos;
6432 if (this->
data->data.front().fe_degree !=
6433 static_cast<unsigned int>(fe_degree))
6437 if (proposed_n_q_points_1d != n_q_points_1d)
6438 correct_pos +=
" ^\n";
6440 correct_pos +=
" \n";
6441 message +=
" " + correct_pos;
6443 Assert(
static_cast<unsigned int>(fe_degree) ==
6444 this->
data->data.front().fe_degree &&
6445 n_q_points == this->n_quadrature_points,
6451 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6462 typename VectorizedArrayType>
6471 Assert(this->matrix_free !=
nullptr,
6472 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6473 " Integer indexing is not possible."));
6481 const unsigned int offsets =
6482 this->mapping_data->data_index_offsets[
cell_index];
6483 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6484 this->J_value = &this->mapping_data->JxW_values[offsets];
6485 if (!this->mapping_data->jacobian_gradients[0].empty())
6487 this->jacobian_gradients =
6488 this->mapping_data->jacobian_gradients[0].data() + offsets;
6489 this->jacobian_gradients_non_inverse =
6490 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6496 for (
unsigned int i = 0; i < n_lanes; ++i)
6497 this->cell_ids[i] =
cell_index * n_lanes + i;
6504 this->cell_ids[i] =
cell_index * n_lanes + i;
6505 for (; i < n_lanes; ++i)
6509 if (this->mapping_data->quadrature_points.empty() ==
false)
6511 &this->mapping_data->quadrature_points
6512 [this->mapping_data->quadrature_point_offsets[this->cell]];
6515 this->is_reinitialized =
true;
6516 this->dof_values_initialized =
false;
6517 this->values_quad_initialized =
false;
6518 this->gradients_quad_initialized =
false;
6519 this->hessians_quad_initialized =
false;
6530 typename VectorizedArrayType>
6537 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6544 this->cell_ids = cell_ids;
6549 for (
unsigned int v = 0; v < n_lanes; ++v)
6563 if (this->mapped_geometry ==
nullptr)
6564 this->mapped_geometry =
6565 std::make_shared<internal::MatrixFreeFunctions::
6566 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6568 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6570 auto &this_jacobian_data = mapping_storage.jacobians[0];
6571 auto &this_J_value_data = mapping_storage.JxW_values;
6572 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6573 auto &this_jacobian_gradients_non_inverse_data =
6574 mapping_storage.jacobian_gradients_non_inverse[0];
6575 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6579 if (this_jacobian_data.size() != 2)
6580 this_jacobian_data.resize_fast(2);
6582 if (this_J_value_data.size() != 1)
6583 this_J_value_data.resize_fast(1);
6585 const auto &update_flags_cells =
6589 this_jacobian_gradients_data.size() != 1)
6591 this_jacobian_gradients_data.resize_fast(1);
6592 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6596 this_quadrature_points_data.size() != 1)
6597 this_quadrature_points_data.resize_fast(1);
6601 if (this_jacobian_data.size() != this->n_quadrature_points)
6602 this_jacobian_data.resize_fast(this->n_quadrature_points);
6604 if (this_J_value_data.size() != this->n_quadrature_points)
6605 this_J_value_data.resize_fast(this->n_quadrature_points);
6607 const auto &update_flags_cells =
6611 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6613 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6614 this_jacobian_gradients_non_inverse_data.resize_fast(
6615 this->n_quadrature_points);
6619 this_quadrature_points_data.size() != this->n_quadrature_points)
6620 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6624 this->jacobian = this_jacobian_data.data();
6625 this->J_value = this_J_value_data.data();
6626 this->jacobian_gradients = this_jacobian_gradients_data.data();
6627 this->jacobian_gradients_non_inverse =
6628 this_jacobian_gradients_non_inverse_data.data();
6632 for (
unsigned int v = 0; v < n_lanes; ++v)
6639 const unsigned int cell_batch_index =
cell_index / n_lanes;
6640 const unsigned int offsets =
6641 this->mapping_data->data_index_offsets[cell_batch_index];
6642 const unsigned int lane =
cell_index % n_lanes;
6644 if (this->cell_type <=
6648 for (
unsigned int q = 0; q < 2; ++q)
6649 for (
unsigned int i = 0; i < dim; ++i)
6650 for (
unsigned int j = 0; j < dim; ++j)
6651 this_jacobian_data[q][i][j][v] =
6652 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6654 const unsigned int q = 0;
6656 this_J_value_data[q][v] =
6657 this->mapping_data->JxW_values[offsets + q][lane];
6659 const auto &update_flags_cells =
6664 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6665 for (
unsigned int j = 0; j < dim; ++j)
6666 this_jacobian_gradients_data[q][i][j][v] =
6668 ->jacobian_gradients[0][offsets + q][i][j][lane];
6670 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6671 for (
unsigned int j = 0; j < dim; ++j)
6672 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6674 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6679 for (
unsigned int i = 0; i < dim; ++i)
6680 this_quadrature_points_data[q][i][v] =
6681 this->mapping_data->quadrature_points
6683 ->quadrature_point_offsets[cell_batch_index] +
6689 const auto cell_type =
6693 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6695 const unsigned int q_src =
6701 this_J_value_data[q][v] =
6702 this->mapping_data->JxW_values[offsets + q_src][lane];
6704 for (
unsigned int i = 0; i < dim; ++i)
6705 for (
unsigned int j = 0; j < dim; ++j)
6706 this_jacobian_data[q][i][j][v] =
6708 ->jacobians[0][offsets + q_src][i][j][lane];
6710 const auto &update_flags_cells =
6715 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6716 for (
unsigned int j = 0; j < dim; ++j)
6717 this_jacobian_gradients_data[q][i][j][v] =
6719 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6721 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6722 for (
unsigned int j = 0; j < dim; ++j)
6723 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6725 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6738 this->mapping_data->quadrature_points
6740 ->quadrature_point_offsets[cell_batch_index] +
6744 this->mapping_data->jacobians[0][offsets + 1];
6746 for (
unsigned int d = 0;
d < dim; ++
d)
6749 static_cast<Number
>(
6750 this->descriptor->quadrature.point(q)[
d]);
6752 for (
unsigned int d = 0;
d < dim; ++
d)
6753 for (
unsigned int e = 0;
e < dim; ++
e)
6756 static_cast<Number
>(
6757 this->descriptor->quadrature.point(q)[
e]);
6759 for (
unsigned int i = 0; i < dim; ++i)
6760 this_quadrature_points_data[q][i][v] = point[i][lane];
6765 for (
unsigned int i = 0; i < dim; ++i)
6766 this_quadrature_points_data[q][i][v] =
6767 this->mapping_data->quadrature_points
6769 ->quadrature_point_offsets[cell_batch_index] +
6778 this->is_reinitialized =
true;
6779 this->dof_values_initialized =
false;
6780 this->values_quad_initialized =
false;
6781 this->gradients_quad_initialized =
false;
6782 this->hessians_quad_initialized =
false;
6793 typename VectorizedArrayType>
6794template <
bool level_dof_access>
6801 VectorizedArrayType>
::
6804 Assert(this->matrix_free ==
nullptr,
6805 ExcMessage(
"Cannot use initialization from cell iterator if "
6806 "initialized from MatrixFree object. Use variant for "
6807 "on the fly computation with arguments as for FEValues "
6810 this->mapped_geometry->reinit(
6812 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6813 if (level_dof_access)
6814 cell->get_mg_dof_indices(this->local_dof_indices);
6816 cell->get_dof_indices(this->local_dof_indices);
6819 this->is_reinitialized =
true;
6830 typename VectorizedArrayType>
6837 VectorizedArrayType>
::
6840 Assert(this->matrix_free == 0,
6841 ExcMessage(
"Cannot use initialization from cell iterator if "
6842 "initialized from MatrixFree object. Use variant for "
6843 "on the fly computation with arguments as for FEValues "
6846 this->mapped_geometry->reinit(cell);
6849 this->is_reinitialized =
true;
6860 typename VectorizedArrayType>
6867 VectorizedArrayType>::
6871 Assert(this->dof_values_initialized ==
true,
6874 evaluate(this->values_dofs, evaluation_flags);
6884 typename VectorizedArrayType>
6891 VectorizedArrayType>::
6892 evaluate(
const VectorizedArrayType *values_array,
6895 const bool hessians_on_general_cells =
6899 if (hessians_on_general_cells)
6902 if (this->
data->element_type ==
6908 if constexpr (fe_degree > -1)
6911 template run<fe_degree, n_q_points_1d>(n_components,
6912 evaluation_flag_actual,
6920 evaluation_flag_actual,
6921 const_cast<VectorizedArrayType *
>(values_array),
6926 this->values_quad_initialized =
6928 this->gradients_quad_initialized =
6930 this->hessians_quad_initialized =
6941 template <
typename Number,
6942 typename VectorizedArrayType,
6944 typename EvaluatorType,
6945 std::enable_if_t<internal::has_begin<VectorType> &&
6948 VectorizedArrayType *
6949 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6955 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6956 const auto &dof_info = fe_eval.get_dof_info();
6963 if (std::is_same_v<typename VectorType::value_type, Number> &&
6967 interleaved_contiguous &&
6968 reinterpret_cast<
std::size_t>(
6970 dof_info.dof_indices_contiguous
6971 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6972 [cell * VectorizedArrayType::
size()]) %
6973 sizeof(VectorizedArrayType) ==
6976 return reinterpret_cast<VectorizedArrayType *
>(
6980 [cell * VectorizedArrayType::size()] +
6982 [fe_eval.get_active_fe_index()]
6983 [fe_eval.get_first_selected_component()] *
6984 VectorizedArrayType::size());
6993 template <
typename Number,
6994 typename VectorizedArrayType,
6996 typename EvaluatorType,
6997 std::enable_if_t<!internal::has_begin<VectorType> ||
7000 VectorizedArrayType *
7001 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7014 typename VectorizedArrayType>
7015template <
typename VectorType>
7022 VectorizedArrayType>::
7023 gather_evaluate(
const VectorType &input_vector,
7026 const VectorizedArrayType *src_ptr =
7027 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7028 *
this, input_vector);
7029 if (src_ptr !=
nullptr)
7030 evaluate(src_ptr, evaluation_flag);
7033 this->read_dof_values(input_vector);
7034 evaluate(this->begin_dof_values(), evaluation_flag);
7045 typename VectorizedArrayType>
7052 VectorizedArrayType>::
7055 integrate(integration_flag, this->values_dofs);
7058 this->dof_values_initialized =
true;
7069 typename VectorizedArrayType>
7076 VectorizedArrayType>::
7078 VectorizedArrayType *values_array,
7079 const bool sum_into_values_array)
7083 Assert(this->values_quad_submitted ==
true,
7086 Assert(this->gradients_quad_submitted ==
true,
7089 Assert(this->hessians_quad_submitted ==
true,
7092 Assert(this->matrix_free !=
nullptr ||
7093 this->mapped_geometry->is_initialized(),
7099 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7100 "EvaluationFlags::hessians are supported."));
7106 unsigned int size = n_components * dim * n_q_points;
7109 for (
unsigned int i = 0; i <
size; ++i)
7110 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7114 for (
unsigned int i = 0; i <
size; ++i)
7115 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7120 if (n_components == dim &&
7121 this->
data->element_type ==
7125 this->divergence_is_requested ==
false)
7127 unsigned int size = n_components * n_q_points;
7130 for (
unsigned int i = 0; i <
size; ++i)
7131 this->values_quad[i] += this->values_from_gradients_quad[i];
7135 for (
unsigned int i = 0; i <
size; ++i)
7136 this->values_quad[i] = this->values_from_gradients_quad[i];
7141 if constexpr (fe_degree > -1)
7144 template run<fe_degree, n_q_points_1d>(n_components,
7145 integration_flag_actual,
7148 sum_into_values_array);
7154 integration_flag_actual,
7157 sum_into_values_array);
7161 this->dof_values_initialized =
true;
7172 typename VectorizedArrayType>
7173template <
typename VectorType>
7180 VectorizedArrayType>::
7182 VectorType &destination)
7184 VectorizedArrayType *dst_ptr =
7185 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7186 *
this, destination);
7187 if (dst_ptr !=
nullptr)
7188 integrate(integration_flag, dst_ptr,
true);
7191 integrate(integration_flag, this->begin_dof_values());
7192 this->distribute_local_to_global(destination);
7203 typename VectorizedArrayType>
7210 VectorizedArrayType>::dof_indices()
const
7212 return {0U, dofs_per_cell};
7226 typename VectorizedArrayType>
7232 VectorizedArrayType>::
7235 const bool is_interior_face,
7236 const unsigned int dof_no,
7237 const unsigned int quad_no,
7238 const unsigned int first_selected_component,
7239 const unsigned int active_fe_index,
7240 const unsigned int active_quad_index,
7241 const unsigned int face_type)
7242 : BaseClass(matrix_free,
7244 first_selected_component,
7252 , dofs_per_component(this->
data->dofs_per_component_on_cell)
7253 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
7254 , n_q_points(this->n_quadrature_points)
7264 typename VectorizedArrayType>
7270 VectorizedArrayType>::
7273 const std::pair<unsigned int, unsigned int> &rang