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_info !=
nullptr &&
3916 this->dof_access_index ==
3918 this->active_fe_index,
3922 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3924 apply_hanging_node_constraints(
false);
3927 this->dof_values_initialized =
true;
3937 typename VectorizedArrayType>
3938template <
typename VectorType>
3942 const unsigned int first_index,
3943 const std::bitset<n_lanes> &
mask)
3945 const auto src_data = internal::get_vector_data<n_components_>(
3948 this->dof_access_index ==
3950 this->active_fe_index,
3954 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3957 this->dof_values_initialized =
true;
3967 typename VectorizedArrayType>
3968template <
typename VectorType>
3972 const unsigned int first_index,
3973 const std::bitset<n_lanes> &
mask)
const
3976 Assert(this->dof_values_initialized ==
true,
3980 apply_hanging_node_constraints(
true);
3982 const auto dst_data = internal::get_vector_data<n_components_>(
3985 this->dof_access_index ==
3987 this->active_fe_index,
3992 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
4001 typename VectorizedArrayType>
4002template <
typename VectorType>
4006 const unsigned int first_index,
4007 const std::bitset<n_lanes> &
mask)
const
4010 Assert(this->dof_values_initialized ==
true,
4014 const auto dst_data = internal::get_vector_data<n_components_>(
4017 this->dof_access_index ==
4019 this->active_fe_index,
4023 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4032 typename VectorizedArrayType>
4033template <
typename VectorType>
4037 const unsigned int first_index,
4038 const std::bitset<n_lanes> &
mask)
const
4041 Assert(this->dof_values_initialized ==
true,
4045 const auto dst_data = internal::get_vector_data<n_components_>(
4048 this->dof_access_index ==
4050 this->active_fe_index,
4054 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4067 typename VectorizedArrayType>
4073 VectorizedArrayType>::value_type
4078 if constexpr (n_components == 1)
4079 return this->values_dofs[dof];
4082 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4083 Tensor<1, n_components_, VectorizedArrayType> return_value;
4084 for (
unsigned int comp = 0; comp < n_components; ++comp)
4085 return_value[comp] = this->values_dofs[comp * dofs + dof];
4086 return return_value;
4096 typename VectorizedArrayType>
4102 VectorizedArrayType>::value_type
4104 get_value(
const unsigned int q_point)
const
4107 Assert(this->values_quad_initialized ==
true,
4112 if constexpr (n_components == 1)
4113 return this->values_quad[q_point];
4116 if (n_components == dim &&
4117 this->
data->element_type ==
4122 Assert(this->values_quad_initialized ==
true,
4127 Assert(this->J_value !=
nullptr,
4130 const std::size_t nqp = this->n_quadrature_points;
4138 const VectorizedArrayType inv_det =
4139 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4140 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4141 this->jacobian[0][2][2];
4144 for (
unsigned int comp = 0; comp < n_components; ++comp)
4145 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4146 jac[comp][comp] * inv_det;
4153 this->jacobian[q_point] :
4162 const VectorizedArrayType inv_det =
4163 (is_face && dim == 2 && this->get_face_no() < 2) ?
4167 for (
unsigned int comp = 0; comp < n_components; ++comp)
4169 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4170 for (
unsigned int e = 1;
e < dim; ++
e)
4172 this->values_quad[e * nqp + q_point] * jac[comp][e];
4173 value_out[comp] *= inv_det;
4180 const std::size_t nqp = this->n_quadrature_points;
4182 for (
unsigned int comp = 0; comp < n_components; ++comp)
4183 return_value[comp] = this->values_quad[comp * nqp + q_point];
4184 return return_value;
4195 typename VectorizedArrayType>
4201 VectorizedArrayType>::gradient_type
4206 Assert(this->gradients_quad_initialized ==
true,
4211 Assert(this->jacobian !=
nullptr,
4213 "update_gradients"));
4214 const std::size_t nqp = this->n_quadrature_points;
4216 if constexpr (n_components == dim && dim > 1)
4218 if (this->
data->element_type ==
4223 Assert(this->gradients_quad_initialized ==
true,
4228 Assert(this->jacobian !=
nullptr,
4230 "update_gradients"));
4231 const std::size_t nqp = this->n_quadrature_points;
4232 const std::size_t nqp_d = nqp * dim;
4235 this->gradients_quad + q_point * dim;
4246 const VectorizedArrayType inv_det =
4247 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4248 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4249 this->jacobian[0][2][2];
4252 for (
unsigned int d = 0;
d < dim; ++
d)
4253 for (
unsigned int comp = 0; comp < n_components; ++comp)
4254 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4256 (jac[comp][comp] * inv_det);
4268 const VectorizedArrayType inv_det =
4269 (is_face && dim == 2 && this->get_face_no() < 2) ?
4273 VectorizedArrayType tmp[dim][dim];
4275 for (
unsigned int d = 0;
d < dim; ++
d)
4276 for (
unsigned int e = 0;
e < dim; ++
e)
4279 for (
unsigned int f = 1; f < dim; ++f)
4280 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4282 for (
unsigned int comp = 0; comp < n_components; ++comp)
4283 for (
unsigned int d = 0;
d < dim; ++
d)
4285 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4286 for (
unsigned int f = 1; f < dim; ++f)
4287 res += jac[comp][f] * tmp[d][f];
4289 grad_out[comp][
d] = res * inv_det;
4300 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4302 "update_hessians"));
4304 const auto jac_grad =
4305 this->jacobian_gradients_non_inverse[q_point];
4307 this->jacobian[q_point];
4311 const VectorizedArrayType inv_det =
4312 (is_face && dim == 2 && this->get_face_no() < 2) ?
4319 VectorizedArrayType tmp[dim][dim];
4320 for (
unsigned int d = 0;
d < dim; ++
d)
4321 for (
unsigned int e = 0;
e < dim; ++
e)
4324 for (
unsigned int f = 1; f < dim; ++f)
4325 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4330 for (
unsigned int d = 0;
d < dim; ++
d)
4332 for (
unsigned int e = 0;
e < dim; ++
e)
4334 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4335 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4336 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4339 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4341 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4346 for (
unsigned int d = 0;
d < dim; ++
d)
4347 for (
unsigned int e = 0;
e < dim; ++
e)
4349 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4350 for (
unsigned int f = 1; f < dim; ++f)
4351 res += tmp[f][d] * inv_t_jac[e][f];
4352 grad_out[
d][
e] = res;
4358 VectorizedArrayType tmp3[dim], tmp4[dim];
4359 for (
unsigned int d = 0;
d < dim; ++
d)
4361 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4362 for (
unsigned int e = 1;
e < dim; ++
e)
4363 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4365 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4366 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4367 for (
unsigned int d = 0;
d < dim; ++
d)
4369 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4370 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4372 for (
unsigned int d = 0;
d < dim; ++
d)
4374 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4375 for (
unsigned int e = 1;
e < dim; ++
e)
4376 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4379 VectorizedArrayType tmp2[dim];
4380 for (
unsigned int d = 0;
d < dim; ++
d)
4382 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4383 for (
unsigned e = 1;
e < dim; ++
e)
4385 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4388 for (
unsigned int d = 0;
d < dim; ++
d)
4389 for (
unsigned int e = 0;
e < dim; ++
e)
4391 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4395 grad_out[
d][
e] *= inv_det;
4406 for (
unsigned int comp = 0; comp < n_components; ++comp)
4407 for (
unsigned int d = 0;
d < dim; ++
d)
4409 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4410 this->jacobian[0][
d][
d];
4419 for (
unsigned int comp = 0; comp < n_components; ++comp)
4420 for (
unsigned int d = 0;
d < dim; ++
d)
4423 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4424 for (
unsigned int e = 1;
e < dim; ++
e)
4425 grad_out[comp][d] +=
4427 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4430 if constexpr (n_components == 1)
4442 typename VectorizedArrayType>
4448 VectorizedArrayType>::value_type
4454 Assert(this->gradients_quad_initialized ==
true,
4458 Assert(this->normal_x_jacobian !=
nullptr,
4460 "update_gradients"));
4462 const std::size_t nqp = this->n_quadrature_points;
4466 for (
unsigned int comp = 0; comp < n_components; ++comp)
4468 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4469 (this->normal_x_jacobian[0][dim - 1]);
4472 const std::size_t
index =
4474 for (
unsigned int comp = 0; comp < n_components; ++comp)
4476 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4477 this->normal_x_jacobian[
index][0];
4478 for (
unsigned int d = 1;
d < dim; ++
d)
4480 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4481 this->normal_x_jacobian[
index][
d];
4484 if constexpr (n_components == 1)
4496 template <
typename VectorizedArrayType>
4499 const VectorizedArrayType *
const hessians,
4501 VectorizedArrayType (&tmp)[1][1])
4503 tmp[0][0] = jac[0][0] *
hessians[0];
4506 template <
typename VectorizedArrayType>
4509 const VectorizedArrayType *
const hessians,
4510 const unsigned int nqp,
4511 VectorizedArrayType (&tmp)[2][2])
4513 for (
unsigned int d = 0;
d < 2; ++
d)
4521 template <
typename VectorizedArrayType>
4524 const VectorizedArrayType *
const hessians,
4525 const unsigned int nqp,
4526 VectorizedArrayType (&tmp)[3][3])
4528 for (
unsigned int d = 0;
d < 3; ++
d)
4549 typename VectorizedArrayType>
4554 VectorizedArrayType>::hessian_type
4559 Assert(this->hessians_quad_initialized ==
true,
4564 Assert(this->jacobian !=
nullptr,
4574 const std::size_t nqp = this->n_quadrature_points;
4575 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4580 for (
unsigned int comp = 0; comp < n_components; ++comp)
4582 for (
unsigned int d = 0;
d < dim; ++
d)
4583 hessian_out[comp][d][d] =
4584 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4585 (jac[
d][
d] * jac[
d][
d]);
4591 hessian_out[comp][0][1] =
4592 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4593 (jac[0][0] * jac[1][1]);
4596 hessian_out[comp][0][1] =
4597 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4598 (jac[0][0] * jac[1][1]);
4599 hessian_out[comp][0][2] =
4600 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4601 (jac[0][0] * jac[2][2]);
4602 hessian_out[comp][1][2] =
4603 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4604 (jac[1][1] * jac[2][2]);
4609 for (
unsigned int d = 0;
d < dim; ++
d)
4610 for (
unsigned int e = d + 1;
e < dim; ++
e)
4611 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4617 for (
unsigned int comp = 0; comp < n_components; ++comp)
4619 VectorizedArrayType tmp[dim][dim];
4620 internal::hessian_unit_times_jac(
4621 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4624 for (
unsigned int d = 0;
d < dim; ++
d)
4625 for (
unsigned int e = d;
e < dim; ++
e)
4627 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4628 for (
unsigned int f = 1; f < dim; ++f)
4629 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4636 for (
unsigned int d = 0;
d < dim; ++
d)
4637 for (
unsigned int e = d + 1;
e < dim; ++
e)
4638 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4644 const auto &jac_grad = this->jacobian_gradients[q_point];
4645 for (
unsigned int comp = 0; comp < n_components; ++comp)
4647 VectorizedArrayType tmp[dim][dim];
4648 internal::hessian_unit_times_jac(
4649 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4652 for (
unsigned int d = 0;
d < dim; ++
d)
4653 for (
unsigned int e = d;
e < dim; ++
e)
4655 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4656 for (
unsigned int f = 1; f < dim; ++f)
4657 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4661 for (
unsigned int d = 0;
d < dim; ++
d)
4662 for (
unsigned int e = 0;
e < dim; ++
e)
4663 hessian_out[comp][d][d] +=
4665 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4668 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4669 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4670 for (
unsigned int f = 0; f < dim; ++f)
4671 hessian_out[comp][d][e] +=
4672 jac_grad[count][f] *
4673 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4676 for (
unsigned int d = 0;
d < dim; ++
d)
4677 for (
unsigned int e = d + 1;
e < dim; ++
e)
4678 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4681 if constexpr (n_components == 1)
4682 return hessian_out[0];
4693 typename VectorizedArrayType>
4698 VectorizedArrayType>::gradient_type
4704 Assert(this->hessians_quad_initialized ==
true,
4715 const std::size_t nqp = this->n_quadrature_points;
4716 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4722 for (
unsigned int comp = 0; comp < n_components; ++comp)
4723 for (
unsigned int d = 0;
d < dim; ++
d)
4724 hessian_out[comp][d] =
4725 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4726 (jac[
d][
d] * jac[
d][
d]);
4731 for (
unsigned int comp = 0; comp < n_components; ++comp)
4735 VectorizedArrayType tmp[dim][dim];
4736 internal::hessian_unit_times_jac(
4737 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4741 for (
unsigned int d = 0;
d < dim; ++
d)
4743 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4744 for (
unsigned int f = 1; f < dim; ++f)
4745 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4752 const auto &jac_grad = this->jacobian_gradients[q_point];
4753 for (
unsigned int comp = 0; comp < n_components; ++comp)
4757 VectorizedArrayType tmp[dim][dim];
4758 internal::hessian_unit_times_jac(
4759 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4763 for (
unsigned int d = 0;
d < dim; ++
d)
4765 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4766 for (
unsigned int f = 1; f < dim; ++f)
4767 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4770 for (
unsigned int d = 0;
d < dim; ++
d)
4771 for (
unsigned int e = 0;
e < dim; ++
e)
4772 hessian_out[comp][d] +=
4774 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4778 if constexpr (n_components == 1)
4779 return hessian_out[0];
4790 typename VectorizedArrayType>
4795 VectorizedArrayType>::value_type
4801 Assert(this->hessians_quad_initialized ==
true,
4806 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4807 if constexpr (n_components == 1)
4809 VectorizedArrayType
sum = hess_diag[0];
4810 for (
unsigned int d = 1;
d < dim; ++
d)
4811 sum += hess_diag[d];
4817 for (
unsigned int comp = 0; comp < n_components; ++comp)
4819 laplacian_out[comp] = hess_diag[comp][0];
4820 for (
unsigned int d = 1;
d < dim; ++
d)
4821 laplacian_out[comp] += hess_diag[comp][d];
4823 return laplacian_out;
4833 typename VectorizedArrayType>
4838 VectorizedArrayType>::value_type
4843 Assert(this->hessians_quad_initialized ==
true,
4848 Assert(this->normal_x_jacobian !=
nullptr,
4850 "update_hessians"));
4854 const std::size_t nqp = this->n_quadrature_points;
4855 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4859 const auto nxj = this->normal_x_jacobian[0];
4861 for (
unsigned int comp = 0; comp < n_components; ++comp)
4863 for (
unsigned int d = 0;
d < dim; ++
d)
4864 hessian_out[comp] +=
4865 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4866 (nxj[
d]) * (nxj[d]);
4873 hessian_out[comp] +=
4874 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4878 hessian_out[comp] +=
4879 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4881 hessian_out[comp] +=
4882 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4884 hessian_out[comp] +=
4885 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4896 const auto normal = this->normal_vector(q_point);
4897 const auto hessian = get_hessian(q_point);
4899 if constexpr (n_components == 1)
4900 hessian_out[0] =
hessian * normal * normal;
4902 for (
unsigned int comp = 0; comp < n_components; ++comp)
4903 hessian_out[comp] =
hessian[comp] * normal * normal;
4905 if constexpr (n_components == 1)
4906 return hessian_out[0];
4917 typename VectorizedArrayType>
4923 this->dof_values_initialized =
true;
4925 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4927 for (
unsigned int comp = 0; comp < n_components; ++comp)
4928 if constexpr (n_components == 1)
4929 this->values_dofs[comp * dofs + dof] = val_in;
4931 this->values_dofs[comp * dofs + dof] = val_in[comp];
4940 typename VectorizedArrayType>
4943 submit_value(
const value_type val_in,
const unsigned int q_point)
4949 Assert(this->J_value !=
nullptr,
4953 this->values_quad_submitted =
true;
4956 const std::size_t nqp = this->n_quadrature_points;
4957 VectorizedArrayType *
values = this->values_quad + q_point;
4959 const VectorizedArrayType JxW =
4961 this->J_value[0] * this->quadrature_weights[q_point] :
4962 this->J_value[q_point];
4963 if constexpr (n_components == 1)
4964 values[0] = val_in * JxW;
4967 if (n_components == dim &&
4968 this->
data->element_type ==
4973 Assert(this->J_value !=
nullptr,
4978 this->values_quad_submitted =
true;
4981 VectorizedArrayType *
values = this->values_quad + q_point;
4982 const std::size_t nqp = this->n_quadrature_points;
4988 const VectorizedArrayType weight =
4989 this->quadrature_weights[q_point];
4991 for (
unsigned int comp = 0; comp < n_components; ++comp)
4992 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4999 this->jacobian[q_point] :
5004 const VectorizedArrayType fac =
5006 this->quadrature_weights[q_point] :
5008 this->J_value[q_point] :
5009 this->J_value[0] * this->quadrature_weights[q_point]) *
5010 ((dim == 2 && this->get_face_no() < 2) ?
5019 for (
unsigned int comp = 0; comp < n_components; ++comp)
5021 values[comp * nqp] = val_in[0] * jac[0][comp];
5022 for (
unsigned int e = 1;
e < dim; ++
e)
5023 values[comp * nqp] += val_in[e] * jac[e][comp];
5024 values[comp * nqp] *= fac;
5029 for (
unsigned int comp = 0; comp < n_components; ++comp)
5030 values[comp * nqp] = val_in[comp] * JxW;
5040 typename VectorizedArrayType>
5041template <
int,
typename>
5045 const unsigned int q_point)
5047 static_assert(n_components == 1,
5048 "Do not try to modify the default template parameters used for"
5049 " selectively enabling this function via std::enable_if!");
5050 submit_value(val_in[0], q_point);
5059 typename VectorizedArrayType>
5062 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
5068 Assert(this->J_value !=
nullptr,
5070 "update_gradients"));
5071 Assert(this->jacobian !=
nullptr,
5073 "update_gradients"));
5075 this->gradients_quad_submitted =
true;
5078 if constexpr (dim > 1 && n_components == dim)
5080 if (this->
data->element_type ==
5089 Assert(this->J_value !=
nullptr,
5091 "update_gradients"));
5092 Assert(this->jacobian !=
nullptr,
5094 "update_gradients"));
5096 this->gradients_quad_submitted =
true;
5099 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5100 VectorizedArrayType *
values =
5101 this->values_from_gradients_quad + q_point;
5102 const std::size_t nqp = this->n_quadrature_points;
5103 const std::size_t nqp_d = nqp * dim;
5113 const VectorizedArrayType weight =
5114 this->quadrature_weights[q_point];
5115 for (
unsigned int d = 0;
d < dim; ++
d)
5116 for (
unsigned int comp = 0; comp < n_components; ++comp)
5117 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5119 (jac[comp][comp] * weight);
5131 const VectorizedArrayType fac =
5133 this->quadrature_weights[q_point] :
5134 this->J_value[0] * this->quadrature_weights[q_point] *
5135 ((dim == 2 && this->get_face_no() < 2) ?
5140 VectorizedArrayType tmp[dim][dim];
5141 for (
unsigned int d = 0;
d < dim; ++
d)
5142 for (
unsigned int e = 0;
e < dim; ++
e)
5144 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5145 for (
unsigned int f = 1; f < dim; ++f)
5146 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5148 for (
unsigned int comp = 0; comp < n_components; ++comp)
5149 for (
unsigned int d = 0;
d < dim; ++
d)
5151 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5152 for (
unsigned int f = 1; f < dim; ++f)
5153 res += jac[f][comp] * tmp[d][f];
5162 const auto jac_grad =
5163 this->jacobian_gradients_non_inverse[q_point];
5165 this->jacobian[q_point];
5169 const VectorizedArrayType fac =
5170 (!is_face) ? this->quadrature_weights[q_point] :
5171 this->J_value[q_point] *
5172 ((dim == 2 && this->get_face_no() < 2) ?
5181 VectorizedArrayType tmp3[dim], tmp4[dim];
5182 for (
unsigned int d = 0;
d < dim; ++
d)
5184 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5185 for (
unsigned int e = 1;
e < dim; ++
e)
5186 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5188 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5189 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5190 for (
unsigned int d = 0;
d < dim; ++
d)
5192 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5193 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5195 for (
unsigned int d = 0;
d < dim; ++
d)
5197 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5198 for (
unsigned int e = 1;
e < dim; ++
e)
5199 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5205 VectorizedArrayType tmp[dim][dim];
5208 for (
unsigned int d = 0;
d < dim; ++
d)
5209 for (
unsigned int e = 0;
e < dim; ++
e)
5211 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5212 for (
unsigned int f = 1; f < dim; ++f)
5213 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5216 for (
unsigned int d = 0;
d < dim; ++
d)
5217 for (
unsigned int e = 0;
e < dim; ++
e)
5219 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5220 for (
unsigned int f = 1; f < dim; ++f)
5221 res += t_jac[d][f] * tmp[e][f];
5228 VectorizedArrayType
value[dim];
5229 for (
unsigned int d = 0;
d < dim; ++
d)
5231 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5232 for (
unsigned int e = 1;
e < dim; ++
e)
5233 value[d] += tmp[d][e] * jac_grad[d][e];
5235 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5236 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5237 for (
unsigned int d = 0;
d < dim; ++
d)
5239 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5240 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5245 for (
unsigned int d = 0;
d < dim; ++
d)
5247 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5248 for (
unsigned int e = 1;
e < dim; ++
e)
5249 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5250 for (
unsigned int e = 0;
e < dim; ++
e)
5251 value[e] -= t_jac[e][d] * tmp2;
5254 for (
unsigned int d = 0;
d < dim; ++
d)
5255 values[d * nqp] =
value[d];
5261 const std::size_t nqp_d = this->n_quadrature_points * dim;
5262 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5266 const VectorizedArrayType JxW =
5267 this->J_value[0] * this->quadrature_weights[q_point];
5273 std::array<VectorizedArrayType, dim> jac;
5274 for (
unsigned int d = 0;
d < dim; ++
d)
5275 jac[d] = this->jacobian[0][d][d];
5277 for (
unsigned int d = 0;
d < dim; ++
d)
5279 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5280 if constexpr (n_components == 1)
5281 gradients[d] = grad_in[d] * factor;
5283 for (
unsigned int comp = 0; comp < n_components; ++comp)
5284 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5291 this->jacobian[q_point] :
5293 const VectorizedArrayType JxW =
5295 this->J_value[q_point] :
5296 this->J_value[0] * this->quadrature_weights[q_point];
5297 if constexpr (n_components == 1)
5298 for (
unsigned int d = 0;
d < dim; ++
d)
5300 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5301 for (
unsigned int e = 1;
e < dim; ++
e)
5302 new_val += (jac[e][d] * grad_in[e]);
5306 for (
unsigned int comp = 0; comp < n_components; ++comp)
5307 for (
unsigned int d = 0;
d < dim; ++
d)
5309 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5310 for (
unsigned int e = 1;
e < dim; ++
e)
5311 new_val += (jac[e][d] * grad_in[comp][e]);
5323 typename VectorizedArrayType>
5324template <
int,
typename>
5328 const unsigned int q_point)
5330 static_assert(n_components == 1 && dim == 1,
5331 "Do not try to modify the default template parameters used for"
5332 " selectively enabling this function via std::enable_if!");
5333 submit_gradient(grad_in[0], q_point);
5342 typename VectorizedArrayType>
5348 Assert(this->normal_x_jacobian !=
nullptr,
5350 "update_gradients"));
5352 this->gradients_quad_submitted =
true;
5355 const std::size_t nqp_d = this->n_quadrature_points * dim;
5356 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5360 const VectorizedArrayType JxW_jac = this->J_value[0] *
5361 this->quadrature_weights[q_point] *
5362 this->normal_x_jacobian[0][dim - 1];
5363 for (
unsigned int comp = 0; comp < n_components; ++comp)
5365 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5366 gradients[comp * nqp_d + d] = VectorizedArrayType();
5367 if constexpr (n_components == 1)
5368 gradients[dim - 1] = grad_in * JxW_jac;
5370 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5375 const unsigned int index =
5378 this->normal_x_jacobian[
index];
5379 const VectorizedArrayType JxW =
5381 this->J_value[
index] * this->quadrature_weights[q_point] :
5382 this->J_value[
index];
5383 for (
unsigned int comp = 0; comp < n_components; ++comp)
5384 for (
unsigned int d = 0;
d < dim; ++
d)
5385 if constexpr (n_components == 1)
5388 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5398 typename VectorizedArrayType>
5401 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5407 Assert(this->J_value !=
nullptr,
5409 "update_hessians"));
5410 Assert(this->jacobian !=
nullptr,
5412 "update_hessians"));
5414 this->hessians_quad_submitted =
true;
5418 const std::size_t nqp = this->n_quadrature_points;
5419 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5422 const VectorizedArrayType JxW =
5423 this->J_value[0] * this->quadrature_weights[q_point];
5426 for (
unsigned int d = 0;
d < dim; ++
d)
5428 const auto jac_d = this->jacobian[0][
d][
d];
5429 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5430 for (
unsigned int comp = 0; comp < n_components; ++comp)
5431 if constexpr (n_components == 1)
5432 this->hessians_quad[
d * nqp + q_point] =
5433 hessian_in[
d][
d] * factor;
5435 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5436 hessian_in[comp][d][d] * factor;
5440 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5441 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5443 const auto jac_d = this->jacobian[0][
d][
d];
5444 const auto jac_e = this->jacobian[0][
e][
e];
5445 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5446 for (
unsigned int comp = 0; comp < n_components; ++comp)
5447 if constexpr (n_components == 1)
5448 this->hessians_quad[off_dia * nqp + q_point] =
5449 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5451 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5452 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5459 const VectorizedArrayType JxW =
5460 this->J_value[0] * this->quadrature_weights[q_point];
5461 for (
unsigned int comp = 0; comp < n_components; ++comp)
5464 if constexpr (n_components == 1)
5465 hessian_c = hessian_in;
5467 hessian_c = hessian_in[comp];
5470 VectorizedArrayType tmp[dim][dim];
5471 for (
unsigned int i = 0; i < dim; ++i)
5472 for (
unsigned int j = 0; j < dim; ++j)
5474 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5475 for (
unsigned int k = 1; k < dim; ++k)
5476 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5480 VectorizedArrayType tmp2[dim][dim];
5481 for (
unsigned int i = 0; i < dim; ++i)
5482 for (
unsigned int j = 0; j < dim; ++j)
5484 tmp2[i][j] = jac[0][i] * tmp[0][j];
5485 for (
unsigned int k = 1; k < dim; ++k)
5486 tmp2[i][j] += jac[k][i] * tmp[k][j];
5490 for (
unsigned int d = 0;
d < dim; ++
d)
5491 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5495 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5496 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5497 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5498 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5504 const VectorizedArrayType JxW = this->J_value[q_point];
5505 const auto &jac_grad = this->jacobian_gradients[q_point];
5506 for (
unsigned int comp = 0; comp < n_components; ++comp)
5509 if constexpr (n_components == 1)
5510 hessian_c = hessian_in;
5512 hessian_c = hessian_in[comp];
5515 VectorizedArrayType tmp[dim][dim];
5516 for (
unsigned int i = 0; i < dim; ++i)
5517 for (
unsigned int j = 0; j < dim; ++j)
5519 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5520 for (
unsigned int k = 1; k < dim; ++k)
5521 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5525 VectorizedArrayType tmp2[dim][dim];
5526 for (
unsigned int i = 0; i < dim; ++i)
5527 for (
unsigned int j = 0; j < dim; ++j)
5529 tmp2[i][j] = jac[0][i] * tmp[0][j];
5530 for (
unsigned int k = 1; k < dim; ++k)
5531 tmp2[i][j] += jac[k][i] * tmp[k][j];
5535 for (
unsigned int d = 0;
d < dim; ++
d)
5536 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5540 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5541 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5542 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5543 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5546 for (
unsigned int d = 0;
d < dim; ++
d)
5548 VectorizedArrayType
sum = 0;
5549 for (
unsigned int e = 0;
e < dim; ++
e)
5550 sum += hessian_c[e][e] * jac_grad[e][d];
5551 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5552 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5554 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5555 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5568 typename VectorizedArrayType>
5572 const unsigned int q_point)
5578 Assert(this->J_value !=
nullptr,
5580 "update_hessians"));
5581 Assert(this->jacobian !=
nullptr,
5583 "update_hessians"));
5585 this->hessians_quad_submitted =
true;
5589 const std::size_t nqp = this->n_quadrature_points;
5590 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5593 const VectorizedArrayType JxW =
5594 this->J_value[0] * this->quadrature_weights[q_point];
5596 const auto nxj = this->normal_x_jacobian[0];
5599 for (
unsigned int d = 0;
d < dim; ++
d)
5601 const auto nxj_d = nxj[
d];
5602 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5603 for (
unsigned int comp = 0; comp < n_components; ++comp)
5604 if constexpr (n_components == 1)
5605 this->hessians_quad[
d * nqp + q_point] =
5606 normal_hessian_in * factor;
5608 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5609 normal_hessian_in[comp] * factor;
5613 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5614 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5616 const auto jac_d = nxj[
d];
5617 const auto jac_e = nxj[
e];
5618 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5619 for (
unsigned int comp = 0; comp < n_components; ++comp)
5620 if constexpr (n_components == 1)
5621 this->hessians_quad[off_dia * nqp + q_point] =
5622 2. * normal_hessian_in * factor;
5624 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5625 2. * normal_hessian_in[comp] * factor;
5630 const auto normal = this->normal_vector(q_point);
5631 const auto normal_projector =
outer_product(normal, normal);
5632 if constexpr (n_components == 1)
5633 submit_hessian(normal_hessian_in * normal_projector, q_point);
5637 for (
unsigned int comp = 0; comp < n_components; ++comp)
5638 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5639 submit_hessian(tmp, q_point);
5650 typename VectorizedArrayType>
5655 VectorizedArrayType>::value_type
5661 Assert(this->values_quad_submitted ==
true,
5666 const std::size_t nqp = this->n_quadrature_points;
5667 for (
unsigned int q = 0; q < nqp; ++q)
5668 for (
unsigned int comp = 0; comp < n_components; ++comp)
5669 return_value[comp] += this->values_quad[comp * nqp + q];
5670 if constexpr (n_components == 1)
5671 return return_value[0];
5673 return return_value;
5682 typename VectorizedArrayType>
5683template <
int,
typename>
5688 static_assert(n_components == dim,
5689 "Do not try to modify the default template parameters used for"
5690 " selectively enabling this function via std::enable_if!");
5693 Assert(this->gradients_quad_initialized ==
true,
5697 Assert(this->jacobian !=
nullptr,
5699 "update_gradients"));
5701 VectorizedArrayType divergence;
5702 const std::size_t nqp = this->n_quadrature_points;
5705 this->
data->element_type ==
5708 VectorizedArrayType inv_det =
5711 this->jacobian[0][0][0] *
5712 ((dim == 2) ? this->jacobian[0][1][1] :
5713 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5721 if (is_face && dim == 2 && this->get_face_no() < 2)
5725 divergence = this->gradients_quad[q_point * dim];
5726 for (
unsigned int d = 1;
d < dim; ++
d)
5727 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5728 divergence *= inv_det;
5737 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5738 for (
unsigned int d = 1;
d < dim; ++
d)
5739 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5740 this->jacobian[0][
d][
d];
5747 this->jacobian[q_point] :
5749 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5750 for (
unsigned int e = 1;
e < dim; ++
e)
5751 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5752 for (
unsigned int d = 1;
d < dim; ++
d)
5753 for (
unsigned int e = 0;
e < dim; ++
e)
5755 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5767 typename VectorizedArrayType>
5768template <
int,
typename>
5773 static_assert(n_components == dim,
5774 "Do not try to modify the default template parameters used for"
5775 " selectively enabling this function via std::enable_if!");
5778 const auto grad = get_gradient(q_point);
5779 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5780 VectorizedArrayType half = Number(0.5);
5781 for (
unsigned int d = 0;
d < dim; ++
d)
5782 symmetrized[d] = grad[d][d];
5788 symmetrized[2] = grad[0][1] + grad[1][0];
5789 symmetrized[2] *= half;
5792 symmetrized[3] = grad[0][1] + grad[1][0];
5793 symmetrized[3] *= half;
5794 symmetrized[4] = grad[0][2] + grad[2][0];
5795 symmetrized[4] *= half;
5796 symmetrized[5] = grad[1][2] + grad[2][1];
5797 symmetrized[5] *= half;
5811 typename VectorizedArrayType>
5812template <
int,
typename>
5814 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5816 get_curl(const unsigned
int q_point) const
5818 static_assert(dim > 1 && n_components == dim,
5819 "Do not try to modify the default template parameters used for"
5820 " selectively enabling this function via std::enable_if!");
5824 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5828 curl[0] = grad[1][0] - grad[0][1];
5831 curl[0] = grad[2][1] - grad[1][2];
5832 curl[1] = grad[0][2] - grad[2][0];
5833 curl[2] = grad[1][0] - grad[0][1];
5847 typename VectorizedArrayType>
5848template <
int,
typename>
5852 const unsigned int q_point)
5854 static_assert(n_components == dim,
5855 "Do not try to modify the default template parameters used for"
5856 " selectively enabling this function via std::enable_if!");
5862 Assert(this->J_value !=
nullptr,
5864 "update_gradients"));
5865 Assert(this->jacobian !=
nullptr,
5867 "update_gradients"));
5869 this->gradients_quad_submitted =
true;
5872 const std::size_t nqp_d = this->n_quadrature_points * dim;
5873 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5875 if (this->
data->element_type ==
5882 const VectorizedArrayType fac =
5884 this->quadrature_weights[q_point] * div_in :
5886 this->J_value[q_point] :
5887 this->J_value[0] * this->quadrature_weights[q_point]) *
5890 this->jacobian[this->cell_type >
5894 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5896 for (
unsigned int d = 0;
d < dim; ++
d)
5898 for (
unsigned int e = 0;
e < dim; ++
e)
5899 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5901 this->divergence_is_requested =
true;
5908 const VectorizedArrayType fac =
5909 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5910 for (
unsigned int d = 0;
d < dim; ++
d)
5912 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5913 for (
unsigned int e = 0;
e < dim; ++
e)
5914 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5921 this->jacobian[q_point] :
5923 const VectorizedArrayType fac =
5925 this->J_value[q_point] :
5926 this->J_value[0] * this->quadrature_weights[q_point]) *
5928 for (
unsigned int d = 0;
d < dim; ++
d)
5930 for (
unsigned int e = 0;
e < dim; ++
e)
5931 gradients[d * nqp_d + e] = jac[d][e] * fac;
5943 typename VectorizedArrayType>
5944template <
int,
typename>
5949 const unsigned int q_point)
5951 static_assert(n_components == dim,
5952 "Do not try to modify the default template parameters used for"
5953 " selectively enabling this function via std::enable_if!");
5956 this->
data->element_type !=
5967 Assert(this->J_value !=
nullptr,
5969 "update_gradients"));
5970 Assert(this->jacobian !=
nullptr,
5972 "update_gradients"));
5974 this->gradients_quad_submitted =
true;
5977 const std::size_t nqp_d = this->n_quadrature_points * dim;
5978 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5981 const VectorizedArrayType JxW =
5982 this->J_value[0] * this->quadrature_weights[q_point];
5984 for (
unsigned int d = 0;
d < dim; ++
d)
5985 gradients[d * nqp_d + d] =
5987 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5988 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5990 const VectorizedArrayType
value =
5999 const VectorizedArrayType JxW =
6001 this->J_value[q_point] :
6002 this->J_value[0] * this->quadrature_weights[q_point];
6005 this->jacobian[q_point] :
6007 VectorizedArrayType weighted[dim][dim];
6008 for (
unsigned int i = 0; i < dim; ++i)
6010 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6011 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6013 const VectorizedArrayType
value =
6015 weighted[i][j] =
value;
6016 weighted[j][i] =
value;
6018 for (
unsigned int comp = 0; comp < dim; ++comp)
6019 for (
unsigned int d = 0;
d < dim; ++
d)
6021 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6022 for (
unsigned int e = 1;
e < dim; ++
e)
6023 new_val += jac[e][d] * weighted[comp][e];
6035 typename VectorizedArrayType>
6036template <
int,
typename>
6040 const unsigned int q_point)
6042 static_assert(n_components == dim,
6043 "Do not try to modify the default template parameters used for"
6044 " selectively enabling this function via std::enable_if!");
6050 grad[1][0] = curl[0];
6051 grad[0][1] = -curl[0];
6054 grad[2][1] = curl[0];
6055 grad[1][2] = -curl[0];
6056 grad[0][2] = curl[1];
6057 grad[2][0] = -curl[1];
6058 grad[1][0] = curl[2];
6059 grad[0][1] = -curl[2];
6064 submit_gradient(grad, q_point);
6077 typename VectorizedArrayType>
6083 VectorizedArrayType>::
6085 const unsigned int fe_no,
6086 const unsigned int quad_no,
6087 const unsigned int first_selected_component,
6088 const unsigned int active_fe_index,
6089 const unsigned int active_quad_index)
6090 : BaseClass(matrix_free,
6092 first_selected_component,
6100 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6101 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6102 , n_q_points(this->
data->n_q_points)
6104 check_template_arguments(fe_no, 0);
6114 typename VectorizedArrayType>
6120 VectorizedArrayType>::
6122 const std::pair<unsigned int, unsigned int> &range,
6123 const unsigned int dof_no,
6124 const unsigned int quad_no,
6125 const unsigned int first_selected_component)
6129 first_selected_component,
6130 matrix_free.get_cell_active_fe_index(range))
6140 typename VectorizedArrayType>
6146 VectorizedArrayType>::
6151 const unsigned int first_selected_component)
6152 : BaseClass(mapping,
6156 first_selected_component,
6158 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6159 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6160 , n_q_points(this->
data->n_q_points)
6172 typename VectorizedArrayType>
6178 VectorizedArrayType>::
6182 const unsigned int first_selected_component)
6187 first_selected_component,
6189 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6190 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6191 , n_q_points(this->
data->n_q_points)
6203 typename VectorizedArrayType>
6209 VectorizedArrayType>::
6212 const unsigned int first_selected_component)
6213 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6215 other.mapped_geometry->get_quadrature(),
6216 other.mapped_geometry->get_fe_values().get_update_flags(),
6217 first_selected_component,
6219 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6220 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6221 , n_q_points(this->
data->n_q_points)
6233 typename VectorizedArrayType>
6242 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6243 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6244 , n_q_points(this->
data->n_q_points)
6256 typename VectorizedArrayType>
6262 VectorizedArrayType> &
6268 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6270 BaseClass::operator=(other);
6282 typename VectorizedArrayType>
6289 VectorizedArrayType>::
6290 check_template_arguments(
const unsigned int dof_no,
6291 const unsigned int first_selected_component)
6294 (void)first_selected_component;
6297 this->
data->dofs_per_component_on_cell > 0,
6299 "There is nothing useful you can do with an FEEvaluation object with "
6300 "FE_Nothing, i.e., without DoFs! If you have passed to "
6301 "MatrixFree::reinit() a collection of finite elements also containing "
6302 "FE_Nothing, please check - before creating FEEvaluation - the category "
6303 "of the current range by calling either "
6304 "MatrixFree::get_cell_range_category(range) or "
6305 "MatrixFree::get_face_range_category(range). The returned category "
6306 "is the index of the active FE, which you can use to exclude "
6313 static_cast<unsigned int>(fe_degree) !=
6314 this->
data->data.front().fe_degree) ||
6315 n_q_points != this->n_quadrature_points)
6317 std::string message =
6318 "-------------------------------------------------------\n";
6319 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6320 message +=
" Called --> FEEvaluation<dim,";
6324 message +=
",Number>(data";
6340 if (
static_cast<unsigned int>(fe_degree) ==
6341 this->
data->data.front().fe_degree)
6343 proposed_dof_comp = dof_no;
6344 proposed_fe_comp = first_selected_component;
6347 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6349 for (
unsigned int nf = 0;
6352 if (this->matrix_free
6353 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6355 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6357 proposed_dof_comp = no;
6358 proposed_fe_comp = nf;
6362 this->mapping_data->descriptor[this->active_quad_index]
6364 proposed_quad_comp = this->quad_no;
6366 for (
unsigned int no = 0;
6371 .descriptor[this->active_quad_index]
6372 .n_q_points == n_q_points)
6374 proposed_quad_comp = no;
6381 if (proposed_dof_comp != first_selected_component)
6382 message +=
"Wrong vector component selection:\n";
6384 message +=
"Wrong quadrature formula selection:\n";
6385 message +=
" Did you mean FEEvaluation<dim,";
6389 message +=
",Number>(data";
6398 std::string correct_pos;
6399 if (proposed_dof_comp != dof_no)
6400 correct_pos =
" ^ ";
6403 if (proposed_quad_comp != this->quad_no)
6404 correct_pos +=
" ^ ";
6407 if (proposed_fe_comp != first_selected_component)
6408 correct_pos +=
" ^\n";
6410 correct_pos +=
" \n";
6416 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6417 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6418 message +=
"Wrong template arguments:\n";
6419 message +=
" Did you mean FEEvaluation<dim,";
6424 message +=
",Number>(data";
6432 std::string correct_pos;
6433 if (this->
data->data.front().fe_degree !=
6434 static_cast<unsigned int>(fe_degree))
6438 if (proposed_n_q_points_1d != n_q_points_1d)
6439 correct_pos +=
" ^\n";
6441 correct_pos +=
" \n";
6442 message +=
" " + correct_pos;
6444 Assert(
static_cast<unsigned int>(fe_degree) ==
6445 this->
data->data.front().fe_degree &&
6446 n_q_points == this->n_quadrature_points,
6452 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6463 typename VectorizedArrayType>
6472 Assert(this->matrix_free !=
nullptr,
6473 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6474 " Integer indexing is not possible."));
6482 const unsigned int offsets =
6483 this->mapping_data->data_index_offsets[
cell_index];
6484 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6485 this->J_value = &this->mapping_data->JxW_values[offsets];
6486 if (!this->mapping_data->jacobian_gradients[0].empty())
6488 this->jacobian_gradients =
6489 this->mapping_data->jacobian_gradients[0].data() + offsets;
6490 this->jacobian_gradients_non_inverse =
6491 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6497 for (
unsigned int i = 0; i < n_lanes; ++i)
6498 this->cell_ids[i] =
cell_index * n_lanes + i;
6505 this->cell_ids[i] =
cell_index * n_lanes + i;
6506 for (; i < n_lanes; ++i)
6510 if (this->mapping_data->quadrature_points.empty() ==
false)
6512 &this->mapping_data->quadrature_points
6513 [this->mapping_data->quadrature_point_offsets[this->cell]];
6516 this->is_reinitialized =
true;
6517 this->dof_values_initialized =
false;
6518 this->values_quad_initialized =
false;
6519 this->gradients_quad_initialized =
false;
6520 this->hessians_quad_initialized =
false;
6531 typename VectorizedArrayType>
6538 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6545 this->cell_ids = cell_ids;
6550 for (
unsigned int v = 0; v < n_lanes; ++v)
6564 if (this->mapped_geometry ==
nullptr)
6565 this->mapped_geometry =
6566 std::make_shared<internal::MatrixFreeFunctions::
6567 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6569 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6571 auto &this_jacobian_data = mapping_storage.jacobians[0];
6572 auto &this_J_value_data = mapping_storage.JxW_values;
6573 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6574 auto &this_jacobian_gradients_non_inverse_data =
6575 mapping_storage.jacobian_gradients_non_inverse[0];
6576 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6580 if (this_jacobian_data.size() != 2)
6581 this_jacobian_data.resize_fast(2);
6583 if (this_J_value_data.size() != 1)
6584 this_J_value_data.resize_fast(1);
6586 const auto &update_flags_cells =
6590 this_jacobian_gradients_data.size() != 1)
6592 this_jacobian_gradients_data.resize_fast(1);
6593 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6597 this_quadrature_points_data.size() != 1)
6598 this_quadrature_points_data.resize_fast(1);
6602 if (this_jacobian_data.size() != this->n_quadrature_points)
6603 this_jacobian_data.resize_fast(this->n_quadrature_points);
6605 if (this_J_value_data.size() != this->n_quadrature_points)
6606 this_J_value_data.resize_fast(this->n_quadrature_points);
6608 const auto &update_flags_cells =
6612 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6614 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6615 this_jacobian_gradients_non_inverse_data.resize_fast(
6616 this->n_quadrature_points);
6620 this_quadrature_points_data.size() != this->n_quadrature_points)
6621 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6625 this->jacobian = this_jacobian_data.data();
6626 this->J_value = this_J_value_data.data();
6627 this->jacobian_gradients = this_jacobian_gradients_data.data();
6628 this->jacobian_gradients_non_inverse =
6629 this_jacobian_gradients_non_inverse_data.data();
6633 for (
unsigned int v = 0; v < n_lanes; ++v)
6640 const unsigned int cell_batch_index =
cell_index / n_lanes;
6641 const unsigned int offsets =
6642 this->mapping_data->data_index_offsets[cell_batch_index];
6643 const unsigned int lane =
cell_index % n_lanes;
6645 if (this->cell_type <=
6649 for (
unsigned int q = 0; q < 2; ++q)
6650 for (
unsigned int i = 0; i < dim; ++i)
6651 for (
unsigned int j = 0; j < dim; ++j)
6652 this_jacobian_data[q][i][j][v] =
6653 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6655 const unsigned int q = 0;
6657 this_J_value_data[q][v] =
6658 this->mapping_data->JxW_values[offsets + q][lane];
6660 const auto &update_flags_cells =
6665 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6666 for (
unsigned int j = 0; j < dim; ++j)
6667 this_jacobian_gradients_data[q][i][j][v] =
6669 ->jacobian_gradients[0][offsets + q][i][j][lane];
6671 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6672 for (
unsigned int j = 0; j < dim; ++j)
6673 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6675 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6680 for (
unsigned int i = 0; i < dim; ++i)
6681 this_quadrature_points_data[q][i][v] =
6682 this->mapping_data->quadrature_points
6684 ->quadrature_point_offsets[cell_batch_index] +
6690 const auto cell_type =
6694 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6696 const unsigned int q_src =
6702 this_J_value_data[q][v] =
6703 this->mapping_data->JxW_values[offsets + q_src][lane];
6705 for (
unsigned int i = 0; i < dim; ++i)
6706 for (
unsigned int j = 0; j < dim; ++j)
6707 this_jacobian_data[q][i][j][v] =
6709 ->jacobians[0][offsets + q_src][i][j][lane];
6711 const auto &update_flags_cells =
6716 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6717 for (
unsigned int j = 0; j < dim; ++j)
6718 this_jacobian_gradients_data[q][i][j][v] =
6720 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6722 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6723 for (
unsigned int j = 0; j < dim; ++j)
6724 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6726 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6739 this->mapping_data->quadrature_points
6741 ->quadrature_point_offsets[cell_batch_index] +
6745 this->mapping_data->jacobians[0][offsets + 1];
6747 for (
unsigned int d = 0;
d < dim; ++
d)
6750 static_cast<Number
>(
6751 this->descriptor->quadrature.point(q)[
d]);
6753 for (
unsigned int d = 0;
d < dim; ++
d)
6754 for (
unsigned int e = 0;
e < dim; ++
e)
6757 static_cast<Number
>(
6758 this->descriptor->quadrature.point(q)[
e]);
6760 for (
unsigned int i = 0; i < dim; ++i)
6761 this_quadrature_points_data[q][i][v] = point[i][lane];
6766 for (
unsigned int i = 0; i < dim; ++i)
6767 this_quadrature_points_data[q][i][v] =
6768 this->mapping_data->quadrature_points
6770 ->quadrature_point_offsets[cell_batch_index] +
6779 this->is_reinitialized =
true;
6780 this->dof_values_initialized =
false;
6781 this->values_quad_initialized =
false;
6782 this->gradients_quad_initialized =
false;
6783 this->hessians_quad_initialized =
false;
6794 typename VectorizedArrayType>
6795template <
bool level_dof_access>
6802 VectorizedArrayType>
::
6805 Assert(this->matrix_free ==
nullptr,
6806 ExcMessage(
"Cannot use initialization from cell iterator if "
6807 "initialized from MatrixFree object. Use variant for "
6808 "on the fly computation with arguments as for FEValues "
6811 this->mapped_geometry->reinit(
6813 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6814 if (level_dof_access)
6815 cell->get_mg_dof_indices(this->local_dof_indices);
6817 cell->get_dof_indices(this->local_dof_indices);
6820 this->is_reinitialized =
true;
6831 typename VectorizedArrayType>
6838 VectorizedArrayType>
::
6841 Assert(this->matrix_free == 0,
6842 ExcMessage(
"Cannot use initialization from cell iterator if "
6843 "initialized from MatrixFree object. Use variant for "
6844 "on the fly computation with arguments as for FEValues "
6847 this->mapped_geometry->reinit(cell);
6850 this->is_reinitialized =
true;
6861 typename VectorizedArrayType>
6868 VectorizedArrayType>::
6872 Assert(this->dof_values_initialized ==
true,
6875 evaluate(this->values_dofs, evaluation_flags);
6885 typename VectorizedArrayType>
6892 VectorizedArrayType>::
6893 evaluate(
const VectorizedArrayType *values_array,
6896 const bool hessians_on_general_cells =
6900 if (hessians_on_general_cells)
6903 if (this->
data->element_type ==
6909 if constexpr (fe_degree > -1)
6912 template run<fe_degree, n_q_points_1d>(n_components,
6913 evaluation_flag_actual,
6921 evaluation_flag_actual,
6922 const_cast<VectorizedArrayType *
>(values_array),
6927 this->values_quad_initialized =
6929 this->gradients_quad_initialized =
6931 this->hessians_quad_initialized =
6942 template <
typename Number,
6943 typename VectorizedArrayType,
6945 typename EvaluatorType,
6946 std::enable_if_t<internal::has_begin<VectorType> &&
6949 VectorizedArrayType *
6950 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6956 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6957 const auto &dof_info = fe_eval.get_dof_info();
6964 if (std::is_same_v<typename VectorType::value_type, Number> &&
6968 interleaved_contiguous &&
6969 reinterpret_cast<
std::size_t>(
6971 dof_info.dof_indices_contiguous
6972 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6973 [cell * VectorizedArrayType::
size()]) %
6974 sizeof(VectorizedArrayType) ==
6977 return reinterpret_cast<VectorizedArrayType *
>(
6981 [cell * VectorizedArrayType::size()] +
6983 [fe_eval.get_active_fe_index()]
6984 [fe_eval.get_first_selected_component()] *
6985 VectorizedArrayType::size());
6994 template <
typename Number,
6995 typename VectorizedArrayType,
6997 typename EvaluatorType,
6998 std::enable_if_t<!internal::has_begin<VectorType> ||
7001 VectorizedArrayType *
7002 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7015 typename VectorizedArrayType>
7016template <
typename VectorType>
7023 VectorizedArrayType>::
7024 gather_evaluate(
const VectorType &input_vector,
7027 const VectorizedArrayType *src_ptr =
7028 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7029 *
this, input_vector);
7030 if (src_ptr !=
nullptr)
7031 evaluate(src_ptr, evaluation_flag);
7034 this->read_dof_values(input_vector);
7035 evaluate(this->begin_dof_values(), evaluation_flag);
7046 typename VectorizedArrayType>
7053 VectorizedArrayType>::
7056 integrate(integration_flag, this->values_dofs);
7059 this->dof_values_initialized =
true;
7070 typename VectorizedArrayType>
7077 VectorizedArrayType>::
7079 VectorizedArrayType *values_array,
7080 const bool sum_into_values_array)
7084 Assert(this->values_quad_submitted ==
true,
7087 Assert(this->gradients_quad_submitted ==
true,
7090 Assert(this->hessians_quad_submitted ==
true,
7093 Assert(this->matrix_free !=
nullptr ||
7094 this->mapped_geometry->is_initialized(),
7100 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7101 "EvaluationFlags::hessians are supported."));
7107 unsigned int size = n_components * dim * n_q_points;
7110 for (
unsigned int i = 0; i <
size; ++i)
7111 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7115 for (
unsigned int i = 0; i <
size; ++i)
7116 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7121 if (n_components == dim &&
7122 this->
data->element_type ==
7126 this->divergence_is_requested ==
false)
7128 unsigned int size = n_components * n_q_points;
7131 for (
unsigned int i = 0; i <
size; ++i)
7132 this->values_quad[i] += this->values_from_gradients_quad[i];
7136 for (
unsigned int i = 0; i <
size; ++i)
7137 this->values_quad[i] = this->values_from_gradients_quad[i];
7142 if constexpr (fe_degree > -1)
7145 template run<fe_degree, n_q_points_1d>(n_components,
7146 integration_flag_actual,
7149 sum_into_values_array);
7155 integration_flag_actual,
7158 sum_into_values_array);
7162 this->dof_values_initialized =
true;
7173 typename VectorizedArrayType>
7174template <
typename VectorType>
7181 VectorizedArrayType>::
7183 VectorType &destination)
7185 VectorizedArrayType *dst_ptr =
7186 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7187 *
this, destination);
7188 if (dst_ptr !=
nullptr)
7189 integrate(integration_flag, dst_ptr,
true);
7192 integrate(integration_flag, this->begin_dof_values());
7193 this->distribute_local_to_global(destination);
7204 typename VectorizedArrayType>
7211 VectorizedArrayType>::dof_indices()
const
7213 return {0U, dofs_per_cell};
7227 typename VectorizedArrayType>
7233 VectorizedArrayType>::
7236 const bool is_interior_face,
7237 const unsigned int dof_no,
7238 const unsigned int quad_no,
7239 const unsigned int first_selected_component,
7240 const unsigned int active_fe_index,
7241 const unsigned int active_quad_index,
7242 const unsigned int face_type)
7243 : BaseClass(matrix_free,
7245 first_selected_component,
7253 , dofs_per_component(this->
data->dofs_per_component_on_cell)
7254 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
7255 , n_q_points(this->n_quadrature_points)
7265 typename VectorizedArrayType>
7271 VectorizedArrayType>::