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>
1704 VectorizedArrayType *values_array,
1705 const bool sum_into_values =
false);
1720 template <
typename VectorType>
1723 VectorType &output_vector);
1807 int n_q_points_1d = fe_degree + 1,
1808 int n_components_ = 1,
1809 typename Number = double,
1815 VectorizedArrayType>
1818 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1819 "Type of Number and of VectorizedArrayType do not match.");
1861 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1963 const unsigned int dof_no = 0,
1964 const unsigned int quad_no = 0,
1979 const std::pair<unsigned int, unsigned int> &range,
1981 const unsigned int dof_no = 0,
1982 const unsigned int quad_no = 0,
1996 reinit(
const unsigned int face_batch_number);
2006 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2013 const unsigned int given_n_q_points_1d);
2077 template <
typename VectorType>
2093 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;
2444 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2446 fe.component_to_base_index(first_selected_component).
first);
2448 this->set_data_pointers(scratch_data_array, n_components_);
2457 typename VectorizedArrayType>
2462 VectorizedArrayType>::
2467 VectorizedArrayType> &other)
2469 , scratch_data_array(other.matrix_free == nullptr ?
2471 other.matrix_free->acquire_scratch_data())
2472 , matrix_free(other.matrix_free)
2474 if (other.matrix_free ==
nullptr)
2482 this->mapped_geometry =
2483 std::make_shared<internal::MatrixFreeFunctions::
2484 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2488 this->mapping_data = &this->mapped_geometry->get_data_storage();
2492 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2494 this->mapped_geometry->get_data_storage().JxW_values.begin();
2495 this->jacobian_gradients =
2496 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2497 this->jacobian_gradients_non_inverse =
2498 this->mapped_geometry->get_data_storage()
2499 .jacobian_gradients_non_inverse[0]
2502 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2505 this->set_data_pointers(scratch_data_array, n_components_);
2514 typename VectorizedArrayType>
2519 VectorizedArrayType> &
2525 VectorizedArrayType> &other)
2528 if (matrix_free ==
nullptr)
2531 delete scratch_data_array;
2540 matrix_free = other.matrix_free;
2542 if (other.matrix_free ==
nullptr)
2551 this->mapped_geometry =
2552 std::make_shared<internal::MatrixFreeFunctions::
2553 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2558 this->mapping_data = &this->mapped_geometry->get_data_storage();
2560 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2562 this->mapped_geometry->get_data_storage().JxW_values.begin();
2563 this->jacobian_gradients =
2564 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2565 this->jacobian_gradients_non_inverse =
2566 this->mapped_geometry->get_data_storage()
2567 .jacobian_gradients_non_inverse[0]
2570 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2577 this->set_data_pointers(scratch_data_array, n_components_);
2588 typename VectorizedArrayType>
2593 VectorizedArrayType>::~FEEvaluationBase()
2595 if (matrix_free !=
nullptr)
2606 delete scratch_data_array;
2617 typename VectorizedArrayType>
2622 Assert(matrix_free !=
nullptr,
2624 "FEEvaluation was not initialized with a MatrixFree object!"));
2625 return *matrix_free;
2634 template <
typename VectorType,
bool>
2635 struct ConstBlockVectorSelector;
2637 template <
typename VectorType>
2638 struct ConstBlockVectorSelector<VectorType, true>
2640 using BaseVectorType =
const typename VectorType::BlockType;
2643 template <
typename VectorType>
2644 struct ConstBlockVectorSelector<VectorType, false>
2646 using BaseVectorType =
typename VectorType::BlockType;
2652 template <
typename VectorType,
bool>
2653 struct BlockVectorSelector;
2655 template <
typename VectorType>
2656 struct BlockVectorSelector<VectorType, true>
2658 using BaseVectorType =
typename ConstBlockVectorSelector<
2660 std::is_const_v<VectorType>>::BaseVectorType;
2662 static BaseVectorType *
2663 get_vector_component(VectorType &vec,
const unsigned int component)
2666 return &vec.block(component);
2670 template <
typename VectorType>
2671 struct BlockVectorSelector<VectorType, false>
2673 using BaseVectorType = VectorType;
2675 static BaseVectorType *
2676 get_vector_component(VectorType &vec,
const unsigned int component)
2695 template <
typename VectorType>
2696 struct BlockVectorSelector<
std::vector<VectorType>, false>
2698 using BaseVectorType = VectorType;
2700 static BaseVectorType *
2701 get_vector_component(std::vector<VectorType> &vec,
2702 const unsigned int component)
2705 return &vec[component];
2709 template <
typename VectorType>
2710 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2712 using BaseVectorType =
const VectorType;
2714 static const BaseVectorType *
2715 get_vector_component(
const std::vector<VectorType> &vec,
2716 const unsigned int component)
2719 return &vec[component];
2723 template <
typename VectorType>
2724 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2726 using BaseVectorType = VectorType;
2728 static BaseVectorType *
2729 get_vector_component(std::vector<VectorType *> &vec,
2730 const unsigned int component)
2733 return vec[component];
2737 template <
typename VectorType>
2738 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2740 using BaseVectorType =
const VectorType;
2742 static const BaseVectorType *
2743 get_vector_component(
const std::vector<VectorType *> &vec,
2744 const unsigned int component)
2747 return vec[component];
2758 typename VectorizedArrayType>
2759template <
typename VectorType,
typename VectorOperation>
2764 const std::array<VectorType *, n_components_> &src,
2767 n_components_> &src_sm,
2768 const std::bitset<n_lanes> &
mask,
2769 const bool apply_constraints)
const
2774 if (this->matrix_free ==
nullptr)
2776 read_write_operation_global(operation, src);
2783 if (this->n_fe_components == 1)
2784 for (
unsigned int comp = 0; comp < n_components; ++comp)
2786 Assert(src[comp] !=
nullptr,
2787 ExcMessage(
"The finite element underlying this FEEvaluation "
2788 "object is scalar, but you requested " +
2789 std::to_string(n_components) +
2790 " components via the template argument in "
2791 "FEEvaluation. In that case, you must pass an "
2792 "std::vector<VectorType> or a BlockVector to " +
2793 "read_dof_values and distribute_local_to_global."));
2805 const bool accesses_exterior_dofs =
2806 this->dof_access_index ==
2808 this->is_interior_face() ==
false;
2818 bool is_contiguous =
true;
2820 if (accesses_exterior_dofs)
2822 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2823 const unsigned int n_filled_lanes =
2828 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2829 if (
mask[v] ==
true &&
2832 [cells[v] / n_lanes] <
2835 is_contiguous = false;
2839 this->dof_access_index :
2840 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2841 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2844 is_contiguous =
false;
2849 read_write_operation_contiguous(operation, src, src_sm,
mask);
2856 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2858 const bool masking_is_active =
mask.count() < n_lanes;
2859 if (masking_is_active)
2860 for (
unsigned int v = 0; v < n_lanes; ++v)
2861 if (
mask[v] ==
false)
2864 bool has_hn_constraints =
false;
2866 if (is_face ==
false)
2872 [this->first_selected_component])
2873 for (
unsigned int v = 0; v < n_lanes; ++v)
2878 has_hn_constraints = true;
2881 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2884 const bool use_vectorized_path =
2885 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2887 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2888 std::array<VectorizedArrayType *, n_components> values_dofs;
2889 for (
unsigned int c = 0; c < n_components; ++c)
2890 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2891 c * dofs_per_component;
2895 [is_face ? this->dof_access_index :
2896 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2897 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2898 IndexStorageVariants::interleaved &&
2899 use_vectorized_path)
2901 const unsigned int *dof_indices =
2903 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2907 [this->first_selected_component] *
2910 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2911 if (n_components == 1 || this->n_fe_components == 1)
2912 for (
unsigned int comp = 0; comp < n_components; ++comp)
2914 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2917 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2919 if (n_components == 1 || this->n_fe_components == 1)
2920 for (
unsigned int i = 0; i < dofs_per_component;
2921 ++i, dof_indices += n_lanes)
2922 for (
unsigned int comp = 0; comp < n_components; ++comp)
2923 operation.process_dof_gather(dof_indices,
2927 values_dofs[comp][i],
2930 for (
unsigned int comp = 0; comp < n_components; ++comp)
2931 for (
unsigned int i = 0; i < dofs_per_component;
2932 ++i, dof_indices += n_lanes)
2933 operation.process_dof_gather(dof_indices,
2937 values_dofs[comp][i],
2944 std::array<const unsigned int *, n_lanes> dof_indices;
2945 dof_indices.fill(
nullptr);
2950 bool has_constraints =
false;
2951 const unsigned int n_components_read =
2952 this->n_fe_components > 1 ? n_components : 1;
2956 for (
unsigned int v = 0; v < n_lanes; ++v)
2962 const std::pair<unsigned int, unsigned int> *my_index_start =
2963 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2964 this->first_selected_component];
2969 if (my_index_start[n_components_read].
second !=
2970 my_index_start[0].
second)
2971 has_constraints =
true;
2974 dof_info.
dof_indices.data() + my_index_start[0].first;
2979 for (
unsigned int v = 0; v < n_lanes; ++v)
2984 const std::pair<unsigned int, unsigned int> *my_index_start =
2985 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2986 this->first_selected_component];
2987 if (my_index_start[n_components_read].
second !=
2988 my_index_start[0].
second)
2989 has_constraints =
true;
2996 dof_info.hanging_node_constraint_masks_comp
2997 [this->active_fe_index][this->first_selected_component])
2998 has_hn_constraints = true;
3001 my_index_start[0].
first ||
3004 my_index_start[0].
first,
3007 dof_info.
dof_indices.data() + my_index_start[0].first;
3011 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3012 return i != numbers::invalid_unsigned_int;
3014 for (
unsigned int comp = 0; comp < n_components; ++comp)
3015 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3016 operation.process_empty(values_dofs[comp][i]);
3020 if (!has_constraints && apply_constraints)
3022 if (n_components == 1 || this->n_fe_components == 1)
3024 for (
unsigned int v = 0; v < n_lanes; ++v)
3029 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3030 for (
unsigned int comp = 0; comp < n_components; ++comp)
3031 operation.process_dof(dof_indices[v][i],
3033 values_dofs[comp][i][v]);
3038 for (
unsigned int comp = 0; comp < n_components; ++comp)
3039 for (
unsigned int v = 0; v < n_lanes; ++v)
3044 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3045 operation.process_dof(
3046 dof_indices[v][comp * dofs_per_component + i],
3048 values_dofs[comp][i][v]);
3060 for (
unsigned int v = 0; v < n_lanes; ++v)
3066 const unsigned int cell_dof_index =
3067 cell_index * this->n_fe_components + this->first_selected_component;
3068 const unsigned int n_components_read =
3069 this->n_fe_components > 1 ? n_components : 1;
3070 unsigned int index_indicators =
3072 unsigned int next_index_indicators =
3073 dof_info.
row_starts[cell_dof_index + 1].second;
3077 if (apply_constraints ==
false &&
3078 (dof_info.
row_starts[cell_dof_index].second !=
3079 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3085 dof_info.hanging_node_constraint_masks_comp
3086 [this->active_fe_index][this->first_selected_component])))
3095 [this->first_selected_component] +
3097 next_index_indicators = index_indicators;
3100 if (n_components == 1 || this->n_fe_components == 1)
3102 unsigned int ind_local = 0;
3103 for (; index_indicators != next_index_indicators; ++index_indicators)
3105 const std::pair<unsigned short, unsigned short> indicator =
3108 for (
unsigned int j = 0; j < indicator.first; ++j)
3109 for (
unsigned int comp = 0; comp < n_components; ++comp)
3110 operation.process_dof(dof_indices[v][j],
3112 values_dofs[comp][ind_local + j][v]);
3114 ind_local += indicator.first;
3115 dof_indices[v] += indicator.first;
3119 Number
value[n_components];
3120 for (
unsigned int comp = 0; comp < n_components; ++comp)
3121 operation.pre_constraints(values_dofs[comp][ind_local][v],
3124 const Number *data_val =
3126 const Number *end_pool =
3128 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3129 for (
unsigned int comp = 0; comp < n_components; ++comp)
3130 operation.process_constraint(*dof_indices[v],
3135 for (
unsigned int comp = 0; comp < n_components; ++comp)
3136 operation.post_constraints(
value[comp],
3137 values_dofs[comp][ind_local][v]);
3143 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3144 for (
unsigned int comp = 0; comp < n_components; ++comp)
3145 operation.process_dof(*dof_indices[v],
3147 values_dofs[comp][ind_local][v]);
3156 for (
unsigned int comp = 0; comp < n_components; ++comp)
3158 unsigned int ind_local = 0;
3161 for (; index_indicators != next_index_indicators;
3164 const std::pair<unsigned short, unsigned short> indicator =
3168 for (
unsigned int j = 0; j < indicator.first; ++j)
3169 operation.process_dof(dof_indices[v][j],
3171 values_dofs[comp][ind_local + j][v]);
3172 ind_local += indicator.first;
3173 dof_indices[v] += indicator.first;
3178 operation.pre_constraints(values_dofs[comp][ind_local][v],
3181 const Number *data_val =
3183 const Number *end_pool =
3186 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3187 operation.process_constraint(*dof_indices[v],
3192 operation.post_constraints(
value,
3193 values_dofs[comp][ind_local][v]);
3200 for (; ind_local < dofs_per_component;
3201 ++dof_indices[v], ++ind_local)
3204 operation.process_dof(*dof_indices[v],
3206 values_dofs[comp][ind_local][v]);
3209 if (apply_constraints ==
true && comp + 1 < n_components)
3210 next_index_indicators =
3211 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3223 typename VectorizedArrayType>
3224template <
typename VectorType,
typename VectorOperation>
3229 const std::array<VectorType *, n_components_> &src)
const
3233 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3234 unsigned int index = this->first_selected_component * dofs_per_component;
3235 for (
unsigned int comp = 0; comp < n_components; ++comp)
3237 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3239 operation.process_empty(
3240 this->values_dofs[comp * dofs_per_component + i]);
3241 operation.process_dof_global(
3242 local_dof_indices[this->data->lexicographic_numbering[
index]],
3244 this->values_dofs[comp * dofs_per_component + i][0]);
3255 typename VectorizedArrayType>
3256template <
typename VectorType,
typename VectorOperation>
3261 const std::array<VectorType *, n_components_> &src,
3264 n_components_> &vectors_sm,
3265 const std::bitset<n_lanes> &
mask)
const
3274 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3277 is_face ? this->dof_access_index :
3279 const unsigned int n_active_lanes =
mask.count();
3282 const std::vector<unsigned int> &dof_indices_cont =
3285 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3286 std::array<VectorizedArrayType *, n_components> values_dofs;
3287 for (
unsigned int c = 0; c < n_components; ++c)
3288 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3289 c * dofs_per_component;
3293 const bool accesses_exterior_dofs =
3294 this->dof_access_index ==
3296 this->is_interior_face() ==
false;
3302 interleaved_contiguous &&
3303 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3305 const unsigned int dof_index =
3306 dof_indices_cont[this->cell * n_lanes] +
3309 [this->first_selected_component] *
3311 if (n_components == 1 || this->n_fe_components == 1)
3312 for (
unsigned int comp = 0; comp < n_components; ++comp)
3313 operation.process_dofs_vectorized(dofs_per_component,
3319 operation.process_dofs_vectorized(dofs_per_component * n_components,
3327 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3331 const unsigned int n_filled_lanes =
3334 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3335 n_active_lanes == n_lanes &&
3336 !accesses_exterior_dofs;
3338 if (vectors_sm[0] !=
nullptr)
3340 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3341 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3343 const auto upper_bound =
3344 std::min<unsigned int>(n_filled_lanes, n_lanes);
3345 for (
unsigned int v = 0; v < upper_bound; ++v)
3347 if (
mask[v] ==
false)
3349 vector_ptrs[v] =
nullptr;
3369 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3370 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3372 [this->active_fe_index][this->first_selected_component]);
3374 vector_ptrs[v] =
nullptr;
3376 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3377 vector_ptrs[v] =
nullptr;
3382 if (use_vectorized_path)
3384 if (n_components == 1 || this->n_fe_components == 1)
3386 for (
unsigned int comp = 0; comp < n_components; ++comp)
3388 auto vector_ptrs = compute_vector_ptrs(comp);
3389 operation.process_dofs_vectorized_transpose(
3398 auto vector_ptrs = compute_vector_ptrs(0);
3399 operation.process_dofs_vectorized_transpose(dofs_per_component *
3407 for (
unsigned int comp = 0; comp < n_components; ++comp)
3409 auto vector_ptrs = compute_vector_ptrs(
3410 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3412 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3413 operation.process_empty(values_dofs[comp][i]);
3415 if (n_components == 1 || this->n_fe_components == 1)
3417 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3418 if (
mask[v] ==
true)
3419 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3420 operation.process_dof(vector_ptrs[v][i],
3421 values_dofs[comp][i][v]);
3425 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3426 if (
mask[v] ==
true)
3427 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3428 operation.process_dof(
3429 vector_ptrs[v][i + comp * dofs_per_component],
3430 values_dofs[comp][i][v]);
3436 std::array<unsigned int, n_lanes> dof_indices;
3437 std::fill(dof_indices.begin(),
3442 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3446 if (
mask[v] ==
true)
3448 dof_indices_cont[cells[v]] +
3451 [this->first_selected_component] *
3457 if (use_vectorized_path)
3463 if (n_components == 1 || this->n_fe_components == 1)
3464 for (
unsigned int comp = 0; comp < n_components; ++comp)
3465 operation.process_dofs_vectorized_transpose(dofs_per_component,
3471 operation.process_dofs_vectorized_transpose(dofs_per_component *
3480 interleaved_contiguous_strided)
3482 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3483 if (n_components == 1 || this->n_fe_components == 1)
3484 for (
unsigned int comp = 0; comp < n_components; ++comp)
3485 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3486 src[comp]->
begin());
3489 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3491 if (n_components == 1 || this->n_fe_components == 1)
3492 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3494 for (
unsigned int comp = 0; comp < n_components; ++comp)
3495 operation.process_dof_gather(dof_indices.data(),
3498 src_ptrs[comp] + i * n_lanes,
3499 values_dofs[comp][i],
3503 for (
unsigned int comp = 0; comp < n_components; ++comp)
3504 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3506 operation.process_dof_gather(
3509 (comp * dofs_per_component + i) * n_lanes,
3510 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3511 values_dofs[comp][i],
3519 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3521 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (
unsigned int comp = 0; comp < n_components; ++comp)
3524 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3525 src[comp]->
begin());
3528 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3530 const unsigned int *offsets =
3532 if (n_components == 1 || this->n_fe_components == 1)
3533 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3535 for (
unsigned int comp = 0; comp < n_components; ++comp)
3536 operation.process_dof_gather(dof_indices.data(),
3540 values_dofs[comp][i],
3543 for (
unsigned int v = 0; v < n_lanes; ++v)
3544 dof_indices[v] += offsets[v];
3547 for (
unsigned int comp = 0; comp < n_components; ++comp)
3548 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3550 operation.process_dof_gather(dof_indices.data(),
3554 values_dofs[comp][i],
3557 for (
unsigned int v = 0; v < n_lanes; ++v)
3558 dof_indices[v] += offsets[v];
3563 for (
unsigned int comp = 0; comp < n_components; ++comp)
3565 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3566 operation.process_empty(values_dofs[comp][i]);
3567 if (accesses_exterior_dofs)
3569 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3570 if (
mask[v] ==
true)
3573 [ind][cells[v] / VectorizedArrayType::size()] ==
3577 if (n_components == 1 || this->n_fe_components == 1)
3579 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3580 operation.process_dof(dof_indices[v] + i,
3582 values_dofs[comp][i][v]);
3586 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3587 operation.process_dof(dof_indices[v] + i +
3588 comp * dofs_per_component,
3590 values_dofs[comp][i][v]);
3595 const unsigned int offset =
3598 if (n_components == 1 || this->n_fe_components == 1)
3600 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3601 operation.process_dof(dof_indices[v] + i * offset,
3603 values_dofs[comp][i][v]);
3607 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3608 operation.process_dof(
3610 (i + comp * dofs_per_component) * offset,
3612 values_dofs[comp][i][v]);
3623 if (n_components == 1 || this->n_fe_components == 1)
3625 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3626 if (
mask[v] ==
true)
3627 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3628 operation.process_dof(dof_indices[v] + i,
3630 values_dofs[comp][i][v]);
3634 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3635 if (
mask[v] ==
true)
3636 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3637 operation.process_dof(dof_indices[v] + i +
3638 comp * dofs_per_component,
3640 values_dofs[comp][i][v]);
3645 const unsigned int *offsets =
3647 [ind][VectorizedArrayType::size() * this->cell];
3648 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3650 if (n_components == 1 || this->n_fe_components == 1)
3651 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3653 if (
mask[v] ==
true)
3654 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3655 operation.process_dof(dof_indices[v] + i * offsets[v],
3657 values_dofs[comp][i][v]);
3661 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3662 if (
mask[v] ==
true)
3663 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3664 operation.process_dof(
3666 (i + comp * dofs_per_component) * offsets[v],
3668 values_dofs[comp][i][v]);
3679 typename VectorType,
3680 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
3681 decltype(std::declval<VectorType>().begin())
3682 get_beginning(VectorType &vec)
3689 typename VectorType,
3690 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
3691 typename VectorType::value_type *
3692 get_beginning(VectorType &)
3697 template <
typename VectorType,
3698 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
3700 const std::vector<ArrayView<const typename VectorType::value_type>> *
3701 get_shared_vector_data(VectorType *vec,
3702 const bool is_valid_mode_for_sm,
3703 const unsigned int active_fe_index,
3707 if (is_valid_mode_for_sm &&
3710 active_fe_index == 0)
3711 return &vec->shared_vector_data();
3716 template <
typename VectorType,
3717 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
3719 const std::vector<ArrayView<const typename VectorType::value_type>> *
3720 get_shared_vector_data(VectorType *,
3728 template <
int n_components,
typename VectorType>
3730 std::array<
typename internal::BlockVectorSelector<
3735 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3739 get_vector_data(VectorType &src,
3740 const unsigned int first_index,
3741 const bool is_valid_mode_for_sm,
3742 const unsigned int active_fe_index,
3748 std::array<
typename internal::BlockVectorSelector<
3754 ArrayView<
const typename internal::BlockVectorSelector<
3760 for (
unsigned int d = 0;
d < n_components; ++
d)
3761 src_data.first[d] = internal::BlockVectorSelector<
3767 for (
unsigned int d = 0;
d < n_components; ++
d)
3768 src_data.second[d] = get_shared_vector_data(
3769 const_cast<typename internal::BlockVectorSelector<
3770 std::remove_const_t<VectorType>,
3772 *>(src_data.first[d]),
3773 is_valid_mode_for_sm,
3787 typename VectorizedArrayType>
3792 if (this->dof_info ==
nullptr ||
3794 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3795 this->dof_info->hanging_node_constraint_masks_comp
3796 [this->active_fe_index][this->first_selected_component] ==
false)
3799 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3802 bool hn_available =
false;
3804 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3806 for (
unsigned int v = 0; v < n_lanes; ++v)
3818 constraint_mask[v] =
mask;
3824 if (hn_available ==
false)
3829 this->data->data.front().fe_degree,
3830 this->get_shape_info(),
3842 typename VectorizedArrayType>
3843template <
typename VectorType>
3847 const unsigned int first_index,
3848 const std::bitset<n_lanes> &
mask)
3850 const auto src_data = internal::get_vector_data<n_components_>(
3853 this->dof_access_index ==
3855 this->active_fe_index,
3859 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3861 apply_hanging_node_constraints(
false);
3864 this->dof_values_initialized =
true;
3874 typename VectorizedArrayType>
3875template <
typename VectorType>
3879 const unsigned int first_index,
3880 const std::bitset<n_lanes> &
mask)
3882 const auto src_data = internal::get_vector_data<n_components_>(
3885 this->dof_access_index ==
3887 this->active_fe_index,
3891 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3894 this->dof_values_initialized =
true;
3904 typename VectorizedArrayType>
3905template <
typename VectorType>
3909 const unsigned int first_index,
3910 const std::bitset<n_lanes> &
mask)
const
3913 Assert(this->dof_values_initialized ==
true,
3917 apply_hanging_node_constraints(
true);
3919 const auto dst_data = internal::get_vector_data<n_components_>(
3922 this->dof_access_index ==
3924 this->active_fe_index,
3929 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
3938 typename VectorizedArrayType>
3939template <
typename VectorType>
3943 const unsigned int first_index,
3944 const std::bitset<n_lanes> &
mask)
const
3947 Assert(this->dof_values_initialized ==
true,
3951 const auto dst_data = internal::get_vector_data<n_components_>(
3954 this->dof_access_index ==
3956 this->active_fe_index,
3960 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
3969 typename VectorizedArrayType>
3970template <
typename VectorType>
3974 const unsigned int first_index,
3975 const std::bitset<n_lanes> &
mask)
const
3978 Assert(this->dof_values_initialized ==
true,
3982 const auto dst_data = internal::get_vector_data<n_components_>(
3985 this->dof_access_index ==
3987 this->active_fe_index,
3991 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4004 typename VectorizedArrayType>
4010 VectorizedArrayType>::value_type
4015 if constexpr (n_components == 1)
4016 return this->values_dofs[dof];
4019 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4021 for (
unsigned int comp = 0; comp < n_components; ++comp)
4022 return_value[comp] = this->values_dofs[comp * dofs + dof];
4023 return return_value;
4033 typename VectorizedArrayType>
4039 VectorizedArrayType>::value_type
4041 get_value(
const unsigned int q_point)
const
4044 Assert(this->values_quad_initialized ==
true,
4049 if constexpr (n_components == 1)
4050 return this->values_quad[q_point];
4053 if (n_components == dim &&
4054 this->data->element_type ==
4059 Assert(this->values_quad_initialized ==
true,
4064 Assert(this->J_value !=
nullptr,
4067 const std::size_t nqp = this->n_quadrature_points;
4075 const VectorizedArrayType inv_det =
4076 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4077 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4078 this->jacobian[0][2][2];
4081 for (
unsigned int comp = 0; comp < n_components; ++comp)
4082 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4083 jac[comp][comp] * inv_det;
4090 this->jacobian[q_point] :
4099 const VectorizedArrayType inv_det =
4100 (is_face && dim == 2 && this->get_face_no() < 2) ?
4104 for (
unsigned int comp = 0; comp < n_components; ++comp)
4106 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4107 for (
unsigned int e = 1;
e < dim; ++
e)
4109 this->values_quad[e * nqp + q_point] * jac[comp][e];
4110 value_out[comp] *= inv_det;
4117 const std::size_t nqp = this->n_quadrature_points;
4119 for (
unsigned int comp = 0; comp < n_components; ++comp)
4120 return_value[comp] = this->values_quad[comp * nqp + q_point];
4121 return return_value;
4132 typename VectorizedArrayType>
4138 VectorizedArrayType>::gradient_type
4143 Assert(this->gradients_quad_initialized ==
true,
4148 Assert(this->jacobian !=
nullptr,
4150 "update_gradients"));
4151 const std::size_t nqp = this->n_quadrature_points;
4153 if constexpr (n_components == dim && dim > 1)
4155 if (this->data->element_type ==
4160 Assert(this->gradients_quad_initialized ==
true,
4165 Assert(this->jacobian !=
nullptr,
4167 "update_gradients"));
4168 const std::size_t nqp = this->n_quadrature_points;
4169 const std::size_t nqp_d = nqp * dim;
4172 this->gradients_quad + q_point * dim;
4183 const VectorizedArrayType inv_det =
4184 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4185 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4186 this->jacobian[0][2][2];
4189 for (
unsigned int d = 0;
d < dim; ++
d)
4190 for (
unsigned int comp = 0; comp < n_components; ++comp)
4191 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4193 (jac[comp][comp] * inv_det);
4205 const VectorizedArrayType inv_det =
4206 (is_face && dim == 2 && this->get_face_no() < 2) ?
4210 VectorizedArrayType tmp[dim][dim];
4212 for (
unsigned int d = 0;
d < dim; ++
d)
4213 for (
unsigned int e = 0;
e < dim; ++
e)
4216 for (
unsigned int f = 1; f < dim; ++f)
4217 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4219 for (
unsigned int comp = 0; comp < n_components; ++comp)
4220 for (
unsigned int d = 0;
d < dim; ++
d)
4222 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4223 for (
unsigned int f = 1; f < dim; ++f)
4224 res += jac[comp][f] * tmp[d][f];
4226 grad_out[comp][
d] = res * inv_det;
4237 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4239 "update_hessians"));
4241 const auto jac_grad =
4242 this->jacobian_gradients_non_inverse[q_point];
4244 this->jacobian[q_point];
4248 const VectorizedArrayType inv_det =
4249 (is_face && dim == 2 && this->get_face_no() < 2) ?
4256 VectorizedArrayType tmp[dim][dim];
4257 for (
unsigned int d = 0;
d < dim; ++
d)
4258 for (
unsigned int e = 0;
e < dim; ++
e)
4261 for (
unsigned int f = 1; f < dim; ++f)
4262 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4267 for (
unsigned int d = 0;
d < dim; ++
d)
4269 for (
unsigned int e = 0;
e < dim; ++
e)
4271 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4272 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4273 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4276 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4278 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4283 for (
unsigned int d = 0;
d < dim; ++
d)
4284 for (
unsigned int e = 0;
e < dim; ++
e)
4286 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4287 for (
unsigned int f = 1; f < dim; ++f)
4288 res += tmp[f][d] * inv_t_jac[e][f];
4289 grad_out[
d][
e] = res;
4295 VectorizedArrayType tmp3[dim], tmp4[dim];
4296 for (
unsigned int d = 0;
d < dim; ++
d)
4298 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4299 for (
unsigned int e = 1;
e < dim; ++
e)
4300 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4302 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4303 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4304 for (
unsigned int d = 0;
d < dim; ++
d)
4306 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4307 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4309 for (
unsigned int d = 0;
d < dim; ++
d)
4311 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4312 for (
unsigned int e = 1;
e < dim; ++
e)
4313 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4316 VectorizedArrayType tmp2[dim];
4317 for (
unsigned int d = 0;
d < dim; ++
d)
4319 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4320 for (
unsigned e = 1;
e < dim; ++
e)
4322 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4325 for (
unsigned int d = 0;
d < dim; ++
d)
4326 for (
unsigned int e = 0;
e < dim; ++
e)
4328 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4332 grad_out[
d][
e] *= inv_det;
4343 for (
unsigned int comp = 0; comp < n_components; ++comp)
4344 for (
unsigned int d = 0;
d < dim; ++
d)
4346 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4347 this->jacobian[0][
d][
d];
4356 for (
unsigned int comp = 0; comp < n_components; ++comp)
4357 for (
unsigned int d = 0;
d < dim; ++
d)
4360 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4361 for (
unsigned int e = 1;
e < dim; ++
e)
4362 grad_out[comp][d] +=
4364 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4367 if constexpr (n_components == 1)
4379 typename VectorizedArrayType>
4385 VectorizedArrayType>::value_type
4391 Assert(this->gradients_quad_initialized ==
true,
4395 Assert(this->normal_x_jacobian !=
nullptr,
4397 "update_gradients"));
4399 const std::size_t nqp = this->n_quadrature_points;
4403 for (
unsigned int comp = 0; comp < n_components; ++comp)
4405 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4406 (this->normal_x_jacobian[0][dim - 1]);
4409 const std::size_t
index =
4411 for (
unsigned int comp = 0; comp < n_components; ++comp)
4413 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4414 this->normal_x_jacobian[
index][0];
4415 for (
unsigned int d = 1;
d < dim; ++
d)
4417 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4418 this->normal_x_jacobian[
index][
d];
4421 if constexpr (n_components == 1)
4433 template <
typename VectorizedArrayType>
4436 const VectorizedArrayType *
const hessians,
4438 VectorizedArrayType (&tmp)[1][1])
4440 tmp[0][0] = jac[0][0] *
hessians[0];
4443 template <
typename VectorizedArrayType>
4446 const VectorizedArrayType *
const hessians,
4447 const unsigned int nqp,
4448 VectorizedArrayType (&tmp)[2][2])
4450 for (
unsigned int d = 0;
d < 2; ++
d)
4458 template <
typename VectorizedArrayType>
4461 const VectorizedArrayType *
const hessians,
4462 const unsigned int nqp,
4463 VectorizedArrayType (&tmp)[3][3])
4465 for (
unsigned int d = 0;
d < 3; ++
d)
4486 typename VectorizedArrayType>
4491 VectorizedArrayType>::hessian_type
4496 Assert(this->hessians_quad_initialized ==
true,
4501 Assert(this->jacobian !=
nullptr,
4511 const std::size_t nqp = this->n_quadrature_points;
4512 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4517 for (
unsigned int comp = 0; comp < n_components; ++comp)
4519 for (
unsigned int d = 0;
d < dim; ++
d)
4520 hessian_out[comp][d][d] =
4521 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4522 (jac[
d][
d] * jac[
d][
d]);
4528 hessian_out[comp][0][1] =
4529 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4530 (jac[0][0] * jac[1][1]);
4533 hessian_out[comp][0][1] =
4534 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4535 (jac[0][0] * jac[1][1]);
4536 hessian_out[comp][0][2] =
4537 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4538 (jac[0][0] * jac[2][2]);
4539 hessian_out[comp][1][2] =
4540 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4541 (jac[1][1] * jac[2][2]);
4546 for (
unsigned int d = 0;
d < dim; ++
d)
4547 for (
unsigned int e = d + 1;
e < dim; ++
e)
4548 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4554 for (
unsigned int comp = 0; comp < n_components; ++comp)
4556 VectorizedArrayType tmp[dim][dim];
4557 internal::hessian_unit_times_jac(
4558 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4561 for (
unsigned int d = 0;
d < dim; ++
d)
4562 for (
unsigned int e = d;
e < dim; ++
e)
4564 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4565 for (
unsigned int f = 1; f < dim; ++f)
4566 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4573 for (
unsigned int d = 0;
d < dim; ++
d)
4574 for (
unsigned int e = d + 1;
e < dim; ++
e)
4575 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4581 const auto &jac_grad = this->jacobian_gradients[q_point];
4582 for (
unsigned int comp = 0; comp < n_components; ++comp)
4584 VectorizedArrayType tmp[dim][dim];
4585 internal::hessian_unit_times_jac(
4586 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4589 for (
unsigned int d = 0;
d < dim; ++
d)
4590 for (
unsigned int e = d;
e < dim; ++
e)
4592 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4593 for (
unsigned int f = 1; f < dim; ++f)
4594 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4598 for (
unsigned int d = 0;
d < dim; ++
d)
4599 for (
unsigned int e = 0;
e < dim; ++
e)
4600 hessian_out[comp][d][d] +=
4602 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4605 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4606 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4607 for (
unsigned int f = 0; f < dim; ++f)
4608 hessian_out[comp][d][e] +=
4609 jac_grad[count][f] *
4610 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4613 for (
unsigned int d = 0;
d < dim; ++
d)
4614 for (
unsigned int e = d + 1;
e < dim; ++
e)
4615 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4618 if constexpr (n_components == 1)
4619 return hessian_out[0];
4630 typename VectorizedArrayType>
4635 VectorizedArrayType>::gradient_type
4641 Assert(this->hessians_quad_initialized ==
true,
4652 const std::size_t nqp = this->n_quadrature_points;
4653 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4659 for (
unsigned int comp = 0; comp < n_components; ++comp)
4660 for (
unsigned int d = 0;
d < dim; ++
d)
4661 hessian_out[comp][d] =
4662 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4663 (jac[
d][
d] * jac[
d][
d]);
4668 for (
unsigned int comp = 0; comp < n_components; ++comp)
4672 VectorizedArrayType tmp[dim][dim];
4673 internal::hessian_unit_times_jac(
4674 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4678 for (
unsigned int d = 0;
d < dim; ++
d)
4680 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4681 for (
unsigned int f = 1; f < dim; ++f)
4682 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4689 const auto &jac_grad = this->jacobian_gradients[q_point];
4690 for (
unsigned int comp = 0; comp < n_components; ++comp)
4694 VectorizedArrayType tmp[dim][dim];
4695 internal::hessian_unit_times_jac(
4696 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4700 for (
unsigned int d = 0;
d < dim; ++
d)
4702 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4703 for (
unsigned int f = 1; f < dim; ++f)
4704 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4707 for (
unsigned int d = 0;
d < dim; ++
d)
4708 for (
unsigned int e = 0;
e < dim; ++
e)
4709 hessian_out[comp][d] +=
4711 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4715 if constexpr (n_components == 1)
4716 return hessian_out[0];
4727 typename VectorizedArrayType>
4732 VectorizedArrayType>::value_type
4738 Assert(this->hessians_quad_initialized ==
true,
4743 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4744 if constexpr (n_components == 1)
4746 VectorizedArrayType
sum = hess_diag[0];
4747 for (
unsigned int d = 1;
d < dim; ++
d)
4748 sum += hess_diag[d];
4754 for (
unsigned int comp = 0; comp < n_components; ++comp)
4756 laplacian_out[comp] = hess_diag[comp][0];
4757 for (
unsigned int d = 1;
d < dim; ++
d)
4758 laplacian_out[comp] += hess_diag[comp][d];
4760 return laplacian_out;
4770 typename VectorizedArrayType>
4776 this->dof_values_initialized =
true;
4778 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4780 for (
unsigned int comp = 0; comp < n_components; ++comp)
4781 if constexpr (n_components == 1)
4782 this->values_dofs[comp * dofs + dof] = val_in;
4784 this->values_dofs[comp * dofs + dof] = val_in[comp];
4793 typename VectorizedArrayType>
4796 submit_value(
const value_type val_in,
const unsigned int q_point)
4802 Assert(this->J_value !=
nullptr,
4806 this->values_quad_submitted =
true;
4809 const std::size_t nqp = this->n_quadrature_points;
4810 VectorizedArrayType *
values = this->values_quad + q_point;
4812 const VectorizedArrayType JxW =
4814 this->J_value[0] * this->quadrature_weights[q_point] :
4815 this->J_value[q_point];
4816 if constexpr (n_components == 1)
4817 values[0] = val_in * JxW;
4820 if (n_components == dim &&
4821 this->data->element_type ==
4826 Assert(this->J_value !=
nullptr,
4831 this->values_quad_submitted =
true;
4834 VectorizedArrayType *
values = this->values_quad + q_point;
4835 const std::size_t nqp = this->n_quadrature_points;
4841 const VectorizedArrayType weight =
4842 this->quadrature_weights[q_point];
4844 for (
unsigned int comp = 0; comp < n_components; ++comp)
4845 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4852 this->jacobian[q_point] :
4857 const VectorizedArrayType fac =
4859 this->quadrature_weights[q_point] :
4861 this->J_value[q_point] :
4862 this->J_value[0] * this->quadrature_weights[q_point]) *
4863 ((dim == 2 && this->get_face_no() < 2) ?
4872 for (
unsigned int comp = 0; comp < n_components; ++comp)
4874 values[comp * nqp] = val_in[0] * jac[0][comp];
4875 for (
unsigned int e = 1;
e < dim; ++
e)
4876 values[comp * nqp] += val_in[e] * jac[e][comp];
4877 values[comp * nqp] *= fac;
4882 for (
unsigned int comp = 0; comp < n_components; ++comp)
4883 values[comp * nqp] = val_in[comp] * JxW;
4893 typename VectorizedArrayType>
4894template <
int,
typename>
4898 const unsigned int q_point)
4900 static_assert(n_components == 1 && dim == 1,
4901 "Do not try to modify the default template parameters used for"
4902 " selectively enabling this function via std::enable_if!");
4903 submit_value(val_in[0], q_point);
4912 typename VectorizedArrayType>
4915 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
4921 Assert(this->J_value !=
nullptr,
4923 "update_gradients"));
4924 Assert(this->jacobian !=
nullptr,
4926 "update_gradients"));
4928 this->gradients_quad_submitted =
true;
4931 if constexpr (dim > 1 && n_components == dim)
4933 if (this->data->element_type ==
4942 Assert(this->J_value !=
nullptr,
4944 "update_gradients"));
4945 Assert(this->jacobian !=
nullptr,
4947 "update_gradients"));
4949 this->gradients_quad_submitted =
true;
4952 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
4953 VectorizedArrayType *
values =
4954 this->values_from_gradients_quad + q_point;
4955 const std::size_t nqp = this->n_quadrature_points;
4956 const std::size_t nqp_d = nqp * dim;
4966 const VectorizedArrayType weight =
4967 this->quadrature_weights[q_point];
4968 for (
unsigned int d = 0;
d < dim; ++
d)
4969 for (
unsigned int comp = 0; comp < n_components; ++comp)
4970 gradients[comp * nqp_d + d] = grad_in[comp][d] *
4972 (jac[comp][comp] * weight);
4984 const VectorizedArrayType fac =
4986 this->quadrature_weights[q_point] :
4987 this->J_value[0] * this->quadrature_weights[q_point] *
4988 ((dim == 2 && this->get_face_no() < 2) ?
4993 VectorizedArrayType tmp[dim][dim];
4994 for (
unsigned int d = 0;
d < dim; ++
d)
4995 for (
unsigned int e = 0;
e < dim; ++
e)
4997 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
4998 for (
unsigned int f = 1; f < dim; ++f)
4999 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5001 for (
unsigned int comp = 0; comp < n_components; ++comp)
5002 for (
unsigned int d = 0;
d < dim; ++
d)
5004 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5005 for (
unsigned int f = 1; f < dim; ++f)
5006 res += jac[f][comp] * tmp[d][f];
5015 const auto jac_grad =
5016 this->jacobian_gradients_non_inverse[q_point];
5018 this->jacobian[q_point];
5022 const VectorizedArrayType fac =
5023 (!is_face) ? this->quadrature_weights[q_point] :
5024 this->J_value[q_point] *
5025 ((dim == 2 && this->get_face_no() < 2) ?
5034 VectorizedArrayType tmp3[dim], tmp4[dim];
5035 for (
unsigned int d = 0;
d < dim; ++
d)
5037 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5038 for (
unsigned int e = 1;
e < dim; ++
e)
5039 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5041 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5042 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5043 for (
unsigned int d = 0;
d < dim; ++
d)
5045 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5046 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5048 for (
unsigned int d = 0;
d < dim; ++
d)
5050 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5051 for (
unsigned int e = 1;
e < dim; ++
e)
5052 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5058 VectorizedArrayType tmp[dim][dim];
5061 for (
unsigned int d = 0;
d < dim; ++
d)
5062 for (
unsigned int e = 0;
e < dim; ++
e)
5064 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5065 for (
unsigned int f = 1; f < dim; ++f)
5066 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5069 for (
unsigned int d = 0;
d < dim; ++
d)
5070 for (
unsigned int e = 0;
e < dim; ++
e)
5072 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5073 for (
unsigned int f = 1; f < dim; ++f)
5074 res += t_jac[d][f] * tmp[e][f];
5081 VectorizedArrayType
value[dim];
5082 for (
unsigned int d = 0;
d < dim; ++
d)
5084 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5085 for (
unsigned int e = 1;
e < dim; ++
e)
5086 value[d] += tmp[d][e] * jac_grad[d][e];
5088 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5089 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5090 for (
unsigned int d = 0;
d < dim; ++
d)
5092 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5093 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5098 for (
unsigned int d = 0;
d < dim; ++
d)
5100 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5101 for (
unsigned int e = 1;
e < dim; ++
e)
5102 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5103 for (
unsigned int e = 0;
e < dim; ++
e)
5104 value[e] -= t_jac[e][d] * tmp2;
5107 for (
unsigned int d = 0;
d < dim; ++
d)
5108 values[d * nqp] =
value[d];
5114 const std::size_t nqp_d = this->n_quadrature_points * dim;
5115 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5119 const VectorizedArrayType JxW =
5120 this->J_value[0] * this->quadrature_weights[q_point];
5126 std::array<VectorizedArrayType, dim> jac;
5127 for (
unsigned int d = 0;
d < dim; ++
d)
5128 jac[d] = this->jacobian[0][d][d];
5130 for (
unsigned int d = 0;
d < dim; ++
d)
5132 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5133 if constexpr (n_components == 1)
5134 gradients[d] = grad_in[d] * factor;
5136 for (
unsigned int comp = 0; comp < n_components; ++comp)
5137 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5144 this->jacobian[q_point] :
5146 const VectorizedArrayType JxW =
5148 this->J_value[q_point] :
5149 this->J_value[0] * this->quadrature_weights[q_point];
5150 if constexpr (n_components == 1)
5151 for (
unsigned int d = 0;
d < dim; ++
d)
5153 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5154 for (
unsigned int e = 1;
e < dim; ++
e)
5155 new_val += (jac[e][d] * grad_in[e]);
5159 for (
unsigned int comp = 0; comp < n_components; ++comp)
5160 for (
unsigned int d = 0;
d < dim; ++
d)
5162 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5163 for (
unsigned int e = 1;
e < dim; ++
e)
5164 new_val += (jac[e][d] * grad_in[comp][e]);
5176 typename VectorizedArrayType>
5177template <
int,
typename>
5181 const unsigned int q_point)
5183 static_assert(n_components == 1 && dim == 1,
5184 "Do not try to modify the default template parameters used for"
5185 " selectively enabling this function via std::enable_if!");
5186 submit_gradient(grad_in[0], q_point);
5195 typename VectorizedArrayType>
5201 Assert(this->normal_x_jacobian !=
nullptr,
5203 "update_gradients"));
5205 this->gradients_quad_submitted =
true;
5208 const std::size_t nqp_d = this->n_quadrature_points * dim;
5209 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5213 const VectorizedArrayType JxW_jac = this->J_value[0] *
5214 this->quadrature_weights[q_point] *
5215 this->normal_x_jacobian[0][dim - 1];
5216 for (
unsigned int comp = 0; comp < n_components; ++comp)
5218 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5219 gradients[comp * nqp_d + d] = VectorizedArrayType();
5220 if constexpr (n_components == 1)
5221 gradients[dim - 1] = grad_in * JxW_jac;
5223 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5228 const unsigned int index =
5231 this->normal_x_jacobian[
index];
5232 const VectorizedArrayType JxW =
5234 this->J_value[
index] * this->quadrature_weights[q_point] :
5235 this->J_value[
index];
5236 for (
unsigned int comp = 0; comp < n_components; ++comp)
5237 for (
unsigned int d = 0;
d < dim; ++
d)
5238 if constexpr (n_components == 1)
5241 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5251 typename VectorizedArrayType>
5254 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5260 Assert(this->J_value !=
nullptr,
5262 "update_hessians"));
5263 Assert(this->jacobian !=
nullptr,
5265 "update_hessians"));
5267 this->hessians_quad_submitted =
true;
5271 const std::size_t nqp = this->n_quadrature_points;
5272 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5275 const VectorizedArrayType JxW =
5276 this->J_value[0] * this->quadrature_weights[q_point];
5279 for (
unsigned int d = 0;
d < dim; ++
d)
5281 const auto jac_d = this->jacobian[0][
d][
d];
5282 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5283 for (
unsigned int comp = 0; comp < n_components; ++comp)
5284 if constexpr (n_components == 1)
5285 this->hessians_quad[
d * nqp + q_point] =
5286 hessian_in[
d][
d] * factor;
5288 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5289 hessian_in[comp][d][d] * factor;
5293 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5294 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5296 const auto jac_d = this->jacobian[0][
d][
d];
5297 const auto jac_e = this->jacobian[0][
e][
e];
5298 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5299 for (
unsigned int comp = 0; comp < n_components; ++comp)
5300 if constexpr (n_components == 1)
5301 this->hessians_quad[off_dia * nqp + q_point] =
5302 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5304 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5305 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5312 const VectorizedArrayType JxW =
5313 this->J_value[0] * this->quadrature_weights[q_point];
5314 for (
unsigned int comp = 0; comp < n_components; ++comp)
5317 if constexpr (n_components == 1)
5318 hessian_c = hessian_in;
5320 hessian_c = hessian_in[comp];
5323 VectorizedArrayType tmp[dim][dim];
5324 for (
unsigned int i = 0; i < dim; ++i)
5325 for (
unsigned int j = 0; j < dim; ++j)
5327 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5328 for (
unsigned int k = 1; k < dim; ++k)
5329 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5333 VectorizedArrayType tmp2[dim][dim];
5334 for (
unsigned int i = 0; i < dim; ++i)
5335 for (
unsigned int j = 0; j < dim; ++j)
5337 tmp2[i][j] = jac[0][i] * tmp[0][j];
5338 for (
unsigned int k = 1; k < dim; ++k)
5339 tmp2[i][j] += jac[k][i] * tmp[k][j];
5343 for (
unsigned int d = 0;
d < dim; ++
d)
5344 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5348 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5349 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5350 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5351 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5357 const VectorizedArrayType JxW = this->J_value[q_point];
5358 const auto &jac_grad = this->jacobian_gradients[q_point];
5359 for (
unsigned int comp = 0; comp < n_components; ++comp)
5362 if constexpr (n_components == 1)
5363 hessian_c = hessian_in;
5365 hessian_c = hessian_in[comp];
5368 VectorizedArrayType tmp[dim][dim];
5369 for (
unsigned int i = 0; i < dim; ++i)
5370 for (
unsigned int j = 0; j < dim; ++j)
5372 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5373 for (
unsigned int k = 1; k < dim; ++k)
5374 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5378 VectorizedArrayType tmp2[dim][dim];
5379 for (
unsigned int i = 0; i < dim; ++i)
5380 for (
unsigned int j = 0; j < dim; ++j)
5382 tmp2[i][j] = jac[0][i] * tmp[0][j];
5383 for (
unsigned int k = 1; k < dim; ++k)
5384 tmp2[i][j] += jac[k][i] * tmp[k][j];
5388 for (
unsigned int d = 0;
d < dim; ++
d)
5389 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5393 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5394 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5395 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5396 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5399 for (
unsigned int d = 0;
d < dim; ++
d)
5401 VectorizedArrayType
sum = 0;
5402 for (
unsigned int e = 0;
e < dim; ++
e)
5403 sum += hessian_c[e][e] * jac_grad[e][d];
5404 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5405 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5407 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5408 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5421 typename VectorizedArrayType>
5426 VectorizedArrayType>::value_type
5432 Assert(this->values_quad_submitted ==
true,
5437 const std::size_t nqp = this->n_quadrature_points;
5438 for (
unsigned int q = 0; q < nqp; ++q)
5439 for (
unsigned int comp = 0; comp < n_components; ++comp)
5440 return_value[comp] += this->values_quad[comp * nqp + q];
5441 if constexpr (n_components == 1)
5442 return return_value[0];
5444 return return_value;
5453 typename VectorizedArrayType>
5454template <
int,
typename>
5459 static_assert(n_components == dim,
5460 "Do not try to modify the default template parameters used for"
5461 " selectively enabling this function via std::enable_if!");
5464 Assert(this->gradients_quad_initialized ==
true,
5468 Assert(this->jacobian !=
nullptr,
5470 "update_gradients"));
5472 VectorizedArrayType divergence;
5473 const std::size_t nqp = this->n_quadrature_points;
5476 this->data->element_type ==
5479 VectorizedArrayType inv_det =
5482 this->jacobian[0][0][0] *
5483 ((dim == 2) ? this->jacobian[0][1][1] :
5484 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5492 if (is_face && dim == 2 && this->get_face_no() < 2)
5496 divergence = this->gradients_quad[q_point * dim];
5497 for (
unsigned int d = 1;
d < dim; ++
d)
5498 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5499 divergence *= inv_det;
5508 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5509 for (
unsigned int d = 1;
d < dim; ++
d)
5510 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5511 this->jacobian[0][
d][
d];
5518 this->jacobian[q_point] :
5520 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5521 for (
unsigned int e = 1;
e < dim; ++
e)
5522 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5523 for (
unsigned int d = 1;
d < dim; ++
d)
5524 for (
unsigned int e = 0;
e < dim; ++
e)
5526 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5538 typename VectorizedArrayType>
5539template <
int,
typename>
5544 static_assert(n_components == dim,
5545 "Do not try to modify the default template parameters used for"
5546 " selectively enabling this function via std::enable_if!");
5549 const auto grad = get_gradient(q_point);
5550 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5551 VectorizedArrayType half = Number(0.5);
5552 for (
unsigned int d = 0;
d < dim; ++
d)
5553 symmetrized[d] = grad[d][d];
5559 symmetrized[2] = grad[0][1] + grad[1][0];
5560 symmetrized[2] *= half;
5563 symmetrized[3] = grad[0][1] + grad[1][0];
5564 symmetrized[3] *= half;
5565 symmetrized[4] = grad[0][2] + grad[2][0];
5566 symmetrized[4] *= half;
5567 symmetrized[5] = grad[1][2] + grad[2][1];
5568 symmetrized[5] *= half;
5582 typename VectorizedArrayType>
5583template <
int,
typename>
5585 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5587 get_curl(const unsigned
int q_point) const
5589 static_assert(dim > 1 && n_components == dim,
5590 "Do not try to modify the default template parameters used for"
5591 " selectively enabling this function via std::enable_if!");
5595 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5599 curl[0] = grad[1][0] - grad[0][1];
5602 curl[0] = grad[2][1] - grad[1][2];
5603 curl[1] = grad[0][2] - grad[2][0];
5604 curl[2] = grad[1][0] - grad[0][1];
5618 typename VectorizedArrayType>
5619template <
int,
typename>
5623 const unsigned int q_point)
5625 static_assert(n_components == dim,
5626 "Do not try to modify the default template parameters used for"
5627 " selectively enabling this function via std::enable_if!");
5633 Assert(this->J_value !=
nullptr,
5635 "update_gradients"));
5636 Assert(this->jacobian !=
nullptr,
5638 "update_gradients"));
5640 this->gradients_quad_submitted =
true;
5643 const std::size_t nqp_d = this->n_quadrature_points * dim;
5644 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5646 if (this->data->element_type ==
5653 const VectorizedArrayType fac =
5655 this->quadrature_weights[q_point] * div_in :
5657 this->J_value[q_point] :
5658 this->J_value[0] * this->quadrature_weights[q_point]) *
5661 this->jacobian[this->cell_type >
5665 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5667 for (
unsigned int d = 0;
d < dim; ++
d)
5669 for (
unsigned int e = 0;
e < dim; ++
e)
5670 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5672 this->divergence_is_requested =
true;
5679 const VectorizedArrayType fac =
5680 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5681 for (
unsigned int d = 0;
d < dim; ++
d)
5683 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5684 for (
unsigned int e = 0;
e < dim; ++
e)
5685 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5692 this->jacobian[q_point] :
5694 const VectorizedArrayType fac =
5696 this->J_value[q_point] :
5697 this->J_value[0] * this->quadrature_weights[q_point]) *
5699 for (
unsigned int d = 0;
d < dim; ++
d)
5701 for (
unsigned int e = 0;
e < dim; ++
e)
5702 gradients[d * nqp_d + e] = jac[d][e] * fac;
5714 typename VectorizedArrayType>
5715template <
int,
typename>
5720 const unsigned int q_point)
5722 static_assert(n_components == dim,
5723 "Do not try to modify the default template parameters used for"
5724 " selectively enabling this function via std::enable_if!");
5727 this->data->element_type !=
5738 Assert(this->J_value !=
nullptr,
5740 "update_gradients"));
5741 Assert(this->jacobian !=
nullptr,
5743 "update_gradients"));
5745 this->gradients_quad_submitted =
true;
5748 const std::size_t nqp_d = this->n_quadrature_points * dim;
5749 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5752 const VectorizedArrayType JxW =
5753 this->J_value[0] * this->quadrature_weights[q_point];
5755 for (
unsigned int d = 0;
d < dim; ++
d)
5756 gradients[d * nqp_d + d] =
5758 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5759 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5761 const VectorizedArrayType
value =
5770 const VectorizedArrayType JxW =
5772 this->J_value[q_point] :
5773 this->J_value[0] * this->quadrature_weights[q_point];
5776 this->jacobian[q_point] :
5778 VectorizedArrayType weighted[dim][dim];
5779 for (
unsigned int i = 0; i < dim; ++i)
5781 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
5782 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
5784 const VectorizedArrayType
value =
5786 weighted[i][j] =
value;
5787 weighted[j][i] =
value;
5789 for (
unsigned int comp = 0; comp < dim; ++comp)
5790 for (
unsigned int d = 0;
d < dim; ++
d)
5792 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
5793 for (
unsigned int e = 1;
e < dim; ++
e)
5794 new_val += jac[e][d] * weighted[comp][e];
5806 typename VectorizedArrayType>
5807template <
int,
typename>
5811 const unsigned int q_point)
5813 static_assert(n_components == dim,
5814 "Do not try to modify the default template parameters used for"
5815 " selectively enabling this function via std::enable_if!");
5821 grad[1][0] = curl[0];
5822 grad[0][1] = -curl[0];
5825 grad[2][1] = curl[0];
5826 grad[1][2] = -curl[0];
5827 grad[0][2] = curl[1];
5828 grad[2][0] = -curl[1];
5829 grad[1][0] = curl[2];
5830 grad[0][1] = -curl[2];
5835 submit_gradient(grad, q_point);
5848 typename VectorizedArrayType>
5854 VectorizedArrayType>::
5856 const unsigned int fe_no,
5857 const unsigned int quad_no,
5858 const unsigned int first_selected_component,
5859 const unsigned int active_fe_index,
5860 const unsigned int active_quad_index)
5861 : BaseClass(matrix_free,
5863 first_selected_component,
5870 numbers::invalid_unsigned_int )
5871 , dofs_per_component(this->data->dofs_per_component_on_cell)
5872 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5873 , n_q_points(this->data->n_q_points)
5875 check_template_arguments(fe_no, 0);
5885 typename VectorizedArrayType>
5891 VectorizedArrayType>::
5893 const std::pair<unsigned int, unsigned int> &range,
5894 const unsigned int dof_no,
5895 const unsigned int quad_no,
5896 const unsigned int first_selected_component)
5900 first_selected_component,
5901 matrix_free.get_cell_active_fe_index(range))
5911 typename VectorizedArrayType>
5917 VectorizedArrayType>::
5922 const unsigned int first_selected_component)
5927 first_selected_component,
5929 , dofs_per_component(this->data->dofs_per_component_on_cell)
5930 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5931 , n_q_points(this->data->n_q_points)
5943 typename VectorizedArrayType>
5949 VectorizedArrayType>::
5953 const unsigned int first_selected_component)
5958 first_selected_component,
5960 , dofs_per_component(this->data->dofs_per_component_on_cell)
5961 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5962 , n_q_points(this->data->n_q_points)
5974 typename VectorizedArrayType>
5980 VectorizedArrayType>::
5983 const unsigned int first_selected_component)
5984 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
5986 other.mapped_geometry->get_quadrature(),
5987 other.mapped_geometry->get_fe_values().get_update_flags(),
5988 first_selected_component,
5990 , dofs_per_component(this->data->dofs_per_component_on_cell)
5991 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5992 , n_q_points(this->data->n_q_points)
6004 typename VectorizedArrayType>
6013 , dofs_per_component(this->data->dofs_per_component_on_cell)
6014 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6015 , n_q_points(this->data->n_q_points)
6027 typename VectorizedArrayType>
6033 VectorizedArrayType> &
6039 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6041 BaseClass::operator=(other);
6053 typename VectorizedArrayType>
6060 VectorizedArrayType>::
6061 check_template_arguments(
const unsigned int dof_no,
6062 const unsigned int first_selected_component)
6065 (void)first_selected_component;
6068 this->data->dofs_per_component_on_cell > 0,
6070 "There is nothing useful you can do with an FEEvaluation object with "
6071 "FE_Nothing, i.e., without DoFs! If you have passed to "
6072 "MatrixFree::reinit() a collection of finite elements also containing "
6073 "FE_Nothing, please check - before creating FEEvaluation - the category "
6074 "of the current range by calling either "
6075 "MatrixFree::get_cell_range_category(range) or "
6076 "MatrixFree::get_face_range_category(range). The returned category "
6077 "is the index of the active FE, which you can use to exclude "
6084 static_cast<unsigned int>(fe_degree) !=
6085 this->data->data.front().fe_degree) ||
6086 n_q_points != this->n_quadrature_points)
6088 std::string message =
6089 "-------------------------------------------------------\n";
6090 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6091 message +=
" Called --> FEEvaluation<dim,";
6095 message +=
",Number>(data";
6111 if (
static_cast<unsigned int>(fe_degree) ==
6112 this->data->data.front().fe_degree)
6114 proposed_dof_comp = dof_no;
6115 proposed_fe_comp = first_selected_component;
6118 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6120 for (
unsigned int nf = 0;
6123 if (this->matrix_free
6124 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6126 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6128 proposed_dof_comp = no;
6129 proposed_fe_comp = nf;
6133 this->mapping_data->descriptor[this->active_quad_index]
6135 proposed_quad_comp = this->quad_no;
6137 for (
unsigned int no = 0;
6142 .descriptor[this->active_quad_index]
6143 .n_q_points == n_q_points)
6145 proposed_quad_comp = no;
6152 if (proposed_dof_comp != first_selected_component)
6153 message +=
"Wrong vector component selection:\n";
6155 message +=
"Wrong quadrature formula selection:\n";
6156 message +=
" Did you mean FEEvaluation<dim,";
6160 message +=
",Number>(data";
6169 std::string correct_pos;
6170 if (proposed_dof_comp != dof_no)
6171 correct_pos =
" ^ ";
6174 if (proposed_quad_comp != this->quad_no)
6175 correct_pos +=
" ^ ";
6178 if (proposed_fe_comp != first_selected_component)
6179 correct_pos +=
" ^\n";
6181 correct_pos +=
" \n";
6187 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6188 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6189 message +=
"Wrong template arguments:\n";
6190 message +=
" Did you mean FEEvaluation<dim,";
6195 message +=
",Number>(data";
6203 std::string correct_pos;
6204 if (this->data->data.front().fe_degree !=
6205 static_cast<unsigned int>(fe_degree))
6209 if (proposed_n_q_points_1d != n_q_points_1d)
6210 correct_pos +=
" ^\n";
6212 correct_pos +=
" \n";
6213 message +=
" " + correct_pos;
6215 Assert(
static_cast<unsigned int>(fe_degree) ==
6216 this->data->data.front().fe_degree &&
6217 n_q_points == this->n_quadrature_points,
6223 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6234 typename VectorizedArrayType>
6243 Assert(this->matrix_free !=
nullptr,
6244 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6245 " Integer indexing is not possible."));
6253 const unsigned int offsets =
6254 this->mapping_data->data_index_offsets[
cell_index];
6255 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6256 this->J_value = &this->mapping_data->JxW_values[offsets];
6257 if (!this->mapping_data->jacobian_gradients[0].empty())
6259 this->jacobian_gradients =
6260 this->mapping_data->jacobian_gradients[0].data() + offsets;
6261 this->jacobian_gradients_non_inverse =
6262 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6268 for (
unsigned int i = 0; i < n_lanes; ++i)
6269 this->cell_ids[i] =
cell_index * n_lanes + i;
6276 this->cell_ids[i] =
cell_index * n_lanes + i;
6277 for (; i < n_lanes; ++i)
6281 if (this->mapping_data->quadrature_points.empty() ==
false)
6283 &this->mapping_data->quadrature_points
6284 [this->mapping_data->quadrature_point_offsets[this->cell]];
6287 this->is_reinitialized =
true;
6288 this->dof_values_initialized =
false;
6289 this->values_quad_initialized =
false;
6290 this->gradients_quad_initialized =
false;
6291 this->hessians_quad_initialized =
false;
6302 typename VectorizedArrayType>
6309 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6316 this->cell_ids = cell_ids;
6321 for (
unsigned int v = 0; v < n_lanes; ++v)
6335 if (this->mapped_geometry ==
nullptr)
6336 this->mapped_geometry =
6337 std::make_shared<internal::MatrixFreeFunctions::
6338 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6340 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6342 auto &this_jacobian_data = mapping_storage.jacobians[0];
6343 auto &this_J_value_data = mapping_storage.JxW_values;
6344 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6345 auto &this_jacobian_gradients_non_inverse_data =
6346 mapping_storage.jacobian_gradients_non_inverse[0];
6347 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6351 if (this_jacobian_data.size() != 2)
6352 this_jacobian_data.resize_fast(2);
6354 if (this_J_value_data.size() != 1)
6355 this_J_value_data.resize_fast(1);
6357 const auto &update_flags_cells =
6361 this_jacobian_gradients_data.size() != 1)
6363 this_jacobian_gradients_data.resize_fast(1);
6364 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6368 this_quadrature_points_data.size() != 1)
6369 this_quadrature_points_data.resize_fast(1);
6373 if (this_jacobian_data.size() != this->n_quadrature_points)
6374 this_jacobian_data.resize_fast(this->n_quadrature_points);
6376 if (this_J_value_data.size() != this->n_quadrature_points)
6377 this_J_value_data.resize_fast(this->n_quadrature_points);
6379 const auto &update_flags_cells =
6383 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6385 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6386 this_jacobian_gradients_non_inverse_data.resize_fast(
6387 this->n_quadrature_points);
6391 this_quadrature_points_data.size() != this->n_quadrature_points)
6392 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6396 this->jacobian = this_jacobian_data.data();
6397 this->J_value = this_J_value_data.data();
6398 this->jacobian_gradients = this_jacobian_gradients_data.data();
6399 this->jacobian_gradients_non_inverse =
6400 this_jacobian_gradients_non_inverse_data.data();
6404 for (
unsigned int v = 0; v < n_lanes; ++v)
6411 const unsigned int cell_batch_index =
cell_index / n_lanes;
6412 const unsigned int offsets =
6413 this->mapping_data->data_index_offsets[cell_batch_index];
6414 const unsigned int lane =
cell_index % n_lanes;
6416 if (this->cell_type <=
6420 for (
unsigned int q = 0; q < 2; ++q)
6421 for (
unsigned int i = 0; i < dim; ++i)
6422 for (
unsigned int j = 0; j < dim; ++j)
6423 this_jacobian_data[q][i][j][v] =
6424 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6426 const unsigned int q = 0;
6428 this_J_value_data[q][v] =
6429 this->mapping_data->JxW_values[offsets + q][lane];
6431 const auto &update_flags_cells =
6436 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6437 for (
unsigned int j = 0; j < dim; ++j)
6438 this_jacobian_gradients_data[q][i][j][v] =
6440 ->jacobian_gradients[0][offsets + q][i][j][lane];
6442 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6443 for (
unsigned int j = 0; j < dim; ++j)
6444 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6446 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6451 for (
unsigned int i = 0; i < dim; ++i)
6452 this_quadrature_points_data[q][i][v] =
6453 this->mapping_data->quadrature_points
6455 ->quadrature_point_offsets[cell_batch_index] +
6461 const auto cell_type =
6465 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6467 const unsigned int q_src =
6473 this_J_value_data[q][v] =
6474 this->mapping_data->JxW_values[offsets + q_src][lane];
6476 for (
unsigned int i = 0; i < dim; ++i)
6477 for (
unsigned int j = 0; j < dim; ++j)
6478 this_jacobian_data[q][i][j][v] =
6480 ->jacobians[0][offsets + q_src][i][j][lane];
6482 const auto &update_flags_cells =
6487 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6488 for (
unsigned int j = 0; j < dim; ++j)
6489 this_jacobian_gradients_data[q][i][j][v] =
6491 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6493 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6494 for (
unsigned int j = 0; j < dim; ++j)
6495 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6497 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6510 this->mapping_data->quadrature_points
6512 ->quadrature_point_offsets[cell_batch_index] +
6516 this->mapping_data->jacobians[0][offsets + 1];
6518 for (
unsigned int d = 0;
d < dim; ++
d)
6521 static_cast<Number
>(
6522 this->descriptor->quadrature.point(q)[
d]);
6524 for (
unsigned int d = 0;
d < dim; ++
d)
6525 for (
unsigned int e = 0;
e < dim; ++
e)
6528 static_cast<Number
>(
6529 this->descriptor->quadrature.point(q)[
e]);
6531 for (
unsigned int i = 0; i < dim; ++i)
6532 this_quadrature_points_data[q][i][v] = point[i][lane];
6537 for (
unsigned int i = 0; i < dim; ++i)
6538 this_quadrature_points_data[q][i][v] =
6539 this->mapping_data->quadrature_points
6541 ->quadrature_point_offsets[cell_batch_index] +
6550 this->is_reinitialized =
true;
6551 this->dof_values_initialized =
false;
6552 this->values_quad_initialized =
false;
6553 this->gradients_quad_initialized =
false;
6554 this->hessians_quad_initialized =
false;
6565 typename VectorizedArrayType>
6566template <
bool level_dof_access>
6573 VectorizedArrayType>
::
6576 Assert(this->matrix_free ==
nullptr,
6577 ExcMessage(
"Cannot use initialization from cell iterator if "
6578 "initialized from MatrixFree object. Use variant for "
6579 "on the fly computation with arguments as for FEValues "
6582 this->mapped_geometry->reinit(
6584 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6585 if (level_dof_access)
6586 cell->get_mg_dof_indices(this->local_dof_indices);
6588 cell->get_dof_indices(this->local_dof_indices);
6591 this->is_reinitialized =
true;
6602 typename VectorizedArrayType>
6609 VectorizedArrayType>
::
6612 Assert(this->matrix_free == 0,
6613 ExcMessage(
"Cannot use initialization from cell iterator if "
6614 "initialized from MatrixFree object. Use variant for "
6615 "on the fly computation with arguments as for FEValues "
6618 this->mapped_geometry->reinit(cell);
6621 this->is_reinitialized =
true;
6632 typename VectorizedArrayType>
6639 VectorizedArrayType>::
6643 Assert(this->dof_values_initialized ==
true,
6646 evaluate(this->values_dofs, evaluation_flags);
6656 typename VectorizedArrayType>
6663 VectorizedArrayType>::
6664 evaluate(
const VectorizedArrayType *values_array,
6667 const bool hessians_on_general_cells =
6671 if (hessians_on_general_cells)
6674 if (this->data->element_type ==
6680 if constexpr (fe_degree > -1)
6683 template run<fe_degree, n_q_points_1d>(n_components,
6684 evaluation_flag_actual,
6692 evaluation_flag_actual,
6693 const_cast<VectorizedArrayType *
>(values_array),
6698 this->values_quad_initialized =
6700 this->gradients_quad_initialized =
6702 this->hessians_quad_initialized =
6713 template <
typename Number,
6714 typename VectorizedArrayType,
6715 typename VectorType,
6716 typename EvaluatorType,
6717 std::enable_if_t<internal::has_begin<VectorType> &&
6719 VectorType> * =
nullptr>
6720 VectorizedArrayType *
6721 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6727 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6728 const auto &dof_info = fe_eval.get_dof_info();
6735 if (std::is_same_v<typename VectorType::value_type, Number> &&
6739 interleaved_contiguous &&
6740 reinterpret_cast<
std::size_t>(
6742 dof_info.dof_indices_contiguous
6743 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6744 [cell * VectorizedArrayType::size()]) %
6745 sizeof(VectorizedArrayType) ==
6748 return reinterpret_cast<VectorizedArrayType *
>(
6752 [cell * VectorizedArrayType::size()] +
6754 [fe_eval.get_active_fe_index()]
6755 [fe_eval.get_first_selected_component()] *
6756 VectorizedArrayType::size());
6765 template <
typename Number,
6766 typename VectorizedArrayType,
6767 typename VectorType,
6768 typename EvaluatorType,
6769 std::enable_if_t<!internal::has_begin<VectorType> ||
6771 VectorType> * =
nullptr>
6772 VectorizedArrayType *
6773 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
6786 typename VectorizedArrayType>
6787template <
typename VectorType>
6794 VectorizedArrayType>::
6795 gather_evaluate(
const VectorType &input_vector,
6798 const VectorizedArrayType *src_ptr =
6799 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
6800 *
this, input_vector);
6801 if (src_ptr !=
nullptr)
6802 evaluate(src_ptr, evaluation_flag);
6805 this->read_dof_values(input_vector);
6806 evaluate(this->begin_dof_values(), evaluation_flag);
6817 typename VectorizedArrayType>
6824 VectorizedArrayType>::
6827 integrate(integration_flag, this->values_dofs);
6830 this->dof_values_initialized =
true;
6841 typename VectorizedArrayType>
6848 VectorizedArrayType>::
6850 VectorizedArrayType *values_array,
6851 const bool sum_into_values_array)
6855 Assert(this->values_quad_submitted ==
true,
6858 Assert(this->gradients_quad_submitted ==
true,
6861 Assert(this->hessians_quad_submitted ==
true,
6864 Assert(this->matrix_free !=
nullptr ||
6865 this->mapped_geometry->is_initialized(),
6871 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
6872 "EvaluationFlags::hessians are supported."));
6878 unsigned int size = n_components * dim * n_q_points;
6881 for (
unsigned int i = 0; i < size; ++i)
6882 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
6886 for (
unsigned int i = 0; i < size; ++i)
6887 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
6892 if (n_components == dim &&
6893 this->data->element_type ==
6897 this->divergence_is_requested ==
false)
6899 unsigned int size = n_components * n_q_points;
6902 for (
unsigned int i = 0; i < size; ++i)
6903 this->values_quad[i] += this->values_from_gradients_quad[i];
6907 for (
unsigned int i = 0; i < size; ++i)
6908 this->values_quad[i] = this->values_from_gradients_quad[i];
6913 if constexpr (fe_degree > -1)
6916 template run<fe_degree, n_q_points_1d>(n_components,
6917 integration_flag_actual,
6920 sum_into_values_array);
6926 integration_flag_actual,
6929 sum_into_values_array);
6933 this->dof_values_initialized =
true;
6944 typename VectorizedArrayType>
6945template <
typename VectorType>
6952 VectorizedArrayType>::
6954 VectorType &destination)
6956 VectorizedArrayType *dst_ptr =
6957 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
6958 *
this, destination);
6959 if (dst_ptr !=
nullptr)
6960 integrate(integration_flag, dst_ptr,
true);
6963 integrate(integration_flag, this->begin_dof_values());
6964 this->distribute_local_to_global(destination);
6975 typename VectorizedArrayType>
6982 VectorizedArrayType>::dof_indices()
const
6984 return {0U, dofs_per_cell};
6998 typename VectorizedArrayType>
7004 VectorizedArrayType>::
7007 const bool is_interior_face,
7008 const unsigned int dof_no,
7009 const unsigned int quad_no,
7010 const unsigned int first_selected_component,
7011 const unsigned int active_fe_index,
7012 const unsigned int active_quad_index,
7013 const unsigned int face_type)
7014 : BaseClass(matrix_free,
7016 first_selected_component,
7024 , dofs_per_component(this->data->dofs_per_component_on_cell)
7025 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7026 , n_q_points(this->n_quadrature_points)
7036 typename VectorizedArrayType>
7042 VectorizedArrayType>::
7045 const std::pair<unsigned int, unsigned int> &range,
7046 const bool is_interior_face,
7047 const unsigned int dof_no,
7048 const unsigned int quad_no,
7049 const unsigned int first_selected_component)
7054 first_selected_component,
7055 matrix_free.get_face_active_fe_index(range,
7057 numbers::invalid_unsigned_int,
7058 matrix_free.get_face_info(range.
first).face_type)
7068 typename VectorizedArrayType>
7075 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7077 Assert(this->mapped_geometry ==
nullptr,
7078 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7079 " Integer indexing is not possible"));
7080 if (this->mapped_geometry !=
nullptr)
7083 this->cell = face_index;
7084 this->dof_access_index =
7085 this->is_interior_face() ?
7093 this->matrix_free->get_task_info().boundary_partition_data.back())
7094 Assert(this->is_interior_face(),
7096 "Boundary faces do not have a neighbor. When looping over "
7097 "boundary faces use FEFaceEvaluation with the parameter "
7098 "is_interior_face set to true. "));
7100 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7105 this->face_ids[i] = face_index * n_lanes + i;
7106 for (; i < n_lanes; ++i)
7109 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7110 const unsigned int offsets =
7111 this->mapping_data->data_index_offsets[face_index];
7112 this->J_value = &this->mapping_data->JxW_values[offsets];
7113 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7115 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7116 this->normal_x_jacobian =
7118 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7119 this->jacobian_gradients =
7120 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7122 this->jacobian_gradients_non_inverse =
7124 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7128 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7131 this->mapping_data->quadrature_point_offsets.size());
7133 this->mapping_data->quadrature_points.data() +
7134 this->mapping_data->quadrature_point_offsets[this->cell];
7138 this->is_reinitialized =
true;
7139 this->dof_values_initialized =
false;
7140 this->values_quad_initialized =
false;
7141 this->gradients_quad_initialized =
false;
7142 this->hessians_quad_initialized =
false;
7153 typename VectorizedArrayType>
7161 const unsigned int face_number)
7167 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7171 Assert(this->mapped_geometry ==
nullptr,
7172 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7173 " Integer indexing is not possible"));
7174 if (this->mapped_geometry !=
nullptr)
7179 .faces_by_cells_type[
cell_index][face_number];
7182 this->dof_access_index =
7185 if (this->is_interior_face() ==
false)
7190 for (
unsigned int i = 0; i < n_lanes; ++i)
7193 const unsigned int cell_this =
cell_index * n_lanes + i;
7195 unsigned int face_index =
7200 this->face_ids[i] = face_index;
7205 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7206 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7213 auto cell_m = faces.cells_interior[face_index % n_lanes];
7214 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7216 const bool face_identifies_as_interior = cell_m != cell_this;
7218 Assert(cell_m == cell_this || cell_p == cell_this,
7222 if (face_identifies_as_interior)
7224 this->cell_ids[i] = cell_m;
7225 this->face_numbers[i] = faces.interior_face_no;
7229 this->cell_ids[i] = cell_p;
7230 this->face_numbers[i] = faces.exterior_face_no;
7233 const bool orientation_interior_face = faces.face_orientation >= 8;
7234 unsigned int face_orientation = faces.face_orientation % 8;
7235 if (face_identifies_as_interior != orientation_interior_face)
7237 constexpr std::array<std::uint8_t, 8> table{
7238 {0, 1, 2, 3, 6, 5, 4, 7}};
7239 face_orientation = table[face_orientation];
7241 this->face_orientations[i] = face_orientation;
7246 this->face_orientations[0] = 0;
7247 this->face_numbers[0] = face_number;
7252 for (
unsigned int i = 0; i < n_lanes; ++i)
7253 this->cell_ids[i] =
cell_index * n_lanes + i;
7261 this->cell_ids[i] =
cell_index * n_lanes + i;
7262 for (; i < n_lanes; ++i)
7265 for (
unsigned int i = 0; i < n_lanes; ++i)
7272 const unsigned int offsets =
7274 .face_data_by_cells[this->quad_no]
7279 .face_data_by_cells[this->quad_no]
7280 .JxW_values.size());
7282 .face_data_by_cells[this->quad_no]
7283 .JxW_values[offsets];
7285 .face_data_by_cells[this->quad_no]
7286 .normal_vectors[offsets];
7288 .face_data_by_cells[this->quad_no]
7289 .jacobians[!this->is_interior_face()][offsets];
7290 this->normal_x_jacobian =
7292 .face_data_by_cells[this->quad_no]
7293 .normals_times_jacobians[!this->is_interior_face()][offsets];
7294 this->jacobian_gradients =
7295 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7297 this->jacobian_gradients_non_inverse =
7299 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7304 .face_data_by_cells[this->quad_no]
7305 .quadrature_point_offsets.empty() ==
false)
7307 const unsigned int index =
7311 .face_data_by_cells[this->quad_no]
7312 .quadrature_point_offsets.size());
7314 .face_data_by_cells[this->quad_no]
7315 .quadrature_points.data() +
7317 .face_data_by_cells[this->quad_no]
7318 .quadrature_point_offsets[
index];
7322 this->is_reinitialized =
true;
7323 this->dof_values_initialized =
false;
7324 this->values_quad_initialized =
false;
7325 this->gradients_quad_initialized =
false;
7326 this->hessians_quad_initialized =
false;
7337 typename VectorizedArrayType>
7344 VectorizedArrayType>::
7351 evaluate(this->values_dofs, evaluation_flag);
7361 typename VectorizedArrayType>
7368 VectorizedArrayType>::
7369 evaluate(
const VectorizedArrayType *values_array,
7372 Assert((evaluation_flag &
7375 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7376 "and EvaluationFlags::hessians are supported."));
7378 const bool hessians_on_general_cells =
7382 if (hessians_on_general_cells)
7385 if (this->data->element_type ==
7391 if constexpr (fe_degree > -1)
7393 template
run<fe_degree, n_q_points_1d>(n_components,
7394 evaluation_flag_actual,
7399 n_components, evaluation_flag_actual, values_array, *
this);
7402 this->values_quad_initialized =
7404 this->gradients_quad_initialized =
7406 this->hessians_quad_initialized =
7418 typename VectorizedArrayType>
7425 VectorizedArrayType>::
7432 project_to_face(this->values_dofs, evaluation_flag);
7442 typename VectorizedArrayType>
7449 VectorizedArrayType>::
7450 project_to_face(
const VectorizedArrayType *values_array,
7453 Assert((evaluation_flag &
7456 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7457 "and EvaluationFlags::hessians are supported."));
7459 const bool hessians_on_general_cells =
7463 if (hessians_on_general_cells)
7466 if (this->data->element_type ==
7472 if constexpr (fe_degree > -1)
7475 VectorizedArrayType>::template run<fe_degree>(n_components,
7476 evaluation_flag_actual,
7482 evaluation_flag_actual,
7496 typename VectorizedArrayType>
7503 VectorizedArrayType>::
7506 Assert((evaluation_flag &
7509 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7510 "and EvaluationFlags::hessians are supported."));
7512 const bool hessians_on_general_cells =
7516 if (hessians_on_general_cells)
7519 if (this->data->element_type ==
7525 if constexpr (fe_degree > -1)
7528 VectorizedArrayType>::template run<fe_degree>(n_components,
7529 evaluation_flag_actual,
7536 this->values_quad_initialized =
7538 this->gradients_quad_initialized =
7540 this->hessians_quad_initialized =
7552 typename VectorizedArrayType>
7559 VectorizedArrayType>::
7561 const bool sum_into_values)
7563 integrate(integration_flag, this->values_dofs, sum_into_values);
7566 this->dof_values_initialized =
true;
7577 typename VectorizedArrayType>
7584 VectorizedArrayType>::
7586 VectorizedArrayType *values_array,
7587 const bool sum_into_values)
7589 Assert((integration_flag &
7592 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7593 "and EvaluationFlags::hessians are supported."));
7599 unsigned int size = n_components * dim * n_q_points;
7602 for (
unsigned int i = 0; i < size; ++i)
7603 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7607 for (
unsigned int i = 0; i < size; ++i)
7608 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7613 if (this->data->element_type ==
7617 this->divergence_is_requested ==
false)
7619 unsigned int size = n_components * n_q_points;
7622 for (
unsigned int i = 0; i < size; ++i)
7623 this->values_quad[i] += this->values_from_gradients_quad[i];
7627 for (
unsigned int i = 0; i < size; ++i)
7628 this->values_quad[i] = this->values_from_gradients_quad[i];
7633 if constexpr (fe_degree > -1)
7635 template
run<fe_degree, n_q_points_1d>(n_components,
7636 integration_flag_actual,
7643 integration_flag_actual,
7656 typename VectorizedArrayType>
7663 VectorizedArrayType>::
7666 Assert((integration_flag &
7669 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7670 "and EvaluationFlags::hessians are supported."));
7676 unsigned int size = n_components * dim * n_q_points;
7679 for (
unsigned int i = 0; i < size; ++i)
7680 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7684 for (
unsigned int i = 0; i < size; ++i)
7685 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7690 if (this->data->element_type ==
7694 this->divergence_is_requested ==
false)
7696 unsigned int size = n_components * n_q_points;
7699 for (
unsigned int i = 0; i < size; ++i)
7700 this->values_quad[i] += this->values_from_gradients_quad[i];
7704 for (
unsigned int i = 0; i < size; ++i)
7705 this->values_quad[i] = this->values_from_gradients_quad[i];
7710 if constexpr (fe_degree > -1)
7713 VectorizedArrayType>::template run<fe_degree>(n_components,
7714 integration_flag_actual,
7730 typename VectorizedArrayType>
7737 VectorizedArrayType>::
7739 const bool sum_into_values)
7741 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
7744 this->dof_values_initialized =
true;
7755 typename VectorizedArrayType>
7762 VectorizedArrayType>::
7764 VectorizedArrayType *values_array,
7765 const bool sum_into_values)
7767 Assert((integration_flag &
7770 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7771 "and EvaluationFlags::hessians are supported."));
7778 if (this->data->element_type ==
7782 this->divergence_is_requested ==
false)
7785 if constexpr (fe_degree > -1)
7788 VectorizedArrayType>::template run<fe_degree>(n_components,
7789 integration_flag_actual,
7796 integration_flag_actual,
7809 typename VectorizedArrayType>
7810template <
typename VectorType>
7817 VectorizedArrayType>::
7818 gather_evaluate(
const VectorType &input_vector,
7821 Assert((evaluation_flag &
7824 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7825 "and EvaluationFlags::hessians are supported."));
7827 const auto shared_vector_data = internal::get_shared_vector_data(
7829 this->dof_access_index ==
7831 this->active_fe_index,
7834 if (this->data->data.front().fe_degree > 0 &&
7835 fast_evaluation_supported(this->data->data.front().fe_degree,
7836 this->data->data.front().n_q_points_1d) &&
7839 typename VectorType::value_type,
7840 VectorizedArrayType>::
7841 supports(evaluation_flag,
7843 internal::get_beginning<typename VectorType::value_type>(
7845 this->dof_info->index_storage_variants[this->dof_access_index]
7848 if constexpr (fe_degree > -1)
7852 typename VectorType::value_type,
7853 VectorizedArrayType>::template
run<fe_degree,
7857 internal::get_beginning<typename VectorType::value_type>(
7866 typename VectorType::value_type,
7867 VectorizedArrayType>::evaluate(n_components,
7869 internal::get_beginning<
7870 typename VectorType::value_type>(
7878 this->read_dof_values(input_vector);
7879 this->evaluate(evaluation_flag);
7884 this->gradients_quad_initialized =
7897 typename VectorizedArrayType>
7898template <
typename VectorType>
7906 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
7907 const bool integrate_gradients,
7908 VectorType &destination)
7915 integrate_scatter(flag, destination);
7925 typename VectorizedArrayType>
7926template <
typename VectorType>
7933 VectorizedArrayType>::
7935 VectorType &destination)
7937 Assert((this->dof_access_index ==
7939 this->is_interior_face() ==
false) ==
false,
7942 const auto shared_vector_data = internal::get_shared_vector_data(
7944 this->dof_access_index ==
7946 this->active_fe_index,
7949 if (this->data->data.front().fe_degree > 0 &&
7950 fast_evaluation_supported(this->data->data.front().fe_degree,
7951 this->data->data.front().n_q_points_1d) &&
7954 typename VectorType::value_type,
7955 VectorizedArrayType>::
7956 supports(integration_flag,
7958 internal::get_beginning<typename VectorType::value_type>(
7960 this->dof_info->index_storage_variants[this->dof_access_index]
7963 if constexpr (fe_degree > -1)
7967 typename VectorType::value_type,
7968 VectorizedArrayType>::template
run<fe_degree,
7972 internal::get_beginning<typename VectorType::value_type>(
7981 typename VectorType::value_type,
7982 VectorizedArrayType>::integrate(n_components,
7984 internal::get_beginning<
7985 typename VectorType::value_type>(
7993 integrate(integration_flag);
7994 this->distribute_local_to_global(destination);
8005 typename VectorizedArrayType>
8012 VectorizedArrayType>::dof_indices()
const
8014 return {0U, dofs_per_cell};
8024 typename VectorizedArrayType>
8031 VectorizedArrayType>::
8032 fast_evaluation_supported(
const unsigned int given_degree,
8033 const unsigned int given_n_q_points_1d)
8035 return fe_degree == -1 ?
8048 typename VectorizedArrayType>
8055 VectorizedArrayType>::
8056 fast_evaluation_supported(
const unsigned int given_degree,
8057 const unsigned int given_n_q_points_1d)
8059 return fe_degree == -1 ?
8072 typename VectorizedArrayType>
8079 VectorizedArrayType>::at_boundary()
const
8081 Assert(this->dof_access_index !=
8085 if (this->is_interior_face() ==
false)
8090 this->matrix_free->n_boundary_face_batches()))
8103 typename VectorizedArrayType>
8110 VectorizedArrayType>::boundary_id()
const
8112 Assert(this->dof_access_index !=
8129 typename VectorizedArrayType>
8137 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8139 return this->data->dofs_per_component_on_face;
8149 typename VectorizedArrayType>
8156 VectorizedArrayType>::get_dofs_projected_to_face()
8158 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
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & 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
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
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.
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)
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< VectorizedArrayType > &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 > &)