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();
157 template <
typename VectorType>
160 const VectorType &src,
161 const unsigned int first_index = 0,
162 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
195 template <
typename VectorType>
198 const VectorType &src,
199 const unsigned int first_index = 0,
200 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
233 template <
typename VectorType>
237 const unsigned int first_index = 0,
238 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
278 template <
typename VectorType>
282 const unsigned int first_index = 0,
283 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
288 template <
typename VectorType>
292 const unsigned int first_index = 0,
293 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
369 typename = std::enable_if_t<n_components == n_components_local>>
372 const unsigned int q_point);
423 template <
int dim_ = dim,
424 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
427 const unsigned int q_point);
446 const unsigned int q_point);
522 const unsigned int q_point);
531 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
550 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
553 const unsigned int q_point);
563 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
582 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
586 const unsigned int q_point);
596 template <
int dim_ = dim,
597 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
598 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
614 template <
int dim_ = dim,
615 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
618 const unsigned int q_point);
659 const unsigned int dof_no,
662 const unsigned int fe_degree,
663 const unsigned int n_q_points,
667 const unsigned int face_type);
741 template <
typename VectorType,
typename VectorOperation>
745 const std::array<VectorType *, n_components_> &vectors,
748 n_components_> &vectors_sm,
749 const std::bitset<n_lanes> &mask,
750 const bool apply_constraints =
true)
const;
759 template <
typename VectorType,
typename VectorOperation>
763 const std::array<VectorType *, n_components_> &vectors,
766 n_components_> &vectors_sm,
767 const std::bitset<n_lanes> &mask)
const;
776 template <
typename VectorType,
typename VectorOperation>
780 const std::array<VectorType *, n_components_> &vectors)
const;
1384 typename VectorizedArrayType>
1389 VectorizedArrayType>
1392 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1393 "Type of Number and of VectorizedArrayType do not match.");
1435 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1520 const unsigned int dof_no = 0,
1521 const unsigned int quad_no = 0,
1534 const std::pair<unsigned int, unsigned int> &range,
1535 const unsigned int dof_no = 0,
1536 const unsigned int quad_no = 0,
1646 template <
bool level_dof_access>
1668 const unsigned int given_n_q_points_1d);
1711 template <
typename VectorType>
1741 VectorizedArrayType *values_array,
1742 const bool sum_into_values =
false);
1757 template <
typename VectorType>
1760 VectorType &output_vector);
1844 int n_q_points_1d = fe_degree + 1,
1845 int n_components_ = 1,
1846 typename Number = double,
1852 VectorizedArrayType>
1855 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1856 "Type of Number and of VectorizedArrayType do not match.");
1898 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
2000 const unsigned int dof_no = 0,
2001 const unsigned int quad_no = 0,
2016 const std::pair<unsigned int, unsigned int> &range,
2018 const unsigned int dof_no = 0,
2019 const unsigned int quad_no = 0,
2033 reinit(
const unsigned int face_batch_number);
2043 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2050 const unsigned int given_n_q_points_1d);
2114 template <
typename VectorType>
2130 const bool sum_into_values =
false);
2143 VectorizedArrayType *values_array,
2144 const bool sum_into_values =
false);
2161 const bool sum_into_values =
false);
2169 VectorizedArrayType *values_array,
2170 const bool sum_into_values =
false);
2183 template <
typename VectorType>
2186 VectorType &output_vector);
2191 template <
typename VectorType>
2194 const bool integrate_gradients,
2195 VectorType &output_vector);
2271 namespace MatrixFreeFunctions
2275 template <
int dim,
int degree>
2284 template <
int degree>
2287 static constexpr unsigned int value = degree + 1;
2303 template <
bool is_face,
2306 typename VectorizedArrayType>
2309 extract_initialization_data(
2311 const unsigned int dof_no,
2312 const unsigned int first_selected_component,
2313 const unsigned int quad_no,
2314 const unsigned int fe_degree,
2315 const unsigned int n_q_points,
2316 const unsigned int active_fe_index_given,
2317 const unsigned int active_quad_index_given,
2318 const unsigned int face_type)
2321 InitializationData init_data;
2325 &internal::MatrixFreeFunctions::
2326 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2334 active_fe_index_given :
2340 active_quad_index_given :
2341 std::min<unsigned int>(
2344 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2373 typename VectorizedArrayType>
2378 VectorizedArrayType>::
2381 const unsigned int dof_no,
2382 const unsigned int first_selected_component,
2383 const unsigned int quad_no,
2384 const unsigned int fe_degree,
2385 const unsigned int n_q_points,
2386 const bool is_interior_face,
2387 const unsigned int active_fe_index,
2388 const unsigned int active_quad_index,
2389 const unsigned int face_type)
2391 internal::extract_initialization_data<is_face>(matrix_free,
2393 first_selected_component,
2402 first_selected_component)
2403 , scratch_data_array(matrix_free.acquire_scratch_data())
2404 , matrix_free(&matrix_free)
2406 this->set_data_pointers(scratch_data_array, n_components_);
2408 this->dof_info->start_components.back() == 1 ||
2409 static_cast<int>(n_components_) <=
2411 this->dof_info->start_components
2412 [this->dof_info->component_to_base_index[first_selected_component] +
2414 first_selected_component,
2416 "You tried to construct a vector-valued evaluator with " +
2417 std::to_string(n_components) +
2418 " components. However, "
2419 "the current base element has only " +
2421 this->dof_info->start_components
2422 [this->dof_info->component_to_base_index[first_selected_component] +
2424 first_selected_component) +
2425 " components left when starting from local element index " +
2427 first_selected_component -
2428 this->dof_info->start_components
2429 [this->dof_info->component_to_base_index[first_selected_component]]) +
2430 " (global index " + std::to_string(first_selected_component) +
")"));
2442 typename VectorizedArrayType>
2447 VectorizedArrayType>::
2453 const unsigned int first_selected_component,
2457 other->mapped_geometry->get_quadrature() == quadrature ?
2458 other->mapped_geometry :
2460 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2465 first_selected_component)
2466 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2467 , matrix_free(nullptr)
2469 const unsigned int base_element_number =
2473 first_selected_component >=
2475 ExcMessage(
"The underlying element must at least contain as many "
2476 "components as requested by this class"));
2477 (void)base_element_number;
2481 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2483 fe.component_to_base_index(first_selected_component).
first);
2485 this->set_data_pointers(scratch_data_array, n_components_);
2494 typename VectorizedArrayType>
2499 VectorizedArrayType>::
2504 VectorizedArrayType> &other)
2506 , scratch_data_array(other.matrix_free == nullptr ?
2508 other.matrix_free->acquire_scratch_data())
2509 , matrix_free(other.matrix_free)
2511 if (other.matrix_free ==
nullptr)
2518 this->mapped_geometry =
2519 std::make_shared<internal::MatrixFreeFunctions::
2520 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2525 if constexpr (is_face ==
false)
2526 this->mapping_data = &this->mapped_geometry->get_data_storage();
2530 "face evaluators is not currently "
2536 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2538 this->mapped_geometry->get_data_storage().JxW_values.begin();
2539 this->jacobian_gradients =
2540 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2541 this->jacobian_gradients_non_inverse =
2542 this->mapped_geometry->get_data_storage()
2543 .jacobian_gradients_non_inverse[0]
2546 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2549 this->set_data_pointers(scratch_data_array, n_components_);
2558 typename VectorizedArrayType>
2563 VectorizedArrayType> &
2569 VectorizedArrayType> &other)
2572 if (matrix_free ==
nullptr)
2575 delete scratch_data_array;
2584 matrix_free = other.matrix_free;
2586 if (other.matrix_free ==
nullptr)
2594 this->mapped_geometry =
2595 std::make_shared<internal::MatrixFreeFunctions::
2596 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2601 if constexpr (is_face ==
false)
2602 this->mapping_data = &this->mapped_geometry->get_data_storage();
2606 "face evaluators is not currently "
2611 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2613 this->mapped_geometry->get_data_storage().JxW_values.begin();
2614 this->jacobian_gradients =
2615 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2616 this->jacobian_gradients_non_inverse =
2617 this->mapped_geometry->get_data_storage()
2618 .jacobian_gradients_non_inverse[0]
2621 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2628 this->set_data_pointers(scratch_data_array, n_components_);
2639 typename VectorizedArrayType>
2644 VectorizedArrayType>::~FEEvaluationBase()
2646 if (matrix_free !=
nullptr)
2657 delete scratch_data_array;
2668 typename VectorizedArrayType>
2673 Assert(matrix_free !=
nullptr,
2675 "FEEvaluation was not initialized with a MatrixFree object!"));
2676 return *matrix_free;
2685 template <
typename VectorType,
bool>
2686 struct ConstBlockVectorSelector;
2688 template <
typename VectorType>
2689 struct ConstBlockVectorSelector<
VectorType, true>
2691 using BaseVectorType =
const typename VectorType::BlockType;
2694 template <
typename VectorType>
2695 struct ConstBlockVectorSelector<
VectorType, false>
2697 using BaseVectorType =
typename VectorType::BlockType;
2703 template <
typename VectorType,
bool>
2704 struct BlockVectorSelector;
2706 template <
typename VectorType>
2709 using BaseVectorType =
typename ConstBlockVectorSelector<
2711 std::is_const_v<VectorType>>::BaseVectorType;
2713 static BaseVectorType *
2714 get_vector_component(VectorType &vec,
const unsigned int component)
2717 return &vec.block(component);
2721 template <
typename VectorType>
2722 struct BlockVectorSelector<
VectorType, false>
2726 static BaseVectorType *
2727 get_vector_component(VectorType &vec,
const unsigned int component)
2746 template <
typename VectorType>
2747 struct BlockVectorSelector<
std::vector<VectorType>, false>
2751 static BaseVectorType *
2752 get_vector_component(std::vector<VectorType> &vec,
2753 const unsigned int component)
2756 return &vec[component];
2760 template <
typename VectorType>
2761 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2765 static const BaseVectorType *
2766 get_vector_component(
const std::vector<VectorType> &vec,
2767 const unsigned int component)
2770 return &vec[component];
2774 template <
typename VectorType>
2775 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2779 static BaseVectorType *
2780 get_vector_component(std::vector<VectorType *> &vec,
2781 const unsigned int component)
2784 return vec[component];
2788 template <
typename VectorType>
2789 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2793 static const BaseVectorType *
2794 get_vector_component(
const std::vector<VectorType *> &vec,
2795 const unsigned int component)
2798 return vec[component];
2802 template <
typename VectorType, std::
size_t N>
2803 struct BlockVectorSelector<
std::array<VectorType *, N>, false>
2807 static BaseVectorType *
2808 get_vector_component(std::array<VectorType *, N> &vec,
2809 const unsigned int component)
2812 return vec[component];
2823 typename VectorizedArrayType>
2824template <
typename VectorType,
typename VectorOperation>
2829 const std::array<VectorType *, n_components_> &src,
2832 n_components_> &src_sm,
2833 const std::bitset<n_lanes> &
mask,
2834 const bool apply_constraints)
const
2839 if (this->matrix_free ==
nullptr)
2841 read_write_operation_global(operation, src);
2848 if (this->n_fe_components == 1)
2849 for (
unsigned int comp = 0; comp < n_components; ++comp)
2851 Assert(src[comp] !=
nullptr,
2852 ExcMessage(
"The finite element underlying this FEEvaluation "
2853 "object is scalar, but you requested " +
2854 std::to_string(n_components) +
2855 " components via the template argument in "
2856 "FEEvaluation. In that case, you must pass an "
2857 "std::vector<VectorType> or a BlockVector to " +
2858 "read_dof_values and distribute_local_to_global."));
2870 const bool accesses_exterior_dofs =
2871 this->dof_access_index ==
2873 this->is_interior_face() ==
false;
2883 bool is_contiguous =
true;
2885 if (accesses_exterior_dofs)
2887 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2888 const unsigned int n_filled_lanes =
2893 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2894 if (
mask[v] ==
true &&
2897 [cells[v] / n_lanes] <
2900 is_contiguous = false;
2904 this->dof_access_index :
2905 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2906 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2909 is_contiguous =
false;
2914 read_write_operation_contiguous(operation, src, src_sm,
mask);
2921 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2923 const bool masking_is_active =
mask.count() < n_lanes;
2924 if (masking_is_active)
2925 for (
unsigned int v = 0; v < n_lanes; ++v)
2926 if (
mask[v] ==
false)
2929 bool has_hn_constraints =
false;
2931 if (is_face ==
false)
2937 [this->first_selected_component])
2938 for (
unsigned int v = 0; v < n_lanes; ++v)
2943 has_hn_constraints = true;
2946 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2949 const bool use_vectorized_path =
2950 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2952 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
2953 std::array<VectorizedArrayType *, n_components> values_dofs;
2954 for (
unsigned int c = 0; c < n_components; ++c)
2955 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2956 c * dofs_per_component;
2960 [is_face ? this->dof_access_index :
2961 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2962 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2963 IndexStorageVariants::interleaved &&
2964 use_vectorized_path)
2966 const unsigned int *dof_indices =
2968 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2972 [this->first_selected_component] *
2975 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2976 if (n_components == 1 || this->n_fe_components == 1)
2977 for (
unsigned int comp = 0; comp < n_components; ++comp)
2979 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2982 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2984 if (n_components == 1 || this->n_fe_components == 1)
2985 for (
unsigned int i = 0; i < dofs_per_component;
2986 ++i, dof_indices += n_lanes)
2987 for (
unsigned int comp = 0; comp < n_components; ++comp)
2988 operation.process_dof_gather(dof_indices,
2992 values_dofs[comp][i],
2995 for (
unsigned int comp = 0; comp < n_components; ++comp)
2996 for (
unsigned int i = 0; i < dofs_per_component;
2997 ++i, dof_indices += n_lanes)
2998 operation.process_dof_gather(dof_indices,
3002 values_dofs[comp][i],
3009 std::array<const unsigned int *, n_lanes> dof_indices;
3010 dof_indices.fill(
nullptr);
3015 bool has_constraints =
false;
3016 const unsigned int n_components_read =
3017 this->n_fe_components > 1 ? n_components : 1;
3021 for (
unsigned int v = 0; v < n_lanes; ++v)
3027 const std::pair<unsigned int, unsigned int> *my_index_start =
3028 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3029 this->first_selected_component];
3034 if (my_index_start[n_components_read].
second !=
3035 my_index_start[0].
second)
3036 has_constraints =
true;
3039 dof_info.
dof_indices.data() + my_index_start[0].first;
3044 for (
unsigned int v = 0; v < n_lanes; ++v)
3049 const std::pair<unsigned int, unsigned int> *my_index_start =
3050 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3051 this->first_selected_component];
3052 if (my_index_start[n_components_read].
second !=
3053 my_index_start[0].
second)
3054 has_constraints =
true;
3061 dof_info.hanging_node_constraint_masks_comp
3062 [this->active_fe_index][this->first_selected_component])
3063 has_hn_constraints = true;
3066 my_index_start[0].
first ||
3069 my_index_start[0].
first,
3072 dof_info.
dof_indices.data() + my_index_start[0].first;
3076 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3077 return i != numbers::invalid_unsigned_int;
3079 for (
unsigned int comp = 0; comp < n_components; ++comp)
3080 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3081 operation.process_empty(values_dofs[comp][i]);
3085 if (!has_constraints && apply_constraints)
3087 if (n_components == 1 || this->n_fe_components == 1)
3089 for (
unsigned int v = 0; v < n_lanes; ++v)
3094 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3095 for (
unsigned int comp = 0; comp < n_components; ++comp)
3096 operation.process_dof(dof_indices[v][i],
3098 values_dofs[comp][i][v]);
3103 for (
unsigned int comp = 0; comp < n_components; ++comp)
3104 for (
unsigned int v = 0; v < n_lanes; ++v)
3109 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3110 operation.process_dof(
3111 dof_indices[v][comp * dofs_per_component + i],
3113 values_dofs[comp][i][v]);
3125 for (
unsigned int v = 0; v < n_lanes; ++v)
3131 const unsigned int cell_dof_index =
3132 cell_index * this->n_fe_components + this->first_selected_component;
3133 const unsigned int n_components_read =
3134 this->n_fe_components > 1 ? n_components : 1;
3135 unsigned int index_indicators =
3137 unsigned int next_index_indicators =
3138 dof_info.
row_starts[cell_dof_index + 1].second;
3142 if (apply_constraints ==
false &&
3143 (dof_info.
row_starts[cell_dof_index].second !=
3144 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3150 dof_info.hanging_node_constraint_masks_comp
3151 [this->active_fe_index][this->first_selected_component])))
3160 [this->first_selected_component] +
3162 next_index_indicators = index_indicators;
3165 if (n_components == 1 || this->n_fe_components == 1)
3167 unsigned int ind_local = 0;
3168 for (; index_indicators != next_index_indicators; ++index_indicators)
3170 const std::pair<unsigned short, unsigned short> indicator =
3173 for (
unsigned int j = 0; j < indicator.first; ++j)
3174 for (
unsigned int comp = 0; comp < n_components; ++comp)
3175 operation.process_dof(dof_indices[v][j],
3177 values_dofs[comp][ind_local + j][v]);
3179 ind_local += indicator.first;
3180 dof_indices[v] += indicator.first;
3184 Number
value[n_components];
3185 for (
unsigned int comp = 0; comp < n_components; ++comp)
3186 operation.pre_constraints(values_dofs[comp][ind_local][v],
3189 const Number *data_val =
3191 const Number *end_pool =
3193 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3194 for (
unsigned int comp = 0; comp < n_components; ++comp)
3195 operation.process_constraint(*dof_indices[v],
3200 for (
unsigned int comp = 0; comp < n_components; ++comp)
3201 operation.post_constraints(
value[comp],
3202 values_dofs[comp][ind_local][v]);
3208 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3209 for (
unsigned int comp = 0; comp < n_components; ++comp)
3210 operation.process_dof(*dof_indices[v],
3212 values_dofs[comp][ind_local][v]);
3221 for (
unsigned int comp = 0; comp < n_components; ++comp)
3223 unsigned int ind_local = 0;
3226 for (; index_indicators != next_index_indicators;
3229 const std::pair<unsigned short, unsigned short> indicator =
3233 for (
unsigned int j = 0; j < indicator.first; ++j)
3234 operation.process_dof(dof_indices[v][j],
3236 values_dofs[comp][ind_local + j][v]);
3237 ind_local += indicator.first;
3238 dof_indices[v] += indicator.first;
3243 operation.pre_constraints(values_dofs[comp][ind_local][v],
3246 const Number *data_val =
3248 const Number *end_pool =
3251 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3252 operation.process_constraint(*dof_indices[v],
3257 operation.post_constraints(
value,
3258 values_dofs[comp][ind_local][v]);
3265 for (; ind_local < dofs_per_component;
3266 ++dof_indices[v], ++ind_local)
3269 operation.process_dof(*dof_indices[v],
3271 values_dofs[comp][ind_local][v]);
3274 if (apply_constraints ==
true && comp + 1 < n_components)
3275 next_index_indicators =
3276 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3288 typename VectorizedArrayType>
3289template <
typename VectorType,
typename VectorOperation>
3294 const std::array<VectorType *, n_components_> &src)
const
3298 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3299 unsigned int index = this->first_selected_component * dofs_per_component;
3300 for (
unsigned int comp = 0; comp < n_components; ++comp)
3302 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3304 operation.process_empty(
3305 this->values_dofs[comp * dofs_per_component + i]);
3306 operation.process_dof_global(
3307 local_dof_indices[this->
data->lexicographic_numbering[index]],
3309 this->values_dofs[comp * dofs_per_component + i][0]);
3320 typename VectorizedArrayType>
3321template <
typename VectorType,
typename VectorOperation>
3326 const std::array<VectorType *, n_components_> &src,
3329 n_components_> &vectors_sm,
3330 const std::bitset<n_lanes> &
mask)
const
3339 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3342 is_face ? this->dof_access_index :
3344 const unsigned int n_active_lanes =
mask.count();
3347 const std::vector<unsigned int> &dof_indices_cont =
3350 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3351 std::array<VectorizedArrayType *, n_components> values_dofs{{
nullptr}};
3352 for (
unsigned int c = 0; c < n_components; ++c)
3353 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3354 c * dofs_per_component;
3358 const bool accesses_exterior_dofs =
3359 this->dof_access_index ==
3361 this->is_interior_face() ==
false;
3367 interleaved_contiguous &&
3368 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3370 const unsigned int dof_index =
3371 dof_indices_cont[this->cell * n_lanes] +
3374 [this->first_selected_component] *
3376 if (n_components == 1 || this->n_fe_components == 1)
3377 for (
unsigned int comp = 0; comp < n_components; ++comp)
3378 operation.process_dofs_vectorized(dofs_per_component,
3384 operation.process_dofs_vectorized(dofs_per_component * n_components,
3392 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3396 const unsigned int n_filled_lanes =
3399 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3400 n_active_lanes == n_lanes &&
3401 !accesses_exterior_dofs;
3403 if (vectors_sm[0] !=
nullptr)
3405 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3406 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3409 const auto upper_bound =
3410 std::min<unsigned int>(n_filled_lanes, n_lanes);
3411 for (
unsigned int v = 0; v < upper_bound; ++v)
3413 if (
mask[v] ==
false)
3415 vector_ptrs[v] =
nullptr;
3435 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3436 vectors_sm[comp]->operator[](temp.first).
data() + temp.second +
3438 [this->active_fe_index][this->first_selected_component]);
3440 vector_ptrs[v] =
nullptr;
3442 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3443 vector_ptrs[v] =
nullptr;
3448 if (use_vectorized_path)
3450 if (n_components == 1 || this->n_fe_components == 1)
3452 for (
unsigned int comp = 0; comp < n_components; ++comp)
3454 auto vector_ptrs = compute_vector_ptrs(comp);
3455 operation.process_dofs_vectorized_transpose(
3464 auto vector_ptrs = compute_vector_ptrs(0);
3465 operation.process_dofs_vectorized_transpose(dofs_per_component *
3473 for (
unsigned int comp = 0; comp < n_components; ++comp)
3475 auto vector_ptrs = compute_vector_ptrs(
3476 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3478 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3479 operation.process_empty(values_dofs[comp][i]);
3481 if (n_components == 1 || this->n_fe_components == 1)
3483 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3484 if (
mask[v] ==
true)
3485 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3486 operation.process_dof(vector_ptrs[v][i],
3487 values_dofs[comp][i][v]);
3491 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3492 if (
mask[v] ==
true)
3493 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3494 operation.process_dof(
3495 vector_ptrs[v][i + comp * dofs_per_component],
3496 values_dofs[comp][i][v]);
3502 std::array<unsigned int, n_lanes> dof_indices{
3505 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3509 if (
mask[v] ==
true)
3511 dof_indices_cont[cells[v]] +
3514 [this->first_selected_component] *
3520 if (use_vectorized_path)
3526 if (n_components == 1 || this->n_fe_components == 1)
3527 for (
unsigned int comp = 0; comp < n_components; ++comp)
3528 operation.process_dofs_vectorized_transpose(dofs_per_component,
3534 operation.process_dofs_vectorized_transpose(dofs_per_component *
3543 interleaved_contiguous_strided)
3545 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3547 if (n_components == 1 || this->n_fe_components == 1)
3548 for (
unsigned int comp = 0; comp < n_components; ++comp)
3549 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3550 src[comp]->
begin());
3553 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3555 if (n_components == 1 || this->n_fe_components == 1)
3556 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3558 for (
unsigned int comp = 0; comp < n_components; ++comp)
3559 operation.process_dof_gather(dof_indices.data(),
3562 src_ptrs[comp] + i * n_lanes,
3563 values_dofs[comp][i],
3567 for (
unsigned int comp = 0; comp < n_components; ++comp)
3568 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3570 operation.process_dof_gather(
3573 (comp * dofs_per_component + i) * n_lanes,
3574 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3575 values_dofs[comp][i],
3583 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3585 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3587 if (n_components == 1 || this->n_fe_components == 1)
3588 for (
unsigned int comp = 0; comp < n_components; ++comp)
3589 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3590 src[comp]->
begin());
3593 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3595 const unsigned int *offsets =
3597 if (n_components == 1 || this->n_fe_components == 1)
3598 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3600 for (
unsigned int comp = 0; comp < n_components; ++comp)
3601 operation.process_dof_gather(dof_indices.data(),
3605 values_dofs[comp][i],
3608 for (
unsigned int v = 0; v < n_lanes; ++v)
3609 dof_indices[v] += offsets[v];
3612 for (
unsigned int comp = 0; comp < n_components; ++comp)
3613 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3615 operation.process_dof_gather(dof_indices.data(),
3619 values_dofs[comp][i],
3622 for (
unsigned int v = 0; v < n_lanes; ++v)
3623 dof_indices[v] += offsets[v];
3628 for (
unsigned int comp = 0; comp < n_components; ++comp)
3630 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3631 operation.process_empty(values_dofs[comp][i]);
3632 if (accesses_exterior_dofs)
3634 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3635 if (
mask[v] ==
true)
3638 [ind][cells[v] / VectorizedArrayType::size()] ==
3642 if (n_components == 1 || this->n_fe_components == 1)
3644 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3645 operation.process_dof(dof_indices[v] + i,
3647 values_dofs[comp][i][v]);
3651 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3652 operation.process_dof(dof_indices[v] + i +
3653 comp * dofs_per_component,
3655 values_dofs[comp][i][v]);
3660 const unsigned int offset =
3663 if (n_components == 1 || this->n_fe_components == 1)
3665 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3666 operation.process_dof(dof_indices[v] + i * offset,
3668 values_dofs[comp][i][v]);
3672 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3673 operation.process_dof(
3675 (i + comp * dofs_per_component) * offset,
3677 values_dofs[comp][i][v]);
3688 if (n_components == 1 || this->n_fe_components == 1)
3690 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3691 if (
mask[v] ==
true)
3692 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3693 operation.process_dof(dof_indices[v] + i,
3695 values_dofs[comp][i][v]);
3699 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3700 if (
mask[v] ==
true)
3701 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3702 operation.process_dof(dof_indices[v] + i +
3703 comp * dofs_per_component,
3705 values_dofs[comp][i][v]);
3710 const unsigned int *offsets =
3712 [ind][VectorizedArrayType::size() * this->cell];
3713 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3715 if (n_components == 1 || this->n_fe_components == 1)
3716 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3718 if (
mask[v] ==
true)
3719 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3720 operation.process_dof(dof_indices[v] + i * offsets[v],
3722 values_dofs[comp][i][v]);
3726 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3727 if (
mask[v] ==
true)
3728 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3729 operation.process_dof(
3731 (i + comp * dofs_per_component) * offsets[v],
3733 values_dofs[comp][i][v]);
3745 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3746 decltype(std::declval<VectorType>().begin())
3747 get_beginning(VectorType &vec)
3755 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3756 typename VectorType::value_type *
3757 get_beginning(VectorType &)
3763 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3765 const std::vector<ArrayView<const typename VectorType::value_type>> *
3766 get_shared_vector_data(VectorType *vec,
3767 const bool is_valid_mode_for_sm,
3768 const unsigned int active_fe_index,
3772 if (is_valid_mode_for_sm &&
3775 active_fe_index == 0)
3776 return &vec->shared_vector_data();
3782 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3784 const std::vector<ArrayView<const typename VectorType::value_type>> *
3785 get_shared_vector_data(VectorType *,
3793 template <
int n_components,
typename VectorType>
3795 std::array<
typename internal::BlockVectorSelector<
3800 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3804 get_vector_data(VectorType &src,
3805 const unsigned int first_index,
3806 const bool is_valid_mode_for_sm,
3807 const unsigned int active_fe_index,
3813 std::array<
typename internal::BlockVectorSelector<
3819 ArrayView<
const typename internal::BlockVectorSelector<
3825 for (
unsigned int d = 0;
d < n_components; ++
d)
3826 src_data.first[d] = internal::BlockVectorSelector<
3832 for (
unsigned int d = 0;
d < n_components; ++
d)
3833 src_data.second[d] = get_shared_vector_data(
3834 const_cast<typename internal::BlockVectorSelector<
3835 std::remove_const_t<VectorType>,
3837 *>(src_data.first[d]),
3838 is_valid_mode_for_sm,
3852 typename VectorizedArrayType>
3857 if (this->dof_info ==
nullptr ||
3859 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3860 this->dof_info->hanging_node_constraint_masks_comp
3861 [this->active_fe_index][this->first_selected_component] ==
false)
3864 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3865 constraint_mask{{internal::MatrixFreeFunctions::
3866 unconstrained_compressed_constraint_kind}};
3868 bool hn_available =
false;
3870 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3872 for (
unsigned int v = 0; v < n_lanes; ++v)
3884 constraint_mask[v] =
mask;
3890 if (hn_available ==
false)
3895 this->
data->data.front().fe_degree,
3896 this->get_shape_info(),
3908 typename VectorizedArrayType>
3909template <
typename VectorType>
3913 const unsigned int first_index,
3914 const std::bitset<n_lanes> &
mask)
3916 const auto src_data = internal::get_vector_data<n_components_>(
3919 this->dof_info !=
nullptr &&
3920 this->dof_access_index ==
3922 this->active_fe_index,
3926 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3928 apply_hanging_node_constraints(
false);
3932 this->dof_values_initialized =
true;
3942 typename VectorizedArrayType>
3943template <
typename VectorType>
3947 const unsigned int first_index,
3948 const std::bitset<n_lanes> &
mask)
3950 const auto src_data = internal::get_vector_data<n_components_>(
3953 this->dof_access_index ==
3955 this->active_fe_index,
3959 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3963 this->dof_values_initialized =
true;
3973 typename VectorizedArrayType>
3974template <
typename VectorType>
3978 const unsigned int first_index,
3979 const std::bitset<n_lanes> &
mask)
const
3983 Assert(this->dof_values_initialized ==
true,
3987 apply_hanging_node_constraints(
true);
3989 const auto dst_data = internal::get_vector_data<n_components_>(
3992 this->dof_access_index ==
3994 this->active_fe_index,
3999 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
4008 typename VectorizedArrayType>
4009template <
typename VectorType>
4013 const unsigned int first_index,
4014 const std::bitset<n_lanes> &
mask)
const
4018 Assert(this->dof_values_initialized ==
true,
4022 const auto dst_data = internal::get_vector_data<n_components_>(
4025 this->dof_access_index ==
4027 this->active_fe_index,
4031 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4040 typename VectorizedArrayType>
4041template <
typename VectorType>
4045 const unsigned int first_index,
4046 const std::bitset<n_lanes> &
mask)
const
4050 Assert(this->dof_values_initialized ==
true,
4054 const auto dst_data = internal::get_vector_data<n_components_>(
4057 this->dof_access_index ==
4059 this->active_fe_index,
4063 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4076 typename VectorizedArrayType>
4082 VectorizedArrayType>::value_type
4087 if constexpr (n_components == 1)
4088 return this->values_dofs[dof];
4091 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4092 Tensor<1, n_components_, VectorizedArrayType> return_value;
4093 for (
unsigned int comp = 0; comp < n_components; ++comp)
4094 return_value[comp] = this->values_dofs[comp * dofs + dof];
4095 return return_value;
4105 typename VectorizedArrayType>
4111 VectorizedArrayType>::value_type
4113 get_value(
const unsigned int q_point)
const
4117 Assert(this->values_quad_initialized ==
true,
4122 if constexpr (n_components == 1)
4123 return this->values_quad[q_point];
4126 if (n_components == dim &&
4127 this->
data->element_type ==
4133 Assert(this->values_quad_initialized ==
true,
4138 Assert(this->J_value !=
nullptr,
4141 const std::size_t nqp = this->n_quadrature_points;
4149 const VectorizedArrayType inv_det =
4150 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4151 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4152 this->jacobian[0][2][2];
4155 for (
unsigned int comp = 0; comp < n_components; ++comp)
4156 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4157 jac[comp][comp] * inv_det;
4164 this->jacobian[q_point] :
4173 const VectorizedArrayType inv_det =
4174 (is_face && dim == 2 && this->get_face_no() < 2) ?
4178 for (
unsigned int comp = 0; comp < n_components; ++comp)
4180 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4181 for (
unsigned int e = 1;
e < dim; ++
e)
4183 this->values_quad[e * nqp + q_point] * jac[comp][e];
4184 value_out[comp] *= inv_det;
4191 const std::size_t nqp = this->n_quadrature_points;
4193 for (
unsigned int comp = 0; comp < n_components; ++comp)
4194 return_value[comp] = this->values_quad[comp * nqp + q_point];
4195 return return_value;
4206 typename VectorizedArrayType>
4212 VectorizedArrayType>::gradient_type
4218 Assert(this->gradients_quad_initialized ==
true,
4223 Assert(this->jacobian !=
nullptr,
4225 "update_gradients"));
4226 const std::size_t nqp = this->n_quadrature_points;
4228 if constexpr (n_components == dim && dim > 1)
4230 if (this->
data->element_type ==
4236 Assert(this->gradients_quad_initialized ==
true,
4241 Assert(this->jacobian !=
nullptr,
4243 "update_gradients"));
4244 const std::size_t nqp = this->n_quadrature_points;
4245 const std::size_t nqp_d = nqp * dim;
4248 this->gradients_quad + q_point * dim;
4259 const VectorizedArrayType inv_det =
4260 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4261 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4262 this->jacobian[0][2][2];
4265 for (
unsigned int d = 0;
d < dim; ++
d)
4266 for (
unsigned int comp = 0; comp < n_components; ++comp)
4267 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4269 (jac[comp][comp] * inv_det);
4281 const VectorizedArrayType inv_det =
4282 (is_face && dim == 2 && this->get_face_no() < 2) ?
4286 VectorizedArrayType tmp[dim][dim];
4288 for (
unsigned int d = 0;
d < dim; ++
d)
4289 for (
unsigned int e = 0;
e < dim; ++
e)
4292 for (
unsigned int f = 1; f < dim; ++f)
4293 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4295 for (
unsigned int comp = 0; comp < n_components; ++comp)
4296 for (
unsigned int d = 0;
d < dim; ++
d)
4298 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4299 for (
unsigned int f = 1; f < dim; ++f)
4300 res += jac[comp][f] * tmp[d][f];
4302 grad_out[comp][
d] = res * inv_det;
4313 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4315 "update_hessians"));
4317 const auto jac_grad =
4318 this->jacobian_gradients_non_inverse[q_point];
4320 this->jacobian[q_point];
4324 const VectorizedArrayType inv_det =
4325 (is_face && dim == 2 && this->get_face_no() < 2) ?
4332 VectorizedArrayType tmp[dim][dim];
4333 for (
unsigned int d = 0;
d < dim; ++
d)
4334 for (
unsigned int e = 0;
e < dim; ++
e)
4337 for (
unsigned int f = 1; f < dim; ++f)
4338 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4343 for (
unsigned int d = 0;
d < dim; ++
d)
4345 for (
unsigned int e = 0;
e < dim; ++
e)
4347 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4348 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4349 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4352 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4354 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4359 for (
unsigned int d = 0;
d < dim; ++
d)
4360 for (
unsigned int e = 0;
e < dim; ++
e)
4362 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4363 for (
unsigned int f = 1; f < dim; ++f)
4364 res += tmp[f][d] * inv_t_jac[e][f];
4365 grad_out[
d][
e] = res;
4371 VectorizedArrayType tmp3[dim], tmp4[dim];
4372 for (
unsigned int d = 0;
d < dim; ++
d)
4374 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4375 for (
unsigned int e = 1;
e < dim; ++
e)
4376 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4378 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4379 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4380 for (
unsigned int d = 0;
d < dim; ++
d)
4382 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4383 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4385 for (
unsigned int d = 0;
d < dim; ++
d)
4387 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4388 for (
unsigned int e = 1;
e < dim; ++
e)
4389 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4392 VectorizedArrayType tmp2[dim];
4393 for (
unsigned int d = 0;
d < dim; ++
d)
4395 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4396 for (
unsigned e = 1;
e < dim; ++
e)
4398 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4401 for (
unsigned int d = 0;
d < dim; ++
d)
4402 for (
unsigned int e = 0;
e < dim; ++
e)
4404 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4408 grad_out[
d][
e] *= inv_det;
4419 for (
unsigned int comp = 0; comp < n_components; ++comp)
4420 for (
unsigned int d = 0;
d < dim; ++
d)
4422 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4423 this->jacobian[0][
d][
d];
4432 for (
unsigned int comp = 0; comp < n_components; ++comp)
4433 for (
unsigned int d = 0;
d < dim; ++
d)
4436 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4437 for (
unsigned int e = 1;
e < dim; ++
e)
4438 grad_out[comp][d] +=
4440 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4443 if constexpr (n_components == 1)
4455 typename VectorizedArrayType>
4461 VectorizedArrayType>::value_type
4468 Assert(this->gradients_quad_initialized ==
true,
4472 Assert(this->normal_x_jacobian !=
nullptr,
4474 "update_gradients"));
4476 const std::size_t nqp = this->n_quadrature_points;
4480 for (
unsigned int comp = 0; comp < n_components; ++comp)
4482 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4483 (this->normal_x_jacobian[0][dim - 1]);
4486 const std::size_t
index =
4488 for (
unsigned int comp = 0; comp < n_components; ++comp)
4490 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4491 this->normal_x_jacobian[
index][0];
4492 for (
unsigned int d = 1;
d < dim; ++
d)
4494 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4495 this->normal_x_jacobian[
index][
d];
4498 if constexpr (n_components == 1)
4510 template <
typename VectorizedArrayType>
4513 const VectorizedArrayType *
const hessians,
4515 VectorizedArrayType (&tmp)[1][1])
4517 tmp[0][0] = jac[0][0] *
hessians[0];
4520 template <
typename VectorizedArrayType>
4523 const VectorizedArrayType *
const hessians,
4524 const unsigned int nqp,
4525 VectorizedArrayType (&tmp)[2][2])
4527 for (
unsigned int d = 0;
d < 2; ++
d)
4535 template <
typename VectorizedArrayType>
4538 const VectorizedArrayType *
const hessians,
4539 const unsigned int nqp,
4540 VectorizedArrayType (&tmp)[3][3])
4542 for (
unsigned int d = 0;
d < 3; ++
d)
4563 typename VectorizedArrayType>
4568 VectorizedArrayType>::hessian_type
4574 Assert(this->hessians_quad_initialized ==
true,
4579 Assert(this->jacobian !=
nullptr,
4589 const std::size_t nqp = this->n_quadrature_points;
4590 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4595 for (
unsigned int comp = 0; comp < n_components; ++comp)
4597 for (
unsigned int d = 0;
d < dim; ++
d)
4598 hessian_out[comp][d][d] =
4599 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4600 (jac[
d][
d] * jac[
d][
d]);
4606 hessian_out[comp][0][1] =
4607 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4608 (jac[0][0] * jac[1][1]);
4611 hessian_out[comp][0][1] =
4612 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4613 (jac[0][0] * jac[1][1]);
4614 hessian_out[comp][0][2] =
4615 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4616 (jac[0][0] * jac[2][2]);
4617 hessian_out[comp][1][2] =
4618 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4619 (jac[1][1] * jac[2][2]);
4624 for (
unsigned int d = 0;
d < dim; ++
d)
4625 for (
unsigned int e = d + 1;
e < dim; ++
e)
4626 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4632 for (
unsigned int comp = 0; comp < n_components; ++comp)
4634 VectorizedArrayType tmp[dim][dim];
4635 internal::hessian_unit_times_jac(
4636 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4639 for (
unsigned int d = 0;
d < dim; ++
d)
4640 for (
unsigned int e = d;
e < dim; ++
e)
4642 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4643 for (
unsigned int f = 1; f < dim; ++f)
4644 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4651 for (
unsigned int d = 0;
d < dim; ++
d)
4652 for (
unsigned int e = d + 1;
e < dim; ++
e)
4653 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4659 const auto &jac_grad = this->jacobian_gradients[q_point];
4660 for (
unsigned int comp = 0; comp < n_components; ++comp)
4662 VectorizedArrayType tmp[dim][dim];
4663 internal::hessian_unit_times_jac(
4664 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4667 for (
unsigned int d = 0;
d < dim; ++
d)
4668 for (
unsigned int e = d;
e < dim; ++
e)
4670 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4671 for (
unsigned int f = 1; f < dim; ++f)
4672 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4676 for (
unsigned int d = 0;
d < dim; ++
d)
4677 for (
unsigned int e = 0;
e < dim; ++
e)
4678 hessian_out[comp][d][d] +=
4680 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4683 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4684 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4685 for (
unsigned int f = 0; f < dim; ++f)
4686 hessian_out[comp][d][e] +=
4687 jac_grad[count][f] *
4688 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4691 for (
unsigned int d = 0;
d < dim; ++
d)
4692 for (
unsigned int e = d + 1;
e < dim; ++
e)
4693 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4696 if constexpr (n_components == 1)
4697 return hessian_out[0];
4708 typename VectorizedArrayType>
4713 VectorizedArrayType>::gradient_type
4720 Assert(this->hessians_quad_initialized ==
true,
4731 const std::size_t nqp = this->n_quadrature_points;
4732 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4738 for (
unsigned int comp = 0; comp < n_components; ++comp)
4739 for (
unsigned int d = 0;
d < dim; ++
d)
4740 hessian_out[comp][d] =
4741 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4742 (jac[
d][
d] * jac[
d][
d]);
4747 for (
unsigned int comp = 0; comp < n_components; ++comp)
4751 VectorizedArrayType tmp[dim][dim];
4752 internal::hessian_unit_times_jac(
4753 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4757 for (
unsigned int d = 0;
d < dim; ++
d)
4759 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4760 for (
unsigned int f = 1; f < dim; ++f)
4761 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4768 const auto &jac_grad = this->jacobian_gradients[q_point];
4769 for (
unsigned int comp = 0; comp < n_components; ++comp)
4773 VectorizedArrayType tmp[dim][dim];
4774 internal::hessian_unit_times_jac(
4775 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4779 for (
unsigned int d = 0;
d < dim; ++
d)
4781 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4782 for (
unsigned int f = 1; f < dim; ++f)
4783 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4786 for (
unsigned int d = 0;
d < dim; ++
d)
4787 for (
unsigned int e = 0;
e < dim; ++
e)
4788 hessian_out[comp][d] +=
4790 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4794 if constexpr (n_components == 1)
4795 return hessian_out[0];
4806 typename VectorizedArrayType>
4811 VectorizedArrayType>::value_type
4818 Assert(this->hessians_quad_initialized ==
true,
4823 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4824 if constexpr (n_components == 1)
4826 VectorizedArrayType
sum = hess_diag[0];
4827 for (
unsigned int d = 1;
d < dim; ++
d)
4828 sum += hess_diag[d];
4834 for (
unsigned int comp = 0; comp < n_components; ++comp)
4836 laplacian_out[comp] = hess_diag[comp][0];
4837 for (
unsigned int d = 1;
d < dim; ++
d)
4838 laplacian_out[comp] += hess_diag[comp][d];
4840 return laplacian_out;
4850 typename VectorizedArrayType>
4855 VectorizedArrayType>::value_type
4861 Assert(this->hessians_quad_initialized ==
true,
4866 Assert(this->normal_x_jacobian !=
nullptr,
4868 "update_hessians"));
4872 const std::size_t nqp = this->n_quadrature_points;
4873 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4877 const auto nxj = this->normal_x_jacobian[0];
4879 for (
unsigned int comp = 0; comp < n_components; ++comp)
4881 for (
unsigned int d = 0;
d < dim; ++
d)
4882 hessian_out[comp] +=
4883 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4884 (nxj[
d]) * (nxj[d]);
4891 hessian_out[comp] +=
4892 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4896 hessian_out[comp] +=
4897 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4899 hessian_out[comp] +=
4900 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4902 hessian_out[comp] +=
4903 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4914 const auto normal = this->normal_vector(q_point);
4915 const auto hessian = get_hessian(q_point);
4917 if constexpr (n_components == 1)
4918 hessian_out[0] =
hessian * normal * normal;
4920 for (
unsigned int comp = 0; comp < n_components; ++comp)
4921 hessian_out[comp] =
hessian[comp] * normal * normal;
4923 if constexpr (n_components == 1)
4924 return hessian_out[0];
4935 typename VectorizedArrayType>
4942 this->dof_values_initialized =
true;
4944 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4946 for (
unsigned int comp = 0; comp < n_components; ++comp)
4947 if constexpr (n_components == 1)
4948 this->values_dofs[comp * dofs + dof] = val_in;
4950 this->values_dofs[comp * dofs + dof] = val_in[comp];
4959 typename VectorizedArrayType>
4962 submit_value(
const value_type val_in,
const unsigned int q_point)
4969 Assert(this->J_value !=
nullptr,
4974 this->values_quad_submitted =
true;
4977 const std::size_t nqp = this->n_quadrature_points;
4978 VectorizedArrayType *
values = this->values_quad + q_point;
4980 const VectorizedArrayType JxW =
4982 this->J_value[0] * this->quadrature_weights[q_point] :
4983 this->J_value[q_point];
4984 if constexpr (n_components == 1)
4985 values[0] = val_in * JxW;
4988 if (n_components == dim &&
4989 this->
data->element_type ==
4994 Assert(this->J_value !=
nullptr,
5000 this->values_quad_submitted =
true;
5003 VectorizedArrayType *
values = this->values_quad + q_point;
5004 const std::size_t nqp = this->n_quadrature_points;
5010 const VectorizedArrayType weight =
5011 this->quadrature_weights[q_point];
5013 for (
unsigned int comp = 0; comp < n_components; ++comp)
5014 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
5021 this->jacobian[q_point] :
5026 const VectorizedArrayType fac =
5028 this->quadrature_weights[q_point] :
5030 this->J_value[q_point] :
5031 this->J_value[0] * this->quadrature_weights[q_point]) *
5032 ((dim == 2 && this->get_face_no() < 2) ?
5041 for (
unsigned int comp = 0; comp < n_components; ++comp)
5043 values[comp * nqp] = val_in[0] * jac[0][comp];
5044 for (
unsigned int e = 1;
e < dim; ++
e)
5045 values[comp * nqp] += val_in[e] * jac[e][comp];
5046 values[comp * nqp] *= fac;
5051 for (
unsigned int comp = 0; comp < n_components; ++comp)
5052 values[comp * nqp] = val_in[comp] * JxW;
5062 typename VectorizedArrayType>
5063template <
int,
typename>
5067 const unsigned int q_point)
5069 static_assert(n_components == 1,
5070 "Do not try to modify the default template parameters used for"
5071 " selectively enabling this function via std::enable_if!");
5072 submit_value(val_in[0], q_point);
5081 typename VectorizedArrayType>
5084 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
5091 Assert(this->J_value !=
nullptr,
5093 "update_gradients"));
5094 Assert(this->jacobian !=
nullptr,
5096 "update_gradients"));
5099 this->gradients_quad_submitted =
true;
5102 if constexpr (dim > 1 && n_components == dim)
5104 if (this->
data->element_type ==
5114 Assert(this->J_value !=
nullptr,
5116 "update_gradients"));
5117 Assert(this->jacobian !=
nullptr,
5119 "update_gradients"));
5122 this->gradients_quad_submitted =
true;
5125 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5126 VectorizedArrayType *
values =
5127 this->values_from_gradients_quad + q_point;
5128 const std::size_t nqp = this->n_quadrature_points;
5129 const std::size_t nqp_d = nqp * dim;
5139 const VectorizedArrayType weight =
5140 this->quadrature_weights[q_point];
5141 for (
unsigned int d = 0;
d < dim; ++
d)
5142 for (
unsigned int comp = 0; comp < n_components; ++comp)
5143 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5145 (jac[comp][comp] * weight);
5157 const VectorizedArrayType fac =
5159 this->quadrature_weights[q_point] :
5160 this->J_value[0] * this->quadrature_weights[q_point] *
5161 ((dim == 2 && this->get_face_no() < 2) ?
5166 VectorizedArrayType tmp[dim][dim];
5167 for (
unsigned int d = 0;
d < dim; ++
d)
5168 for (
unsigned int e = 0;
e < dim; ++
e)
5170 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5171 for (
unsigned int f = 1; f < dim; ++f)
5172 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5174 for (
unsigned int comp = 0; comp < n_components; ++comp)
5175 for (
unsigned int d = 0;
d < dim; ++
d)
5177 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5178 for (
unsigned int f = 1; f < dim; ++f)
5179 res += jac[f][comp] * tmp[d][f];
5188 const auto jac_grad =
5189 this->jacobian_gradients_non_inverse[q_point];
5191 this->jacobian[q_point];
5195 const VectorizedArrayType fac =
5196 (!is_face) ? this->quadrature_weights[q_point] :
5197 this->J_value[q_point] *
5198 ((dim == 2 && this->get_face_no() < 2) ?
5207 VectorizedArrayType tmp3[dim], tmp4[dim];
5208 for (
unsigned int d = 0;
d < dim; ++
d)
5210 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5211 for (
unsigned int e = 1;
e < dim; ++
e)
5212 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5214 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5215 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5216 for (
unsigned int d = 0;
d < dim; ++
d)
5218 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5219 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5221 for (
unsigned int d = 0;
d < dim; ++
d)
5223 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5224 for (
unsigned int e = 1;
e < dim; ++
e)
5225 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5231 VectorizedArrayType tmp[dim][dim];
5234 for (
unsigned int d = 0;
d < dim; ++
d)
5235 for (
unsigned int e = 0;
e < dim; ++
e)
5237 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5238 for (
unsigned int f = 1; f < dim; ++f)
5239 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5242 for (
unsigned int d = 0;
d < dim; ++
d)
5243 for (
unsigned int e = 0;
e < dim; ++
e)
5245 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5246 for (
unsigned int f = 1; f < dim; ++f)
5247 res += t_jac[d][f] * tmp[e][f];
5254 VectorizedArrayType
value[dim];
5255 for (
unsigned int d = 0;
d < dim; ++
d)
5257 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5258 for (
unsigned int e = 1;
e < dim; ++
e)
5259 value[d] += tmp[d][e] * jac_grad[d][e];
5261 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5262 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5263 for (
unsigned int d = 0;
d < dim; ++
d)
5265 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5266 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5271 for (
unsigned int d = 0;
d < dim; ++
d)
5273 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5274 for (
unsigned int e = 1;
e < dim; ++
e)
5275 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5276 for (
unsigned int e = 0;
e < dim; ++
e)
5277 value[e] -= t_jac[e][d] * tmp2;
5280 for (
unsigned int d = 0;
d < dim; ++
d)
5281 values[d * nqp] =
value[d];
5287 const std::size_t nqp_d = this->n_quadrature_points * dim;
5288 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5292 const VectorizedArrayType JxW =
5293 this->J_value[0] * this->quadrature_weights[q_point];
5299 std::array<VectorizedArrayType, dim> jac;
5300 for (
unsigned int d = 0;
d < dim; ++
d)
5301 jac[d] = this->jacobian[0][d][d];
5303 for (
unsigned int d = 0;
d < dim; ++
d)
5305 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5306 if constexpr (n_components == 1)
5307 gradients[d] = grad_in[d] * factor;
5309 for (
unsigned int comp = 0; comp < n_components; ++comp)
5310 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5317 this->jacobian[q_point] :
5319 const VectorizedArrayType JxW =
5321 this->J_value[q_point] :
5322 this->J_value[0] * this->quadrature_weights[q_point];
5323 if constexpr (n_components == 1)
5324 for (
unsigned int d = 0;
d < dim; ++
d)
5326 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5327 for (
unsigned int e = 1;
e < dim; ++
e)
5328 new_val += (jac[e][d] * grad_in[e]);
5332 for (
unsigned int comp = 0; comp < n_components; ++comp)
5333 for (
unsigned int d = 0;
d < dim; ++
d)
5335 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5336 for (
unsigned int e = 1;
e < dim; ++
e)
5337 new_val += (jac[e][d] * grad_in[comp][e]);
5349 typename VectorizedArrayType>
5350template <
int,
typename>
5354 const unsigned int q_point)
5356 static_assert(n_components == 1 && dim == 1,
5357 "Do not try to modify the default template parameters used for"
5358 " selectively enabling this function via std::enable_if!");
5359 submit_gradient(grad_in[0], q_point);
5368 typename VectorizedArrayType>
5374 Assert(this->normal_x_jacobian !=
nullptr,
5376 "update_gradients"));
5379 this->gradients_quad_submitted =
true;
5382 const std::size_t nqp_d = this->n_quadrature_points * dim;
5383 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5387 const VectorizedArrayType JxW_jac = this->J_value[0] *
5388 this->quadrature_weights[q_point] *
5389 this->normal_x_jacobian[0][dim - 1];
5390 for (
unsigned int comp = 0; comp < n_components; ++comp)
5392 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5393 gradients[comp * nqp_d + d] = VectorizedArrayType();
5394 if constexpr (n_components == 1)
5395 gradients[dim - 1] = grad_in * JxW_jac;
5397 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5402 const unsigned int index =
5405 this->normal_x_jacobian[
index];
5406 const VectorizedArrayType JxW =
5408 this->J_value[
index] * this->quadrature_weights[q_point] :
5409 this->J_value[
index];
5410 for (
unsigned int comp = 0; comp < n_components; ++comp)
5411 for (
unsigned int d = 0;
d < dim; ++
d)
5412 if constexpr (n_components == 1)
5415 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5425 typename VectorizedArrayType>
5428 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5435 Assert(this->J_value !=
nullptr,
5437 "update_hessians"));
5438 Assert(this->jacobian !=
nullptr,
5440 "update_hessians"));
5443 this->hessians_quad_submitted =
true;
5447 const std::size_t nqp = this->n_quadrature_points;
5448 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5451 const VectorizedArrayType JxW =
5452 this->J_value[0] * this->quadrature_weights[q_point];
5455 for (
unsigned int d = 0;
d < dim; ++
d)
5457 const auto jac_d = this->jacobian[0][
d][
d];
5458 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5459 for (
unsigned int comp = 0; comp < n_components; ++comp)
5460 if constexpr (n_components == 1)
5461 this->hessians_quad[
d * nqp + q_point] =
5462 hessian_in[
d][
d] * factor;
5464 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5465 hessian_in[comp][d][d] * factor;
5469 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5470 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5472 const auto jac_d = this->jacobian[0][
d][
d];
5473 const auto jac_e = this->jacobian[0][
e][
e];
5474 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5475 for (
unsigned int comp = 0; comp < n_components; ++comp)
5476 if constexpr (n_components == 1)
5477 this->hessians_quad[off_dia * nqp + q_point] =
5478 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5480 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5481 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5488 const VectorizedArrayType JxW =
5489 this->J_value[0] * this->quadrature_weights[q_point];
5490 for (
unsigned int comp = 0; comp < n_components; ++comp)
5493 if constexpr (n_components == 1)
5494 hessian_c = hessian_in;
5496 hessian_c = hessian_in[comp];
5499 VectorizedArrayType tmp[dim][dim];
5500 for (
unsigned int i = 0; i < dim; ++i)
5501 for (
unsigned int j = 0; j < dim; ++j)
5503 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5504 for (
unsigned int k = 1; k < dim; ++k)
5505 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5509 VectorizedArrayType tmp2[dim][dim];
5510 for (
unsigned int i = 0; i < dim; ++i)
5511 for (
unsigned int j = 0; j < dim; ++j)
5513 tmp2[i][j] = jac[0][i] * tmp[0][j];
5514 for (
unsigned int k = 1; k < dim; ++k)
5515 tmp2[i][j] += jac[k][i] * tmp[k][j];
5519 for (
unsigned int d = 0;
d < dim; ++
d)
5520 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5524 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5525 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5526 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5527 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5533 const VectorizedArrayType JxW = this->J_value[q_point];
5534 const auto &jac_grad = this->jacobian_gradients[q_point];
5535 for (
unsigned int comp = 0; comp < n_components; ++comp)
5538 if constexpr (n_components == 1)
5539 hessian_c = hessian_in;
5541 hessian_c = hessian_in[comp];
5544 VectorizedArrayType tmp[dim][dim];
5545 for (
unsigned int i = 0; i < dim; ++i)
5546 for (
unsigned int j = 0; j < dim; ++j)
5548 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5549 for (
unsigned int k = 1; k < dim; ++k)
5550 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5554 VectorizedArrayType tmp2[dim][dim];
5555 for (
unsigned int i = 0; i < dim; ++i)
5556 for (
unsigned int j = 0; j < dim; ++j)
5558 tmp2[i][j] = jac[0][i] * tmp[0][j];
5559 for (
unsigned int k = 1; k < dim; ++k)
5560 tmp2[i][j] += jac[k][i] * tmp[k][j];
5564 for (
unsigned int d = 0;
d < dim; ++
d)
5565 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5569 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5570 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5571 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5572 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5575 for (
unsigned int d = 0;
d < dim; ++
d)
5577 VectorizedArrayType
sum = 0;
5578 for (
unsigned int e = 0;
e < dim; ++
e)
5579 sum += hessian_c[e][e] * jac_grad[e][d];
5580 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5581 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5583 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5584 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5597 typename VectorizedArrayType>
5601 const unsigned int q_point)
5608 Assert(this->J_value !=
nullptr,
5610 "update_hessians"));
5611 Assert(this->jacobian !=
nullptr,
5613 "update_hessians"));
5616 this->hessians_quad_submitted =
true;
5620 const std::size_t nqp = this->n_quadrature_points;
5621 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5624 const VectorizedArrayType JxW =
5625 this->J_value[0] * this->quadrature_weights[q_point];
5627 const auto nxj = this->normal_x_jacobian[0];
5630 for (
unsigned int d = 0;
d < dim; ++
d)
5632 const auto nxj_d = nxj[
d];
5633 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5634 for (
unsigned int comp = 0; comp < n_components; ++comp)
5635 if constexpr (n_components == 1)
5636 this->hessians_quad[
d * nqp + q_point] =
5637 normal_hessian_in * factor;
5639 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5640 normal_hessian_in[comp] * factor;
5644 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5645 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5647 const auto jac_d = nxj[
d];
5648 const auto jac_e = nxj[
e];
5649 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5650 for (
unsigned int comp = 0; comp < n_components; ++comp)
5651 if constexpr (n_components == 1)
5652 this->hessians_quad[off_dia * nqp + q_point] =
5653 2. * normal_hessian_in * factor;
5655 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5656 2. * normal_hessian_in[comp] * factor;
5661 const auto normal = this->normal_vector(q_point);
5662 const auto normal_projector =
outer_product(normal, normal);
5663 if constexpr (n_components == 1)
5664 submit_hessian(normal_hessian_in * normal_projector, q_point);
5668 for (
unsigned int comp = 0; comp < n_components; ++comp)
5669 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5670 submit_hessian(tmp, q_point);
5681 typename VectorizedArrayType>
5686 VectorizedArrayType>::value_type
5693 Assert(this->values_quad_submitted ==
true,
5698 const std::size_t nqp = this->n_quadrature_points;
5699 for (
unsigned int q = 0; q < nqp; ++q)
5700 for (
unsigned int comp = 0; comp < n_components; ++comp)
5701 return_value[comp] += this->values_quad[comp * nqp + q];
5702 if constexpr (n_components == 1)
5703 return return_value[0];
5705 return return_value;
5714 typename VectorizedArrayType>
5715template <
int,
typename>
5720 static_assert(n_components == dim,
5721 "Do not try to modify the default template parameters used for"
5722 " selectively enabling this function via std::enable_if!");
5726 Assert(this->gradients_quad_initialized ==
true,
5730 Assert(this->jacobian !=
nullptr,
5732 "update_gradients"));
5734 VectorizedArrayType divergence;
5735 const std::size_t nqp = this->n_quadrature_points;
5738 this->
data->element_type ==
5741 VectorizedArrayType inv_det =
5744 this->jacobian[0][0][0] *
5745 ((dim == 2) ? this->jacobian[0][1][1] :
5746 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5754 if (is_face && dim == 2 && this->get_face_no() < 2)
5758 divergence = this->gradients_quad[q_point * dim];
5759 for (
unsigned int d = 1;
d < dim; ++
d)
5760 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5761 divergence *= inv_det;
5770 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5771 for (
unsigned int d = 1;
d < dim; ++
d)
5772 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5773 this->jacobian[0][
d][
d];
5780 this->jacobian[q_point] :
5782 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5783 for (
unsigned int e = 1;
e < dim; ++
e)
5784 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5785 for (
unsigned int d = 1;
d < dim; ++
d)
5786 for (
unsigned int e = 0;
e < dim; ++
e)
5788 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5800 typename VectorizedArrayType>
5801template <
int,
typename>
5806 static_assert(n_components == dim,
5807 "Do not try to modify the default template parameters used for"
5808 " selectively enabling this function via std::enable_if!");
5811 const auto grad = get_gradient(q_point);
5812 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5813 VectorizedArrayType half = Number(0.5);
5814 for (
unsigned int d = 0;
d < dim; ++
d)
5815 symmetrized[d] = grad[d][d];
5821 symmetrized[2] = grad[0][1] + grad[1][0];
5822 symmetrized[2] *= half;
5825 symmetrized[3] = grad[0][1] + grad[1][0];
5826 symmetrized[3] *= half;
5827 symmetrized[4] = grad[0][2] + grad[2][0];
5828 symmetrized[4] *= half;
5829 symmetrized[5] = grad[1][2] + grad[2][1];
5830 symmetrized[5] *= half;
5844 typename VectorizedArrayType>
5845template <
int,
typename>
5847 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5849 get_curl(const unsigned
int q_point) const
5851 static_assert(dim > 1 && n_components == dim,
5852 "Do not try to modify the default template parameters used for"
5853 " selectively enabling this function via std::enable_if!");
5857 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5861 curl[0] = grad[1][0] - grad[0][1];
5864 curl[0] = grad[2][1] - grad[1][2];
5865 curl[1] = grad[0][2] - grad[2][0];
5866 curl[2] = grad[1][0] - grad[0][1];
5880 typename VectorizedArrayType>
5881template <
int,
typename>
5885 const unsigned int q_point)
5887 static_assert(n_components == dim,
5888 "Do not try to modify the default template parameters used for"
5889 " selectively enabling this function via std::enable_if!");
5896 Assert(this->J_value !=
nullptr,
5898 "update_gradients"));
5899 Assert(this->jacobian !=
nullptr,
5901 "update_gradients"));
5904 this->gradients_quad_submitted =
true;
5907 const std::size_t nqp_d = this->n_quadrature_points * dim;
5908 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5910 if (this->
data->element_type ==
5917 const VectorizedArrayType fac =
5919 this->quadrature_weights[q_point] * div_in :
5921 this->J_value[q_point] :
5922 this->J_value[0] * this->quadrature_weights[q_point]) *
5925 this->jacobian[this->cell_type >
5929 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5931 for (
unsigned int d = 0;
d < dim; ++
d)
5933 for (
unsigned int e = 0;
e < dim; ++
e)
5934 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5936 this->divergence_is_requested =
true;
5943 const VectorizedArrayType fac =
5944 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5945 for (
unsigned int d = 0;
d < dim; ++
d)
5947 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5948 for (
unsigned int e = 0;
e < dim; ++
e)
5949 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5956 this->jacobian[q_point] :
5958 const VectorizedArrayType fac =
5960 this->J_value[q_point] :
5961 this->J_value[0] * this->quadrature_weights[q_point]) *
5963 for (
unsigned int d = 0;
d < dim; ++
d)
5965 for (
unsigned int e = 0;
e < dim; ++
e)
5966 gradients[d * nqp_d + e] = jac[d][e] * fac;
5978 typename VectorizedArrayType>
5979template <
int,
typename>
5984 const unsigned int q_point)
5986 static_assert(n_components == dim,
5987 "Do not try to modify the default template parameters used for"
5988 " selectively enabling this function via std::enable_if!");
5991 this->
data->element_type !=
6003 Assert(this->J_value !=
nullptr,
6005 "update_gradients"));
6006 Assert(this->jacobian !=
nullptr,
6008 "update_gradients"));
6011 this->gradients_quad_submitted =
true;
6014 const std::size_t nqp_d = this->n_quadrature_points * dim;
6015 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
6018 const VectorizedArrayType JxW =
6019 this->J_value[0] * this->quadrature_weights[q_point];
6021 for (
unsigned int d = 0;
d < dim; ++
d)
6022 gradients[d * nqp_d + d] =
6024 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6025 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
6027 const VectorizedArrayType
value =
6036 const VectorizedArrayType JxW =
6038 this->J_value[q_point] :
6039 this->J_value[0] * this->quadrature_weights[q_point];
6042 this->jacobian[q_point] :
6044 VectorizedArrayType weighted[dim][dim];
6045 for (
unsigned int i = 0; i < dim; ++i)
6047 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6048 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6050 const VectorizedArrayType
value =
6052 weighted[i][j] =
value;
6053 weighted[j][i] =
value;
6055 for (
unsigned int comp = 0; comp < dim; ++comp)
6056 for (
unsigned int d = 0;
d < dim; ++
d)
6058 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6059 for (
unsigned int e = 1;
e < dim; ++
e)
6060 new_val += jac[e][d] * weighted[comp][e];
6072 typename VectorizedArrayType>
6073template <
int,
typename>
6077 const unsigned int q_point)
6079 static_assert(n_components == dim,
6080 "Do not try to modify the default template parameters used for"
6081 " selectively enabling this function via std::enable_if!");
6087 grad[1][0] = curl[0];
6088 grad[0][1] = -curl[0];
6091 grad[2][1] = curl[0];
6092 grad[1][2] = -curl[0];
6093 grad[0][2] = curl[1];
6094 grad[2][0] = -curl[1];
6095 grad[1][0] = curl[2];
6096 grad[0][1] = -curl[2];
6101 submit_gradient(grad, q_point);
6114 typename VectorizedArrayType>
6120 VectorizedArrayType>::
6122 const unsigned int fe_no,
6123 const unsigned int quad_no,
6124 const unsigned int first_selected_component,
6125 const unsigned int active_fe_index,
6126 const unsigned int active_quad_index)
6127 : BaseClass(matrix_free,
6129 first_selected_component,
6137 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6138 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6139 , n_q_points(this->
data->n_q_points)
6141 check_template_arguments(fe_no, 0);
6151 typename VectorizedArrayType>
6157 VectorizedArrayType>::
6159 const std::pair<unsigned int, unsigned int> &range,
6160 const unsigned int dof_no,
6161 const unsigned int quad_no,
6162 const unsigned int first_selected_component)
6166 first_selected_component,
6167 matrix_free.get_cell_active_fe_index(range, dof_no))
6177 typename VectorizedArrayType>
6183 VectorizedArrayType>::
6188 const unsigned int first_selected_component)
6189 : BaseClass(mapping,
6193 first_selected_component,
6195 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6196 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6197 , n_q_points(this->
data->n_q_points)
6209 typename VectorizedArrayType>
6215 VectorizedArrayType>::
6219 const unsigned int first_selected_component)
6224 first_selected_component,
6226 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6227 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6228 , n_q_points(this->
data->n_q_points)
6240 typename VectorizedArrayType>
6246 VectorizedArrayType>::
6249 const unsigned int first_selected_component)
6250 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6252 other.mapped_geometry->get_quadrature(),
6253 other.mapped_geometry->get_fe_values().get_update_flags(),
6254 first_selected_component,
6256 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6257 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6258 , n_q_points(this->
data->n_q_points)
6270 typename VectorizedArrayType>
6279 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6280 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6281 , n_q_points(this->
data->n_q_points)
6293 typename VectorizedArrayType>
6299 VectorizedArrayType> &
6305 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6307 BaseClass::operator=(other);
6319 typename VectorizedArrayType>
6326 VectorizedArrayType>::
6327 check_template_arguments(
const unsigned int dof_no,
6328 const unsigned int first_selected_component)
6331 (void)first_selected_component;
6334 this->
data->dofs_per_component_on_cell > 0,
6336 "There is nothing useful you can do with an FEEvaluation object with "
6337 "FE_Nothing, i.e., without DoFs! If you have passed to "
6338 "MatrixFree::reinit() a collection of finite elements also containing "
6339 "FE_Nothing, please check - before creating FEEvaluation - the category "
6340 "of the current range by calling either "
6341 "MatrixFree::get_cell_range_category(range) or "
6342 "MatrixFree::get_face_range_category(range). The returned category "
6343 "is the index of the active FE, which you can use to exclude "
6350 if ((
static_cast<unsigned int>(fe_degree) !=
6352 static_cast<unsigned int>(fe_degree) !=
6353 this->
data->data.front().fe_degree) ||
6354 n_q_points != this->n_quadrature_points)
6356 std::string message =
6357 "-------------------------------------------------------\n";
6359 "Illegal arguments in constructor/wrong template arguments!\n";
6360 message +=
" Called --> FEEvaluation<dim,";
6364 message +=
",Number>(data";
6380 if (
static_cast<unsigned int>(fe_degree) ==
6381 this->
data->data.front().fe_degree)
6383 proposed_dof_comp = dof_no;
6384 proposed_fe_comp = first_selected_component;
6387 for (
unsigned int no = 0;
6390 for (
unsigned int nf = 0;
6393 if (this->matrix_free
6394 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6396 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6398 proposed_dof_comp = no;
6399 proposed_fe_comp = nf;
6403 this->mapping_data->descriptor[this->active_quad_index]
6405 proposed_quad_comp = this->quad_no;
6407 for (
unsigned int no = 0;
6413 .descriptor[this->active_quad_index]
6414 .n_q_points == n_q_points)
6416 proposed_quad_comp = no;
6423 if (proposed_dof_comp != first_selected_component)
6424 message +=
"Wrong vector component selection:\n";
6426 message +=
"Wrong quadrature formula selection:\n";
6427 message +=
" Did you mean FEEvaluation<dim,";
6431 message +=
",Number>(data";
6441 std::string correct_pos;
6442 if (proposed_dof_comp != dof_no)
6443 correct_pos =
" ^ ";
6446 if (proposed_quad_comp != this->quad_no)
6447 correct_pos +=
" ^ ";
6450 if (proposed_fe_comp != first_selected_component)
6451 correct_pos +=
" ^\n";
6453 correct_pos +=
" \n";
6460 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6461 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6462 message +=
"Wrong template arguments:\n";
6463 message +=
" Did you mean FEEvaluation<dim,";
6468 message +=
",Number>(data";
6477 std::string correct_pos;
6478 if (this->
data->data.front().fe_degree !=
6479 static_cast<unsigned int>(fe_degree))
6483 if (proposed_n_q_points_1d != n_q_points_1d)
6484 correct_pos +=
" ^\n";
6486 correct_pos +=
" \n";
6487 message +=
" " + correct_pos;
6489 Assert(
static_cast<unsigned int>(fe_degree) ==
6490 this->
data->data.front().fe_degree &&
6491 n_q_points == this->n_quadrature_points,
6497 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6508 typename VectorizedArrayType>
6517 Assert(this->matrix_free !=
nullptr,
6518 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6519 " Integer indexing is not possible."));
6527 const unsigned int offsets =
6528 this->mapping_data->data_index_offsets[
cell_index];
6529 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6530 this->J_value = &this->mapping_data->JxW_values[offsets];
6531 if (!this->mapping_data->jacobian_gradients[0].empty())
6533 this->jacobian_gradients =
6534 this->mapping_data->jacobian_gradients[0].data() + offsets;
6535 this->jacobian_gradients_non_inverse =
6536 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6542 for (
unsigned int i = 0; i < n_lanes; ++i)
6543 this->cell_ids[i] =
cell_index * n_lanes + i;
6550 this->cell_ids[i] =
cell_index * n_lanes + i;
6551 for (; i < n_lanes; ++i)
6555 if (this->mapping_data->quadrature_points.empty() ==
false)
6557 &this->mapping_data->quadrature_points
6558 [this->mapping_data->quadrature_point_offsets[this->cell]];
6562 this->is_reinitialized =
true;
6563 this->dof_values_initialized =
false;
6564 this->values_quad_initialized =
false;
6565 this->gradients_quad_initialized =
false;
6566 this->hessians_quad_initialized =
false;
6577 typename VectorizedArrayType>
6584 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6591 this->cell_ids = cell_ids;
6596 for (
unsigned int v = 0; v < n_lanes; ++v)
6610 if (this->mapped_geometry ==
nullptr)
6611 this->mapped_geometry =
6612 std::make_shared<internal::MatrixFreeFunctions::
6613 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6615 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6617 auto &this_jacobian_data = mapping_storage.jacobians[0];
6618 auto &this_J_value_data = mapping_storage.JxW_values;
6619 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6620 auto &this_jacobian_gradients_non_inverse_data =
6621 mapping_storage.jacobian_gradients_non_inverse[0];
6622 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6626 if (this_jacobian_data.size() != 2)
6627 this_jacobian_data.resize_fast(2);
6629 if (this_J_value_data.size() != 1)
6630 this_J_value_data.resize_fast(1);
6632 const auto &update_flags_cells =
6636 this_jacobian_gradients_data.size() != 1)
6638 this_jacobian_gradients_data.resize_fast(1);
6639 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6643 this_quadrature_points_data.size() != 1)
6644 this_quadrature_points_data.resize_fast(1);
6648 if (this_jacobian_data.size() != this->n_quadrature_points)
6649 this_jacobian_data.resize_fast(this->n_quadrature_points);
6651 if (this_J_value_data.size() != this->n_quadrature_points)
6652 this_J_value_data.resize_fast(this->n_quadrature_points);
6654 const auto &update_flags_cells =
6658 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6660 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6661 this_jacobian_gradients_non_inverse_data.resize_fast(
6662 this->n_quadrature_points);
6666 this_quadrature_points_data.size() != this->n_quadrature_points)
6667 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6671 this->jacobian = this_jacobian_data.data();
6672 this->J_value = this_J_value_data.data();
6673 this->jacobian_gradients = this_jacobian_gradients_data.data();
6674 this->jacobian_gradients_non_inverse =
6675 this_jacobian_gradients_non_inverse_data.data();
6679 for (
unsigned int v = 0; v < n_lanes; ++v)
6686 const unsigned int cell_batch_index =
cell_index / n_lanes;
6687 const unsigned int offsets =
6688 this->mapping_data->data_index_offsets[cell_batch_index];
6689 const unsigned int lane =
cell_index % n_lanes;
6691 if (this->cell_type <=
6695 for (
unsigned int q = 0; q < 2; ++q)
6696 for (
unsigned int i = 0; i < dim; ++i)
6697 for (
unsigned int j = 0; j < dim; ++j)
6698 this_jacobian_data[q][i][j][v] =
6699 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6701 const unsigned int q = 0;
6703 this_J_value_data[q][v] =
6704 this->mapping_data->JxW_values[offsets + q][lane];
6706 const auto &update_flags_cells =
6711 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6712 for (
unsigned int j = 0; j < dim; ++j)
6713 this_jacobian_gradients_data[q][i][j][v] =
6715 ->jacobian_gradients[0][offsets + q][i][j][lane];
6717 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6718 for (
unsigned int j = 0; j < dim; ++j)
6719 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6721 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6726 for (
unsigned int i = 0; i < dim; ++i)
6727 this_quadrature_points_data[q][i][v] =
6728 this->mapping_data->quadrature_points
6730 ->quadrature_point_offsets[cell_batch_index] +
6736 const auto cell_type =
6740 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6742 const unsigned int q_src =
6748 this_J_value_data[q][v] =
6749 this->mapping_data->JxW_values[offsets + q_src][lane];
6751 for (
unsigned int i = 0; i < dim; ++i)
6752 for (
unsigned int j = 0; j < dim; ++j)
6753 this_jacobian_data[q][i][j][v] =
6755 ->jacobians[0][offsets + q_src][i][j][lane];
6757 const auto &update_flags_cells =
6762 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6763 for (
unsigned int j = 0; j < dim; ++j)
6764 this_jacobian_gradients_data[q][i][j][v] =
6766 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6768 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6769 for (
unsigned int j = 0; j < dim; ++j)
6770 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6772 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6785 this->mapping_data->quadrature_points
6787 ->quadrature_point_offsets[cell_batch_index] +
6791 this->mapping_data->jacobians[0][offsets + 1];
6793 for (
unsigned int d = 0;
d < dim; ++
d)
6796 static_cast<Number
>(
6797 this->descriptor->quadrature.point(q)[
d]);
6799 for (
unsigned int d = 0;
d < dim; ++
d)
6800 for (
unsigned int e = 0;
e < dim; ++
e)
6803 static_cast<Number
>(
6804 this->descriptor->quadrature.point(q)[
e]);
6806 for (
unsigned int i = 0; i < dim; ++i)
6807 this_quadrature_points_data[q][i][v] = point[i][lane];
6812 for (
unsigned int i = 0; i < dim; ++i)
6813 this_quadrature_points_data[q][i][v] =
6814 this->mapping_data->quadrature_points
6816 ->quadrature_point_offsets[cell_batch_index] +
6826 this->is_reinitialized =
true;
6827 this->dof_values_initialized =
false;
6828 this->values_quad_initialized =
false;
6829 this->gradients_quad_initialized =
false;
6830 this->hessians_quad_initialized =
false;
6841 typename VectorizedArrayType>
6842template <
bool level_dof_access>
6849 VectorizedArrayType>
::
6852 Assert(this->matrix_free ==
nullptr,
6853 ExcMessage(
"Cannot use initialization from cell iterator if "
6854 "initialized from MatrixFree object. Use variant for "
6855 "on the fly computation with arguments as for FEValues "
6858 this->mapped_geometry->reinit(
6860 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6861 if (level_dof_access)
6862 cell->get_mg_dof_indices(this->local_dof_indices);
6864 cell->get_dof_indices(this->local_dof_indices);
6868 this->is_reinitialized =
true;
6879 typename VectorizedArrayType>
6886 VectorizedArrayType>
::
6889 Assert(this->matrix_free == 0,
6890 ExcMessage(
"Cannot use initialization from cell iterator if "
6891 "initialized from MatrixFree object. Use variant for "
6892 "on the fly computation with arguments as for FEValues "
6895 this->mapped_geometry->reinit(cell);
6899 this->is_reinitialized =
true;
6910 typename VectorizedArrayType>
6917 VectorizedArrayType>::
6922 Assert(this->dof_values_initialized ==
true,
6925 evaluate(this->values_dofs, evaluation_flags);
6935 typename VectorizedArrayType>
6942 VectorizedArrayType>::
6943 evaluate(
const VectorizedArrayType *values_array,
6946 const bool hessians_on_general_cells =
6950 if (hessians_on_general_cells)
6953 if (this->
data->element_type ==
6959 if constexpr (fe_degree > -1)
6962 template run<fe_degree, n_q_points_1d>(n_components,
6963 evaluation_flag_actual,
6971 evaluation_flag_actual,
6972 const_cast<VectorizedArrayType *
>(values_array),
6978 this->values_quad_initialized =
6980 this->gradients_quad_initialized =
6982 this->hessians_quad_initialized =
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 &fe_eval, VectorType &vector)
7007 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7008 const auto &dof_info = fe_eval.get_dof_info();
7015 if (std::is_same_v<typename VectorType::value_type, Number> &&
7019 interleaved_contiguous &&
7020 reinterpret_cast<
std::size_t>(
7022 dof_info.dof_indices_contiguous
7023 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
7024 [cell * VectorizedArrayType::
size()]) %
7025 sizeof(VectorizedArrayType) ==
7028 return reinterpret_cast<VectorizedArrayType *
>(
7032 [cell * VectorizedArrayType::size()] +
7034 [fe_eval.get_active_fe_index()]
7035 [fe_eval.get_first_selected_component()] *
7036 VectorizedArrayType::size());
7045 template <
typename Number,
7046 typename VectorizedArrayType,
7048 typename EvaluatorType,
7049 std::enable_if_t<!internal::has_begin<VectorType> ||
7052 VectorizedArrayType *
7053 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7066 typename VectorizedArrayType>
7067template <
typename VectorType>
7074 VectorizedArrayType>::
7075 gather_evaluate(
const VectorType &input_vector,
7078 const VectorizedArrayType *src_ptr =
7079 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7080 *
this, input_vector);
7081 if (src_ptr !=
nullptr)
7082 evaluate(src_ptr, evaluation_flag);
7085 this->read_dof_values(input_vector);
7086 evaluate(this->begin_dof_values(), evaluation_flag);
7097 typename VectorizedArrayType>
7104 VectorizedArrayType>::
7107 integrate(integration_flag, this->values_dofs);
7111 this->dof_values_initialized =
true;
7122 typename VectorizedArrayType>
7129 VectorizedArrayType>::
7131 VectorizedArrayType *values_array,
7132 const bool sum_into_values_array)
7137 Assert(this->values_quad_submitted ==
true,
7140 Assert(this->gradients_quad_submitted ==
true,
7143 Assert(this->hessians_quad_submitted ==
true,
7146 Assert(this->matrix_free !=
nullptr ||
7147 this->mapped_geometry->is_initialized(),
7153 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7154 "EvaluationFlags::hessians are supported."));
7160 unsigned int size = n_components * dim * n_q_points;
7163 for (
unsigned int i = 0; i <
size; ++i)
7164 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7168 for (
unsigned int i = 0; i <
size; ++i)
7169 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7174 if (n_components == dim &&
7175 this->
data->element_type ==
7179 this->divergence_is_requested ==
false)
7181 unsigned int size = n_components * n_q_points;
7184 for (
unsigned int i = 0; i <
size; ++i)
7185 this->values_quad[i] += this->values_from_gradients_quad[i];
7189 for (
unsigned int i = 0; i <
size; ++i)
7190 this->values_quad[i] = this->values_from_gradients_quad[i];
7195 if constexpr (fe_degree > -1)
7198 template run<fe_degree, n_q_points_1d>(n_components,
7199 integration_flag_actual,
7202 sum_into_values_array);
7208 integration_flag_actual,
7211 sum_into_values_array);
7216 this->dof_values_initialized =
true;
7227 typename VectorizedArrayType>
7228template <
typename VectorType>
7235 VectorizedArrayType>::
7237 VectorType &destination)
7239 VectorizedArrayType *dst_ptr =
7240 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7241 *
this, destination);
7242 if (dst_ptr !=
nullptr)
7243 integrate(integration_flag, dst_ptr,
true);
7246 integrate(integration_flag, this->begin_dof_values());
7247 this->distribute_local_to_global(destination);
7258 typename VectorizedArrayType>
7265 VectorizedArrayType>::dof_indices()
const
7282 typename VectorizedArrayType>
7288 VectorizedArrayType>::
7291 const bool is_interior_face,
7292 const unsigned int dof_no,
7293 const unsigned int quad_no,
7294 const unsigned int first_selected_component,
7295 const unsigned int active_fe_index,
7296 const unsigned int active_quad_index,
7297 const unsigned int face_type)
7298 : BaseClass(matrix_free,
7300 first_selected_component,
7308 , dofs_per_component(this->
data->dofs_per_component_on_cell)
7309 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
7310 , n_q_points(this->n_quadrature_points)
7320 typename VectorizedArrayType>
7326 VectorizedArrayType>::
7329 const std::pair<unsigned int, unsigned int> &range,
7330 const bool is_interior_face,
7331 const unsigned int dof_no,
7332 const unsigned int quad_no,
7333 const unsigned int first_selected_component)
7338 first_selected_component,
7339 matrix_free.get_face_active_fe_index(range,
7342 matrix_free.get_face_info(range.
first).face_type)
7352 typename VectorizedArrayType>
7359 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7361 Assert(this->mapped_geometry ==
nullptr,
7362 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7363 " Integer indexing is not possible"));
7364 if (this->mapped_geometry !=
nullptr)
7367 this->cell = face_index;
7368 this->dof_access_index =
7369 this->is_interior_face() ?
7377 this->matrix_free->get_task_info().boundary_partition_data.back())
7378 Assert(this->is_interior_face(),
7380 "Boundary faces do not have a neighbor. When looping over "
7381 "boundary faces use FEFaceEvaluation with the parameter "
7382 "is_interior_face set to true. "));
7384 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7389 this->face_ids[i] = face_index * n_lanes + i;
7390 for (; i < n_lanes; ++i)
7393 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7394 const unsigned int offsets =
7395 this->mapping_data->data_index_offsets[face_index];
7396 this->J_value = &this->mapping_data->JxW_values[offsets];
7397 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7399 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7400 this->normal_x_jacobian =
7402 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7403 this->jacobian_gradients =
7404 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7406 this->jacobian_gradients_non_inverse =
7408 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7412 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7415 this->mapping_data->quadrature_point_offsets.size());
7417 this->mapping_data->quadrature_points.data() +
7418 this->mapping_data->quadrature_point_offsets[this->cell];
7423 this->is_reinitialized =
true;
7424 this->dof_values_initialized =
false;
7425 this->values_quad_initialized =
false;
7426 this->gradients_quad_initialized =
false;
7427 this->hessians_quad_initialized =
false;
7438 typename VectorizedArrayType>
7446 const unsigned int face_number)
7452 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7456 Assert(this->mapped_geometry ==
nullptr,
7457 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7458 " Integer indexing is not possible"));
7459 if (this->mapped_geometry !=
nullptr)
7464 .faces_by_cells_type[
cell_index][face_number];
7467 this->dof_access_index =
7470 if (this->is_interior_face() ==
false)
7475 for (
unsigned int i = 0; i < n_lanes; ++i)
7478 const unsigned int cell_this =
cell_index * n_lanes + i;
7480 unsigned int face_index =
7485 this->face_ids[i] = face_index;
7490 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7491 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7498 auto cell_m = faces.cells_interior[face_index % n_lanes];
7499 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7501 const bool face_identifies_as_interior = cell_m != cell_this;
7503 Assert(cell_m == cell_this || cell_p == cell_this,
7507 if (face_identifies_as_interior)
7509 this->cell_ids[i] = cell_m;
7510 this->face_numbers[i] = faces.interior_face_no;
7514 this->cell_ids[i] = cell_p;
7515 this->face_numbers[i] = faces.exterior_face_no;
7518 const bool orientation_interior_face = faces.face_orientation >= 8;
7519 auto face_orientation = faces.face_orientation % 8;
7520 if (face_identifies_as_interior != orientation_interior_face)
7523 ->reference_cell() ==
7524 ReferenceCells::get_hypercube<dim>(),
7528 .get_inverse_combined_orientation(face_orientation);
7530 this->face_orientations[i] = face_orientation;
7535 this->face_orientations[0] = 0;
7536 this->face_numbers[0] = face_number;
7541 for (
unsigned int i = 0; i < n_lanes; ++i)
7542 this->cell_ids[i] =
cell_index * n_lanes + i;
7550 this->cell_ids[i] =
cell_index * n_lanes + i;
7551 for (; i < n_lanes; ++i)
7554 for (
unsigned int i = 0; i < n_lanes; ++i)
7561 const unsigned int offsets =
7563 .face_data_by_cells[this->quad_no]
7568 .face_data_by_cells[this->quad_no]
7569 .JxW_values.size());
7571 .face_data_by_cells[this->quad_no]
7572 .JxW_values[offsets];
7574 .face_data_by_cells[this->quad_no]
7575 .normal_vectors[offsets];
7577 .face_data_by_cells[this->quad_no]
7578 .jacobians[!this->is_interior_face()][offsets];
7579 this->normal_x_jacobian =
7581 .face_data_by_cells[this->quad_no]
7582 .normals_times_jacobians[!this->is_interior_face()][offsets];
7583 this->jacobian_gradients =
7584 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7586 this->jacobian_gradients_non_inverse =
7588 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7593 .face_data_by_cells[this->quad_no]
7594 .quadrature_point_offsets.empty() ==
false)
7596 const unsigned int index =
7600 .face_data_by_cells[this->quad_no]
7601 .quadrature_point_offsets.size());
7603 .face_data_by_cells[this->quad_no]
7604 .quadrature_points.data() +
7606 .face_data_by_cells[this->quad_no]
7607 .quadrature_point_offsets[
index];
7612 this->is_reinitialized =
true;
7613 this->dof_values_initialized =
false;
7614 this->values_quad_initialized =
false;
7615 this->gradients_quad_initialized =
false;
7616 this->hessians_quad_initialized =
false;
7627 typename VectorizedArrayType>
7634 VectorizedArrayType>::
7642 evaluate(this->values_dofs, evaluation_flag);
7652 typename VectorizedArrayType>
7659 VectorizedArrayType>::
7660 evaluate(
const VectorizedArrayType *values_array,
7663 Assert((evaluation_flag &
7666 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7667 "and EvaluationFlags::hessians are supported."));
7669 const bool hessians_on_general_cells =
7673 if (hessians_on_general_cells)
7676 if (this->
data->element_type ==
7682 if constexpr (fe_degree > -1)
7684 template
run<fe_degree, n_q_points_1d>(n_components,
7685 evaluation_flag_actual,
7690 n_components, evaluation_flag_actual, values_array, *
this);
7694 this->values_quad_initialized =
7696 this->gradients_quad_initialized =
7698 this->hessians_quad_initialized =
7710 typename VectorizedArrayType>
7717 VectorizedArrayType>::
7725 project_to_face(this->values_dofs, evaluation_flag);
7735 typename VectorizedArrayType>
7742 VectorizedArrayType>::
7743 project_to_face(
const VectorizedArrayType *values_array,
7746 Assert((evaluation_flag &
7749 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7750 "and EvaluationFlags::hessians are supported."));
7752 const bool hessians_on_general_cells =
7756 if (hessians_on_general_cells)
7759 if (this->
data->element_type ==
7765 if constexpr (fe_degree > -1)
7768 VectorizedArrayType>::template run<fe_degree>(n_components,
7769 evaluation_flag_actual,
7775 evaluation_flag_actual,
7789 typename VectorizedArrayType>
7796 VectorizedArrayType>::
7799 Assert((evaluation_flag &
7802 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7803 "and EvaluationFlags::hessians are supported."));
7805 const bool hessians_on_general_cells =
7809 if (hessians_on_general_cells)
7812 if (this->
data->element_type ==
7818 if constexpr (fe_degree > -1)
7821 VectorizedArrayType>::template run<fe_degree>(n_components,
7822 evaluation_flag_actual,
7830 this->values_quad_initialized =
7832 this->gradients_quad_initialized =
7834 this->hessians_quad_initialized =
7846 typename VectorizedArrayType>
7853 VectorizedArrayType>::
7855 const bool sum_into_values)
7857 integrate(integration_flag, this->values_dofs, sum_into_values);
7861 this->dof_values_initialized =
true;
7872 typename VectorizedArrayType>
7879 VectorizedArrayType>::
7881 VectorizedArrayType *values_array,
7882 const bool sum_into_values)
7884 Assert((integration_flag &
7887 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7888 "and EvaluationFlags::hessians are supported."));
7894 unsigned int size = n_components * dim * n_q_points;
7897 for (
unsigned int i = 0; i <
size; ++i)
7898 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7902 for (
unsigned int i = 0; i <
size; ++i)
7903 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7908 if (this->
data->element_type ==
7912 this->divergence_is_requested ==
false)
7914 unsigned int size = n_components * n_q_points;
7917 for (
unsigned int i = 0; i <
size; ++i)
7918 this->values_quad[i] += this->values_from_gradients_quad[i];
7922 for (
unsigned int i = 0; i <
size; ++i)
7923 this->values_quad[i] = this->values_from_gradients_quad[i];
7928 if constexpr (fe_degree > -1)
7930 template
run<fe_degree, n_q_points_1d>(n_components,
7931 integration_flag_actual,
7938 integration_flag_actual,
7951 typename VectorizedArrayType>
7958 VectorizedArrayType>::
7961 Assert((integration_flag &
7964 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7965 "and EvaluationFlags::hessians are supported."));
7971 unsigned int size = n_components * dim * n_q_points;
7974 for (
unsigned int i = 0; i <
size; ++i)
7975 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7979 for (
unsigned int i = 0; i <
size; ++i)
7980 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7985 if (this->
data->element_type ==
7989 this->divergence_is_requested ==
false)
7991 unsigned int size = n_components * n_q_points;
7994 for (
unsigned int i = 0; i <
size; ++i)
7995 this->values_quad[i] += this->values_from_gradients_quad[i];
7999 for (
unsigned int i = 0; i <
size; ++i)
8000 this->values_quad[i] = this->values_from_gradients_quad[i];
8005 if constexpr (fe_degree > -1)
8008 VectorizedArrayType>::template run<fe_degree>(n_components,
8009 integration_flag_actual,
8025 typename VectorizedArrayType>
8032 VectorizedArrayType>::
8034 const bool sum_into_values)
8036 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
8040 this->dof_values_initialized =
true;
8051 typename VectorizedArrayType>
8058 VectorizedArrayType>::
8060 VectorizedArrayType *values_array,
8061 const bool sum_into_values)
8063 Assert((integration_flag &
8066 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8067 "and EvaluationFlags::hessians are supported."));
8074 if (this->
data->element_type ==
8078 this->divergence_is_requested ==
false)
8081 if constexpr (fe_degree > -1)
8084 VectorizedArrayType>::template run<fe_degree>(n_components,
8085 integration_flag_actual,
8092 integration_flag_actual,
8105 typename VectorizedArrayType>
8106template <
typename VectorType>
8113 VectorizedArrayType>::
8114 gather_evaluate(
const VectorType &input_vector,
8117 Assert((evaluation_flag &
8120 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8121 "and EvaluationFlags::hessians are supported."));
8123 const auto shared_vector_data = internal::get_shared_vector_data(
8125 this->dof_access_index ==
8127 this->active_fe_index,
8130 if (this->
data->data.front().fe_degree > 0 &&
8131 fast_evaluation_supported(this->
data->data.front().fe_degree,
8132 this->data->data.front().n_q_points_1d) &&
8135 typename VectorType::value_type,
8136 VectorizedArrayType>::
8137 supports(evaluation_flag,
8141 this->dof_info->index_storage_variants[this->dof_access_index]
8144 if constexpr (fe_degree > -1)
8148 typename VectorType::value_type,
8149 VectorizedArrayType>::template
run<fe_degree,
8153 internal::get_beginning<typename VectorType::value_type>(
8162 typename VectorType::value_type,
8163 VectorizedArrayType>::evaluate(n_components,
8165 internal::get_beginning<
8166 typename VectorType::value_type>(
8174 this->read_dof_values(input_vector);
8175 this->evaluate(evaluation_flag);
8181 this->gradients_quad_initialized =
8183 this->hessians_quad_initialized =
8195 typename VectorizedArrayType>
8196template <
typename VectorType>
8204 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8205 const bool integrate_gradients,
8206 VectorType &destination)
8213 integrate_scatter(flag, destination);
8223 typename VectorizedArrayType>
8224template <
typename VectorType>
8231 VectorizedArrayType>::
8233 VectorType &destination)
8235 Assert((this->dof_access_index ==
8237 this->is_interior_face() ==
false) ==
false,
8240 const auto shared_vector_data = internal::get_shared_vector_data(
8242 this->dof_access_index ==
8244 this->active_fe_index,
8247 if (this->
data->data.front().fe_degree > 0 &&
8248 fast_evaluation_supported(this->
data->data.front().fe_degree,
8249 this->data->data.front().n_q_points_1d) &&
8252 typename VectorType::value_type,
8253 VectorizedArrayType>::
8254 supports(integration_flag,
8258 this->dof_info->index_storage_variants[this->dof_access_index]
8261 if constexpr (fe_degree > -1)
8265 typename VectorType::value_type,
8266 VectorizedArrayType>::template
run<fe_degree,
8270 internal::get_beginning<typename VectorType::value_type>(
8279 typename VectorType::value_type,
8280 VectorizedArrayType>::integrate(n_components,
8282 internal::get_beginning<
8283 typename VectorType::value_type>(
8291 integrate(integration_flag);
8292 this->distribute_local_to_global(destination);
8303 typename VectorizedArrayType>
8310 VectorizedArrayType>::dof_indices()
const
8323 typename VectorizedArrayType>
8330 VectorizedArrayType>::
8331 fast_evaluation_supported(
const unsigned int given_degree,
8332 const unsigned int given_n_q_points_1d)
8334 return fe_degree == -1 ?
8347 typename VectorizedArrayType>
8354 VectorizedArrayType>::
8355 fast_evaluation_supported(
const unsigned int given_degree,
8356 const unsigned int given_n_q_points_1d)
8358 return fe_degree == -1 ?
8371 typename VectorizedArrayType>
8378 VectorizedArrayType>::at_boundary()
const
8380 Assert(this->dof_access_index !=
8384 if (this->is_interior_face() ==
false)
8389 this->matrix_free->n_boundary_face_batches()))
8402 typename VectorizedArrayType>
8409 VectorizedArrayType>::boundary_id()
const
8411 Assert(this->dof_access_index !=
8428 typename VectorizedArrayType>
8436 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8438 return this->
data->dofs_per_component_on_face;
8448 typename VectorizedArrayType>
8455 VectorizedArrayType>::get_dofs_projected_to_face()
8457 return this->
data->dofs_per_component_on_face * n_components_;
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > val_in, const unsigned int q_point)
void submit_normal_hessian(const value_type normal_hessian_in, const unsigned int q_point)
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
std::conditional_t< n_components_==1, Tensor< 2, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 3, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > > > hessian_type
std::conditional_t< n_components_==1, Tensor< 1, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 2, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > > > gradient_type
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
VectorizedArrayType get_divergence(const unsigned int q_point) const
hessian_type get_hessian(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
static constexpr unsigned int n_components
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
std::conditional_t< n_components_==1, VectorizedArrayType, Tensor< 1, n_components_, VectorizedArrayType > > value_type
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
static constexpr unsigned int n_lanes
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
value_type get_normal_hessian(const unsigned int q_point) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
const unsigned int quad_no
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const unsigned int active_fe_index
std::array< unsigned int, n_lanes > cell_ids
const unsigned int active_quad_index
bool is_interior_face() const
FEEvaluationData & operator=(const FEEvaluationData &other)
const unsigned int first_selected_component
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
void reinit(const std::array< unsigned int, n_lanes > &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
static constexpr unsigned int n_lanes
FEEvaluation(const FEEvaluation &other)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
void project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
unsigned int get_dofs_per_component_projected_to_face()
void project_to_face(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
static constexpr unsigned int static_dofs_per_component
static constexpr unsigned int n_lanes
static constexpr unsigned int static_n_q_points
unsigned int get_dofs_projected_to_face()
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
types::boundary_id boundary_id() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
void evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
DoFHandler< dim >::cell_iterator get_cell_iterator(const unsigned int cell_batch_index, const unsigned int lane_index, const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_jacobian_grads
Gradient of volume element.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
constexpr const ReferenceCell & get_hypercube()
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static constexpr unsigned int value
@ dof_access_face_exterior
@ dof_access_face_interior
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
std::vector< std::pair< unsigned int, unsigned int > > row_starts
std::vector< std::vector< unsigned int > > component_dof_indices_offset
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< unsigned int > dof_indices
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
std::vector< unsigned int > row_starts_plain_indices
std::vector< unsigned int > component_to_base_index
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
std::vector< unsigned int > plain_dof_indices
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
std::vector< unsigned int > dof_indices_interleaved
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
std::vector< QuadratureDescriptor > descriptor
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)