16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
92 typename VectorizedArrayType>
99 std::conditional_t<n_components_ == 1,
106 n_components_ == dim,
113 n_components_ == dim,
118 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
156 template <
typename VectorType>
159 const VectorType &src,
160 const unsigned int first_index = 0,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
191 template <
typename VectorType>
194 const VectorType &src,
195 const unsigned int first_index = 0,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
229 template <
typename VectorType>
233 const unsigned int first_index = 0,
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
274 template <
typename VectorType>
278 const unsigned int first_index = 0,
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
284 template <
typename VectorType>
288 const unsigned int first_index = 0,
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
365 typename = std::enable_if_t<n_components == n_components_local>>
368 const unsigned int q_point);
419 template <
int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
423 const unsigned int q_point);
442 const unsigned int q_point);
518 const unsigned int q_point);
527 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
546 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
549 const unsigned int q_point);
559 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
578 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
582 const unsigned int q_point);
592 template <
int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
610 template <
int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
614 const unsigned int q_point);
655 const unsigned int dof_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
663 const unsigned int face_type);
737 template <
typename VectorType,
typename VectorOperation>
741 const std::array<VectorType *, n_components_> &vectors,
744 n_components_> &vectors_sm,
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints =
true)
const;
755 template <
typename VectorType,
typename VectorOperation>
759 const std::array<VectorType *, n_components_> &vectors,
762 n_components_> &vectors_sm,
763 const std::bitset<n_lanes> &mask)
const;
772 template <
typename VectorType,
typename VectorOperation>
776 const std::array<VectorType *, n_components_> &vectors)
const;
1380 typename VectorizedArrayType>
1385 VectorizedArrayType>
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1530 const std::pair<unsigned int, unsigned int> &range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1642 template <
bool level_dof_access>
1664 const unsigned int given_n_q_points_1d);
1707 template <
typename VectorType>
1737 VectorizedArrayType *values_array,
1738 const bool sum_into_values =
false);
1753 template <
typename VectorType>
1756 VectorType &output_vector);
1840 int n_q_points_1d = fe_degree + 1,
1841 int n_components_ = 1,
1842 typename Number = double,
1848 VectorizedArrayType>
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
2012 const std::pair<unsigned int, unsigned int> &range,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2029 reinit(
const unsigned int face_batch_number);
2039 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2046 const unsigned int given_n_q_points_1d);
2110 template <
typename VectorType>
2126 const bool sum_into_values =
false);
2139 VectorizedArrayType *values_array,
2140 const bool sum_into_values =
false);
2157 const bool sum_into_values =
false);
2165 VectorizedArrayType *values_array,
2166 const bool sum_into_values =
false);
2179 template <
typename VectorType>
2182 VectorType &output_vector);
2187 template <
typename VectorType>
2190 const bool integrate_gradients,
2191 VectorType &output_vector);
2267 namespace MatrixFreeFunctions
2271 template <
int dim,
int degree>
2280 template <
int degree>
2283 static constexpr unsigned int value = degree + 1;
2299 template <
bool is_face,
2302 typename VectorizedArrayType>
2305 extract_initialization_data(
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2312 const unsigned int active_fe_index_given,
2313 const unsigned int active_quad_index_given,
2314 const unsigned int face_type)
2317 InitializationData init_data;
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2330 active_fe_index_given :
2336 active_quad_index_given :
2337 std::min<unsigned int>(
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2369 typename VectorizedArrayType>
2374 VectorizedArrayType>::
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2387 internal::extract_initialization_data<is_face>(matrix_free,
2389 first_selected_component,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2402 this->set_data_pointers(scratch_data_array, n_components_);
2404 this->dof_info->start_components.back() == 1 ||
2405 static_cast<int>(n_components_) <=
2407 this->dof_info->start_components
2408 [this->dof_info->component_to_base_index[first_selected_component] +
2410 first_selected_component,
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2417 this->dof_info->start_components
2418 [this->dof_info->component_to_base_index[first_selected_component] +
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) +
")"));
2438 typename VectorizedArrayType>
2443 VectorizedArrayType>::
2449 const unsigned int first_selected_component,
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2461 first_selected_component)
2462 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2463 , matrix_free(nullptr)
2465 const unsigned int base_element_number =
2469 first_selected_component >=
2471 ExcMessage(
"The underlying element must at least contain as many "
2472 "components as requested by this class"));
2473 (void)base_element_number;
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2479 fe.component_to_base_index(first_selected_component).
first);
2481 this->set_data_pointers(scratch_data_array, n_components_);
2490 typename VectorizedArrayType>
2495 VectorizedArrayType>::
2500 VectorizedArrayType> &other)
2502 , scratch_data_array(other.matrix_free == nullptr ?
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(other.matrix_free)
2507 if (other.matrix_free ==
nullptr)
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2520 this->mapping_data = &this->mapped_geometry->get_data_storage();
2524 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2526 this->mapped_geometry->get_data_storage().JxW_values.begin();
2527 this->jacobian_gradients =
2528 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2529 this->jacobian_gradients_non_inverse =
2530 this->mapped_geometry->get_data_storage()
2531 .jacobian_gradients_non_inverse[0]
2534 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2537 this->set_data_pointers(scratch_data_array, n_components_);
2546 typename VectorizedArrayType>
2551 VectorizedArrayType> &
2557 VectorizedArrayType> &other)
2560 if (matrix_free ==
nullptr)
2563 delete scratch_data_array;
2572 matrix_free = other.matrix_free;
2574 if (other.matrix_free ==
nullptr)
2582 this->mapped_geometry =
2583 std::make_shared<internal::MatrixFreeFunctions::
2584 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2589 this->mapping_data = &this->mapped_geometry->get_data_storage();
2591 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2593 this->mapped_geometry->get_data_storage().JxW_values.begin();
2594 this->jacobian_gradients =
2595 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2596 this->jacobian_gradients_non_inverse =
2597 this->mapped_geometry->get_data_storage()
2598 .jacobian_gradients_non_inverse[0]
2601 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2608 this->set_data_pointers(scratch_data_array, n_components_);
2619 typename VectorizedArrayType>
2624 VectorizedArrayType>::~FEEvaluationBase()
2626 if (matrix_free !=
nullptr)
2637 delete scratch_data_array;
2648 typename VectorizedArrayType>
2653 Assert(matrix_free !=
nullptr,
2655 "FEEvaluation was not initialized with a MatrixFree object!"));
2656 return *matrix_free;
2665 template <
typename VectorType,
bool>
2666 struct ConstBlockVectorSelector;
2668 template <
typename VectorType>
2669 struct ConstBlockVectorSelector<
VectorType, true>
2671 using BaseVectorType =
const typename VectorType::BlockType;
2674 template <
typename VectorType>
2675 struct ConstBlockVectorSelector<
VectorType, false>
2677 using BaseVectorType =
typename VectorType::BlockType;
2683 template <
typename VectorType,
bool>
2684 struct BlockVectorSelector;
2686 template <
typename VectorType>
2689 using BaseVectorType =
typename ConstBlockVectorSelector<
2691 std::is_const_v<VectorType>>::BaseVectorType;
2693 static BaseVectorType *
2694 get_vector_component(VectorType &vec,
const unsigned int component)
2697 return &vec.block(component);
2701 template <
typename VectorType>
2702 struct BlockVectorSelector<
VectorType, false>
2706 static BaseVectorType *
2707 get_vector_component(VectorType &vec,
const unsigned int component)
2726 template <
typename VectorType>
2727 struct BlockVectorSelector<
std::vector<VectorType>, false>
2731 static BaseVectorType *
2732 get_vector_component(std::vector<VectorType> &vec,
2733 const unsigned int component)
2736 return &vec[component];
2740 template <
typename VectorType>
2741 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2745 static const BaseVectorType *
2746 get_vector_component(
const std::vector<VectorType> &vec,
2747 const unsigned int component)
2750 return &vec[component];
2754 template <
typename VectorType>
2755 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2759 static BaseVectorType *
2760 get_vector_component(std::vector<VectorType *> &vec,
2761 const unsigned int component)
2764 return vec[component];
2768 template <
typename VectorType>
2769 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2773 static const BaseVectorType *
2774 get_vector_component(
const std::vector<VectorType *> &vec,
2775 const unsigned int component)
2778 return vec[component];
2782 template <
typename VectorType, std::
size_t N>
2783 struct BlockVectorSelector<
std::array<VectorType *, N>, false>
2787 static BaseVectorType *
2788 get_vector_component(std::array<VectorType *, N> &vec,
2789 const unsigned int component)
2792 return vec[component];
2803 typename VectorizedArrayType>
2804template <
typename VectorType,
typename VectorOperation>
2809 const std::array<VectorType *, n_components_> &src,
2812 n_components_> &src_sm,
2813 const std::bitset<n_lanes> &
mask,
2814 const bool apply_constraints)
const
2819 if (this->matrix_free ==
nullptr)
2821 read_write_operation_global(operation, src);
2828 if (this->n_fe_components == 1)
2829 for (
unsigned int comp = 0; comp < n_components; ++comp)
2831 Assert(src[comp] !=
nullptr,
2832 ExcMessage(
"The finite element underlying this FEEvaluation "
2833 "object is scalar, but you requested " +
2834 std::to_string(n_components) +
2835 " components via the template argument in "
2836 "FEEvaluation. In that case, you must pass an "
2837 "std::vector<VectorType> or a BlockVector to " +
2838 "read_dof_values and distribute_local_to_global."));
2850 const bool accesses_exterior_dofs =
2851 this->dof_access_index ==
2853 this->is_interior_face() ==
false;
2863 bool is_contiguous =
true;
2865 if (accesses_exterior_dofs)
2867 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2868 const unsigned int n_filled_lanes =
2873 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2874 if (
mask[v] ==
true &&
2877 [cells[v] / n_lanes] <
2880 is_contiguous = false;
2884 this->dof_access_index :
2885 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2886 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2889 is_contiguous =
false;
2894 read_write_operation_contiguous(operation, src, src_sm,
mask);
2901 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2903 const bool masking_is_active =
mask.count() < n_lanes;
2904 if (masking_is_active)
2905 for (
unsigned int v = 0; v < n_lanes; ++v)
2906 if (
mask[v] ==
false)
2909 bool has_hn_constraints =
false;
2911 if (is_face ==
false)
2917 [this->first_selected_component])
2918 for (
unsigned int v = 0; v < n_lanes; ++v)
2923 has_hn_constraints = true;
2926 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2929 const bool use_vectorized_path =
2930 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2932 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2933 std::array<VectorizedArrayType *, n_components> values_dofs;
2934 for (
unsigned int c = 0; c < n_components; ++c)
2935 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2936 c * dofs_per_component;
2940 [is_face ? this->dof_access_index :
2941 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2942 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2943 IndexStorageVariants::interleaved &&
2944 use_vectorized_path)
2946 const unsigned int *dof_indices =
2948 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2952 [this->first_selected_component] *
2955 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2956 if (n_components == 1 || this->n_fe_components == 1)
2957 for (
unsigned int comp = 0; comp < n_components; ++comp)
2959 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2962 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2964 if (n_components == 1 || this->n_fe_components == 1)
2965 for (
unsigned int i = 0; i < dofs_per_component;
2966 ++i, dof_indices += n_lanes)
2967 for (
unsigned int comp = 0; comp < n_components; ++comp)
2968 operation.process_dof_gather(dof_indices,
2972 values_dofs[comp][i],
2975 for (
unsigned int comp = 0; comp < n_components; ++comp)
2976 for (
unsigned int i = 0; i < dofs_per_component;
2977 ++i, dof_indices += n_lanes)
2978 operation.process_dof_gather(dof_indices,
2982 values_dofs[comp][i],
2989 std::array<const unsigned int *, n_lanes> dof_indices;
2990 dof_indices.fill(
nullptr);
2995 bool has_constraints =
false;
2996 const unsigned int n_components_read =
2997 this->n_fe_components > 1 ? n_components : 1;
3001 for (
unsigned int v = 0; v < n_lanes; ++v)
3007 const std::pair<unsigned int, unsigned int> *my_index_start =
3008 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3009 this->first_selected_component];
3014 if (my_index_start[n_components_read].
second !=
3015 my_index_start[0].
second)
3016 has_constraints =
true;
3019 dof_info.
dof_indices.data() + my_index_start[0].first;
3024 for (
unsigned int v = 0; v < n_lanes; ++v)
3029 const std::pair<unsigned int, unsigned int> *my_index_start =
3030 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3031 this->first_selected_component];
3032 if (my_index_start[n_components_read].
second !=
3033 my_index_start[0].
second)
3034 has_constraints =
true;
3041 dof_info.hanging_node_constraint_masks_comp
3042 [this->active_fe_index][this->first_selected_component])
3043 has_hn_constraints = true;
3046 my_index_start[0].
first ||
3049 my_index_start[0].
first,
3052 dof_info.
dof_indices.data() + my_index_start[0].first;
3056 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3057 return i != numbers::invalid_unsigned_int;
3059 for (
unsigned int comp = 0; comp < n_components; ++comp)
3060 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3061 operation.process_empty(values_dofs[comp][i]);
3065 if (!has_constraints && apply_constraints)
3067 if (n_components == 1 || this->n_fe_components == 1)
3069 for (
unsigned int v = 0; v < n_lanes; ++v)
3074 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3075 for (
unsigned int comp = 0; comp < n_components; ++comp)
3076 operation.process_dof(dof_indices[v][i],
3078 values_dofs[comp][i][v]);
3083 for (
unsigned int comp = 0; comp < n_components; ++comp)
3084 for (
unsigned int v = 0; v < n_lanes; ++v)
3089 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3090 operation.process_dof(
3091 dof_indices[v][comp * dofs_per_component + i],
3093 values_dofs[comp][i][v]);
3105 for (
unsigned int v = 0; v < n_lanes; ++v)
3111 const unsigned int cell_dof_index =
3112 cell_index * this->n_fe_components + this->first_selected_component;
3113 const unsigned int n_components_read =
3114 this->n_fe_components > 1 ? n_components : 1;
3115 unsigned int index_indicators =
3117 unsigned int next_index_indicators =
3118 dof_info.
row_starts[cell_dof_index + 1].second;
3122 if (apply_constraints ==
false &&
3123 (dof_info.
row_starts[cell_dof_index].second !=
3124 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3130 dof_info.hanging_node_constraint_masks_comp
3131 [this->active_fe_index][this->first_selected_component])))
3140 [this->first_selected_component] +
3142 next_index_indicators = index_indicators;
3145 if (n_components == 1 || this->n_fe_components == 1)
3147 unsigned int ind_local = 0;
3148 for (; index_indicators != next_index_indicators; ++index_indicators)
3150 const std::pair<unsigned short, unsigned short> indicator =
3153 for (
unsigned int j = 0; j < indicator.first; ++j)
3154 for (
unsigned int comp = 0; comp < n_components; ++comp)
3155 operation.process_dof(dof_indices[v][j],
3157 values_dofs[comp][ind_local + j][v]);
3159 ind_local += indicator.first;
3160 dof_indices[v] += indicator.first;
3164 Number
value[n_components];
3165 for (
unsigned int comp = 0; comp < n_components; ++comp)
3166 operation.pre_constraints(values_dofs[comp][ind_local][v],
3169 const Number *data_val =
3171 const Number *end_pool =
3173 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3174 for (
unsigned int comp = 0; comp < n_components; ++comp)
3175 operation.process_constraint(*dof_indices[v],
3180 for (
unsigned int comp = 0; comp < n_components; ++comp)
3181 operation.post_constraints(
value[comp],
3182 values_dofs[comp][ind_local][v]);
3188 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3189 for (
unsigned int comp = 0; comp < n_components; ++comp)
3190 operation.process_dof(*dof_indices[v],
3192 values_dofs[comp][ind_local][v]);
3201 for (
unsigned int comp = 0; comp < n_components; ++comp)
3203 unsigned int ind_local = 0;
3206 for (; index_indicators != next_index_indicators;
3209 const std::pair<unsigned short, unsigned short> indicator =
3213 for (
unsigned int j = 0; j < indicator.first; ++j)
3214 operation.process_dof(dof_indices[v][j],
3216 values_dofs[comp][ind_local + j][v]);
3217 ind_local += indicator.first;
3218 dof_indices[v] += indicator.first;
3223 operation.pre_constraints(values_dofs[comp][ind_local][v],
3226 const Number *data_val =
3228 const Number *end_pool =
3231 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3232 operation.process_constraint(*dof_indices[v],
3237 operation.post_constraints(
value,
3238 values_dofs[comp][ind_local][v]);
3245 for (; ind_local < dofs_per_component;
3246 ++dof_indices[v], ++ind_local)
3249 operation.process_dof(*dof_indices[v],
3251 values_dofs[comp][ind_local][v]);
3254 if (apply_constraints ==
true && comp + 1 < n_components)
3255 next_index_indicators =
3256 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3268 typename VectorizedArrayType>
3269template <
typename VectorType,
typename VectorOperation>
3274 const std::array<VectorType *, n_components_> &src)
const
3278 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3279 unsigned int index = this->first_selected_component * dofs_per_component;
3280 for (
unsigned int comp = 0; comp < n_components; ++comp)
3282 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3284 operation.process_empty(
3285 this->values_dofs[comp * dofs_per_component + i]);
3286 operation.process_dof_global(
3287 local_dof_indices[this->data->lexicographic_numbering[
index]],
3289 this->values_dofs[comp * dofs_per_component + i][0]);
3300 typename VectorizedArrayType>
3301template <
typename VectorType,
typename VectorOperation>
3306 const std::array<VectorType *, n_components_> &src,
3309 n_components_> &vectors_sm,
3310 const std::bitset<n_lanes> &
mask)
const
3319 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3322 is_face ? this->dof_access_index :
3324 const unsigned int n_active_lanes =
mask.count();
3327 const std::vector<unsigned int> &dof_indices_cont =
3330 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3331 std::array<VectorizedArrayType *, n_components> values_dofs;
3332 for (
unsigned int c = 0; c < n_components; ++c)
3333 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3334 c * dofs_per_component;
3338 const bool accesses_exterior_dofs =
3339 this->dof_access_index ==
3341 this->is_interior_face() ==
false;
3347 interleaved_contiguous &&
3348 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3350 const unsigned int dof_index =
3351 dof_indices_cont[this->cell * n_lanes] +
3354 [this->first_selected_component] *
3356 if (n_components == 1 || this->n_fe_components == 1)
3357 for (
unsigned int comp = 0; comp < n_components; ++comp)
3358 operation.process_dofs_vectorized(dofs_per_component,
3364 operation.process_dofs_vectorized(dofs_per_component * n_components,
3372 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3376 const unsigned int n_filled_lanes =
3379 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3380 n_active_lanes == n_lanes &&
3381 !accesses_exterior_dofs;
3383 if (vectors_sm[0] !=
nullptr)
3385 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3386 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3388 const auto upper_bound =
3389 std::min<unsigned int>(n_filled_lanes, n_lanes);
3390 for (
unsigned int v = 0; v < upper_bound; ++v)
3392 if (
mask[v] ==
false)
3394 vector_ptrs[v] =
nullptr;
3414 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3415 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3417 [this->active_fe_index][this->first_selected_component]);
3419 vector_ptrs[v] =
nullptr;
3421 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3422 vector_ptrs[v] =
nullptr;
3427 if (use_vectorized_path)
3429 if (n_components == 1 || this->n_fe_components == 1)
3431 for (
unsigned int comp = 0; comp < n_components; ++comp)
3433 auto vector_ptrs = compute_vector_ptrs(comp);
3434 operation.process_dofs_vectorized_transpose(
3443 auto vector_ptrs = compute_vector_ptrs(0);
3444 operation.process_dofs_vectorized_transpose(dofs_per_component *
3452 for (
unsigned int comp = 0; comp < n_components; ++comp)
3454 auto vector_ptrs = compute_vector_ptrs(
3455 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3457 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3458 operation.process_empty(values_dofs[comp][i]);
3460 if (n_components == 1 || this->n_fe_components == 1)
3462 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3463 if (
mask[v] ==
true)
3464 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3465 operation.process_dof(vector_ptrs[v][i],
3466 values_dofs[comp][i][v]);
3470 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3471 if (
mask[v] ==
true)
3472 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3473 operation.process_dof(
3474 vector_ptrs[v][i + comp * dofs_per_component],
3475 values_dofs[comp][i][v]);
3481 std::array<unsigned int, n_lanes> dof_indices;
3482 std::fill(dof_indices.begin(),
3487 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3491 if (
mask[v] ==
true)
3493 dof_indices_cont[cells[v]] +
3496 [this->first_selected_component] *
3502 if (use_vectorized_path)
3508 if (n_components == 1 || this->n_fe_components == 1)
3509 for (
unsigned int comp = 0; comp < n_components; ++comp)
3510 operation.process_dofs_vectorized_transpose(dofs_per_component,
3516 operation.process_dofs_vectorized_transpose(dofs_per_component *
3525 interleaved_contiguous_strided)
3527 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3528 if (n_components == 1 || this->n_fe_components == 1)
3529 for (
unsigned int comp = 0; comp < n_components; ++comp)
3530 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3531 src[comp]->
begin());
3534 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3536 if (n_components == 1 || this->n_fe_components == 1)
3537 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3539 for (
unsigned int comp = 0; comp < n_components; ++comp)
3540 operation.process_dof_gather(dof_indices.data(),
3543 src_ptrs[comp] + i * n_lanes,
3544 values_dofs[comp][i],
3548 for (
unsigned int comp = 0; comp < n_components; ++comp)
3549 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3551 operation.process_dof_gather(
3554 (comp * dofs_per_component + i) * n_lanes,
3555 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3556 values_dofs[comp][i],
3564 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3566 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3567 if (n_components == 1 || this->n_fe_components == 1)
3568 for (
unsigned int comp = 0; comp < n_components; ++comp)
3569 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3570 src[comp]->
begin());
3573 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3575 const unsigned int *offsets =
3577 if (n_components == 1 || this->n_fe_components == 1)
3578 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3580 for (
unsigned int comp = 0; comp < n_components; ++comp)
3581 operation.process_dof_gather(dof_indices.data(),
3585 values_dofs[comp][i],
3588 for (
unsigned int v = 0; v < n_lanes; ++v)
3589 dof_indices[v] += offsets[v];
3592 for (
unsigned int comp = 0; comp < n_components; ++comp)
3593 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3595 operation.process_dof_gather(dof_indices.data(),
3599 values_dofs[comp][i],
3602 for (
unsigned int v = 0; v < n_lanes; ++v)
3603 dof_indices[v] += offsets[v];
3608 for (
unsigned int comp = 0; comp < n_components; ++comp)
3610 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3611 operation.process_empty(values_dofs[comp][i]);
3612 if (accesses_exterior_dofs)
3614 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3615 if (
mask[v] ==
true)
3618 [ind][cells[v] / VectorizedArrayType::size()] ==
3622 if (n_components == 1 || this->n_fe_components == 1)
3624 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3625 operation.process_dof(dof_indices[v] + i,
3627 values_dofs[comp][i][v]);
3631 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3632 operation.process_dof(dof_indices[v] + i +
3633 comp * dofs_per_component,
3635 values_dofs[comp][i][v]);
3640 const unsigned int offset =
3643 if (n_components == 1 || this->n_fe_components == 1)
3645 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3646 operation.process_dof(dof_indices[v] + i * offset,
3648 values_dofs[comp][i][v]);
3652 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3653 operation.process_dof(
3655 (i + comp * dofs_per_component) * offset,
3657 values_dofs[comp][i][v]);
3668 if (n_components == 1 || this->n_fe_components == 1)
3670 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3671 if (
mask[v] ==
true)
3672 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3673 operation.process_dof(dof_indices[v] + i,
3675 values_dofs[comp][i][v]);
3679 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3680 if (
mask[v] ==
true)
3681 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3682 operation.process_dof(dof_indices[v] + i +
3683 comp * dofs_per_component,
3685 values_dofs[comp][i][v]);
3690 const unsigned int *offsets =
3692 [ind][VectorizedArrayType::size() * this->cell];
3693 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3695 if (n_components == 1 || this->n_fe_components == 1)
3696 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3698 if (
mask[v] ==
true)
3699 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3700 operation.process_dof(dof_indices[v] + i * offsets[v],
3702 values_dofs[comp][i][v]);
3706 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3707 if (
mask[v] ==
true)
3708 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3709 operation.process_dof(
3711 (i + comp * dofs_per_component) * offsets[v],
3713 values_dofs[comp][i][v]);
3725 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3726 decltype(std::declval<VectorType>().begin())
3727 get_beginning(VectorType &vec)
3735 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3736 typename VectorType::value_type *
3737 get_beginning(VectorType &)
3743 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3745 const std::vector<ArrayView<const typename VectorType::value_type>> *
3746 get_shared_vector_data(VectorType *vec,
3747 const bool is_valid_mode_for_sm,
3748 const unsigned int active_fe_index,
3752 if (is_valid_mode_for_sm &&
3755 active_fe_index == 0)
3756 return &vec->shared_vector_data();
3762 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3764 const std::vector<ArrayView<const typename VectorType::value_type>> *
3765 get_shared_vector_data(VectorType *,
3773 template <
int n_components,
typename VectorType>
3775 std::array<
typename internal::BlockVectorSelector<
3780 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3784 get_vector_data(VectorType &src,
3785 const unsigned int first_index,
3786 const bool is_valid_mode_for_sm,
3787 const unsigned int active_fe_index,
3793 std::array<
typename internal::BlockVectorSelector<
3799 ArrayView<
const typename internal::BlockVectorSelector<
3805 for (
unsigned int d = 0;
d < n_components; ++
d)
3806 src_data.first[d] = internal::BlockVectorSelector<
3812 for (
unsigned int d = 0;
d < n_components; ++
d)
3813 src_data.second[d] = get_shared_vector_data(
3814 const_cast<typename internal::BlockVectorSelector<
3815 std::remove_const_t<VectorType>,
3817 *>(src_data.first[d]),
3818 is_valid_mode_for_sm,
3832 typename VectorizedArrayType>
3837 if (this->dof_info ==
nullptr ||
3839 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3840 this->dof_info->hanging_node_constraint_masks_comp
3841 [this->active_fe_index][this->first_selected_component] ==
false)
3844 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3847 bool hn_available =
false;
3849 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3851 for (
unsigned int v = 0; v < n_lanes; ++v)
3863 constraint_mask[v] =
mask;
3869 if (hn_available ==
false)
3874 this->data->data.front().fe_degree,
3875 this->get_shape_info(),
3887 typename VectorizedArrayType>
3888template <
typename VectorType>
3892 const unsigned int first_index,
3893 const std::bitset<n_lanes> &
mask)
3895 const auto src_data = internal::get_vector_data<n_components_>(
3898 this->dof_access_index ==
3900 this->active_fe_index,
3904 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3906 apply_hanging_node_constraints(
false);
3909 this->dof_values_initialized =
true;
3919 typename VectorizedArrayType>
3920template <
typename VectorType>
3924 const unsigned int first_index,
3925 const std::bitset<n_lanes> &
mask)
3927 const auto src_data = internal::get_vector_data<n_components_>(
3930 this->dof_access_index ==
3932 this->active_fe_index,
3936 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3939 this->dof_values_initialized =
true;
3949 typename VectorizedArrayType>
3950template <
typename VectorType>
3954 const unsigned int first_index,
3955 const std::bitset<n_lanes> &
mask)
const
3958 Assert(this->dof_values_initialized ==
true,
3962 apply_hanging_node_constraints(
true);
3964 const auto dst_data = internal::get_vector_data<n_components_>(
3967 this->dof_access_index ==
3969 this->active_fe_index,
3974 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
3983 typename VectorizedArrayType>
3984template <
typename VectorType>
3988 const unsigned int first_index,
3989 const std::bitset<n_lanes> &
mask)
const
3992 Assert(this->dof_values_initialized ==
true,
3996 const auto dst_data = internal::get_vector_data<n_components_>(
3999 this->dof_access_index ==
4001 this->active_fe_index,
4005 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4014 typename VectorizedArrayType>
4015template <
typename VectorType>
4019 const unsigned int first_index,
4020 const std::bitset<n_lanes> &
mask)
const
4023 Assert(this->dof_values_initialized ==
true,
4027 const auto dst_data = internal::get_vector_data<n_components_>(
4030 this->dof_access_index ==
4032 this->active_fe_index,
4036 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4049 typename VectorizedArrayType>
4055 VectorizedArrayType>::value_type
4060 if constexpr (n_components == 1)
4061 return this->values_dofs[dof];
4064 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4066 for (
unsigned int comp = 0; comp < n_components; ++comp)
4067 return_value[comp] = this->values_dofs[comp * dofs + dof];
4068 return return_value;
4078 typename VectorizedArrayType>
4084 VectorizedArrayType>::value_type
4086 get_value(
const unsigned int q_point)
const
4089 Assert(this->values_quad_initialized ==
true,
4094 if constexpr (n_components == 1)
4095 return this->values_quad[q_point];
4098 if (n_components == dim &&
4099 this->data->element_type ==
4104 Assert(this->values_quad_initialized ==
true,
4109 Assert(this->J_value !=
nullptr,
4112 const std::size_t nqp = this->n_quadrature_points;
4120 const VectorizedArrayType inv_det =
4121 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4122 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4123 this->jacobian[0][2][2];
4126 for (
unsigned int comp = 0; comp < n_components; ++comp)
4127 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4128 jac[comp][comp] * inv_det;
4135 this->jacobian[q_point] :
4144 const VectorizedArrayType inv_det =
4145 (is_face && dim == 2 && this->get_face_no() < 2) ?
4149 for (
unsigned int comp = 0; comp < n_components; ++comp)
4151 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4152 for (
unsigned int e = 1;
e < dim; ++
e)
4154 this->values_quad[e * nqp + q_point] * jac[comp][e];
4155 value_out[comp] *= inv_det;
4162 const std::size_t nqp = this->n_quadrature_points;
4164 for (
unsigned int comp = 0; comp < n_components; ++comp)
4165 return_value[comp] = this->values_quad[comp * nqp + q_point];
4166 return return_value;
4177 typename VectorizedArrayType>
4183 VectorizedArrayType>::gradient_type
4188 Assert(this->gradients_quad_initialized ==
true,
4193 Assert(this->jacobian !=
nullptr,
4195 "update_gradients"));
4196 const std::size_t nqp = this->n_quadrature_points;
4198 if constexpr (n_components == dim && dim > 1)
4200 if (this->data->element_type ==
4205 Assert(this->gradients_quad_initialized ==
true,
4210 Assert(this->jacobian !=
nullptr,
4212 "update_gradients"));
4213 const std::size_t nqp = this->n_quadrature_points;
4214 const std::size_t nqp_d = nqp * dim;
4217 this->gradients_quad + q_point * dim;
4228 const VectorizedArrayType inv_det =
4229 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4230 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4231 this->jacobian[0][2][2];
4234 for (
unsigned int d = 0;
d < dim; ++
d)
4235 for (
unsigned int comp = 0; comp < n_components; ++comp)
4236 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4238 (jac[comp][comp] * inv_det);
4250 const VectorizedArrayType inv_det =
4251 (is_face && dim == 2 && this->get_face_no() < 2) ?
4255 VectorizedArrayType tmp[dim][dim];
4257 for (
unsigned int d = 0;
d < dim; ++
d)
4258 for (
unsigned int e = 0;
e < dim; ++
e)
4261 for (
unsigned int f = 1; f < dim; ++f)
4262 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4264 for (
unsigned int comp = 0; comp < n_components; ++comp)
4265 for (
unsigned int d = 0;
d < dim; ++
d)
4267 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4268 for (
unsigned int f = 1; f < dim; ++f)
4269 res += jac[comp][f] * tmp[d][f];
4271 grad_out[comp][
d] = res * inv_det;
4282 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4284 "update_hessians"));
4286 const auto jac_grad =
4287 this->jacobian_gradients_non_inverse[q_point];
4289 this->jacobian[q_point];
4293 const VectorizedArrayType inv_det =
4294 (is_face && dim == 2 && this->get_face_no() < 2) ?
4301 VectorizedArrayType tmp[dim][dim];
4302 for (
unsigned int d = 0;
d < dim; ++
d)
4303 for (
unsigned int e = 0;
e < dim; ++
e)
4306 for (
unsigned int f = 1; f < dim; ++f)
4307 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4312 for (
unsigned int d = 0;
d < dim; ++
d)
4314 for (
unsigned int e = 0;
e < dim; ++
e)
4316 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4317 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4318 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4321 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4323 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4328 for (
unsigned int d = 0;
d < dim; ++
d)
4329 for (
unsigned int e = 0;
e < dim; ++
e)
4331 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4332 for (
unsigned int f = 1; f < dim; ++f)
4333 res += tmp[f][d] * inv_t_jac[e][f];
4334 grad_out[
d][
e] = res;
4340 VectorizedArrayType tmp3[dim], tmp4[dim];
4341 for (
unsigned int d = 0;
d < dim; ++
d)
4343 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4344 for (
unsigned int e = 1;
e < dim; ++
e)
4345 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4347 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4348 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4349 for (
unsigned int d = 0;
d < dim; ++
d)
4351 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4352 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4354 for (
unsigned int d = 0;
d < dim; ++
d)
4356 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4357 for (
unsigned int e = 1;
e < dim; ++
e)
4358 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4361 VectorizedArrayType tmp2[dim];
4362 for (
unsigned int d = 0;
d < dim; ++
d)
4364 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4365 for (
unsigned e = 1;
e < dim; ++
e)
4367 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4370 for (
unsigned int d = 0;
d < dim; ++
d)
4371 for (
unsigned int e = 0;
e < dim; ++
e)
4373 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4377 grad_out[
d][
e] *= inv_det;
4388 for (
unsigned int comp = 0; comp < n_components; ++comp)
4389 for (
unsigned int d = 0;
d < dim; ++
d)
4391 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4392 this->jacobian[0][
d][
d];
4401 for (
unsigned int comp = 0; comp < n_components; ++comp)
4402 for (
unsigned int d = 0;
d < dim; ++
d)
4405 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4406 for (
unsigned int e = 1;
e < dim; ++
e)
4407 grad_out[comp][d] +=
4409 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4412 if constexpr (n_components == 1)
4424 typename VectorizedArrayType>
4430 VectorizedArrayType>::value_type
4436 Assert(this->gradients_quad_initialized ==
true,
4440 Assert(this->normal_x_jacobian !=
nullptr,
4442 "update_gradients"));
4444 const std::size_t nqp = this->n_quadrature_points;
4448 for (
unsigned int comp = 0; comp < n_components; ++comp)
4450 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4451 (this->normal_x_jacobian[0][dim - 1]);
4454 const std::size_t
index =
4456 for (
unsigned int comp = 0; comp < n_components; ++comp)
4458 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4459 this->normal_x_jacobian[
index][0];
4460 for (
unsigned int d = 1;
d < dim; ++
d)
4462 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4463 this->normal_x_jacobian[
index][
d];
4466 if constexpr (n_components == 1)
4478 template <
typename VectorizedArrayType>
4481 const VectorizedArrayType *
const hessians,
4483 VectorizedArrayType (&tmp)[1][1])
4485 tmp[0][0] = jac[0][0] *
hessians[0];
4488 template <
typename VectorizedArrayType>
4491 const VectorizedArrayType *
const hessians,
4492 const unsigned int nqp,
4493 VectorizedArrayType (&tmp)[2][2])
4495 for (
unsigned int d = 0;
d < 2; ++
d)
4503 template <
typename VectorizedArrayType>
4506 const VectorizedArrayType *
const hessians,
4507 const unsigned int nqp,
4508 VectorizedArrayType (&tmp)[3][3])
4510 for (
unsigned int d = 0;
d < 3; ++
d)
4531 typename VectorizedArrayType>
4536 VectorizedArrayType>::hessian_type
4541 Assert(this->hessians_quad_initialized ==
true,
4546 Assert(this->jacobian !=
nullptr,
4556 const std::size_t nqp = this->n_quadrature_points;
4557 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4562 for (
unsigned int comp = 0; comp < n_components; ++comp)
4564 for (
unsigned int d = 0;
d < dim; ++
d)
4565 hessian_out[comp][d][d] =
4566 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4567 (jac[
d][
d] * jac[
d][
d]);
4573 hessian_out[comp][0][1] =
4574 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4575 (jac[0][0] * jac[1][1]);
4578 hessian_out[comp][0][1] =
4579 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4580 (jac[0][0] * jac[1][1]);
4581 hessian_out[comp][0][2] =
4582 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4583 (jac[0][0] * jac[2][2]);
4584 hessian_out[comp][1][2] =
4585 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4586 (jac[1][1] * jac[2][2]);
4591 for (
unsigned int d = 0;
d < dim; ++
d)
4592 for (
unsigned int e = d + 1;
e < dim; ++
e)
4593 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4599 for (
unsigned int comp = 0; comp < n_components; ++comp)
4601 VectorizedArrayType tmp[dim][dim];
4602 internal::hessian_unit_times_jac(
4603 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4606 for (
unsigned int d = 0;
d < dim; ++
d)
4607 for (
unsigned int e = d;
e < dim; ++
e)
4609 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4610 for (
unsigned int f = 1; f < dim; ++f)
4611 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4618 for (
unsigned int d = 0;
d < dim; ++
d)
4619 for (
unsigned int e = d + 1;
e < dim; ++
e)
4620 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4626 const auto &jac_grad = this->jacobian_gradients[q_point];
4627 for (
unsigned int comp = 0; comp < n_components; ++comp)
4629 VectorizedArrayType tmp[dim][dim];
4630 internal::hessian_unit_times_jac(
4631 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4634 for (
unsigned int d = 0;
d < dim; ++
d)
4635 for (
unsigned int e = d;
e < dim; ++
e)
4637 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4638 for (
unsigned int f = 1; f < dim; ++f)
4639 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4643 for (
unsigned int d = 0;
d < dim; ++
d)
4644 for (
unsigned int e = 0;
e < dim; ++
e)
4645 hessian_out[comp][d][d] +=
4647 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4650 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4651 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4652 for (
unsigned int f = 0; f < dim; ++f)
4653 hessian_out[comp][d][e] +=
4654 jac_grad[count][f] *
4655 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4658 for (
unsigned int d = 0;
d < dim; ++
d)
4659 for (
unsigned int e = d + 1;
e < dim; ++
e)
4660 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4663 if constexpr (n_components == 1)
4664 return hessian_out[0];
4675 typename VectorizedArrayType>
4680 VectorizedArrayType>::gradient_type
4686 Assert(this->hessians_quad_initialized ==
true,
4697 const std::size_t nqp = this->n_quadrature_points;
4698 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4704 for (
unsigned int comp = 0; comp < n_components; ++comp)
4705 for (
unsigned int d = 0;
d < dim; ++
d)
4706 hessian_out[comp][d] =
4707 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4708 (jac[
d][
d] * jac[
d][
d]);
4713 for (
unsigned int comp = 0; comp < n_components; ++comp)
4717 VectorizedArrayType tmp[dim][dim];
4718 internal::hessian_unit_times_jac(
4719 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4723 for (
unsigned int d = 0;
d < dim; ++
d)
4725 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4726 for (
unsigned int f = 1; f < dim; ++f)
4727 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4734 const auto &jac_grad = this->jacobian_gradients[q_point];
4735 for (
unsigned int comp = 0; comp < n_components; ++comp)
4739 VectorizedArrayType tmp[dim][dim];
4740 internal::hessian_unit_times_jac(
4741 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4745 for (
unsigned int d = 0;
d < dim; ++
d)
4747 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4748 for (
unsigned int f = 1; f < dim; ++f)
4749 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4752 for (
unsigned int d = 0;
d < dim; ++
d)
4753 for (
unsigned int e = 0;
e < dim; ++
e)
4754 hessian_out[comp][d] +=
4756 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4760 if constexpr (n_components == 1)
4761 return hessian_out[0];
4772 typename VectorizedArrayType>
4777 VectorizedArrayType>::value_type
4783 Assert(this->hessians_quad_initialized ==
true,
4788 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4789 if constexpr (n_components == 1)
4791 VectorizedArrayType
sum = hess_diag[0];
4792 for (
unsigned int d = 1;
d < dim; ++
d)
4793 sum += hess_diag[d];
4799 for (
unsigned int comp = 0; comp < n_components; ++comp)
4801 laplacian_out[comp] = hess_diag[comp][0];
4802 for (
unsigned int d = 1;
d < dim; ++
d)
4803 laplacian_out[comp] += hess_diag[comp][d];
4805 return laplacian_out;
4815 typename VectorizedArrayType>
4820 VectorizedArrayType>::value_type
4825 Assert(this->hessians_quad_initialized ==
true,
4830 Assert(this->normal_x_jacobian !=
nullptr,
4832 "update_hessians"));
4836 const std::size_t nqp = this->n_quadrature_points;
4837 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4841 const auto nxj = this->normal_x_jacobian[0];
4843 for (
unsigned int comp = 0; comp < n_components; ++comp)
4845 for (
unsigned int d = 0;
d < dim; ++
d)
4846 hessian_out[comp] +=
4847 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4848 (nxj[
d]) * (nxj[d]);
4855 hessian_out[comp] +=
4856 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4860 hessian_out[comp] +=
4861 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4863 hessian_out[comp] +=
4864 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4866 hessian_out[comp] +=
4867 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4875 if constexpr (n_components == 1)
4876 return hessian_out[0];
4883 const auto normal = this->normal_vector(q_point);
4884 return get_hessian(q_point) * normal * normal;
4894 typename VectorizedArrayType>
4900 this->dof_values_initialized =
true;
4902 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4904 for (
unsigned int comp = 0; comp < n_components; ++comp)
4905 if constexpr (n_components == 1)
4906 this->values_dofs[comp * dofs + dof] = val_in;
4908 this->values_dofs[comp * dofs + dof] = val_in[comp];
4917 typename VectorizedArrayType>
4920 submit_value(
const value_type val_in,
const unsigned int q_point)
4926 Assert(this->J_value !=
nullptr,
4930 this->values_quad_submitted =
true;
4933 const std::size_t nqp = this->n_quadrature_points;
4934 VectorizedArrayType *
values = this->values_quad + q_point;
4936 const VectorizedArrayType JxW =
4938 this->J_value[0] * this->quadrature_weights[q_point] :
4939 this->J_value[q_point];
4940 if constexpr (n_components == 1)
4941 values[0] = val_in * JxW;
4944 if (n_components == dim &&
4945 this->data->element_type ==
4950 Assert(this->J_value !=
nullptr,
4955 this->values_quad_submitted =
true;
4958 VectorizedArrayType *
values = this->values_quad + q_point;
4959 const std::size_t nqp = this->n_quadrature_points;
4965 const VectorizedArrayType weight =
4966 this->quadrature_weights[q_point];
4968 for (
unsigned int comp = 0; comp < n_components; ++comp)
4969 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4976 this->jacobian[q_point] :
4981 const VectorizedArrayType fac =
4983 this->quadrature_weights[q_point] :
4985 this->J_value[q_point] :
4986 this->J_value[0] * this->quadrature_weights[q_point]) *
4987 ((dim == 2 && this->get_face_no() < 2) ?
4996 for (
unsigned int comp = 0; comp < n_components; ++comp)
4998 values[comp * nqp] = val_in[0] * jac[0][comp];
4999 for (
unsigned int e = 1;
e < dim; ++
e)
5000 values[comp * nqp] += val_in[e] * jac[e][comp];
5001 values[comp * nqp] *= fac;
5006 for (
unsigned int comp = 0; comp < n_components; ++comp)
5007 values[comp * nqp] = val_in[comp] * JxW;
5017 typename VectorizedArrayType>
5018template <
int,
typename>
5022 const unsigned int q_point)
5024 static_assert(n_components == 1,
5025 "Do not try to modify the default template parameters used for"
5026 " selectively enabling this function via std::enable_if!");
5027 submit_value(val_in[0], q_point);
5036 typename VectorizedArrayType>
5039 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
5045 Assert(this->J_value !=
nullptr,
5047 "update_gradients"));
5048 Assert(this->jacobian !=
nullptr,
5050 "update_gradients"));
5052 this->gradients_quad_submitted =
true;
5055 if constexpr (dim > 1 && n_components == dim)
5057 if (this->data->element_type ==
5066 Assert(this->J_value !=
nullptr,
5068 "update_gradients"));
5069 Assert(this->jacobian !=
nullptr,
5071 "update_gradients"));
5073 this->gradients_quad_submitted =
true;
5076 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5077 VectorizedArrayType *
values =
5078 this->values_from_gradients_quad + q_point;
5079 const std::size_t nqp = this->n_quadrature_points;
5080 const std::size_t nqp_d = nqp * dim;
5090 const VectorizedArrayType weight =
5091 this->quadrature_weights[q_point];
5092 for (
unsigned int d = 0;
d < dim; ++
d)
5093 for (
unsigned int comp = 0; comp < n_components; ++comp)
5094 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5096 (jac[comp][comp] * weight);
5108 const VectorizedArrayType fac =
5110 this->quadrature_weights[q_point] :
5111 this->J_value[0] * this->quadrature_weights[q_point] *
5112 ((dim == 2 && this->get_face_no() < 2) ?
5117 VectorizedArrayType tmp[dim][dim];
5118 for (
unsigned int d = 0;
d < dim; ++
d)
5119 for (
unsigned int e = 0;
e < dim; ++
e)
5121 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5122 for (
unsigned int f = 1; f < dim; ++f)
5123 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5125 for (
unsigned int comp = 0; comp < n_components; ++comp)
5126 for (
unsigned int d = 0;
d < dim; ++
d)
5128 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5129 for (
unsigned int f = 1; f < dim; ++f)
5130 res += jac[f][comp] * tmp[d][f];
5139 const auto jac_grad =
5140 this->jacobian_gradients_non_inverse[q_point];
5142 this->jacobian[q_point];
5146 const VectorizedArrayType fac =
5147 (!is_face) ? this->quadrature_weights[q_point] :
5148 this->J_value[q_point] *
5149 ((dim == 2 && this->get_face_no() < 2) ?
5158 VectorizedArrayType tmp3[dim], tmp4[dim];
5159 for (
unsigned int d = 0;
d < dim; ++
d)
5161 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5162 for (
unsigned int e = 1;
e < dim; ++
e)
5163 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5165 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5166 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5167 for (
unsigned int d = 0;
d < dim; ++
d)
5169 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5170 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5172 for (
unsigned int d = 0;
d < dim; ++
d)
5174 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5175 for (
unsigned int e = 1;
e < dim; ++
e)
5176 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5182 VectorizedArrayType tmp[dim][dim];
5185 for (
unsigned int d = 0;
d < dim; ++
d)
5186 for (
unsigned int e = 0;
e < dim; ++
e)
5188 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5189 for (
unsigned int f = 1; f < dim; ++f)
5190 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5193 for (
unsigned int d = 0;
d < dim; ++
d)
5194 for (
unsigned int e = 0;
e < dim; ++
e)
5196 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5197 for (
unsigned int f = 1; f < dim; ++f)
5198 res += t_jac[d][f] * tmp[e][f];
5205 VectorizedArrayType
value[dim];
5206 for (
unsigned int d = 0;
d < dim; ++
d)
5208 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5209 for (
unsigned int e = 1;
e < dim; ++
e)
5210 value[d] += tmp[d][e] * jac_grad[d][e];
5212 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5213 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5214 for (
unsigned int d = 0;
d < dim; ++
d)
5216 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5217 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5222 for (
unsigned int d = 0;
d < dim; ++
d)
5224 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5225 for (
unsigned int e = 1;
e < dim; ++
e)
5226 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5227 for (
unsigned int e = 0;
e < dim; ++
e)
5228 value[e] -= t_jac[e][d] * tmp2;
5231 for (
unsigned int d = 0;
d < dim; ++
d)
5232 values[d * nqp] =
value[d];
5238 const std::size_t nqp_d = this->n_quadrature_points * dim;
5239 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5243 const VectorizedArrayType JxW =
5244 this->J_value[0] * this->quadrature_weights[q_point];
5250 std::array<VectorizedArrayType, dim> jac;
5251 for (
unsigned int d = 0;
d < dim; ++
d)
5252 jac[d] = this->jacobian[0][d][d];
5254 for (
unsigned int d = 0;
d < dim; ++
d)
5256 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5257 if constexpr (n_components == 1)
5258 gradients[d] = grad_in[d] * factor;
5260 for (
unsigned int comp = 0; comp < n_components; ++comp)
5261 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5268 this->jacobian[q_point] :
5270 const VectorizedArrayType JxW =
5272 this->J_value[q_point] :
5273 this->J_value[0] * this->quadrature_weights[q_point];
5274 if constexpr (n_components == 1)
5275 for (
unsigned int d = 0;
d < dim; ++
d)
5277 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5278 for (
unsigned int e = 1;
e < dim; ++
e)
5279 new_val += (jac[e][d] * grad_in[e]);
5283 for (
unsigned int comp = 0; comp < n_components; ++comp)
5284 for (
unsigned int d = 0;
d < dim; ++
d)
5286 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5287 for (
unsigned int e = 1;
e < dim; ++
e)
5288 new_val += (jac[e][d] * grad_in[comp][e]);
5300 typename VectorizedArrayType>
5301template <
int,
typename>
5305 const unsigned int q_point)
5307 static_assert(n_components == 1 && dim == 1,
5308 "Do not try to modify the default template parameters used for"
5309 " selectively enabling this function via std::enable_if!");
5310 submit_gradient(grad_in[0], q_point);
5319 typename VectorizedArrayType>
5325 Assert(this->normal_x_jacobian !=
nullptr,
5327 "update_gradients"));
5329 this->gradients_quad_submitted =
true;
5332 const std::size_t nqp_d = this->n_quadrature_points * dim;
5333 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5337 const VectorizedArrayType JxW_jac = this->J_value[0] *
5338 this->quadrature_weights[q_point] *
5339 this->normal_x_jacobian[0][dim - 1];
5340 for (
unsigned int comp = 0; comp < n_components; ++comp)
5342 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5343 gradients[comp * nqp_d + d] = VectorizedArrayType();
5344 if constexpr (n_components == 1)
5345 gradients[dim - 1] = grad_in * JxW_jac;
5347 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5352 const unsigned int index =
5355 this->normal_x_jacobian[
index];
5356 const VectorizedArrayType JxW =
5358 this->J_value[
index] * this->quadrature_weights[q_point] :
5359 this->J_value[
index];
5360 for (
unsigned int comp = 0; comp < n_components; ++comp)
5361 for (
unsigned int d = 0;
d < dim; ++
d)
5362 if constexpr (n_components == 1)
5365 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5375 typename VectorizedArrayType>
5378 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5384 Assert(this->J_value !=
nullptr,
5386 "update_hessians"));
5387 Assert(this->jacobian !=
nullptr,
5389 "update_hessians"));
5391 this->hessians_quad_submitted =
true;
5395 const std::size_t nqp = this->n_quadrature_points;
5396 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5399 const VectorizedArrayType JxW =
5400 this->J_value[0] * this->quadrature_weights[q_point];
5403 for (
unsigned int d = 0;
d < dim; ++
d)
5405 const auto jac_d = this->jacobian[0][
d][
d];
5406 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5407 for (
unsigned int comp = 0; comp < n_components; ++comp)
5408 if constexpr (n_components == 1)
5409 this->hessians_quad[
d * nqp + q_point] =
5410 hessian_in[
d][
d] * factor;
5412 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5413 hessian_in[comp][d][d] * factor;
5417 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5418 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5420 const auto jac_d = this->jacobian[0][
d][
d];
5421 const auto jac_e = this->jacobian[0][
e][
e];
5422 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5423 for (
unsigned int comp = 0; comp < n_components; ++comp)
5424 if constexpr (n_components == 1)
5425 this->hessians_quad[off_dia * nqp + q_point] =
5426 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5428 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5429 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5436 const VectorizedArrayType JxW =
5437 this->J_value[0] * this->quadrature_weights[q_point];
5438 for (
unsigned int comp = 0; comp < n_components; ++comp)
5441 if constexpr (n_components == 1)
5442 hessian_c = hessian_in;
5444 hessian_c = hessian_in[comp];
5447 VectorizedArrayType tmp[dim][dim];
5448 for (
unsigned int i = 0; i < dim; ++i)
5449 for (
unsigned int j = 0; j < dim; ++j)
5451 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5452 for (
unsigned int k = 1; k < dim; ++k)
5453 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5457 VectorizedArrayType tmp2[dim][dim];
5458 for (
unsigned int i = 0; i < dim; ++i)
5459 for (
unsigned int j = 0; j < dim; ++j)
5461 tmp2[i][j] = jac[0][i] * tmp[0][j];
5462 for (
unsigned int k = 1; k < dim; ++k)
5463 tmp2[i][j] += jac[k][i] * tmp[k][j];
5467 for (
unsigned int d = 0;
d < dim; ++
d)
5468 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5472 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5473 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5474 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5475 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5481 const VectorizedArrayType JxW = this->J_value[q_point];
5482 const auto &jac_grad = this->jacobian_gradients[q_point];
5483 for (
unsigned int comp = 0; comp < n_components; ++comp)
5486 if constexpr (n_components == 1)
5487 hessian_c = hessian_in;
5489 hessian_c = hessian_in[comp];
5492 VectorizedArrayType tmp[dim][dim];
5493 for (
unsigned int i = 0; i < dim; ++i)
5494 for (
unsigned int j = 0; j < dim; ++j)
5496 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5497 for (
unsigned int k = 1; k < dim; ++k)
5498 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5502 VectorizedArrayType tmp2[dim][dim];
5503 for (
unsigned int i = 0; i < dim; ++i)
5504 for (
unsigned int j = 0; j < dim; ++j)
5506 tmp2[i][j] = jac[0][i] * tmp[0][j];
5507 for (
unsigned int k = 1; k < dim; ++k)
5508 tmp2[i][j] += jac[k][i] * tmp[k][j];
5512 for (
unsigned int d = 0;
d < dim; ++
d)
5513 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5517 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5518 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5519 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5520 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5523 for (
unsigned int d = 0;
d < dim; ++
d)
5525 VectorizedArrayType
sum = 0;
5526 for (
unsigned int e = 0;
e < dim; ++
e)
5527 sum += hessian_c[e][e] * jac_grad[e][d];
5528 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5529 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5531 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5532 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5545 typename VectorizedArrayType>
5549 const unsigned int q_point)
5555 Assert(this->J_value !=
nullptr,
5557 "update_hessians"));
5558 Assert(this->jacobian !=
nullptr,
5560 "update_hessians"));
5562 this->hessians_quad_submitted =
true;
5566 const std::size_t nqp = this->n_quadrature_points;
5567 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5570 const VectorizedArrayType JxW =
5571 this->J_value[0] * this->quadrature_weights[q_point];
5573 const auto nxj = this->normal_x_jacobian[0];
5576 for (
unsigned int d = 0;
d < dim; ++
d)
5578 const auto nxj_d = nxj[
d];
5579 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5580 for (
unsigned int comp = 0; comp < n_components; ++comp)
5581 if constexpr (n_components == 1)
5582 this->hessians_quad[
d * nqp + q_point] =
5583 normal_hessian_in * factor;
5585 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5586 normal_hessian_in[comp] * factor;
5590 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5591 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5593 const auto jac_d = nxj[
d];
5594 const auto jac_e = nxj[
e];
5595 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5596 for (
unsigned int comp = 0; comp < n_components; ++comp)
5597 if constexpr (n_components == 1)
5598 this->hessians_quad[off_dia * nqp + q_point] =
5599 2. * normal_hessian_in * factor;
5601 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5602 2. * normal_hessian_in[comp] * factor;
5607 const auto normal = this->normal_vector(q_point);
5608 const auto normal_projector =
outer_product(normal, normal);
5609 if constexpr (n_components == 1)
5610 submit_hessian(normal_hessian_in * normal_projector, q_point);
5612 submit_hessian(
outer_product(normal_hessian_in, normal_projector),
5623 typename VectorizedArrayType>
5628 VectorizedArrayType>::value_type
5634 Assert(this->values_quad_submitted ==
true,
5639 const std::size_t nqp = this->n_quadrature_points;
5640 for (
unsigned int q = 0; q < nqp; ++q)
5641 for (
unsigned int comp = 0; comp < n_components; ++comp)
5642 return_value[comp] += this->values_quad[comp * nqp + q];
5643 if constexpr (n_components == 1)
5644 return return_value[0];
5646 return return_value;
5655 typename VectorizedArrayType>
5656template <
int,
typename>
5661 static_assert(n_components == dim,
5662 "Do not try to modify the default template parameters used for"
5663 " selectively enabling this function via std::enable_if!");
5666 Assert(this->gradients_quad_initialized ==
true,
5670 Assert(this->jacobian !=
nullptr,
5672 "update_gradients"));
5674 VectorizedArrayType divergence;
5675 const std::size_t nqp = this->n_quadrature_points;
5678 this->data->element_type ==
5681 VectorizedArrayType inv_det =
5684 this->jacobian[0][0][0] *
5685 ((dim == 2) ? this->jacobian[0][1][1] :
5686 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5694 if (is_face && dim == 2 && this->get_face_no() < 2)
5698 divergence = this->gradients_quad[q_point * dim];
5699 for (
unsigned int d = 1;
d < dim; ++
d)
5700 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5701 divergence *= inv_det;
5710 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5711 for (
unsigned int d = 1;
d < dim; ++
d)
5712 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5713 this->jacobian[0][
d][
d];
5720 this->jacobian[q_point] :
5722 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5723 for (
unsigned int e = 1;
e < dim; ++
e)
5724 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5725 for (
unsigned int d = 1;
d < dim; ++
d)
5726 for (
unsigned int e = 0;
e < dim; ++
e)
5728 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5740 typename VectorizedArrayType>
5741template <
int,
typename>
5746 static_assert(n_components == dim,
5747 "Do not try to modify the default template parameters used for"
5748 " selectively enabling this function via std::enable_if!");
5751 const auto grad = get_gradient(q_point);
5752 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5753 VectorizedArrayType half = Number(0.5);
5754 for (
unsigned int d = 0;
d < dim; ++
d)
5755 symmetrized[d] = grad[d][d];
5761 symmetrized[2] = grad[0][1] + grad[1][0];
5762 symmetrized[2] *= half;
5765 symmetrized[3] = grad[0][1] + grad[1][0];
5766 symmetrized[3] *= half;
5767 symmetrized[4] = grad[0][2] + grad[2][0];
5768 symmetrized[4] *= half;
5769 symmetrized[5] = grad[1][2] + grad[2][1];
5770 symmetrized[5] *= half;
5784 typename VectorizedArrayType>
5785template <
int,
typename>
5787 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5789 get_curl(const unsigned
int q_point) const
5791 static_assert(dim > 1 && n_components == dim,
5792 "Do not try to modify the default template parameters used for"
5793 " selectively enabling this function via std::enable_if!");
5797 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5801 curl[0] = grad[1][0] - grad[0][1];
5804 curl[0] = grad[2][1] - grad[1][2];
5805 curl[1] = grad[0][2] - grad[2][0];
5806 curl[2] = grad[1][0] - grad[0][1];
5820 typename VectorizedArrayType>
5821template <
int,
typename>
5825 const unsigned int q_point)
5827 static_assert(n_components == dim,
5828 "Do not try to modify the default template parameters used for"
5829 " selectively enabling this function via std::enable_if!");
5835 Assert(this->J_value !=
nullptr,
5837 "update_gradients"));
5838 Assert(this->jacobian !=
nullptr,
5840 "update_gradients"));
5842 this->gradients_quad_submitted =
true;
5845 const std::size_t nqp_d = this->n_quadrature_points * dim;
5846 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5848 if (this->data->element_type ==
5855 const VectorizedArrayType fac =
5857 this->quadrature_weights[q_point] * div_in :
5859 this->J_value[q_point] :
5860 this->J_value[0] * this->quadrature_weights[q_point]) *
5863 this->jacobian[this->cell_type >
5867 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5869 for (
unsigned int d = 0;
d < dim; ++
d)
5871 for (
unsigned int e = 0;
e < dim; ++
e)
5872 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5874 this->divergence_is_requested =
true;
5881 const VectorizedArrayType fac =
5882 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5883 for (
unsigned int d = 0;
d < dim; ++
d)
5885 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5886 for (
unsigned int e = 0;
e < dim; ++
e)
5887 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5894 this->jacobian[q_point] :
5896 const VectorizedArrayType fac =
5898 this->J_value[q_point] :
5899 this->J_value[0] * this->quadrature_weights[q_point]) *
5901 for (
unsigned int d = 0;
d < dim; ++
d)
5903 for (
unsigned int e = 0;
e < dim; ++
e)
5904 gradients[d * nqp_d + e] = jac[d][e] * fac;
5916 typename VectorizedArrayType>
5917template <
int,
typename>
5922 const unsigned int q_point)
5924 static_assert(n_components == dim,
5925 "Do not try to modify the default template parameters used for"
5926 " selectively enabling this function via std::enable_if!");
5929 this->data->element_type !=
5940 Assert(this->J_value !=
nullptr,
5942 "update_gradients"));
5943 Assert(this->jacobian !=
nullptr,
5945 "update_gradients"));
5947 this->gradients_quad_submitted =
true;
5950 const std::size_t nqp_d = this->n_quadrature_points * dim;
5951 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5954 const VectorizedArrayType JxW =
5955 this->J_value[0] * this->quadrature_weights[q_point];
5957 for (
unsigned int d = 0;
d < dim; ++
d)
5958 gradients[d * nqp_d + d] =
5960 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5961 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5963 const VectorizedArrayType
value =
5972 const VectorizedArrayType JxW =
5974 this->J_value[q_point] :
5975 this->J_value[0] * this->quadrature_weights[q_point];
5978 this->jacobian[q_point] :
5980 VectorizedArrayType weighted[dim][dim];
5981 for (
unsigned int i = 0; i < dim; ++i)
5983 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
5984 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
5986 const VectorizedArrayType
value =
5988 weighted[i][j] =
value;
5989 weighted[j][i] =
value;
5991 for (
unsigned int comp = 0; comp < dim; ++comp)
5992 for (
unsigned int d = 0;
d < dim; ++
d)
5994 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
5995 for (
unsigned int e = 1;
e < dim; ++
e)
5996 new_val += jac[e][d] * weighted[comp][e];
6008 typename VectorizedArrayType>
6009template <
int,
typename>
6013 const unsigned int q_point)
6015 static_assert(n_components == dim,
6016 "Do not try to modify the default template parameters used for"
6017 " selectively enabling this function via std::enable_if!");
6023 grad[1][0] = curl[0];
6024 grad[0][1] = -curl[0];
6027 grad[2][1] = curl[0];
6028 grad[1][2] = -curl[0];
6029 grad[0][2] = curl[1];
6030 grad[2][0] = -curl[1];
6031 grad[1][0] = curl[2];
6032 grad[0][1] = -curl[2];
6037 submit_gradient(grad, q_point);
6050 typename VectorizedArrayType>
6056 VectorizedArrayType>::
6058 const unsigned int fe_no,
6059 const unsigned int quad_no,
6060 const unsigned int first_selected_component,
6061 const unsigned int active_fe_index,
6062 const unsigned int active_quad_index)
6063 : BaseClass(matrix_free,
6065 first_selected_component,
6072 numbers::invalid_unsigned_int )
6073 , dofs_per_component(this->data->dofs_per_component_on_cell)
6074 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6075 , n_q_points(this->data->n_q_points)
6077 check_template_arguments(fe_no, 0);
6087 typename VectorizedArrayType>
6093 VectorizedArrayType>::
6095 const std::pair<unsigned int, unsigned int> &range,
6096 const unsigned int dof_no,
6097 const unsigned int quad_no,
6098 const unsigned int first_selected_component)
6102 first_selected_component,
6103 matrix_free.get_cell_active_fe_index(range))
6113 typename VectorizedArrayType>
6119 VectorizedArrayType>::
6124 const unsigned int first_selected_component)
6125 : BaseClass(mapping,
6129 first_selected_component,
6131 , dofs_per_component(this->data->dofs_per_component_on_cell)
6132 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6133 , n_q_points(this->data->n_q_points)
6145 typename VectorizedArrayType>
6151 VectorizedArrayType>::
6155 const unsigned int first_selected_component)
6160 first_selected_component,
6162 , dofs_per_component(this->data->dofs_per_component_on_cell)
6163 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6164 , n_q_points(this->data->n_q_points)
6176 typename VectorizedArrayType>
6182 VectorizedArrayType>::
6185 const unsigned int first_selected_component)
6186 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6188 other.mapped_geometry->get_quadrature(),
6189 other.mapped_geometry->get_fe_values().get_update_flags(),
6190 first_selected_component,
6192 , dofs_per_component(this->data->dofs_per_component_on_cell)
6193 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6194 , n_q_points(this->data->n_q_points)
6206 typename VectorizedArrayType>
6215 , dofs_per_component(this->data->dofs_per_component_on_cell)
6216 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6217 , n_q_points(this->data->n_q_points)
6229 typename VectorizedArrayType>
6235 VectorizedArrayType> &
6241 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6243 BaseClass::operator=(other);
6255 typename VectorizedArrayType>
6262 VectorizedArrayType>::
6263 check_template_arguments(
const unsigned int dof_no,
6264 const unsigned int first_selected_component)
6267 (void)first_selected_component;
6270 this->data->dofs_per_component_on_cell > 0,
6272 "There is nothing useful you can do with an FEEvaluation object with "
6273 "FE_Nothing, i.e., without DoFs! If you have passed to "
6274 "MatrixFree::reinit() a collection of finite elements also containing "
6275 "FE_Nothing, please check - before creating FEEvaluation - the category "
6276 "of the current range by calling either "
6277 "MatrixFree::get_cell_range_category(range) or "
6278 "MatrixFree::get_face_range_category(range). The returned category "
6279 "is the index of the active FE, which you can use to exclude "
6286 static_cast<unsigned int>(fe_degree) !=
6287 this->data->data.front().fe_degree) ||
6288 n_q_points != this->n_quadrature_points)
6290 std::string message =
6291 "-------------------------------------------------------\n";
6292 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6293 message +=
" Called --> FEEvaluation<dim,";
6297 message +=
",Number>(data";
6313 if (
static_cast<unsigned int>(fe_degree) ==
6314 this->data->data.front().fe_degree)
6316 proposed_dof_comp = dof_no;
6317 proposed_fe_comp = first_selected_component;
6320 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6322 for (
unsigned int nf = 0;
6325 if (this->matrix_free
6326 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6328 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6330 proposed_dof_comp = no;
6331 proposed_fe_comp = nf;
6335 this->mapping_data->descriptor[this->active_quad_index]
6337 proposed_quad_comp = this->quad_no;
6339 for (
unsigned int no = 0;
6344 .descriptor[this->active_quad_index]
6345 .n_q_points == n_q_points)
6347 proposed_quad_comp = no;
6354 if (proposed_dof_comp != first_selected_component)
6355 message +=
"Wrong vector component selection:\n";
6357 message +=
"Wrong quadrature formula selection:\n";
6358 message +=
" Did you mean FEEvaluation<dim,";
6362 message +=
",Number>(data";
6371 std::string correct_pos;
6372 if (proposed_dof_comp != dof_no)
6373 correct_pos =
" ^ ";
6376 if (proposed_quad_comp != this->quad_no)
6377 correct_pos +=
" ^ ";
6380 if (proposed_fe_comp != first_selected_component)
6381 correct_pos +=
" ^\n";
6383 correct_pos +=
" \n";
6389 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6390 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6391 message +=
"Wrong template arguments:\n";
6392 message +=
" Did you mean FEEvaluation<dim,";
6397 message +=
",Number>(data";
6405 std::string correct_pos;
6406 if (this->data->data.front().fe_degree !=
6407 static_cast<unsigned int>(fe_degree))
6411 if (proposed_n_q_points_1d != n_q_points_1d)
6412 correct_pos +=
" ^\n";
6414 correct_pos +=
" \n";
6415 message +=
" " + correct_pos;
6417 Assert(
static_cast<unsigned int>(fe_degree) ==
6418 this->data->data.front().fe_degree &&
6419 n_q_points == this->n_quadrature_points,
6425 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6436 typename VectorizedArrayType>
6445 Assert(this->matrix_free !=
nullptr,
6446 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6447 " Integer indexing is not possible."));
6455 const unsigned int offsets =
6456 this->mapping_data->data_index_offsets[
cell_index];
6457 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6458 this->J_value = &this->mapping_data->JxW_values[offsets];
6459 if (!this->mapping_data->jacobian_gradients[0].empty())
6461 this->jacobian_gradients =
6462 this->mapping_data->jacobian_gradients[0].data() + offsets;
6463 this->jacobian_gradients_non_inverse =
6464 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6470 for (
unsigned int i = 0; i < n_lanes; ++i)
6471 this->cell_ids[i] =
cell_index * n_lanes + i;
6478 this->cell_ids[i] =
cell_index * n_lanes + i;
6479 for (; i < n_lanes; ++i)
6483 if (this->mapping_data->quadrature_points.empty() ==
false)
6485 &this->mapping_data->quadrature_points
6486 [this->mapping_data->quadrature_point_offsets[this->cell]];
6489 this->is_reinitialized =
true;
6490 this->dof_values_initialized =
false;
6491 this->values_quad_initialized =
false;
6492 this->gradients_quad_initialized =
false;
6493 this->hessians_quad_initialized =
false;
6504 typename VectorizedArrayType>
6511 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6518 this->cell_ids = cell_ids;
6523 for (
unsigned int v = 0; v < n_lanes; ++v)
6537 if (this->mapped_geometry ==
nullptr)
6538 this->mapped_geometry =
6539 std::make_shared<internal::MatrixFreeFunctions::
6540 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6542 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6544 auto &this_jacobian_data = mapping_storage.jacobians[0];
6545 auto &this_J_value_data = mapping_storage.JxW_values;
6546 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6547 auto &this_jacobian_gradients_non_inverse_data =
6548 mapping_storage.jacobian_gradients_non_inverse[0];
6549 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6553 if (this_jacobian_data.size() != 2)
6554 this_jacobian_data.resize_fast(2);
6556 if (this_J_value_data.size() != 1)
6557 this_J_value_data.resize_fast(1);
6559 const auto &update_flags_cells =
6563 this_jacobian_gradients_data.size() != 1)
6565 this_jacobian_gradients_data.resize_fast(1);
6566 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6570 this_quadrature_points_data.size() != 1)
6571 this_quadrature_points_data.resize_fast(1);
6575 if (this_jacobian_data.size() != this->n_quadrature_points)
6576 this_jacobian_data.resize_fast(this->n_quadrature_points);
6578 if (this_J_value_data.size() != this->n_quadrature_points)
6579 this_J_value_data.resize_fast(this->n_quadrature_points);
6581 const auto &update_flags_cells =
6585 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6587 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6588 this_jacobian_gradients_non_inverse_data.resize_fast(
6589 this->n_quadrature_points);
6593 this_quadrature_points_data.size() != this->n_quadrature_points)
6594 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6598 this->jacobian = this_jacobian_data.data();
6599 this->J_value = this_J_value_data.data();
6600 this->jacobian_gradients = this_jacobian_gradients_data.data();
6601 this->jacobian_gradients_non_inverse =
6602 this_jacobian_gradients_non_inverse_data.data();
6606 for (
unsigned int v = 0; v < n_lanes; ++v)
6613 const unsigned int cell_batch_index =
cell_index / n_lanes;
6614 const unsigned int offsets =
6615 this->mapping_data->data_index_offsets[cell_batch_index];
6616 const unsigned int lane =
cell_index % n_lanes;
6618 if (this->cell_type <=
6622 for (
unsigned int q = 0; q < 2; ++q)
6623 for (
unsigned int i = 0; i < dim; ++i)
6624 for (
unsigned int j = 0; j < dim; ++j)
6625 this_jacobian_data[q][i][j][v] =
6626 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6628 const unsigned int q = 0;
6630 this_J_value_data[q][v] =
6631 this->mapping_data->JxW_values[offsets + q][lane];
6633 const auto &update_flags_cells =
6638 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6639 for (
unsigned int j = 0; j < dim; ++j)
6640 this_jacobian_gradients_data[q][i][j][v] =
6642 ->jacobian_gradients[0][offsets + q][i][j][lane];
6644 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6645 for (
unsigned int j = 0; j < dim; ++j)
6646 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6648 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6653 for (
unsigned int i = 0; i < dim; ++i)
6654 this_quadrature_points_data[q][i][v] =
6655 this->mapping_data->quadrature_points
6657 ->quadrature_point_offsets[cell_batch_index] +
6663 const auto cell_type =
6667 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6669 const unsigned int q_src =
6675 this_J_value_data[q][v] =
6676 this->mapping_data->JxW_values[offsets + q_src][lane];
6678 for (
unsigned int i = 0; i < dim; ++i)
6679 for (
unsigned int j = 0; j < dim; ++j)
6680 this_jacobian_data[q][i][j][v] =
6682 ->jacobians[0][offsets + q_src][i][j][lane];
6684 const auto &update_flags_cells =
6689 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6690 for (
unsigned int j = 0; j < dim; ++j)
6691 this_jacobian_gradients_data[q][i][j][v] =
6693 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6695 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6696 for (
unsigned int j = 0; j < dim; ++j)
6697 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6699 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6712 this->mapping_data->quadrature_points
6714 ->quadrature_point_offsets[cell_batch_index] +
6718 this->mapping_data->jacobians[0][offsets + 1];
6720 for (
unsigned int d = 0;
d < dim; ++
d)
6723 static_cast<Number
>(
6724 this->descriptor->quadrature.point(q)[
d]);
6726 for (
unsigned int d = 0;
d < dim; ++
d)
6727 for (
unsigned int e = 0;
e < dim; ++
e)
6730 static_cast<Number
>(
6731 this->descriptor->quadrature.point(q)[
e]);
6733 for (
unsigned int i = 0; i < dim; ++i)
6734 this_quadrature_points_data[q][i][v] = point[i][lane];
6739 for (
unsigned int i = 0; i < dim; ++i)
6740 this_quadrature_points_data[q][i][v] =
6741 this->mapping_data->quadrature_points
6743 ->quadrature_point_offsets[cell_batch_index] +
6752 this->is_reinitialized =
true;
6753 this->dof_values_initialized =
false;
6754 this->values_quad_initialized =
false;
6755 this->gradients_quad_initialized =
false;
6756 this->hessians_quad_initialized =
false;
6767 typename VectorizedArrayType>
6768template <
bool level_dof_access>
6775 VectorizedArrayType>
::
6778 Assert(this->matrix_free ==
nullptr,
6779 ExcMessage(
"Cannot use initialization from cell iterator if "
6780 "initialized from MatrixFree object. Use variant for "
6781 "on the fly computation with arguments as for FEValues "
6784 this->mapped_geometry->reinit(
6786 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6787 if (level_dof_access)
6788 cell->get_mg_dof_indices(this->local_dof_indices);
6790 cell->get_dof_indices(this->local_dof_indices);
6793 this->is_reinitialized =
true;
6804 typename VectorizedArrayType>
6811 VectorizedArrayType>
::
6814 Assert(this->matrix_free == 0,
6815 ExcMessage(
"Cannot use initialization from cell iterator if "
6816 "initialized from MatrixFree object. Use variant for "
6817 "on the fly computation with arguments as for FEValues "
6820 this->mapped_geometry->reinit(cell);
6823 this->is_reinitialized =
true;
6834 typename VectorizedArrayType>
6841 VectorizedArrayType>::
6845 Assert(this->dof_values_initialized ==
true,
6848 evaluate(this->values_dofs, evaluation_flags);
6858 typename VectorizedArrayType>
6865 VectorizedArrayType>::
6866 evaluate(
const VectorizedArrayType *values_array,
6869 const bool hessians_on_general_cells =
6873 if (hessians_on_general_cells)
6876 if (this->data->element_type ==
6882 if constexpr (fe_degree > -1)
6885 template run<fe_degree, n_q_points_1d>(n_components,
6886 evaluation_flag_actual,
6894 evaluation_flag_actual,
6895 const_cast<VectorizedArrayType *
>(values_array),
6900 this->values_quad_initialized =
6902 this->gradients_quad_initialized =
6904 this->hessians_quad_initialized =
6915 template <
typename Number,
6916 typename VectorizedArrayType,
6918 typename EvaluatorType,
6919 std::enable_if_t<internal::has_begin<VectorType> &&
6922 VectorizedArrayType *
6923 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6929 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6930 const auto &dof_info = fe_eval.get_dof_info();
6937 if (std::is_same_v<typename VectorType::value_type, Number> &&
6941 interleaved_contiguous &&
6942 reinterpret_cast<
std::size_t>(
6944 dof_info.dof_indices_contiguous
6945 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6946 [cell * VectorizedArrayType::size()]) %
6947 sizeof(VectorizedArrayType) ==
6950 return reinterpret_cast<VectorizedArrayType *
>(
6954 [cell * VectorizedArrayType::size()] +
6956 [fe_eval.get_active_fe_index()]
6957 [fe_eval.get_first_selected_component()] *
6958 VectorizedArrayType::size());
6967 template <
typename Number,
6968 typename VectorizedArrayType,
6970 typename EvaluatorType,
6971 std::enable_if_t<!internal::has_begin<VectorType> ||
6974 VectorizedArrayType *
6975 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
6988 typename VectorizedArrayType>
6989template <
typename VectorType>
6996 VectorizedArrayType>::
6997 gather_evaluate(
const VectorType &input_vector,
7000 const VectorizedArrayType *src_ptr =
7001 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7002 *
this, input_vector);
7003 if (src_ptr !=
nullptr)
7004 evaluate(src_ptr, evaluation_flag);
7007 this->read_dof_values(input_vector);
7008 evaluate(this->begin_dof_values(), evaluation_flag);
7019 typename VectorizedArrayType>
7026 VectorizedArrayType>::
7029 integrate(integration_flag, this->values_dofs);
7032 this->dof_values_initialized =
true;
7043 typename VectorizedArrayType>
7050 VectorizedArrayType>::
7052 VectorizedArrayType *values_array,
7053 const bool sum_into_values_array)
7057 Assert(this->values_quad_submitted ==
true,
7060 Assert(this->gradients_quad_submitted ==
true,
7063 Assert(this->hessians_quad_submitted ==
true,
7066 Assert(this->matrix_free !=
nullptr ||
7067 this->mapped_geometry->is_initialized(),
7073 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7074 "EvaluationFlags::hessians are supported."));
7080 unsigned int size = n_components * dim * n_q_points;
7083 for (
unsigned int i = 0; i < size; ++i)
7084 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7088 for (
unsigned int i = 0; i < size; ++i)
7089 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7094 if (n_components == dim &&
7095 this->data->element_type ==
7099 this->divergence_is_requested ==
false)
7101 unsigned int size = n_components * n_q_points;
7104 for (
unsigned int i = 0; i < size; ++i)
7105 this->values_quad[i] += this->values_from_gradients_quad[i];
7109 for (
unsigned int i = 0; i < size; ++i)
7110 this->values_quad[i] = this->values_from_gradients_quad[i];
7115 if constexpr (fe_degree > -1)
7118 template run<fe_degree, n_q_points_1d>(n_components,
7119 integration_flag_actual,
7122 sum_into_values_array);
7128 integration_flag_actual,
7131 sum_into_values_array);
7135 this->dof_values_initialized =
true;
7146 typename VectorizedArrayType>
7147template <
typename VectorType>
7154 VectorizedArrayType>::
7156 VectorType &destination)
7158 VectorizedArrayType *dst_ptr =
7159 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7160 *
this, destination);
7161 if (dst_ptr !=
nullptr)
7162 integrate(integration_flag, dst_ptr,
true);
7165 integrate(integration_flag, this->begin_dof_values());
7166 this->distribute_local_to_global(destination);
7177 typename VectorizedArrayType>
7184 VectorizedArrayType>::dof_indices()
const
7186 return {0U, dofs_per_cell};
7200 typename VectorizedArrayType>
7206 VectorizedArrayType>::
7209 const bool is_interior_face,
7210 const unsigned int dof_no,
7211 const unsigned int quad_no,
7212 const unsigned int first_selected_component,
7213 const unsigned int active_fe_index,
7214 const unsigned int active_quad_index,
7215 const unsigned int face_type)
7216 : BaseClass(matrix_free,
7218 first_selected_component,
7226 , dofs_per_component(this->data->dofs_per_component_on_cell)
7227 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7228 , n_q_points(this->n_quadrature_points)
7238 typename VectorizedArrayType>
7244 VectorizedArrayType>::
7247 const std::pair<unsigned int, unsigned int> &range,
7248 const bool is_interior_face,
7249 const unsigned int dof_no,
7250 const unsigned int quad_no,
7251 const unsigned int first_selected_component)
7256 first_selected_component,
7257 matrix_free.get_face_active_fe_index(range,
7259 numbers::invalid_unsigned_int,
7260 matrix_free.get_face_info(range.
first).face_type)
7270 typename VectorizedArrayType>
7277 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7279 Assert(this->mapped_geometry ==
nullptr,
7280 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7281 " Integer indexing is not possible"));
7282 if (this->mapped_geometry !=
nullptr)
7285 this->cell = face_index;
7286 this->dof_access_index =
7287 this->is_interior_face() ?
7295 this->matrix_free->get_task_info().boundary_partition_data.back())
7296 Assert(this->is_interior_face(),
7298 "Boundary faces do not have a neighbor. When looping over "
7299 "boundary faces use FEFaceEvaluation with the parameter "
7300 "is_interior_face set to true. "));
7302 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7307 this->face_ids[i] = face_index * n_lanes + i;
7308 for (; i < n_lanes; ++i)
7311 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7312 const unsigned int offsets =
7313 this->mapping_data->data_index_offsets[face_index];
7314 this->J_value = &this->mapping_data->JxW_values[offsets];
7315 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7317 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7318 this->normal_x_jacobian =
7320 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7321 this->jacobian_gradients =
7322 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7324 this->jacobian_gradients_non_inverse =
7326 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7330 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7333 this->mapping_data->quadrature_point_offsets.size());
7335 this->mapping_data->quadrature_points.data() +
7336 this->mapping_data->quadrature_point_offsets[this->cell];
7340 this->is_reinitialized =
true;
7341 this->dof_values_initialized =
false;
7342 this->values_quad_initialized =
false;
7343 this->gradients_quad_initialized =
false;
7344 this->hessians_quad_initialized =
false;
7355 typename VectorizedArrayType>
7363 const unsigned int face_number)
7369 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7373 Assert(this->mapped_geometry ==
nullptr,
7374 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7375 " Integer indexing is not possible"));
7376 if (this->mapped_geometry !=
nullptr)
7381 .faces_by_cells_type[
cell_index][face_number];
7384 this->dof_access_index =
7387 if (this->is_interior_face() ==
false)
7392 for (
unsigned int i = 0; i < n_lanes; ++i)
7395 const unsigned int cell_this =
cell_index * n_lanes + i;
7397 unsigned int face_index =
7402 this->face_ids[i] = face_index;
7407 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7408 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7415 auto cell_m = faces.cells_interior[face_index % n_lanes];
7416 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7418 const bool face_identifies_as_interior = cell_m != cell_this;
7420 Assert(cell_m == cell_this || cell_p == cell_this,
7424 if (face_identifies_as_interior)
7426 this->cell_ids[i] = cell_m;
7427 this->face_numbers[i] = faces.interior_face_no;
7431 this->cell_ids[i] = cell_p;
7432 this->face_numbers[i] = faces.exterior_face_no;
7435 const bool orientation_interior_face = faces.face_orientation >= 8;
7436 unsigned int face_orientation = faces.face_orientation % 8;
7437 if (face_identifies_as_interior != orientation_interior_face)
7439 constexpr std::array<std::uint8_t, 8> table{
7440 {0, 1, 2, 3, 6, 5, 4, 7}};
7441 face_orientation = table[face_orientation];
7443 this->face_orientations[i] = face_orientation;
7448 this->face_orientations[0] = 0;
7449 this->face_numbers[0] = face_number;
7454 for (
unsigned int i = 0; i < n_lanes; ++i)
7455 this->cell_ids[i] =
cell_index * n_lanes + i;
7463 this->cell_ids[i] =
cell_index * n_lanes + i;
7464 for (; i < n_lanes; ++i)
7467 for (
unsigned int i = 0; i < n_lanes; ++i)
7474 const unsigned int offsets =
7476 .face_data_by_cells[this->quad_no]
7481 .face_data_by_cells[this->quad_no]
7482 .JxW_values.size());
7484 .face_data_by_cells[this->quad_no]
7485 .JxW_values[offsets];
7487 .face_data_by_cells[this->quad_no]
7488 .normal_vectors[offsets];
7490 .face_data_by_cells[this->quad_no]
7491 .jacobians[!this->is_interior_face()][offsets];
7492 this->normal_x_jacobian =
7494 .face_data_by_cells[this->quad_no]
7495 .normals_times_jacobians[!this->is_interior_face()][offsets];
7496 this->jacobian_gradients =
7497 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7499 this->jacobian_gradients_non_inverse =
7501 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7506 .face_data_by_cells[this->quad_no]
7507 .quadrature_point_offsets.empty() ==
false)
7509 const unsigned int index =
7513 .face_data_by_cells[this->quad_no]
7514 .quadrature_point_offsets.size());
7516 .face_data_by_cells[this->quad_no]
7517 .quadrature_points.data() +
7519 .face_data_by_cells[this->quad_no]
7520 .quadrature_point_offsets[
index];
7524 this->is_reinitialized =
true;
7525 this->dof_values_initialized =
false;
7526 this->values_quad_initialized =
false;
7527 this->gradients_quad_initialized =
false;
7528 this->hessians_quad_initialized =
false;
7539 typename VectorizedArrayType>
7546 VectorizedArrayType>::
7553 evaluate(this->values_dofs, evaluation_flag);
7563 typename VectorizedArrayType>
7570 VectorizedArrayType>::
7571 evaluate(
const VectorizedArrayType *values_array,
7574 Assert((evaluation_flag &
7577 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7578 "and EvaluationFlags::hessians are supported."));
7580 const bool hessians_on_general_cells =
7584 if (hessians_on_general_cells)
7587 if (this->data->element_type ==
7593 if constexpr (fe_degree > -1)
7595 template
run<fe_degree, n_q_points_1d>(n_components,
7596 evaluation_flag_actual,
7601 n_components, evaluation_flag_actual, values_array, *
this);
7604 this->values_quad_initialized =
7606 this->gradients_quad_initialized =
7608 this->hessians_quad_initialized =
7620 typename VectorizedArrayType>
7627 VectorizedArrayType>::
7634 project_to_face(this->values_dofs, evaluation_flag);
7644 typename VectorizedArrayType>
7651 VectorizedArrayType>::
7652 project_to_face(
const VectorizedArrayType *values_array,
7655 Assert((evaluation_flag &
7658 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7659 "and EvaluationFlags::hessians are supported."));
7661 const bool hessians_on_general_cells =
7665 if (hessians_on_general_cells)
7668 if (this->data->element_type ==
7674 if constexpr (fe_degree > -1)
7677 VectorizedArrayType>::template run<fe_degree>(n_components,
7678 evaluation_flag_actual,
7684 evaluation_flag_actual,
7698 typename VectorizedArrayType>
7705 VectorizedArrayType>::
7708 Assert((evaluation_flag &
7711 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7712 "and EvaluationFlags::hessians are supported."));
7714 const bool hessians_on_general_cells =
7718 if (hessians_on_general_cells)
7721 if (this->data->element_type ==
7727 if constexpr (fe_degree > -1)
7730 VectorizedArrayType>::template run<fe_degree>(n_components,
7731 evaluation_flag_actual,
7738 this->values_quad_initialized =
7740 this->gradients_quad_initialized =
7742 this->hessians_quad_initialized =
7754 typename VectorizedArrayType>
7761 VectorizedArrayType>::
7763 const bool sum_into_values)
7765 integrate(integration_flag, this->values_dofs, sum_into_values);
7768 this->dof_values_initialized =
true;
7779 typename VectorizedArrayType>
7786 VectorizedArrayType>::
7788 VectorizedArrayType *values_array,
7789 const bool sum_into_values)
7791 Assert((integration_flag &
7794 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7795 "and EvaluationFlags::hessians are supported."));
7801 unsigned int size = n_components * dim * n_q_points;
7804 for (
unsigned int i = 0; i < size; ++i)
7805 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7809 for (
unsigned int i = 0; i < size; ++i)
7810 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7815 if (this->data->element_type ==
7819 this->divergence_is_requested ==
false)
7821 unsigned int size = n_components * n_q_points;
7824 for (
unsigned int i = 0; i < size; ++i)
7825 this->values_quad[i] += this->values_from_gradients_quad[i];
7829 for (
unsigned int i = 0; i < size; ++i)
7830 this->values_quad[i] = this->values_from_gradients_quad[i];
7835 if constexpr (fe_degree > -1)
7837 template
run<fe_degree, n_q_points_1d>(n_components,
7838 integration_flag_actual,
7845 integration_flag_actual,
7858 typename VectorizedArrayType>
7865 VectorizedArrayType>::
7868 Assert((integration_flag &
7871 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7872 "and EvaluationFlags::hessians are supported."));
7878 unsigned int size = n_components * dim * n_q_points;
7881 for (
unsigned int i = 0; i < size; ++i)
7882 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7886 for (
unsigned int i = 0; i < size; ++i)
7887 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7892 if (this->data->element_type ==
7896 this->divergence_is_requested ==
false)
7898 unsigned int size = n_components * n_q_points;
7901 for (
unsigned int i = 0; i < size; ++i)
7902 this->values_quad[i] += this->values_from_gradients_quad[i];
7906 for (
unsigned int i = 0; i < size; ++i)
7907 this->values_quad[i] = this->values_from_gradients_quad[i];
7912 if constexpr (fe_degree > -1)
7915 VectorizedArrayType>::template run<fe_degree>(n_components,
7916 integration_flag_actual,
7932 typename VectorizedArrayType>
7939 VectorizedArrayType>::
7941 const bool sum_into_values)
7943 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
7946 this->dof_values_initialized =
true;
7957 typename VectorizedArrayType>
7964 VectorizedArrayType>::
7966 VectorizedArrayType *values_array,
7967 const bool sum_into_values)
7969 Assert((integration_flag &
7972 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7973 "and EvaluationFlags::hessians are supported."));
7980 if (this->data->element_type ==
7984 this->divergence_is_requested ==
false)
7987 if constexpr (fe_degree > -1)
7990 VectorizedArrayType>::template run<fe_degree>(n_components,
7991 integration_flag_actual,
7998 integration_flag_actual,
8011 typename VectorizedArrayType>
8012template <
typename VectorType>
8019 VectorizedArrayType>::
8020 gather_evaluate(
const VectorType &input_vector,
8023 Assert((evaluation_flag &
8026 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8027 "and EvaluationFlags::hessians are supported."));
8029 const auto shared_vector_data = internal::get_shared_vector_data(
8031 this->dof_access_index ==
8033 this->active_fe_index,
8036 if (this->data->data.front().fe_degree > 0 &&
8037 fast_evaluation_supported(this->data->data.front().fe_degree,
8038 this->data->data.front().n_q_points_1d) &&
8041 typename VectorType::value_type,
8042 VectorizedArrayType>::
8043 supports(evaluation_flag,
8047 this->dof_info->index_storage_variants[this->dof_access_index]
8050 if constexpr (fe_degree > -1)
8054 typename VectorType::value_type,
8055 VectorizedArrayType>::template
run<fe_degree,
8059 internal::get_beginning<typename VectorType::value_type>(
8068 typename VectorType::value_type,
8069 VectorizedArrayType>::evaluate(n_components,
8071 internal::get_beginning<
8072 typename VectorType::value_type>(
8080 this->read_dof_values(input_vector);
8081 this->evaluate(evaluation_flag);
8086 this->gradients_quad_initialized =
8099 typename VectorizedArrayType>
8100template <
typename VectorType>
8108 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8109 const bool integrate_gradients,
8110 VectorType &destination)
8117 integrate_scatter(flag, destination);
8127 typename VectorizedArrayType>
8128template <
typename VectorType>
8135 VectorizedArrayType>::
8137 VectorType &destination)
8139 Assert((this->dof_access_index ==
8141 this->is_interior_face() ==
false) ==
false,
8144 const auto shared_vector_data = internal::get_shared_vector_data(
8146 this->dof_access_index ==
8148 this->active_fe_index,
8151 if (this->data->data.front().fe_degree > 0 &&
8152 fast_evaluation_supported(this->data->data.front().fe_degree,
8153 this->data->data.front().n_q_points_1d) &&
8156 typename VectorType::value_type,
8157 VectorizedArrayType>::
8158 supports(integration_flag,
8162 this->dof_info->index_storage_variants[this->dof_access_index]
8165 if constexpr (fe_degree > -1)
8169 typename VectorType::value_type,
8170 VectorizedArrayType>::template
run<fe_degree,
8174 internal::get_beginning<typename VectorType::value_type>(
8183 typename VectorType::value_type,
8184 VectorizedArrayType>::integrate(n_components,
8186 internal::get_beginning<
8187 typename VectorType::value_type>(
8195 integrate(integration_flag);
8196 this->distribute_local_to_global(destination);
8207 typename VectorizedArrayType>
8214 VectorizedArrayType>::dof_indices()
const
8216 return {0U, dofs_per_cell};
8226 typename VectorizedArrayType>
8233 VectorizedArrayType>::
8234 fast_evaluation_supported(
const unsigned int given_degree,
8235 const unsigned int given_n_q_points_1d)
8237 return fe_degree == -1 ?
8250 typename VectorizedArrayType>
8257 VectorizedArrayType>::
8258 fast_evaluation_supported(
const unsigned int given_degree,
8259 const unsigned int given_n_q_points_1d)
8261 return fe_degree == -1 ?
8274 typename VectorizedArrayType>
8281 VectorizedArrayType>::at_boundary()
const
8283 Assert(this->dof_access_index !=
8287 if (this->is_interior_face() ==
false)
8292 this->matrix_free->n_boundary_face_batches()))
8305 typename VectorizedArrayType>
8312 VectorizedArrayType>::boundary_id()
const
8314 Assert(this->dof_access_index !=
8331 typename VectorizedArrayType>
8339 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8341 return this->data->dofs_per_component_on_face;
8351 typename VectorizedArrayType>
8358 VectorizedArrayType>::get_dofs_projected_to_face()
8360 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
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
#define DEAL_II_NAMESPACE_CLOSE
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.
#define DEAL_II_NOT_IMPLEMENTED()
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)
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 check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const types::boundary_id internal_face_boundary_id
static const unsigned int invalid_unsigned_int
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)