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;
364 template <
int dim_ = dim,
365 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
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);
441 const unsigned int q_point);
495 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
513 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
516 const unsigned int q_point);
526 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
545 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
549 const unsigned int q_point);
559 template <
int dim_ = dim,
560 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
561 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
576 template <
int dim_ = dim,
577 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
580 const unsigned int q_point);
621 const unsigned int dof_no,
624 const unsigned int fe_degree,
625 const unsigned int n_q_points,
629 const unsigned int face_type);
703 template <
typename VectorType,
typename VectorOperation>
707 const std::array<VectorType *, n_components_> &vectors,
710 n_components_> &vectors_sm,
711 const std::bitset<n_lanes> &mask,
712 const bool apply_constraints =
true)
const;
721 template <
typename VectorType,
typename VectorOperation>
725 const std::array<VectorType *, n_components_> &vectors,
728 n_components_> &vectors_sm,
729 const std::bitset<n_lanes> &mask)
const;
738 template <
typename VectorType,
typename VectorOperation>
742 const std::array<VectorType *, n_components_> &vectors)
const;
1346 typename VectorizedArrayType>
1351 VectorizedArrayType>
1354 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1355 "Type of Number and of VectorizedArrayType do not match.");
1397 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1482 const unsigned int dof_no = 0,
1483 const unsigned int quad_no = 0,
1496 const std::pair<unsigned int, unsigned int> &range,
1497 const unsigned int dof_no = 0,
1498 const unsigned int quad_no = 0,
1608 template <
bool level_dof_access>
1630 const unsigned int given_n_q_points_1d);
1673 template <
typename VectorType>
1703 VectorizedArrayType *values_array,
1704 const bool sum_into_values =
false);
1719 template <
typename VectorType>
1722 VectorType &output_vector);
1806 int n_q_points_1d = fe_degree + 1,
1807 int n_components_ = 1,
1808 typename Number = double,
1814 VectorizedArrayType>
1817 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1818 "Type of Number and of VectorizedArrayType do not match.");
1860 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1962 const unsigned int dof_no = 0,
1963 const unsigned int quad_no = 0,
1978 const std::pair<unsigned int, unsigned int> &range,
1980 const unsigned int dof_no = 0,
1981 const unsigned int quad_no = 0,
1995 reinit(
const unsigned int face_batch_number);
2005 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2012 const unsigned int given_n_q_points_1d);
2076 template <
typename VectorType>
2092 const bool sum_into_values =
false);
2105 VectorizedArrayType *values_array,
2106 const bool sum_into_values =
false);
2123 const bool sum_into_values =
false);
2131 VectorizedArrayType *values_array,
2132 const bool sum_into_values =
false);
2145 template <
typename VectorType>
2148 VectorType &output_vector);
2153 template <
typename VectorType>
2156 const bool integrate_gradients,
2157 VectorType &output_vector);
2233 namespace MatrixFreeFunctions
2237 template <
int dim,
int degree>
2246 template <
int degree>
2249 static constexpr unsigned int value = degree + 1;
2265 template <
bool is_face,
2268 typename VectorizedArrayType>
2271 extract_initialization_data(
2273 const unsigned int dof_no,
2274 const unsigned int first_selected_component,
2275 const unsigned int quad_no,
2276 const unsigned int fe_degree,
2277 const unsigned int n_q_points,
2278 const unsigned int active_fe_index_given,
2279 const unsigned int active_quad_index_given,
2280 const unsigned int face_type)
2283 InitializationData init_data;
2287 &internal::MatrixFreeFunctions::
2288 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2296 active_fe_index_given :
2302 active_quad_index_given :
2303 std::min<unsigned int>(
2306 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2335 typename VectorizedArrayType>
2340 VectorizedArrayType>::
2343 const unsigned int dof_no,
2344 const unsigned int first_selected_component,
2345 const unsigned int quad_no,
2346 const unsigned int fe_degree,
2347 const unsigned int n_q_points,
2348 const bool is_interior_face,
2349 const unsigned int active_fe_index,
2350 const unsigned int active_quad_index,
2351 const unsigned int face_type)
2353 internal::extract_initialization_data<is_face>(matrix_free,
2355 first_selected_component,
2364 first_selected_component)
2365 , scratch_data_array(matrix_free.acquire_scratch_data())
2366 , matrix_free(&matrix_free)
2368 this->set_data_pointers(scratch_data_array, n_components_);
2370 this->dof_info->start_components.back() == 1 ||
2371 static_cast<int>(n_components_) <=
2373 this->dof_info->start_components
2374 [this->dof_info->component_to_base_index[first_selected_component] +
2376 first_selected_component,
2378 "You tried to construct a vector-valued evaluator with " +
2379 std::to_string(n_components) +
2380 " components. However, "
2381 "the current base element has only " +
2383 this->dof_info->start_components
2384 [this->dof_info->component_to_base_index[first_selected_component] +
2386 first_selected_component) +
2387 " components left when starting from local element index " +
2389 first_selected_component -
2390 this->dof_info->start_components
2391 [this->dof_info->component_to_base_index[first_selected_component]]) +
2392 " (global index " + std::to_string(first_selected_component) +
")"));
2404 typename VectorizedArrayType>
2409 VectorizedArrayType>::
2415 const unsigned int first_selected_component,
2419 other->mapped_geometry->get_quadrature() == quadrature ?
2420 other->mapped_geometry :
2422 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2427 first_selected_component)
2428 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2429 , matrix_free(nullptr)
2431 const unsigned int base_element_number =
2435 first_selected_component >=
2437 ExcMessage(
"The underlying element must at least contain as many "
2438 "components as requested by this class"));
2439 (void)base_element_number;
2443 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2445 fe.component_to_base_index(first_selected_component).
first);
2447 this->set_data_pointers(scratch_data_array, n_components_);
2456 typename VectorizedArrayType>
2461 VectorizedArrayType>::
2466 VectorizedArrayType> &other)
2468 , scratch_data_array(other.matrix_free == nullptr ?
2470 other.matrix_free->acquire_scratch_data())
2471 , matrix_free(other.matrix_free)
2473 if (other.matrix_free ==
nullptr)
2480 this->mapped_geometry =
2481 std::make_shared<internal::MatrixFreeFunctions::
2482 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2486 this->mapping_data = &this->mapped_geometry->get_data_storage();
2490 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2492 this->mapped_geometry->get_data_storage().JxW_values.begin();
2493 this->jacobian_gradients =
2494 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2495 this->jacobian_gradients_non_inverse =
2496 this->mapped_geometry->get_data_storage()
2497 .jacobian_gradients_non_inverse[0]
2500 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2503 this->set_data_pointers(scratch_data_array, n_components_);
2512 typename VectorizedArrayType>
2517 VectorizedArrayType> &
2523 VectorizedArrayType> &other)
2526 if (matrix_free ==
nullptr)
2529 delete scratch_data_array;
2538 matrix_free = other.matrix_free;
2540 if (other.matrix_free ==
nullptr)
2548 this->mapped_geometry =
2549 std::make_shared<internal::MatrixFreeFunctions::
2550 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2555 this->mapping_data = &this->mapped_geometry->get_data_storage();
2557 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2559 this->mapped_geometry->get_data_storage().JxW_values.begin();
2560 this->jacobian_gradients =
2561 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2562 this->jacobian_gradients_non_inverse =
2563 this->mapped_geometry->get_data_storage()
2564 .jacobian_gradients_non_inverse[0]
2567 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2574 this->set_data_pointers(scratch_data_array, n_components_);
2585 typename VectorizedArrayType>
2590 VectorizedArrayType>::~FEEvaluationBase()
2592 if (matrix_free !=
nullptr)
2603 delete scratch_data_array;
2614 typename VectorizedArrayType>
2619 Assert(matrix_free !=
nullptr,
2621 "FEEvaluation was not initialized with a MatrixFree object!"));
2622 return *matrix_free;
2631 template <
typename VectorType,
bool>
2632 struct ConstBlockVectorSelector;
2634 template <
typename VectorType>
2635 struct ConstBlockVectorSelector<
VectorType, true>
2637 using BaseVectorType =
const typename VectorType::BlockType;
2640 template <
typename VectorType>
2641 struct ConstBlockVectorSelector<
VectorType, false>
2643 using BaseVectorType =
typename VectorType::BlockType;
2649 template <
typename VectorType,
bool>
2650 struct BlockVectorSelector;
2652 template <
typename VectorType>
2655 using BaseVectorType =
typename ConstBlockVectorSelector<
2657 std::is_const_v<VectorType>>::BaseVectorType;
2659 static BaseVectorType *
2660 get_vector_component(VectorType &vec,
const unsigned int component)
2663 return &vec.block(component);
2667 template <
typename VectorType>
2668 struct BlockVectorSelector<
VectorType, false>
2672 static BaseVectorType *
2673 get_vector_component(VectorType &vec,
const unsigned int component)
2692 template <
typename VectorType>
2693 struct BlockVectorSelector<
std::vector<VectorType>, false>
2697 static BaseVectorType *
2698 get_vector_component(std::vector<VectorType> &vec,
2699 const unsigned int component)
2702 return &vec[component];
2706 template <
typename VectorType>
2707 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2711 static const BaseVectorType *
2712 get_vector_component(
const std::vector<VectorType> &vec,
2713 const unsigned int component)
2716 return &vec[component];
2720 template <
typename VectorType>
2721 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2725 static BaseVectorType *
2726 get_vector_component(std::vector<VectorType *> &vec,
2727 const unsigned int component)
2730 return vec[component];
2734 template <
typename VectorType>
2735 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2739 static const BaseVectorType *
2740 get_vector_component(
const std::vector<VectorType *> &vec,
2741 const unsigned int component)
2744 return vec[component];
2748 template <
typename VectorType, std::
size_t N>
2749 struct BlockVectorSelector<
std::array<VectorType *, N>, false>
2753 static BaseVectorType *
2754 get_vector_component(std::array<VectorType *, N> &vec,
2755 const unsigned int component)
2758 return vec[component];
2769 typename VectorizedArrayType>
2770template <
typename VectorType,
typename VectorOperation>
2775 const std::array<VectorType *, n_components_> &src,
2778 n_components_> &src_sm,
2779 const std::bitset<n_lanes> &
mask,
2780 const bool apply_constraints)
const
2785 if (this->matrix_free ==
nullptr)
2787 read_write_operation_global(operation, src);
2794 if (this->n_fe_components == 1)
2795 for (
unsigned int comp = 0; comp < n_components; ++comp)
2797 Assert(src[comp] !=
nullptr,
2798 ExcMessage(
"The finite element underlying this FEEvaluation "
2799 "object is scalar, but you requested " +
2800 std::to_string(n_components) +
2801 " components via the template argument in "
2802 "FEEvaluation. In that case, you must pass an "
2803 "std::vector<VectorType> or a BlockVector to " +
2804 "read_dof_values and distribute_local_to_global."));
2816 const bool accesses_exterior_dofs =
2817 this->dof_access_index ==
2819 this->is_interior_face() ==
false;
2829 bool is_contiguous =
true;
2831 if (accesses_exterior_dofs)
2833 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2834 const unsigned int n_filled_lanes =
2839 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2840 if (
mask[v] ==
true &&
2843 [cells[v] / n_lanes] <
2846 is_contiguous = false;
2850 this->dof_access_index :
2851 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2852 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2855 is_contiguous =
false;
2860 read_write_operation_contiguous(operation, src, src_sm,
mask);
2867 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2869 const bool masking_is_active =
mask.count() < n_lanes;
2870 if (masking_is_active)
2871 for (
unsigned int v = 0; v < n_lanes; ++v)
2872 if (
mask[v] ==
false)
2875 bool has_hn_constraints =
false;
2877 if (is_face ==
false)
2883 [this->first_selected_component])
2884 for (
unsigned int v = 0; v < n_lanes; ++v)
2889 has_hn_constraints = true;
2892 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2895 const bool use_vectorized_path =
2896 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2898 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2899 std::array<VectorizedArrayType *, n_components> values_dofs;
2900 for (
unsigned int c = 0; c < n_components; ++c)
2901 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2902 c * dofs_per_component;
2906 [is_face ? this->dof_access_index :
2907 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2908 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2909 IndexStorageVariants::interleaved &&
2910 use_vectorized_path)
2912 const unsigned int *dof_indices =
2914 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2918 [this->first_selected_component] *
2921 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2922 if (n_components == 1 || this->n_fe_components == 1)
2923 for (
unsigned int comp = 0; comp < n_components; ++comp)
2925 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2928 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2930 if (n_components == 1 || this->n_fe_components == 1)
2931 for (
unsigned int i = 0; i < dofs_per_component;
2932 ++i, dof_indices += n_lanes)
2933 for (
unsigned int comp = 0; comp < n_components; ++comp)
2934 operation.process_dof_gather(dof_indices,
2938 values_dofs[comp][i],
2941 for (
unsigned int comp = 0; comp < n_components; ++comp)
2942 for (
unsigned int i = 0; i < dofs_per_component;
2943 ++i, dof_indices += n_lanes)
2944 operation.process_dof_gather(dof_indices,
2948 values_dofs[comp][i],
2955 std::array<const unsigned int *, n_lanes> dof_indices;
2956 dof_indices.fill(
nullptr);
2961 bool has_constraints =
false;
2962 const unsigned int n_components_read =
2963 this->n_fe_components > 1 ? n_components : 1;
2967 for (
unsigned int v = 0; v < n_lanes; ++v)
2973 const std::pair<unsigned int, unsigned int> *my_index_start =
2974 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2975 this->first_selected_component];
2980 if (my_index_start[n_components_read].
second !=
2981 my_index_start[0].
second)
2982 has_constraints =
true;
2985 dof_info.
dof_indices.data() + my_index_start[0].first;
2990 for (
unsigned int v = 0; v < n_lanes; ++v)
2995 const std::pair<unsigned int, unsigned int> *my_index_start =
2996 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2997 this->first_selected_component];
2998 if (my_index_start[n_components_read].
second !=
2999 my_index_start[0].
second)
3000 has_constraints =
true;
3007 dof_info.hanging_node_constraint_masks_comp
3008 [this->active_fe_index][this->first_selected_component])
3009 has_hn_constraints = true;
3012 my_index_start[0].
first ||
3015 my_index_start[0].
first,
3018 dof_info.
dof_indices.data() + my_index_start[0].first;
3022 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3023 return i != numbers::invalid_unsigned_int;
3025 for (
unsigned int comp = 0; comp < n_components; ++comp)
3026 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3027 operation.process_empty(values_dofs[comp][i]);
3031 if (!has_constraints && apply_constraints)
3033 if (n_components == 1 || this->n_fe_components == 1)
3035 for (
unsigned int v = 0; v < n_lanes; ++v)
3040 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3041 for (
unsigned int comp = 0; comp < n_components; ++comp)
3042 operation.process_dof(dof_indices[v][i],
3044 values_dofs[comp][i][v]);
3049 for (
unsigned int comp = 0; comp < n_components; ++comp)
3050 for (
unsigned int v = 0; v < n_lanes; ++v)
3055 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3056 operation.process_dof(
3057 dof_indices[v][comp * dofs_per_component + i],
3059 values_dofs[comp][i][v]);
3071 for (
unsigned int v = 0; v < n_lanes; ++v)
3077 const unsigned int cell_dof_index =
3078 cell_index * this->n_fe_components + this->first_selected_component;
3079 const unsigned int n_components_read =
3080 this->n_fe_components > 1 ? n_components : 1;
3081 unsigned int index_indicators =
3083 unsigned int next_index_indicators =
3084 dof_info.
row_starts[cell_dof_index + 1].second;
3088 if (apply_constraints ==
false &&
3089 (dof_info.
row_starts[cell_dof_index].second !=
3090 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3096 dof_info.hanging_node_constraint_masks_comp
3097 [this->active_fe_index][this->first_selected_component])))
3106 [this->first_selected_component] +
3108 next_index_indicators = index_indicators;
3111 if (n_components == 1 || this->n_fe_components == 1)
3113 unsigned int ind_local = 0;
3114 for (; index_indicators != next_index_indicators; ++index_indicators)
3116 const std::pair<unsigned short, unsigned short> indicator =
3119 for (
unsigned int j = 0; j < indicator.first; ++j)
3120 for (
unsigned int comp = 0; comp < n_components; ++comp)
3121 operation.process_dof(dof_indices[v][j],
3123 values_dofs[comp][ind_local + j][v]);
3125 ind_local += indicator.first;
3126 dof_indices[v] += indicator.first;
3130 Number
value[n_components];
3131 for (
unsigned int comp = 0; comp < n_components; ++comp)
3132 operation.pre_constraints(values_dofs[comp][ind_local][v],
3135 const Number *data_val =
3137 const Number *end_pool =
3139 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3140 for (
unsigned int comp = 0; comp < n_components; ++comp)
3141 operation.process_constraint(*dof_indices[v],
3146 for (
unsigned int comp = 0; comp < n_components; ++comp)
3147 operation.post_constraints(
value[comp],
3148 values_dofs[comp][ind_local][v]);
3154 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3155 for (
unsigned int comp = 0; comp < n_components; ++comp)
3156 operation.process_dof(*dof_indices[v],
3158 values_dofs[comp][ind_local][v]);
3167 for (
unsigned int comp = 0; comp < n_components; ++comp)
3169 unsigned int ind_local = 0;
3172 for (; index_indicators != next_index_indicators;
3175 const std::pair<unsigned short, unsigned short> indicator =
3179 for (
unsigned int j = 0; j < indicator.first; ++j)
3180 operation.process_dof(dof_indices[v][j],
3182 values_dofs[comp][ind_local + j][v]);
3183 ind_local += indicator.first;
3184 dof_indices[v] += indicator.first;
3189 operation.pre_constraints(values_dofs[comp][ind_local][v],
3192 const Number *data_val =
3194 const Number *end_pool =
3197 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3198 operation.process_constraint(*dof_indices[v],
3203 operation.post_constraints(
value,
3204 values_dofs[comp][ind_local][v]);
3211 for (; ind_local < dofs_per_component;
3212 ++dof_indices[v], ++ind_local)
3215 operation.process_dof(*dof_indices[v],
3217 values_dofs[comp][ind_local][v]);
3220 if (apply_constraints ==
true && comp + 1 < n_components)
3221 next_index_indicators =
3222 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3234 typename VectorizedArrayType>
3235template <
typename VectorType,
typename VectorOperation>
3240 const std::array<VectorType *, n_components_> &src)
const
3244 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3245 unsigned int index = this->first_selected_component * dofs_per_component;
3246 for (
unsigned int comp = 0; comp < n_components; ++comp)
3248 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3250 operation.process_empty(
3251 this->values_dofs[comp * dofs_per_component + i]);
3252 operation.process_dof_global(
3253 local_dof_indices[this->data->lexicographic_numbering[
index]],
3255 this->values_dofs[comp * dofs_per_component + i][0]);
3266 typename VectorizedArrayType>
3267template <
typename VectorType,
typename VectorOperation>
3272 const std::array<VectorType *, n_components_> &src,
3275 n_components_> &vectors_sm,
3276 const std::bitset<n_lanes> &
mask)
const
3285 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3288 is_face ? this->dof_access_index :
3290 const unsigned int n_active_lanes =
mask.count();
3293 const std::vector<unsigned int> &dof_indices_cont =
3296 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3297 std::array<VectorizedArrayType *, n_components> values_dofs;
3298 for (
unsigned int c = 0; c < n_components; ++c)
3299 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3300 c * dofs_per_component;
3304 const bool accesses_exterior_dofs =
3305 this->dof_access_index ==
3307 this->is_interior_face() ==
false;
3313 interleaved_contiguous &&
3314 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3316 const unsigned int dof_index =
3317 dof_indices_cont[this->cell * n_lanes] +
3320 [this->first_selected_component] *
3322 if (n_components == 1 || this->n_fe_components == 1)
3323 for (
unsigned int comp = 0; comp < n_components; ++comp)
3324 operation.process_dofs_vectorized(dofs_per_component,
3330 operation.process_dofs_vectorized(dofs_per_component * n_components,
3338 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3342 const unsigned int n_filled_lanes =
3345 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3346 n_active_lanes == n_lanes &&
3347 !accesses_exterior_dofs;
3349 if (vectors_sm[0] !=
nullptr)
3351 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3352 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3354 const auto upper_bound =
3355 std::min<unsigned int>(n_filled_lanes, n_lanes);
3356 for (
unsigned int v = 0; v < upper_bound; ++v)
3358 if (
mask[v] ==
false)
3360 vector_ptrs[v] =
nullptr;
3380 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3381 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3383 [this->active_fe_index][this->first_selected_component]);
3385 vector_ptrs[v] =
nullptr;
3387 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3388 vector_ptrs[v] =
nullptr;
3393 if (use_vectorized_path)
3395 if (n_components == 1 || this->n_fe_components == 1)
3397 for (
unsigned int comp = 0; comp < n_components; ++comp)
3399 auto vector_ptrs = compute_vector_ptrs(comp);
3400 operation.process_dofs_vectorized_transpose(
3409 auto vector_ptrs = compute_vector_ptrs(0);
3410 operation.process_dofs_vectorized_transpose(dofs_per_component *
3418 for (
unsigned int comp = 0; comp < n_components; ++comp)
3420 auto vector_ptrs = compute_vector_ptrs(
3421 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3423 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3424 operation.process_empty(values_dofs[comp][i]);
3426 if (n_components == 1 || this->n_fe_components == 1)
3428 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3429 if (
mask[v] ==
true)
3430 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3431 operation.process_dof(vector_ptrs[v][i],
3432 values_dofs[comp][i][v]);
3436 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3437 if (
mask[v] ==
true)
3438 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3439 operation.process_dof(
3440 vector_ptrs[v][i + comp * dofs_per_component],
3441 values_dofs[comp][i][v]);
3447 std::array<unsigned int, n_lanes> dof_indices;
3448 std::fill(dof_indices.begin(),
3453 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3457 if (
mask[v] ==
true)
3459 dof_indices_cont[cells[v]] +
3462 [this->first_selected_component] *
3468 if (use_vectorized_path)
3474 if (n_components == 1 || this->n_fe_components == 1)
3475 for (
unsigned int comp = 0; comp < n_components; ++comp)
3476 operation.process_dofs_vectorized_transpose(dofs_per_component,
3482 operation.process_dofs_vectorized_transpose(dofs_per_component *
3491 interleaved_contiguous_strided)
3493 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3494 if (n_components == 1 || this->n_fe_components == 1)
3495 for (
unsigned int comp = 0; comp < n_components; ++comp)
3496 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3497 src[comp]->
begin());
3500 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3502 if (n_components == 1 || this->n_fe_components == 1)
3503 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3505 for (
unsigned int comp = 0; comp < n_components; ++comp)
3506 operation.process_dof_gather(dof_indices.data(),
3509 src_ptrs[comp] + i * n_lanes,
3510 values_dofs[comp][i],
3514 for (
unsigned int comp = 0; comp < n_components; ++comp)
3515 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3517 operation.process_dof_gather(
3520 (comp * dofs_per_component + i) * n_lanes,
3521 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3522 values_dofs[comp][i],
3530 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3532 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3533 if (n_components == 1 || this->n_fe_components == 1)
3534 for (
unsigned int comp = 0; comp < n_components; ++comp)
3535 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3536 src[comp]->
begin());
3539 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3541 const unsigned int *offsets =
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3546 for (
unsigned int comp = 0; comp < n_components; ++comp)
3547 operation.process_dof_gather(dof_indices.data(),
3551 values_dofs[comp][i],
3554 for (
unsigned int v = 0; v < n_lanes; ++v)
3555 dof_indices[v] += offsets[v];
3558 for (
unsigned int comp = 0; comp < n_components; ++comp)
3559 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3561 operation.process_dof_gather(dof_indices.data(),
3565 values_dofs[comp][i],
3568 for (
unsigned int v = 0; v < n_lanes; ++v)
3569 dof_indices[v] += offsets[v];
3574 for (
unsigned int comp = 0; comp < n_components; ++comp)
3576 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3577 operation.process_empty(values_dofs[comp][i]);
3578 if (accesses_exterior_dofs)
3580 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3581 if (
mask[v] ==
true)
3584 [ind][cells[v] / VectorizedArrayType::size()] ==
3588 if (n_components == 1 || this->n_fe_components == 1)
3590 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3591 operation.process_dof(dof_indices[v] + i,
3593 values_dofs[comp][i][v]);
3597 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3598 operation.process_dof(dof_indices[v] + i +
3599 comp * dofs_per_component,
3601 values_dofs[comp][i][v]);
3606 const unsigned int offset =
3609 if (n_components == 1 || this->n_fe_components == 1)
3611 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3612 operation.process_dof(dof_indices[v] + i * offset,
3614 values_dofs[comp][i][v]);
3618 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3619 operation.process_dof(
3621 (i + comp * dofs_per_component) * offset,
3623 values_dofs[comp][i][v]);
3634 if (n_components == 1 || this->n_fe_components == 1)
3636 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3637 if (
mask[v] ==
true)
3638 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3639 operation.process_dof(dof_indices[v] + i,
3641 values_dofs[comp][i][v]);
3645 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3646 if (
mask[v] ==
true)
3647 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3651 values_dofs[comp][i][v]);
3656 const unsigned int *offsets =
3658 [ind][VectorizedArrayType::size() * this->cell];
3659 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3661 if (n_components == 1 || this->n_fe_components == 1)
3662 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3664 if (
mask[v] ==
true)
3665 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3666 operation.process_dof(dof_indices[v] + i * offsets[v],
3668 values_dofs[comp][i][v]);
3672 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3673 if (
mask[v] ==
true)
3674 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3675 operation.process_dof(
3677 (i + comp * dofs_per_component) * offsets[v],
3679 values_dofs[comp][i][v]);
3691 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3692 decltype(std::declval<VectorType>().begin())
3693 get_beginning(VectorType &vec)
3701 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3702 typename VectorType::value_type *
3703 get_beginning(VectorType &)
3709 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3711 const std::vector<ArrayView<const typename VectorType::value_type>> *
3712 get_shared_vector_data(VectorType *vec,
3713 const bool is_valid_mode_for_sm,
3714 const unsigned int active_fe_index,
3718 if (is_valid_mode_for_sm &&
3721 active_fe_index == 0)
3722 return &vec->shared_vector_data();
3728 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3730 const std::vector<ArrayView<const typename VectorType::value_type>> *
3731 get_shared_vector_data(VectorType *,
3739 template <
int n_components,
typename VectorType>
3741 std::array<
typename internal::BlockVectorSelector<
3746 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3750 get_vector_data(VectorType &src,
3751 const unsigned int first_index,
3752 const bool is_valid_mode_for_sm,
3753 const unsigned int active_fe_index,
3759 std::array<
typename internal::BlockVectorSelector<
3765 ArrayView<
const typename internal::BlockVectorSelector<
3771 for (
unsigned int d = 0;
d < n_components; ++
d)
3772 src_data.first[d] = internal::BlockVectorSelector<
3778 for (
unsigned int d = 0;
d < n_components; ++
d)
3779 src_data.second[d] = get_shared_vector_data(
3780 const_cast<typename internal::BlockVectorSelector<
3781 std::remove_const_t<VectorType>,
3783 *>(src_data.first[d]),
3784 is_valid_mode_for_sm,
3798 typename VectorizedArrayType>
3803 if (this->dof_info ==
nullptr ||
3805 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3806 this->dof_info->hanging_node_constraint_masks_comp
3807 [this->active_fe_index][this->first_selected_component] ==
false)
3810 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3813 bool hn_available =
false;
3815 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3817 for (
unsigned int v = 0; v < n_lanes; ++v)
3829 constraint_mask[v] =
mask;
3835 if (hn_available ==
false)
3840 this->data->data.front().fe_degree,
3841 this->get_shape_info(),
3853 typename VectorizedArrayType>
3854template <
typename VectorType>
3858 const unsigned int first_index,
3859 const std::bitset<n_lanes> &
mask)
3861 const auto src_data = internal::get_vector_data<n_components_>(
3864 this->dof_access_index ==
3866 this->active_fe_index,
3870 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3872 apply_hanging_node_constraints(
false);
3875 this->dof_values_initialized =
true;
3885 typename VectorizedArrayType>
3886template <
typename VectorType>
3890 const unsigned int first_index,
3891 const std::bitset<n_lanes> &
mask)
3893 const auto src_data = internal::get_vector_data<n_components_>(
3896 this->dof_access_index ==
3898 this->active_fe_index,
3902 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3905 this->dof_values_initialized =
true;
3915 typename VectorizedArrayType>
3916template <
typename VectorType>
3920 const unsigned int first_index,
3921 const std::bitset<n_lanes> &
mask)
const
3924 Assert(this->dof_values_initialized ==
true,
3928 apply_hanging_node_constraints(
true);
3930 const auto dst_data = internal::get_vector_data<n_components_>(
3933 this->dof_access_index ==
3935 this->active_fe_index,
3940 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
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 const auto dst_data = internal::get_vector_data<n_components_>(
3965 this->dof_access_index ==
3967 this->active_fe_index,
3971 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
3980 typename VectorizedArrayType>
3981template <
typename VectorType>
3985 const unsigned int first_index,
3986 const std::bitset<n_lanes> &
mask)
const
3989 Assert(this->dof_values_initialized ==
true,
3993 const auto dst_data = internal::get_vector_data<n_components_>(
3996 this->dof_access_index ==
3998 this->active_fe_index,
4002 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4015 typename VectorizedArrayType>
4021 VectorizedArrayType>::value_type
4026 if constexpr (n_components == 1)
4027 return this->values_dofs[dof];
4030 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4032 for (
unsigned int comp = 0; comp < n_components; ++comp)
4033 return_value[comp] = this->values_dofs[comp * dofs + dof];
4034 return return_value;
4044 typename VectorizedArrayType>
4050 VectorizedArrayType>::value_type
4052 get_value(
const unsigned int q_point)
const
4055 Assert(this->values_quad_initialized ==
true,
4060 if constexpr (n_components == 1)
4061 return this->values_quad[q_point];
4064 if (n_components == dim &&
4065 this->data->element_type ==
4070 Assert(this->values_quad_initialized ==
true,
4075 Assert(this->J_value !=
nullptr,
4078 const std::size_t nqp = this->n_quadrature_points;
4086 const VectorizedArrayType inv_det =
4087 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4088 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4089 this->jacobian[0][2][2];
4092 for (
unsigned int comp = 0; comp < n_components; ++comp)
4093 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4094 jac[comp][comp] * inv_det;
4101 this->jacobian[q_point] :
4110 const VectorizedArrayType inv_det =
4111 (is_face && dim == 2 && this->get_face_no() < 2) ?
4115 for (
unsigned int comp = 0; comp < n_components; ++comp)
4117 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4118 for (
unsigned int e = 1;
e < dim; ++
e)
4120 this->values_quad[e * nqp + q_point] * jac[comp][e];
4121 value_out[comp] *= inv_det;
4128 const std::size_t nqp = this->n_quadrature_points;
4130 for (
unsigned int comp = 0; comp < n_components; ++comp)
4131 return_value[comp] = this->values_quad[comp * nqp + q_point];
4132 return return_value;
4143 typename VectorizedArrayType>
4149 VectorizedArrayType>::gradient_type
4154 Assert(this->gradients_quad_initialized ==
true,
4159 Assert(this->jacobian !=
nullptr,
4161 "update_gradients"));
4162 const std::size_t nqp = this->n_quadrature_points;
4164 if constexpr (n_components == dim && dim > 1)
4166 if (this->data->element_type ==
4171 Assert(this->gradients_quad_initialized ==
true,
4176 Assert(this->jacobian !=
nullptr,
4178 "update_gradients"));
4179 const std::size_t nqp = this->n_quadrature_points;
4180 const std::size_t nqp_d = nqp * dim;
4183 this->gradients_quad + q_point * dim;
4194 const VectorizedArrayType inv_det =
4195 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4196 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4197 this->jacobian[0][2][2];
4200 for (
unsigned int d = 0;
d < dim; ++
d)
4201 for (
unsigned int comp = 0; comp < n_components; ++comp)
4202 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4204 (jac[comp][comp] * inv_det);
4216 const VectorizedArrayType inv_det =
4217 (is_face && dim == 2 && this->get_face_no() < 2) ?
4221 VectorizedArrayType tmp[dim][dim];
4223 for (
unsigned int d = 0;
d < dim; ++
d)
4224 for (
unsigned int e = 0;
e < dim; ++
e)
4227 for (
unsigned int f = 1; f < dim; ++f)
4228 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4230 for (
unsigned int comp = 0; comp < n_components; ++comp)
4231 for (
unsigned int d = 0;
d < dim; ++
d)
4233 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4234 for (
unsigned int f = 1; f < dim; ++f)
4235 res += jac[comp][f] * tmp[d][f];
4237 grad_out[comp][
d] = res * inv_det;
4248 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4250 "update_hessians"));
4252 const auto jac_grad =
4253 this->jacobian_gradients_non_inverse[q_point];
4255 this->jacobian[q_point];
4259 const VectorizedArrayType inv_det =
4260 (is_face && dim == 2 && this->get_face_no() < 2) ?
4267 VectorizedArrayType tmp[dim][dim];
4268 for (
unsigned int d = 0;
d < dim; ++
d)
4269 for (
unsigned int e = 0;
e < dim; ++
e)
4272 for (
unsigned int f = 1; f < dim; ++f)
4273 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4278 for (
unsigned int d = 0;
d < dim; ++
d)
4280 for (
unsigned int e = 0;
e < dim; ++
e)
4282 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4283 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4284 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4287 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4289 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4294 for (
unsigned int d = 0;
d < dim; ++
d)
4295 for (
unsigned int e = 0;
e < dim; ++
e)
4297 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4298 for (
unsigned int f = 1; f < dim; ++f)
4299 res += tmp[f][d] * inv_t_jac[e][f];
4300 grad_out[
d][
e] = res;
4306 VectorizedArrayType tmp3[dim], tmp4[dim];
4307 for (
unsigned int d = 0;
d < dim; ++
d)
4309 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4310 for (
unsigned int e = 1;
e < dim; ++
e)
4311 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4313 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4314 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4315 for (
unsigned int d = 0;
d < dim; ++
d)
4317 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4318 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4320 for (
unsigned int d = 0;
d < dim; ++
d)
4322 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4323 for (
unsigned int e = 1;
e < dim; ++
e)
4324 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4327 VectorizedArrayType tmp2[dim];
4328 for (
unsigned int d = 0;
d < dim; ++
d)
4330 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4331 for (
unsigned e = 1;
e < dim; ++
e)
4333 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4336 for (
unsigned int d = 0;
d < dim; ++
d)
4337 for (
unsigned int e = 0;
e < dim; ++
e)
4339 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4343 grad_out[
d][
e] *= inv_det;
4354 for (
unsigned int comp = 0; comp < n_components; ++comp)
4355 for (
unsigned int d = 0;
d < dim; ++
d)
4357 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4358 this->jacobian[0][
d][
d];
4367 for (
unsigned int comp = 0; comp < n_components; ++comp)
4368 for (
unsigned int d = 0;
d < dim; ++
d)
4371 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4372 for (
unsigned int e = 1;
e < dim; ++
e)
4373 grad_out[comp][d] +=
4375 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4378 if constexpr (n_components == 1)
4390 typename VectorizedArrayType>
4396 VectorizedArrayType>::value_type
4402 Assert(this->gradients_quad_initialized ==
true,
4406 Assert(this->normal_x_jacobian !=
nullptr,
4408 "update_gradients"));
4410 const std::size_t nqp = this->n_quadrature_points;
4414 for (
unsigned int comp = 0; comp < n_components; ++comp)
4416 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4417 (this->normal_x_jacobian[0][dim - 1]);
4420 const std::size_t
index =
4422 for (
unsigned int comp = 0; comp < n_components; ++comp)
4424 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4425 this->normal_x_jacobian[
index][0];
4426 for (
unsigned int d = 1;
d < dim; ++
d)
4428 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4429 this->normal_x_jacobian[
index][
d];
4432 if constexpr (n_components == 1)
4444 template <
typename VectorizedArrayType>
4447 const VectorizedArrayType *
const hessians,
4449 VectorizedArrayType (&tmp)[1][1])
4451 tmp[0][0] = jac[0][0] *
hessians[0];
4454 template <
typename VectorizedArrayType>
4457 const VectorizedArrayType *
const hessians,
4458 const unsigned int nqp,
4459 VectorizedArrayType (&tmp)[2][2])
4461 for (
unsigned int d = 0;
d < 2; ++
d)
4469 template <
typename VectorizedArrayType>
4472 const VectorizedArrayType *
const hessians,
4473 const unsigned int nqp,
4474 VectorizedArrayType (&tmp)[3][3])
4476 for (
unsigned int d = 0;
d < 3; ++
d)
4497 typename VectorizedArrayType>
4502 VectorizedArrayType>::hessian_type
4507 Assert(this->hessians_quad_initialized ==
true,
4512 Assert(this->jacobian !=
nullptr,
4522 const std::size_t nqp = this->n_quadrature_points;
4523 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4528 for (
unsigned int comp = 0; comp < n_components; ++comp)
4530 for (
unsigned int d = 0;
d < dim; ++
d)
4531 hessian_out[comp][d][d] =
4532 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4533 (jac[
d][
d] * jac[
d][
d]);
4539 hessian_out[comp][0][1] =
4540 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4541 (jac[0][0] * jac[1][1]);
4544 hessian_out[comp][0][1] =
4545 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4546 (jac[0][0] * jac[1][1]);
4547 hessian_out[comp][0][2] =
4548 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4549 (jac[0][0] * jac[2][2]);
4550 hessian_out[comp][1][2] =
4551 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4552 (jac[1][1] * jac[2][2]);
4557 for (
unsigned int d = 0;
d < dim; ++
d)
4558 for (
unsigned int e = d + 1;
e < dim; ++
e)
4559 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4565 for (
unsigned int comp = 0; comp < n_components; ++comp)
4567 VectorizedArrayType tmp[dim][dim];
4568 internal::hessian_unit_times_jac(
4569 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4572 for (
unsigned int d = 0;
d < dim; ++
d)
4573 for (
unsigned int e = d;
e < dim; ++
e)
4575 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4576 for (
unsigned int f = 1; f < dim; ++f)
4577 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4584 for (
unsigned int d = 0;
d < dim; ++
d)
4585 for (
unsigned int e = d + 1;
e < dim; ++
e)
4586 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4592 const auto &jac_grad = this->jacobian_gradients[q_point];
4593 for (
unsigned int comp = 0; comp < n_components; ++comp)
4595 VectorizedArrayType tmp[dim][dim];
4596 internal::hessian_unit_times_jac(
4597 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4600 for (
unsigned int d = 0;
d < dim; ++
d)
4601 for (
unsigned int e = d;
e < dim; ++
e)
4603 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4604 for (
unsigned int f = 1; f < dim; ++f)
4605 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4609 for (
unsigned int d = 0;
d < dim; ++
d)
4610 for (
unsigned int e = 0;
e < dim; ++
e)
4611 hessian_out[comp][d][d] +=
4613 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4616 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4617 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4618 for (
unsigned int f = 0; f < dim; ++f)
4619 hessian_out[comp][d][e] +=
4620 jac_grad[count][f] *
4621 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4624 for (
unsigned int d = 0;
d < dim; ++
d)
4625 for (
unsigned int e = d + 1;
e < dim; ++
e)
4626 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4629 if constexpr (n_components == 1)
4630 return hessian_out[0];
4641 typename VectorizedArrayType>
4646 VectorizedArrayType>::gradient_type
4652 Assert(this->hessians_quad_initialized ==
true,
4663 const std::size_t nqp = this->n_quadrature_points;
4664 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4670 for (
unsigned int comp = 0; comp < n_components; ++comp)
4671 for (
unsigned int d = 0;
d < dim; ++
d)
4672 hessian_out[comp][d] =
4673 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4674 (jac[
d][
d] * jac[
d][
d]);
4679 for (
unsigned int comp = 0; comp < n_components; ++comp)
4683 VectorizedArrayType tmp[dim][dim];
4684 internal::hessian_unit_times_jac(
4685 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4689 for (
unsigned int d = 0;
d < dim; ++
d)
4691 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4692 for (
unsigned int f = 1; f < dim; ++f)
4693 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4700 const auto &jac_grad = this->jacobian_gradients[q_point];
4701 for (
unsigned int comp = 0; comp < n_components; ++comp)
4705 VectorizedArrayType tmp[dim][dim];
4706 internal::hessian_unit_times_jac(
4707 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4711 for (
unsigned int d = 0;
d < dim; ++
d)
4713 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4714 for (
unsigned int f = 1; f < dim; ++f)
4715 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4718 for (
unsigned int d = 0;
d < dim; ++
d)
4719 for (
unsigned int e = 0;
e < dim; ++
e)
4720 hessian_out[comp][d] +=
4722 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4726 if constexpr (n_components == 1)
4727 return hessian_out[0];
4738 typename VectorizedArrayType>
4743 VectorizedArrayType>::value_type
4749 Assert(this->hessians_quad_initialized ==
true,
4754 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4755 if constexpr (n_components == 1)
4757 VectorizedArrayType
sum = hess_diag[0];
4758 for (
unsigned int d = 1;
d < dim; ++
d)
4759 sum += hess_diag[d];
4765 for (
unsigned int comp = 0; comp < n_components; ++comp)
4767 laplacian_out[comp] = hess_diag[comp][0];
4768 for (
unsigned int d = 1;
d < dim; ++
d)
4769 laplacian_out[comp] += hess_diag[comp][d];
4771 return laplacian_out;
4781 typename VectorizedArrayType>
4787 this->dof_values_initialized =
true;
4789 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4791 for (
unsigned int comp = 0; comp < n_components; ++comp)
4792 if constexpr (n_components == 1)
4793 this->values_dofs[comp * dofs + dof] = val_in;
4795 this->values_dofs[comp * dofs + dof] = val_in[comp];
4804 typename VectorizedArrayType>
4807 submit_value(
const value_type val_in,
const unsigned int q_point)
4813 Assert(this->J_value !=
nullptr,
4817 this->values_quad_submitted =
true;
4820 const std::size_t nqp = this->n_quadrature_points;
4821 VectorizedArrayType *
values = this->values_quad + q_point;
4823 const VectorizedArrayType JxW =
4825 this->J_value[0] * this->quadrature_weights[q_point] :
4826 this->J_value[q_point];
4827 if constexpr (n_components == 1)
4828 values[0] = val_in * JxW;
4831 if (n_components == dim &&
4832 this->data->element_type ==
4837 Assert(this->J_value !=
nullptr,
4842 this->values_quad_submitted =
true;
4845 VectorizedArrayType *
values = this->values_quad + q_point;
4846 const std::size_t nqp = this->n_quadrature_points;
4852 const VectorizedArrayType weight =
4853 this->quadrature_weights[q_point];
4855 for (
unsigned int comp = 0; comp < n_components; ++comp)
4856 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4863 this->jacobian[q_point] :
4868 const VectorizedArrayType fac =
4870 this->quadrature_weights[q_point] :
4872 this->J_value[q_point] :
4873 this->J_value[0] * this->quadrature_weights[q_point]) *
4874 ((dim == 2 && this->get_face_no() < 2) ?
4883 for (
unsigned int comp = 0; comp < n_components; ++comp)
4885 values[comp * nqp] = val_in[0] * jac[0][comp];
4886 for (
unsigned int e = 1;
e < dim; ++
e)
4887 values[comp * nqp] += val_in[e] * jac[e][comp];
4888 values[comp * nqp] *= fac;
4893 for (
unsigned int comp = 0; comp < n_components; ++comp)
4894 values[comp * nqp] = val_in[comp] * JxW;
4904 typename VectorizedArrayType>
4905template <
int,
typename>
4909 const unsigned int q_point)
4911 static_assert(n_components == 1 && dim == 1,
4912 "Do not try to modify the default template parameters used for"
4913 " selectively enabling this function via std::enable_if!");
4914 submit_value(val_in[0], q_point);
4923 typename VectorizedArrayType>
4926 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
4932 Assert(this->J_value !=
nullptr,
4934 "update_gradients"));
4935 Assert(this->jacobian !=
nullptr,
4937 "update_gradients"));
4939 this->gradients_quad_submitted =
true;
4942 if constexpr (dim > 1 && n_components == dim)
4944 if (this->data->element_type ==
4953 Assert(this->J_value !=
nullptr,
4955 "update_gradients"));
4956 Assert(this->jacobian !=
nullptr,
4958 "update_gradients"));
4960 this->gradients_quad_submitted =
true;
4963 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
4964 VectorizedArrayType *
values =
4965 this->values_from_gradients_quad + q_point;
4966 const std::size_t nqp = this->n_quadrature_points;
4967 const std::size_t nqp_d = nqp * dim;
4977 const VectorizedArrayType weight =
4978 this->quadrature_weights[q_point];
4979 for (
unsigned int d = 0;
d < dim; ++
d)
4980 for (
unsigned int comp = 0; comp < n_components; ++comp)
4981 gradients[comp * nqp_d + d] = grad_in[comp][d] *
4983 (jac[comp][comp] * weight);
4995 const VectorizedArrayType fac =
4997 this->quadrature_weights[q_point] :
4998 this->J_value[0] * this->quadrature_weights[q_point] *
4999 ((dim == 2 && this->get_face_no() < 2) ?
5004 VectorizedArrayType tmp[dim][dim];
5005 for (
unsigned int d = 0;
d < dim; ++
d)
5006 for (
unsigned int e = 0;
e < dim; ++
e)
5008 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5009 for (
unsigned int f = 1; f < dim; ++f)
5010 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5012 for (
unsigned int comp = 0; comp < n_components; ++comp)
5013 for (
unsigned int d = 0;
d < dim; ++
d)
5015 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5016 for (
unsigned int f = 1; f < dim; ++f)
5017 res += jac[f][comp] * tmp[d][f];
5026 const auto jac_grad =
5027 this->jacobian_gradients_non_inverse[q_point];
5029 this->jacobian[q_point];
5033 const VectorizedArrayType fac =
5034 (!is_face) ? this->quadrature_weights[q_point] :
5035 this->J_value[q_point] *
5036 ((dim == 2 && this->get_face_no() < 2) ?
5045 VectorizedArrayType tmp3[dim], tmp4[dim];
5046 for (
unsigned int d = 0;
d < dim; ++
d)
5048 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5049 for (
unsigned int e = 1;
e < dim; ++
e)
5050 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5052 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5053 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5054 for (
unsigned int d = 0;
d < dim; ++
d)
5056 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5057 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5059 for (
unsigned int d = 0;
d < dim; ++
d)
5061 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5062 for (
unsigned int e = 1;
e < dim; ++
e)
5063 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5069 VectorizedArrayType tmp[dim][dim];
5072 for (
unsigned int d = 0;
d < dim; ++
d)
5073 for (
unsigned int e = 0;
e < dim; ++
e)
5075 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5076 for (
unsigned int f = 1; f < dim; ++f)
5077 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5080 for (
unsigned int d = 0;
d < dim; ++
d)
5081 for (
unsigned int e = 0;
e < dim; ++
e)
5083 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5084 for (
unsigned int f = 1; f < dim; ++f)
5085 res += t_jac[d][f] * tmp[e][f];
5092 VectorizedArrayType
value[dim];
5093 for (
unsigned int d = 0;
d < dim; ++
d)
5095 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5096 for (
unsigned int e = 1;
e < dim; ++
e)
5097 value[d] += tmp[d][e] * jac_grad[d][e];
5099 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5100 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5101 for (
unsigned int d = 0;
d < dim; ++
d)
5103 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5104 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5109 for (
unsigned int d = 0;
d < dim; ++
d)
5111 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5112 for (
unsigned int e = 1;
e < dim; ++
e)
5113 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5114 for (
unsigned int e = 0;
e < dim; ++
e)
5115 value[e] -= t_jac[e][d] * tmp2;
5118 for (
unsigned int d = 0;
d < dim; ++
d)
5119 values[d * nqp] =
value[d];
5125 const std::size_t nqp_d = this->n_quadrature_points * dim;
5126 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5130 const VectorizedArrayType JxW =
5131 this->J_value[0] * this->quadrature_weights[q_point];
5137 std::array<VectorizedArrayType, dim> jac;
5138 for (
unsigned int d = 0;
d < dim; ++
d)
5139 jac[d] = this->jacobian[0][d][d];
5141 for (
unsigned int d = 0;
d < dim; ++
d)
5143 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5144 if constexpr (n_components == 1)
5145 gradients[d] = grad_in[d] * factor;
5147 for (
unsigned int comp = 0; comp < n_components; ++comp)
5148 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5155 this->jacobian[q_point] :
5157 const VectorizedArrayType JxW =
5159 this->J_value[q_point] :
5160 this->J_value[0] * this->quadrature_weights[q_point];
5161 if constexpr (n_components == 1)
5162 for (
unsigned int d = 0;
d < dim; ++
d)
5164 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5165 for (
unsigned int e = 1;
e < dim; ++
e)
5166 new_val += (jac[e][d] * grad_in[e]);
5170 for (
unsigned int comp = 0; comp < n_components; ++comp)
5171 for (
unsigned int d = 0;
d < dim; ++
d)
5173 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5174 for (
unsigned int e = 1;
e < dim; ++
e)
5175 new_val += (jac[e][d] * grad_in[comp][e]);
5187 typename VectorizedArrayType>
5188template <
int,
typename>
5192 const unsigned int q_point)
5194 static_assert(n_components == 1 && dim == 1,
5195 "Do not try to modify the default template parameters used for"
5196 " selectively enabling this function via std::enable_if!");
5197 submit_gradient(grad_in[0], q_point);
5206 typename VectorizedArrayType>
5212 Assert(this->normal_x_jacobian !=
nullptr,
5214 "update_gradients"));
5216 this->gradients_quad_submitted =
true;
5219 const std::size_t nqp_d = this->n_quadrature_points * dim;
5220 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5224 const VectorizedArrayType JxW_jac = this->J_value[0] *
5225 this->quadrature_weights[q_point] *
5226 this->normal_x_jacobian[0][dim - 1];
5227 for (
unsigned int comp = 0; comp < n_components; ++comp)
5229 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5230 gradients[comp * nqp_d + d] = VectorizedArrayType();
5231 if constexpr (n_components == 1)
5232 gradients[dim - 1] = grad_in * JxW_jac;
5234 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5239 const unsigned int index =
5242 this->normal_x_jacobian[
index];
5243 const VectorizedArrayType JxW =
5245 this->J_value[
index] * this->quadrature_weights[q_point] :
5246 this->J_value[
index];
5247 for (
unsigned int comp = 0; comp < n_components; ++comp)
5248 for (
unsigned int d = 0;
d < dim; ++
d)
5249 if constexpr (n_components == 1)
5252 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5262 typename VectorizedArrayType>
5265 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5271 Assert(this->J_value !=
nullptr,
5273 "update_hessians"));
5274 Assert(this->jacobian !=
nullptr,
5276 "update_hessians"));
5278 this->hessians_quad_submitted =
true;
5282 const std::size_t nqp = this->n_quadrature_points;
5283 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5286 const VectorizedArrayType JxW =
5287 this->J_value[0] * this->quadrature_weights[q_point];
5290 for (
unsigned int d = 0;
d < dim; ++
d)
5292 const auto jac_d = this->jacobian[0][
d][
d];
5293 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5294 for (
unsigned int comp = 0; comp < n_components; ++comp)
5295 if constexpr (n_components == 1)
5296 this->hessians_quad[
d * nqp + q_point] =
5297 hessian_in[
d][
d] * factor;
5299 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5300 hessian_in[comp][d][d] * factor;
5304 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5305 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5307 const auto jac_d = this->jacobian[0][
d][
d];
5308 const auto jac_e = this->jacobian[0][
e][
e];
5309 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5310 for (
unsigned int comp = 0; comp < n_components; ++comp)
5311 if constexpr (n_components == 1)
5312 this->hessians_quad[off_dia * nqp + q_point] =
5313 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5315 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5316 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5323 const VectorizedArrayType JxW =
5324 this->J_value[0] * this->quadrature_weights[q_point];
5325 for (
unsigned int comp = 0; comp < n_components; ++comp)
5328 if constexpr (n_components == 1)
5329 hessian_c = hessian_in;
5331 hessian_c = hessian_in[comp];
5334 VectorizedArrayType tmp[dim][dim];
5335 for (
unsigned int i = 0; i < dim; ++i)
5336 for (
unsigned int j = 0; j < dim; ++j)
5338 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5339 for (
unsigned int k = 1; k < dim; ++k)
5340 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5344 VectorizedArrayType tmp2[dim][dim];
5345 for (
unsigned int i = 0; i < dim; ++i)
5346 for (
unsigned int j = 0; j < dim; ++j)
5348 tmp2[i][j] = jac[0][i] * tmp[0][j];
5349 for (
unsigned int k = 1; k < dim; ++k)
5350 tmp2[i][j] += jac[k][i] * tmp[k][j];
5354 for (
unsigned int d = 0;
d < dim; ++
d)
5355 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5359 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5360 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5361 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5362 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5368 const VectorizedArrayType JxW = this->J_value[q_point];
5369 const auto &jac_grad = this->jacobian_gradients[q_point];
5370 for (
unsigned int comp = 0; comp < n_components; ++comp)
5373 if constexpr (n_components == 1)
5374 hessian_c = hessian_in;
5376 hessian_c = hessian_in[comp];
5379 VectorizedArrayType tmp[dim][dim];
5380 for (
unsigned int i = 0; i < dim; ++i)
5381 for (
unsigned int j = 0; j < dim; ++j)
5383 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5384 for (
unsigned int k = 1; k < dim; ++k)
5385 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5389 VectorizedArrayType tmp2[dim][dim];
5390 for (
unsigned int i = 0; i < dim; ++i)
5391 for (
unsigned int j = 0; j < dim; ++j)
5393 tmp2[i][j] = jac[0][i] * tmp[0][j];
5394 for (
unsigned int k = 1; k < dim; ++k)
5395 tmp2[i][j] += jac[k][i] * tmp[k][j];
5399 for (
unsigned int d = 0;
d < dim; ++
d)
5400 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5404 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5405 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5406 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5407 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5410 for (
unsigned int d = 0;
d < dim; ++
d)
5412 VectorizedArrayType
sum = 0;
5413 for (
unsigned int e = 0;
e < dim; ++
e)
5414 sum += hessian_c[e][e] * jac_grad[e][d];
5415 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5416 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5418 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5419 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5432 typename VectorizedArrayType>
5437 VectorizedArrayType>::value_type
5443 Assert(this->values_quad_submitted ==
true,
5448 const std::size_t nqp = this->n_quadrature_points;
5449 for (
unsigned int q = 0; q < nqp; ++q)
5450 for (
unsigned int comp = 0; comp < n_components; ++comp)
5451 return_value[comp] += this->values_quad[comp * nqp + q];
5452 if constexpr (n_components == 1)
5453 return return_value[0];
5455 return return_value;
5464 typename VectorizedArrayType>
5465template <
int,
typename>
5470 static_assert(n_components == dim,
5471 "Do not try to modify the default template parameters used for"
5472 " selectively enabling this function via std::enable_if!");
5475 Assert(this->gradients_quad_initialized ==
true,
5479 Assert(this->jacobian !=
nullptr,
5481 "update_gradients"));
5483 VectorizedArrayType divergence;
5484 const std::size_t nqp = this->n_quadrature_points;
5487 this->data->element_type ==
5490 VectorizedArrayType inv_det =
5493 this->jacobian[0][0][0] *
5494 ((dim == 2) ? this->jacobian[0][1][1] :
5495 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5503 if (is_face && dim == 2 && this->get_face_no() < 2)
5507 divergence = this->gradients_quad[q_point * dim];
5508 for (
unsigned int d = 1;
d < dim; ++
d)
5509 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5510 divergence *= inv_det;
5519 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5520 for (
unsigned int d = 1;
d < dim; ++
d)
5521 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5522 this->jacobian[0][
d][
d];
5529 this->jacobian[q_point] :
5531 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5532 for (
unsigned int e = 1;
e < dim; ++
e)
5533 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5534 for (
unsigned int d = 1;
d < dim; ++
d)
5535 for (
unsigned int e = 0;
e < dim; ++
e)
5537 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5549 typename VectorizedArrayType>
5550template <
int,
typename>
5555 static_assert(n_components == dim,
5556 "Do not try to modify the default template parameters used for"
5557 " selectively enabling this function via std::enable_if!");
5560 const auto grad = get_gradient(q_point);
5561 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5562 VectorizedArrayType half = Number(0.5);
5563 for (
unsigned int d = 0;
d < dim; ++
d)
5564 symmetrized[d] = grad[d][d];
5570 symmetrized[2] = grad[0][1] + grad[1][0];
5571 symmetrized[2] *= half;
5574 symmetrized[3] = grad[0][1] + grad[1][0];
5575 symmetrized[3] *= half;
5576 symmetrized[4] = grad[0][2] + grad[2][0];
5577 symmetrized[4] *= half;
5578 symmetrized[5] = grad[1][2] + grad[2][1];
5579 symmetrized[5] *= half;
5593 typename VectorizedArrayType>
5594template <
int,
typename>
5596 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5598 get_curl(const unsigned
int q_point) const
5600 static_assert(dim > 1 && n_components == dim,
5601 "Do not try to modify the default template parameters used for"
5602 " selectively enabling this function via std::enable_if!");
5606 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5610 curl[0] = grad[1][0] - grad[0][1];
5613 curl[0] = grad[2][1] - grad[1][2];
5614 curl[1] = grad[0][2] - grad[2][0];
5615 curl[2] = grad[1][0] - grad[0][1];
5629 typename VectorizedArrayType>
5630template <
int,
typename>
5634 const unsigned int q_point)
5636 static_assert(n_components == dim,
5637 "Do not try to modify the default template parameters used for"
5638 " selectively enabling this function via std::enable_if!");
5644 Assert(this->J_value !=
nullptr,
5646 "update_gradients"));
5647 Assert(this->jacobian !=
nullptr,
5649 "update_gradients"));
5651 this->gradients_quad_submitted =
true;
5654 const std::size_t nqp_d = this->n_quadrature_points * dim;
5655 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5657 if (this->data->element_type ==
5664 const VectorizedArrayType fac =
5666 this->quadrature_weights[q_point] * div_in :
5668 this->J_value[q_point] :
5669 this->J_value[0] * this->quadrature_weights[q_point]) *
5672 this->jacobian[this->cell_type >
5676 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5678 for (
unsigned int d = 0;
d < dim; ++
d)
5680 for (
unsigned int e = 0;
e < dim; ++
e)
5681 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5683 this->divergence_is_requested =
true;
5690 const VectorizedArrayType fac =
5691 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5692 for (
unsigned int d = 0;
d < dim; ++
d)
5694 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5695 for (
unsigned int e = 0;
e < dim; ++
e)
5696 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5703 this->jacobian[q_point] :
5705 const VectorizedArrayType fac =
5707 this->J_value[q_point] :
5708 this->J_value[0] * this->quadrature_weights[q_point]) *
5710 for (
unsigned int d = 0;
d < dim; ++
d)
5712 for (
unsigned int e = 0;
e < dim; ++
e)
5713 gradients[d * nqp_d + e] = jac[d][e] * fac;
5725 typename VectorizedArrayType>
5726template <
int,
typename>
5731 const unsigned int q_point)
5733 static_assert(n_components == dim,
5734 "Do not try to modify the default template parameters used for"
5735 " selectively enabling this function via std::enable_if!");
5738 this->data->element_type !=
5749 Assert(this->J_value !=
nullptr,
5751 "update_gradients"));
5752 Assert(this->jacobian !=
nullptr,
5754 "update_gradients"));
5756 this->gradients_quad_submitted =
true;
5759 const std::size_t nqp_d = this->n_quadrature_points * dim;
5760 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5763 const VectorizedArrayType JxW =
5764 this->J_value[0] * this->quadrature_weights[q_point];
5766 for (
unsigned int d = 0;
d < dim; ++
d)
5767 gradients[d * nqp_d + d] =
5769 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5770 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5772 const VectorizedArrayType
value =
5781 const VectorizedArrayType JxW =
5783 this->J_value[q_point] :
5784 this->J_value[0] * this->quadrature_weights[q_point];
5787 this->jacobian[q_point] :
5789 VectorizedArrayType weighted[dim][dim];
5790 for (
unsigned int i = 0; i < dim; ++i)
5792 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
5793 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
5795 const VectorizedArrayType
value =
5797 weighted[i][j] =
value;
5798 weighted[j][i] =
value;
5800 for (
unsigned int comp = 0; comp < dim; ++comp)
5801 for (
unsigned int d = 0;
d < dim; ++
d)
5803 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
5804 for (
unsigned int e = 1;
e < dim; ++
e)
5805 new_val += jac[e][d] * weighted[comp][e];
5817 typename VectorizedArrayType>
5818template <
int,
typename>
5822 const unsigned int q_point)
5824 static_assert(n_components == dim,
5825 "Do not try to modify the default template parameters used for"
5826 " selectively enabling this function via std::enable_if!");
5832 grad[1][0] = curl[0];
5833 grad[0][1] = -curl[0];
5836 grad[2][1] = curl[0];
5837 grad[1][2] = -curl[0];
5838 grad[0][2] = curl[1];
5839 grad[2][0] = -curl[1];
5840 grad[1][0] = curl[2];
5841 grad[0][1] = -curl[2];
5846 submit_gradient(grad, q_point);
5859 typename VectorizedArrayType>
5865 VectorizedArrayType>::
5867 const unsigned int fe_no,
5868 const unsigned int quad_no,
5869 const unsigned int first_selected_component,
5870 const unsigned int active_fe_index,
5871 const unsigned int active_quad_index)
5872 : BaseClass(matrix_free,
5874 first_selected_component,
5881 numbers::invalid_unsigned_int )
5882 , dofs_per_component(this->data->dofs_per_component_on_cell)
5883 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5884 , n_q_points(this->data->n_q_points)
5886 check_template_arguments(fe_no, 0);
5896 typename VectorizedArrayType>
5902 VectorizedArrayType>::
5904 const std::pair<unsigned int, unsigned int> &range,
5905 const unsigned int dof_no,
5906 const unsigned int quad_no,
5907 const unsigned int first_selected_component)
5911 first_selected_component,
5912 matrix_free.get_cell_active_fe_index(range))
5922 typename VectorizedArrayType>
5928 VectorizedArrayType>::
5933 const unsigned int first_selected_component)
5934 : BaseClass(mapping,
5938 first_selected_component,
5940 , dofs_per_component(this->data->dofs_per_component_on_cell)
5941 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5942 , n_q_points(this->data->n_q_points)
5954 typename VectorizedArrayType>
5960 VectorizedArrayType>::
5964 const unsigned int first_selected_component)
5969 first_selected_component,
5971 , dofs_per_component(this->data->dofs_per_component_on_cell)
5972 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5973 , n_q_points(this->data->n_q_points)
5985 typename VectorizedArrayType>
5991 VectorizedArrayType>::
5994 const unsigned int first_selected_component)
5995 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
5997 other.mapped_geometry->get_quadrature(),
5998 other.mapped_geometry->get_fe_values().get_update_flags(),
5999 first_selected_component,
6001 , dofs_per_component(this->data->dofs_per_component_on_cell)
6002 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6003 , n_q_points(this->data->n_q_points)
6015 typename VectorizedArrayType>
6024 , dofs_per_component(this->data->dofs_per_component_on_cell)
6025 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6026 , n_q_points(this->data->n_q_points)
6038 typename VectorizedArrayType>
6044 VectorizedArrayType> &
6050 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6052 BaseClass::operator=(other);
6064 typename VectorizedArrayType>
6071 VectorizedArrayType>::
6072 check_template_arguments(
const unsigned int dof_no,
6073 const unsigned int first_selected_component)
6076 (void)first_selected_component;
6079 this->data->dofs_per_component_on_cell > 0,
6081 "There is nothing useful you can do with an FEEvaluation object with "
6082 "FE_Nothing, i.e., without DoFs! If you have passed to "
6083 "MatrixFree::reinit() a collection of finite elements also containing "
6084 "FE_Nothing, please check - before creating FEEvaluation - the category "
6085 "of the current range by calling either "
6086 "MatrixFree::get_cell_range_category(range) or "
6087 "MatrixFree::get_face_range_category(range). The returned category "
6088 "is the index of the active FE, which you can use to exclude "
6095 static_cast<unsigned int>(fe_degree) !=
6096 this->data->data.front().fe_degree) ||
6097 n_q_points != this->n_quadrature_points)
6099 std::string message =
6100 "-------------------------------------------------------\n";
6101 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6102 message +=
" Called --> FEEvaluation<dim,";
6106 message +=
",Number>(data";
6122 if (
static_cast<unsigned int>(fe_degree) ==
6123 this->data->data.front().fe_degree)
6125 proposed_dof_comp = dof_no;
6126 proposed_fe_comp = first_selected_component;
6129 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6131 for (
unsigned int nf = 0;
6134 if (this->matrix_free
6135 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6137 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6139 proposed_dof_comp = no;
6140 proposed_fe_comp = nf;
6144 this->mapping_data->descriptor[this->active_quad_index]
6146 proposed_quad_comp = this->quad_no;
6148 for (
unsigned int no = 0;
6153 .descriptor[this->active_quad_index]
6154 .n_q_points == n_q_points)
6156 proposed_quad_comp = no;
6163 if (proposed_dof_comp != first_selected_component)
6164 message +=
"Wrong vector component selection:\n";
6166 message +=
"Wrong quadrature formula selection:\n";
6167 message +=
" Did you mean FEEvaluation<dim,";
6171 message +=
",Number>(data";
6180 std::string correct_pos;
6181 if (proposed_dof_comp != dof_no)
6182 correct_pos =
" ^ ";
6185 if (proposed_quad_comp != this->quad_no)
6186 correct_pos +=
" ^ ";
6189 if (proposed_fe_comp != first_selected_component)
6190 correct_pos +=
" ^\n";
6192 correct_pos +=
" \n";
6198 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6199 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6200 message +=
"Wrong template arguments:\n";
6201 message +=
" Did you mean FEEvaluation<dim,";
6206 message +=
",Number>(data";
6214 std::string correct_pos;
6215 if (this->data->data.front().fe_degree !=
6216 static_cast<unsigned int>(fe_degree))
6220 if (proposed_n_q_points_1d != n_q_points_1d)
6221 correct_pos +=
" ^\n";
6223 correct_pos +=
" \n";
6224 message +=
" " + correct_pos;
6226 Assert(
static_cast<unsigned int>(fe_degree) ==
6227 this->data->data.front().fe_degree &&
6228 n_q_points == this->n_quadrature_points,
6234 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6245 typename VectorizedArrayType>
6254 Assert(this->matrix_free !=
nullptr,
6255 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6256 " Integer indexing is not possible."));
6264 const unsigned int offsets =
6265 this->mapping_data->data_index_offsets[
cell_index];
6266 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6267 this->J_value = &this->mapping_data->JxW_values[offsets];
6268 if (!this->mapping_data->jacobian_gradients[0].empty())
6270 this->jacobian_gradients =
6271 this->mapping_data->jacobian_gradients[0].data() + offsets;
6272 this->jacobian_gradients_non_inverse =
6273 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6279 for (
unsigned int i = 0; i < n_lanes; ++i)
6280 this->cell_ids[i] =
cell_index * n_lanes + i;
6287 this->cell_ids[i] =
cell_index * n_lanes + i;
6288 for (; i < n_lanes; ++i)
6292 if (this->mapping_data->quadrature_points.empty() ==
false)
6294 &this->mapping_data->quadrature_points
6295 [this->mapping_data->quadrature_point_offsets[this->cell]];
6298 this->is_reinitialized =
true;
6299 this->dof_values_initialized =
false;
6300 this->values_quad_initialized =
false;
6301 this->gradients_quad_initialized =
false;
6302 this->hessians_quad_initialized =
false;
6313 typename VectorizedArrayType>
6320 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6327 this->cell_ids = cell_ids;
6332 for (
unsigned int v = 0; v < n_lanes; ++v)
6346 if (this->mapped_geometry ==
nullptr)
6347 this->mapped_geometry =
6348 std::make_shared<internal::MatrixFreeFunctions::
6349 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6351 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6353 auto &this_jacobian_data = mapping_storage.jacobians[0];
6354 auto &this_J_value_data = mapping_storage.JxW_values;
6355 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6356 auto &this_jacobian_gradients_non_inverse_data =
6357 mapping_storage.jacobian_gradients_non_inverse[0];
6358 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6362 if (this_jacobian_data.size() != 2)
6363 this_jacobian_data.resize_fast(2);
6365 if (this_J_value_data.size() != 1)
6366 this_J_value_data.resize_fast(1);
6368 const auto &update_flags_cells =
6372 this_jacobian_gradients_data.size() != 1)
6374 this_jacobian_gradients_data.resize_fast(1);
6375 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6379 this_quadrature_points_data.size() != 1)
6380 this_quadrature_points_data.resize_fast(1);
6384 if (this_jacobian_data.size() != this->n_quadrature_points)
6385 this_jacobian_data.resize_fast(this->n_quadrature_points);
6387 if (this_J_value_data.size() != this->n_quadrature_points)
6388 this_J_value_data.resize_fast(this->n_quadrature_points);
6390 const auto &update_flags_cells =
6394 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6396 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6397 this_jacobian_gradients_non_inverse_data.resize_fast(
6398 this->n_quadrature_points);
6402 this_quadrature_points_data.size() != this->n_quadrature_points)
6403 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6407 this->jacobian = this_jacobian_data.data();
6408 this->J_value = this_J_value_data.data();
6409 this->jacobian_gradients = this_jacobian_gradients_data.data();
6410 this->jacobian_gradients_non_inverse =
6411 this_jacobian_gradients_non_inverse_data.data();
6415 for (
unsigned int v = 0; v < n_lanes; ++v)
6422 const unsigned int cell_batch_index =
cell_index / n_lanes;
6423 const unsigned int offsets =
6424 this->mapping_data->data_index_offsets[cell_batch_index];
6425 const unsigned int lane =
cell_index % n_lanes;
6427 if (this->cell_type <=
6431 for (
unsigned int q = 0; q < 2; ++q)
6432 for (
unsigned int i = 0; i < dim; ++i)
6433 for (
unsigned int j = 0; j < dim; ++j)
6434 this_jacobian_data[q][i][j][v] =
6435 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6437 const unsigned int q = 0;
6439 this_J_value_data[q][v] =
6440 this->mapping_data->JxW_values[offsets + q][lane];
6442 const auto &update_flags_cells =
6447 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6448 for (
unsigned int j = 0; j < dim; ++j)
6449 this_jacobian_gradients_data[q][i][j][v] =
6451 ->jacobian_gradients[0][offsets + q][i][j][lane];
6453 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6454 for (
unsigned int j = 0; j < dim; ++j)
6455 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6457 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6462 for (
unsigned int i = 0; i < dim; ++i)
6463 this_quadrature_points_data[q][i][v] =
6464 this->mapping_data->quadrature_points
6466 ->quadrature_point_offsets[cell_batch_index] +
6472 const auto cell_type =
6476 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6478 const unsigned int q_src =
6484 this_J_value_data[q][v] =
6485 this->mapping_data->JxW_values[offsets + q_src][lane];
6487 for (
unsigned int i = 0; i < dim; ++i)
6488 for (
unsigned int j = 0; j < dim; ++j)
6489 this_jacobian_data[q][i][j][v] =
6491 ->jacobians[0][offsets + q_src][i][j][lane];
6493 const auto &update_flags_cells =
6498 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6499 for (
unsigned int j = 0; j < dim; ++j)
6500 this_jacobian_gradients_data[q][i][j][v] =
6502 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6504 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6505 for (
unsigned int j = 0; j < dim; ++j)
6506 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6508 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6521 this->mapping_data->quadrature_points
6523 ->quadrature_point_offsets[cell_batch_index] +
6527 this->mapping_data->jacobians[0][offsets + 1];
6529 for (
unsigned int d = 0;
d < dim; ++
d)
6532 static_cast<Number
>(
6533 this->descriptor->quadrature.point(q)[
d]);
6535 for (
unsigned int d = 0;
d < dim; ++
d)
6536 for (
unsigned int e = 0;
e < dim; ++
e)
6539 static_cast<Number
>(
6540 this->descriptor->quadrature.point(q)[
e]);
6542 for (
unsigned int i = 0; i < dim; ++i)
6543 this_quadrature_points_data[q][i][v] = point[i][lane];
6548 for (
unsigned int i = 0; i < dim; ++i)
6549 this_quadrature_points_data[q][i][v] =
6550 this->mapping_data->quadrature_points
6552 ->quadrature_point_offsets[cell_batch_index] +
6561 this->is_reinitialized =
true;
6562 this->dof_values_initialized =
false;
6563 this->values_quad_initialized =
false;
6564 this->gradients_quad_initialized =
false;
6565 this->hessians_quad_initialized =
false;
6576 typename VectorizedArrayType>
6577template <
bool level_dof_access>
6584 VectorizedArrayType>
::
6587 Assert(this->matrix_free ==
nullptr,
6588 ExcMessage(
"Cannot use initialization from cell iterator if "
6589 "initialized from MatrixFree object. Use variant for "
6590 "on the fly computation with arguments as for FEValues "
6593 this->mapped_geometry->reinit(
6595 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6596 if (level_dof_access)
6597 cell->get_mg_dof_indices(this->local_dof_indices);
6599 cell->get_dof_indices(this->local_dof_indices);
6602 this->is_reinitialized =
true;
6613 typename VectorizedArrayType>
6620 VectorizedArrayType>
::
6623 Assert(this->matrix_free == 0,
6624 ExcMessage(
"Cannot use initialization from cell iterator if "
6625 "initialized from MatrixFree object. Use variant for "
6626 "on the fly computation with arguments as for FEValues "
6629 this->mapped_geometry->reinit(cell);
6632 this->is_reinitialized =
true;
6643 typename VectorizedArrayType>
6650 VectorizedArrayType>::
6654 Assert(this->dof_values_initialized ==
true,
6657 evaluate(this->values_dofs, evaluation_flags);
6667 typename VectorizedArrayType>
6674 VectorizedArrayType>::
6675 evaluate(
const VectorizedArrayType *values_array,
6678 const bool hessians_on_general_cells =
6682 if (hessians_on_general_cells)
6685 if (this->data->element_type ==
6691 if constexpr (fe_degree > -1)
6694 template run<fe_degree, n_q_points_1d>(n_components,
6695 evaluation_flag_actual,
6703 evaluation_flag_actual,
6704 const_cast<VectorizedArrayType *
>(values_array),
6709 this->values_quad_initialized =
6711 this->gradients_quad_initialized =
6713 this->hessians_quad_initialized =
6724 template <
typename Number,
6725 typename VectorizedArrayType,
6727 typename EvaluatorType,
6728 std::enable_if_t<internal::has_begin<VectorType> &&
6731 VectorizedArrayType *
6732 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6738 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6739 const auto &dof_info = fe_eval.get_dof_info();
6746 if (std::is_same_v<typename VectorType::value_type, Number> &&
6750 interleaved_contiguous &&
6751 reinterpret_cast<
std::size_t>(
6753 dof_info.dof_indices_contiguous
6754 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6755 [cell * VectorizedArrayType::size()]) %
6756 sizeof(VectorizedArrayType) ==
6759 return reinterpret_cast<VectorizedArrayType *
>(
6763 [cell * VectorizedArrayType::size()] +
6765 [fe_eval.get_active_fe_index()]
6766 [fe_eval.get_first_selected_component()] *
6767 VectorizedArrayType::size());
6776 template <
typename Number,
6777 typename VectorizedArrayType,
6779 typename EvaluatorType,
6780 std::enable_if_t<!internal::has_begin<VectorType> ||
6783 VectorizedArrayType *
6784 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
6797 typename VectorizedArrayType>
6798template <
typename VectorType>
6805 VectorizedArrayType>::
6806 gather_evaluate(
const VectorType &input_vector,
6809 const VectorizedArrayType *src_ptr =
6810 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
6811 *
this, input_vector);
6812 if (src_ptr !=
nullptr)
6813 evaluate(src_ptr, evaluation_flag);
6816 this->read_dof_values(input_vector);
6817 evaluate(this->begin_dof_values(), evaluation_flag);
6828 typename VectorizedArrayType>
6835 VectorizedArrayType>::
6838 integrate(integration_flag, this->values_dofs);
6841 this->dof_values_initialized =
true;
6852 typename VectorizedArrayType>
6859 VectorizedArrayType>::
6861 VectorizedArrayType *values_array,
6862 const bool sum_into_values_array)
6866 Assert(this->values_quad_submitted ==
true,
6869 Assert(this->gradients_quad_submitted ==
true,
6872 Assert(this->hessians_quad_submitted ==
true,
6875 Assert(this->matrix_free !=
nullptr ||
6876 this->mapped_geometry->is_initialized(),
6882 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
6883 "EvaluationFlags::hessians are supported."));
6889 unsigned int size = n_components * dim * n_q_points;
6892 for (
unsigned int i = 0; i < size; ++i)
6893 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
6897 for (
unsigned int i = 0; i < size; ++i)
6898 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
6903 if (n_components == dim &&
6904 this->data->element_type ==
6908 this->divergence_is_requested ==
false)
6910 unsigned int size = n_components * n_q_points;
6913 for (
unsigned int i = 0; i < size; ++i)
6914 this->values_quad[i] += this->values_from_gradients_quad[i];
6918 for (
unsigned int i = 0; i < size; ++i)
6919 this->values_quad[i] = this->values_from_gradients_quad[i];
6924 if constexpr (fe_degree > -1)
6927 template run<fe_degree, n_q_points_1d>(n_components,
6928 integration_flag_actual,
6931 sum_into_values_array);
6937 integration_flag_actual,
6940 sum_into_values_array);
6944 this->dof_values_initialized =
true;
6955 typename VectorizedArrayType>
6956template <
typename VectorType>
6963 VectorizedArrayType>::
6965 VectorType &destination)
6967 VectorizedArrayType *dst_ptr =
6968 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
6969 *
this, destination);
6970 if (dst_ptr !=
nullptr)
6971 integrate(integration_flag, dst_ptr,
true);
6974 integrate(integration_flag, this->begin_dof_values());
6975 this->distribute_local_to_global(destination);
6986 typename VectorizedArrayType>
6993 VectorizedArrayType>::dof_indices()
const
6995 return {0U, dofs_per_cell};
7009 typename VectorizedArrayType>
7015 VectorizedArrayType>::
7018 const bool is_interior_face,
7019 const unsigned int dof_no,
7020 const unsigned int quad_no,
7021 const unsigned int first_selected_component,
7022 const unsigned int active_fe_index,
7023 const unsigned int active_quad_index,
7024 const unsigned int face_type)
7025 : BaseClass(matrix_free,
7027 first_selected_component,
7035 , dofs_per_component(this->data->dofs_per_component_on_cell)
7036 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7037 , n_q_points(this->n_quadrature_points)
7047 typename VectorizedArrayType>
7053 VectorizedArrayType>::
7056 const std::pair<unsigned int, unsigned int> &range,
7057 const bool is_interior_face,
7058 const unsigned int dof_no,
7059 const unsigned int quad_no,
7060 const unsigned int first_selected_component)
7065 first_selected_component,
7066 matrix_free.get_face_active_fe_index(range,
7068 numbers::invalid_unsigned_int,
7069 matrix_free.get_face_info(range.
first).face_type)
7079 typename VectorizedArrayType>
7086 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7088 Assert(this->mapped_geometry ==
nullptr,
7089 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7090 " Integer indexing is not possible"));
7091 if (this->mapped_geometry !=
nullptr)
7094 this->cell = face_index;
7095 this->dof_access_index =
7096 this->is_interior_face() ?
7104 this->matrix_free->get_task_info().boundary_partition_data.back())
7105 Assert(this->is_interior_face(),
7107 "Boundary faces do not have a neighbor. When looping over "
7108 "boundary faces use FEFaceEvaluation with the parameter "
7109 "is_interior_face set to true. "));
7111 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7116 this->face_ids[i] = face_index * n_lanes + i;
7117 for (; i < n_lanes; ++i)
7120 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7121 const unsigned int offsets =
7122 this->mapping_data->data_index_offsets[face_index];
7123 this->J_value = &this->mapping_data->JxW_values[offsets];
7124 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7126 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7127 this->normal_x_jacobian =
7129 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7130 this->jacobian_gradients =
7131 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7133 this->jacobian_gradients_non_inverse =
7135 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7139 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7142 this->mapping_data->quadrature_point_offsets.size());
7144 this->mapping_data->quadrature_points.data() +
7145 this->mapping_data->quadrature_point_offsets[this->cell];
7149 this->is_reinitialized =
true;
7150 this->dof_values_initialized =
false;
7151 this->values_quad_initialized =
false;
7152 this->gradients_quad_initialized =
false;
7153 this->hessians_quad_initialized =
false;
7164 typename VectorizedArrayType>
7172 const unsigned int face_number)
7178 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7182 Assert(this->mapped_geometry ==
nullptr,
7183 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7184 " Integer indexing is not possible"));
7185 if (this->mapped_geometry !=
nullptr)
7190 .faces_by_cells_type[
cell_index][face_number];
7193 this->dof_access_index =
7196 if (this->is_interior_face() ==
false)
7201 for (
unsigned int i = 0; i < n_lanes; ++i)
7204 const unsigned int cell_this =
cell_index * n_lanes + i;
7206 unsigned int face_index =
7211 this->face_ids[i] = face_index;
7216 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7217 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7224 auto cell_m = faces.cells_interior[face_index % n_lanes];
7225 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7227 const bool face_identifies_as_interior = cell_m != cell_this;
7229 Assert(cell_m == cell_this || cell_p == cell_this,
7233 if (face_identifies_as_interior)
7235 this->cell_ids[i] = cell_m;
7236 this->face_numbers[i] = faces.interior_face_no;
7240 this->cell_ids[i] = cell_p;
7241 this->face_numbers[i] = faces.exterior_face_no;
7244 const bool orientation_interior_face = faces.face_orientation >= 8;
7245 unsigned int face_orientation = faces.face_orientation % 8;
7246 if (face_identifies_as_interior != orientation_interior_face)
7248 constexpr std::array<std::uint8_t, 8> table{
7249 {0, 1, 2, 3, 6, 5, 4, 7}};
7250 face_orientation = table[face_orientation];
7252 this->face_orientations[i] = face_orientation;
7257 this->face_orientations[0] = 0;
7258 this->face_numbers[0] = face_number;
7263 for (
unsigned int i = 0; i < n_lanes; ++i)
7264 this->cell_ids[i] =
cell_index * n_lanes + i;
7272 this->cell_ids[i] =
cell_index * n_lanes + i;
7273 for (; i < n_lanes; ++i)
7276 for (
unsigned int i = 0; i < n_lanes; ++i)
7283 const unsigned int offsets =
7285 .face_data_by_cells[this->quad_no]
7290 .face_data_by_cells[this->quad_no]
7291 .JxW_values.size());
7293 .face_data_by_cells[this->quad_no]
7294 .JxW_values[offsets];
7296 .face_data_by_cells[this->quad_no]
7297 .normal_vectors[offsets];
7299 .face_data_by_cells[this->quad_no]
7300 .jacobians[!this->is_interior_face()][offsets];
7301 this->normal_x_jacobian =
7303 .face_data_by_cells[this->quad_no]
7304 .normals_times_jacobians[!this->is_interior_face()][offsets];
7305 this->jacobian_gradients =
7306 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7308 this->jacobian_gradients_non_inverse =
7310 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7315 .face_data_by_cells[this->quad_no]
7316 .quadrature_point_offsets.empty() ==
false)
7318 const unsigned int index =
7322 .face_data_by_cells[this->quad_no]
7323 .quadrature_point_offsets.size());
7325 .face_data_by_cells[this->quad_no]
7326 .quadrature_points.data() +
7328 .face_data_by_cells[this->quad_no]
7329 .quadrature_point_offsets[
index];
7333 this->is_reinitialized =
true;
7334 this->dof_values_initialized =
false;
7335 this->values_quad_initialized =
false;
7336 this->gradients_quad_initialized =
false;
7337 this->hessians_quad_initialized =
false;
7348 typename VectorizedArrayType>
7355 VectorizedArrayType>::
7362 evaluate(this->values_dofs, evaluation_flag);
7372 typename VectorizedArrayType>
7379 VectorizedArrayType>::
7380 evaluate(
const VectorizedArrayType *values_array,
7383 Assert((evaluation_flag &
7386 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7387 "and EvaluationFlags::hessians are supported."));
7389 const bool hessians_on_general_cells =
7393 if (hessians_on_general_cells)
7396 if (this->data->element_type ==
7402 if constexpr (fe_degree > -1)
7404 template
run<fe_degree, n_q_points_1d>(n_components,
7405 evaluation_flag_actual,
7410 n_components, evaluation_flag_actual, values_array, *
this);
7413 this->values_quad_initialized =
7415 this->gradients_quad_initialized =
7417 this->hessians_quad_initialized =
7429 typename VectorizedArrayType>
7436 VectorizedArrayType>::
7443 project_to_face(this->values_dofs, evaluation_flag);
7453 typename VectorizedArrayType>
7460 VectorizedArrayType>::
7461 project_to_face(
const VectorizedArrayType *values_array,
7464 Assert((evaluation_flag &
7467 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7468 "and EvaluationFlags::hessians are supported."));
7470 const bool hessians_on_general_cells =
7474 if (hessians_on_general_cells)
7477 if (this->data->element_type ==
7483 if constexpr (fe_degree > -1)
7486 VectorizedArrayType>::template run<fe_degree>(n_components,
7487 evaluation_flag_actual,
7493 evaluation_flag_actual,
7507 typename VectorizedArrayType>
7514 VectorizedArrayType>::
7517 Assert((evaluation_flag &
7520 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7521 "and EvaluationFlags::hessians are supported."));
7523 const bool hessians_on_general_cells =
7527 if (hessians_on_general_cells)
7530 if (this->data->element_type ==
7536 if constexpr (fe_degree > -1)
7539 VectorizedArrayType>::template run<fe_degree>(n_components,
7540 evaluation_flag_actual,
7547 this->values_quad_initialized =
7549 this->gradients_quad_initialized =
7551 this->hessians_quad_initialized =
7563 typename VectorizedArrayType>
7570 VectorizedArrayType>::
7572 const bool sum_into_values)
7574 integrate(integration_flag, this->values_dofs, sum_into_values);
7577 this->dof_values_initialized =
true;
7588 typename VectorizedArrayType>
7595 VectorizedArrayType>::
7597 VectorizedArrayType *values_array,
7598 const bool sum_into_values)
7600 Assert((integration_flag &
7603 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7604 "and EvaluationFlags::hessians are supported."));
7610 unsigned int size = n_components * dim * n_q_points;
7613 for (
unsigned int i = 0; i < size; ++i)
7614 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7618 for (
unsigned int i = 0; i < size; ++i)
7619 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7624 if (this->data->element_type ==
7628 this->divergence_is_requested ==
false)
7630 unsigned int size = n_components * n_q_points;
7633 for (
unsigned int i = 0; i < size; ++i)
7634 this->values_quad[i] += this->values_from_gradients_quad[i];
7638 for (
unsigned int i = 0; i < size; ++i)
7639 this->values_quad[i] = this->values_from_gradients_quad[i];
7644 if constexpr (fe_degree > -1)
7646 template
run<fe_degree, n_q_points_1d>(n_components,
7647 integration_flag_actual,
7654 integration_flag_actual,
7667 typename VectorizedArrayType>
7674 VectorizedArrayType>::
7677 Assert((integration_flag &
7680 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7681 "and EvaluationFlags::hessians are supported."));
7687 unsigned int size = n_components * dim * n_q_points;
7690 for (
unsigned int i = 0; i < size; ++i)
7691 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7695 for (
unsigned int i = 0; i < size; ++i)
7696 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7701 if (this->data->element_type ==
7705 this->divergence_is_requested ==
false)
7707 unsigned int size = n_components * n_q_points;
7710 for (
unsigned int i = 0; i < size; ++i)
7711 this->values_quad[i] += this->values_from_gradients_quad[i];
7715 for (
unsigned int i = 0; i < size; ++i)
7716 this->values_quad[i] = this->values_from_gradients_quad[i];
7721 if constexpr (fe_degree > -1)
7724 VectorizedArrayType>::template run<fe_degree>(n_components,
7725 integration_flag_actual,
7741 typename VectorizedArrayType>
7748 VectorizedArrayType>::
7750 const bool sum_into_values)
7752 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
7755 this->dof_values_initialized =
true;
7766 typename VectorizedArrayType>
7773 VectorizedArrayType>::
7775 VectorizedArrayType *values_array,
7776 const bool sum_into_values)
7778 Assert((integration_flag &
7781 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7782 "and EvaluationFlags::hessians are supported."));
7789 if (this->data->element_type ==
7793 this->divergence_is_requested ==
false)
7796 if constexpr (fe_degree > -1)
7799 VectorizedArrayType>::template run<fe_degree>(n_components,
7800 integration_flag_actual,
7807 integration_flag_actual,
7820 typename VectorizedArrayType>
7821template <
typename VectorType>
7828 VectorizedArrayType>::
7829 gather_evaluate(
const VectorType &input_vector,
7832 Assert((evaluation_flag &
7835 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7836 "and EvaluationFlags::hessians are supported."));
7838 const auto shared_vector_data = internal::get_shared_vector_data(
7840 this->dof_access_index ==
7842 this->active_fe_index,
7845 if (this->data->data.front().fe_degree > 0 &&
7846 fast_evaluation_supported(this->data->data.front().fe_degree,
7847 this->data->data.front().n_q_points_1d) &&
7850 typename VectorType::value_type,
7851 VectorizedArrayType>::
7852 supports(evaluation_flag,
7856 this->dof_info->index_storage_variants[this->dof_access_index]
7859 if constexpr (fe_degree > -1)
7863 typename VectorType::value_type,
7864 VectorizedArrayType>::template
run<fe_degree,
7868 internal::get_beginning<typename VectorType::value_type>(
7877 typename VectorType::value_type,
7878 VectorizedArrayType>::evaluate(n_components,
7880 internal::get_beginning<
7881 typename VectorType::value_type>(
7889 this->read_dof_values(input_vector);
7890 this->evaluate(evaluation_flag);
7895 this->gradients_quad_initialized =
7908 typename VectorizedArrayType>
7909template <
typename VectorType>
7917 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
7918 const bool integrate_gradients,
7919 VectorType &destination)
7926 integrate_scatter(flag, destination);
7936 typename VectorizedArrayType>
7937template <
typename VectorType>
7944 VectorizedArrayType>::
7946 VectorType &destination)
7948 Assert((this->dof_access_index ==
7950 this->is_interior_face() ==
false) ==
false,
7953 const auto shared_vector_data = internal::get_shared_vector_data(
7955 this->dof_access_index ==
7957 this->active_fe_index,
7960 if (this->data->data.front().fe_degree > 0 &&
7961 fast_evaluation_supported(this->data->data.front().fe_degree,
7962 this->data->data.front().n_q_points_1d) &&
7965 typename VectorType::value_type,
7966 VectorizedArrayType>::
7967 supports(integration_flag,
7971 this->dof_info->index_storage_variants[this->dof_access_index]
7974 if constexpr (fe_degree > -1)
7978 typename VectorType::value_type,
7979 VectorizedArrayType>::template
run<fe_degree,
7983 internal::get_beginning<typename VectorType::value_type>(
7992 typename VectorType::value_type,
7993 VectorizedArrayType>::integrate(n_components,
7995 internal::get_beginning<
7996 typename VectorType::value_type>(
8004 integrate(integration_flag);
8005 this->distribute_local_to_global(destination);
8016 typename VectorizedArrayType>
8023 VectorizedArrayType>::dof_indices()
const
8025 return {0U, dofs_per_cell};
8035 typename VectorizedArrayType>
8042 VectorizedArrayType>::
8043 fast_evaluation_supported(
const unsigned int given_degree,
8044 const unsigned int given_n_q_points_1d)
8046 return fe_degree == -1 ?
8059 typename VectorizedArrayType>
8066 VectorizedArrayType>::
8067 fast_evaluation_supported(
const unsigned int given_degree,
8068 const unsigned int given_n_q_points_1d)
8070 return fe_degree == -1 ?
8083 typename VectorizedArrayType>
8090 VectorizedArrayType>::at_boundary()
const
8092 Assert(this->dof_access_index !=
8096 if (this->is_interior_face() ==
false)
8101 this->matrix_free->n_boundary_face_batches()))
8114 typename VectorizedArrayType>
8121 VectorizedArrayType>::boundary_id()
const
8123 Assert(this->dof_access_index !=
8140 typename VectorizedArrayType>
8148 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8150 return this->data->dofs_per_component_on_face;
8160 typename VectorizedArrayType>
8167 VectorizedArrayType>::get_dofs_projected_to_face()
8169 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)
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
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 > &)