16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
92 typename VectorizedArrayType>
99 std::conditional_t<n_components_ == 1,
106 n_components_ == dim,
113 n_components_ == dim,
118 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
156 template <
typename VectorType>
159 const VectorType &src,
160 const unsigned int first_index = 0,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
191 template <
typename VectorType>
194 const VectorType &src,
195 const unsigned int first_index = 0,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
229 template <
typename VectorType>
233 const unsigned int first_index = 0,
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
274 template <
typename VectorType>
278 const unsigned int first_index = 0,
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
284 template <
typename VectorType>
288 const unsigned int first_index = 0,
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip())
const;
365 typename = std::enable_if_t<n_components == n_components_local>>
368 const unsigned int q_point);
419 template <
int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
423 const unsigned int q_point);
442 const unsigned int q_point);
518 const unsigned int q_point);
527 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
546 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
549 const unsigned int q_point);
559 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
578 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
582 const unsigned int q_point);
592 template <
int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
610 template <
int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
614 const unsigned int q_point);
655 const unsigned int dof_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
663 const unsigned int face_type);
737 template <
typename VectorType,
typename VectorOperation>
741 const std::array<VectorType *, n_components_> &vectors,
744 n_components_> &vectors_sm,
745 const std::bitset<n_lanes> &mask,
746 const bool apply_constraints =
true)
const;
755 template <
typename VectorType,
typename VectorOperation>
759 const std::array<VectorType *, n_components_> &vectors,
762 n_components_> &vectors_sm,
763 const std::bitset<n_lanes> &mask)
const;
772 template <
typename VectorType,
typename VectorOperation>
776 const std::array<VectorType *, n_components_> &vectors)
const;
1380 typename VectorizedArrayType>
1385 VectorizedArrayType>
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1530 const std::pair<unsigned int, unsigned int> &range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1642 template <
bool level_dof_access>
1664 const unsigned int given_n_q_points_1d);
1707 template <
typename VectorType>
1737 VectorizedArrayType *values_array,
1738 const bool sum_into_values =
false);
1753 template <
typename VectorType>
1756 VectorType &output_vector);
1840 int n_q_points_1d = fe_degree + 1,
1841 int n_components_ = 1,
1842 typename Number = double,
1848 VectorizedArrayType>
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
2012 const std::pair<unsigned int, unsigned int> &range,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2029 reinit(
const unsigned int face_batch_number);
2039 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2046 const unsigned int given_n_q_points_1d);
2110 template <
typename VectorType>
2126 const bool sum_into_values =
false);
2139 VectorizedArrayType *values_array,
2140 const bool sum_into_values =
false);
2157 const bool sum_into_values =
false);
2165 VectorizedArrayType *values_array,
2166 const bool sum_into_values =
false);
2179 template <
typename VectorType>
2182 VectorType &output_vector);
2187 template <
typename VectorType>
2190 const bool integrate_gradients,
2191 VectorType &output_vector);
2267 namespace MatrixFreeFunctions
2271 template <
int dim,
int degree>
2280 template <
int degree>
2283 static constexpr unsigned int value = degree + 1;
2299 template <
bool is_face,
2302 typename VectorizedArrayType>
2305 extract_initialization_data(
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2312 const unsigned int active_fe_index_given,
2313 const unsigned int active_quad_index_given,
2314 const unsigned int face_type)
2317 InitializationData init_data;
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2330 active_fe_index_given :
2336 active_quad_index_given :
2337 std::min<unsigned int>(
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2369 typename VectorizedArrayType>
2374 VectorizedArrayType>::
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2387 internal::extract_initialization_data<is_face>(matrix_free,
2389 first_selected_component,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2402 this->set_data_pointers(scratch_data_array, n_components_);
2404 this->dof_info->start_components.back() == 1 ||
2405 static_cast<int>(n_components_) <=
2407 this->dof_info->start_components
2408 [this->dof_info->component_to_base_index[first_selected_component] +
2410 first_selected_component,
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2417 this->dof_info->start_components
2418 [this->dof_info->component_to_base_index[first_selected_component] +
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) +
")"));
2438 typename VectorizedArrayType>
2443 VectorizedArrayType>::
2449 const unsigned int first_selected_component,
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2461 first_selected_component)
2462 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2463 , matrix_free(nullptr)
2465 const unsigned int base_element_number =
2469 first_selected_component >=
2471 ExcMessage(
"The underlying element must at least contain as many "
2472 "components as requested by this class"));
2473 (void)base_element_number;
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2479 fe.component_to_base_index(first_selected_component).
first);
2481 this->set_data_pointers(scratch_data_array, n_components_);
2490 typename VectorizedArrayType>
2495 VectorizedArrayType>::
2500 VectorizedArrayType> &other)
2502 , scratch_data_array(other.matrix_free == nullptr ?
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(other.matrix_free)
2507 if (other.matrix_free ==
nullptr)
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2521 if constexpr (is_face ==
false)
2522 this->mapping_data = &this->mapped_geometry->get_data_storage();
2526 "face evaluators is not currently "
2532 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2534 this->mapped_geometry->get_data_storage().JxW_values.begin();
2535 this->jacobian_gradients =
2536 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2537 this->jacobian_gradients_non_inverse =
2538 this->mapped_geometry->get_data_storage()
2539 .jacobian_gradients_non_inverse[0]
2542 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2545 this->set_data_pointers(scratch_data_array, n_components_);
2554 typename VectorizedArrayType>
2559 VectorizedArrayType> &
2565 VectorizedArrayType> &other)
2568 if (matrix_free ==
nullptr)
2571 delete scratch_data_array;
2580 matrix_free = other.matrix_free;
2582 if (other.matrix_free ==
nullptr)
2590 this->mapped_geometry =
2591 std::make_shared<internal::MatrixFreeFunctions::
2592 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2597 if constexpr (is_face ==
false)
2598 this->mapping_data = &this->mapped_geometry->get_data_storage();
2602 "face evaluators is not currently "
2607 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2609 this->mapped_geometry->get_data_storage().JxW_values.begin();
2610 this->jacobian_gradients =
2611 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2612 this->jacobian_gradients_non_inverse =
2613 this->mapped_geometry->get_data_storage()
2614 .jacobian_gradients_non_inverse[0]
2617 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2624 this->set_data_pointers(scratch_data_array, n_components_);
2635 typename VectorizedArrayType>
2640 VectorizedArrayType>::~FEEvaluationBase()
2642 if (matrix_free !=
nullptr)
2653 delete scratch_data_array;
2664 typename VectorizedArrayType>
2669 Assert(matrix_free !=
nullptr,
2671 "FEEvaluation was not initialized with a MatrixFree object!"));
2672 return *matrix_free;
2681 template <
typename VectorType,
bool>
2682 struct ConstBlockVectorSelector;
2684 template <
typename VectorType>
2685 struct ConstBlockVectorSelector<
VectorType, true>
2687 using BaseVectorType =
const typename VectorType::BlockType;
2690 template <
typename VectorType>
2691 struct ConstBlockVectorSelector<
VectorType, false>
2693 using BaseVectorType =
typename VectorType::BlockType;
2699 template <
typename VectorType,
bool>
2700 struct BlockVectorSelector;
2702 template <
typename VectorType>
2705 using BaseVectorType =
typename ConstBlockVectorSelector<
2707 std::is_const_v<VectorType>>::BaseVectorType;
2709 static BaseVectorType *
2710 get_vector_component(VectorType &vec,
const unsigned int component)
2713 return &vec.block(component);
2717 template <
typename VectorType>
2718 struct BlockVectorSelector<
VectorType, false>
2722 static BaseVectorType *
2723 get_vector_component(VectorType &vec,
const unsigned int component)
2742 template <
typename VectorType>
2743 struct BlockVectorSelector<
std::vector<VectorType>, false>
2747 static BaseVectorType *
2748 get_vector_component(std::vector<VectorType> &vec,
2749 const unsigned int component)
2752 return &vec[component];
2756 template <
typename VectorType>
2757 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2761 static const BaseVectorType *
2762 get_vector_component(
const std::vector<VectorType> &vec,
2763 const unsigned int component)
2766 return &vec[component];
2770 template <
typename VectorType>
2771 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2775 static BaseVectorType *
2776 get_vector_component(std::vector<VectorType *> &vec,
2777 const unsigned int component)
2780 return vec[component];
2784 template <
typename VectorType>
2785 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2789 static const BaseVectorType *
2790 get_vector_component(
const std::vector<VectorType *> &vec,
2791 const unsigned int component)
2794 return vec[component];
2798 template <
typename VectorType, std::
size_t N>
2799 struct BlockVectorSelector<
std::array<VectorType *, N>, false>
2803 static BaseVectorType *
2804 get_vector_component(std::array<VectorType *, N> &vec,
2805 const unsigned int component)
2808 return vec[component];
2819 typename VectorizedArrayType>
2820template <
typename VectorType,
typename VectorOperation>
2825 const std::array<VectorType *, n_components_> &src,
2828 n_components_> &src_sm,
2829 const std::bitset<n_lanes> &
mask,
2830 const bool apply_constraints)
const
2835 if (this->matrix_free ==
nullptr)
2837 read_write_operation_global(operation, src);
2844 if (this->n_fe_components == 1)
2845 for (
unsigned int comp = 0; comp < n_components; ++comp)
2847 Assert(src[comp] !=
nullptr,
2848 ExcMessage(
"The finite element underlying this FEEvaluation "
2849 "object is scalar, but you requested " +
2850 std::to_string(n_components) +
2851 " components via the template argument in "
2852 "FEEvaluation. In that case, you must pass an "
2853 "std::vector<VectorType> or a BlockVector to " +
2854 "read_dof_values and distribute_local_to_global."));
2866 const bool accesses_exterior_dofs =
2867 this->dof_access_index ==
2869 this->is_interior_face() ==
false;
2879 bool is_contiguous =
true;
2881 if (accesses_exterior_dofs)
2883 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2884 const unsigned int n_filled_lanes =
2889 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2890 if (
mask[v] ==
true &&
2893 [cells[v] / n_lanes] <
2896 is_contiguous = false;
2900 this->dof_access_index :
2901 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2902 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2905 is_contiguous =
false;
2910 read_write_operation_contiguous(operation, src, src_sm,
mask);
2917 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2919 const bool masking_is_active =
mask.count() < n_lanes;
2920 if (masking_is_active)
2921 for (
unsigned int v = 0; v < n_lanes; ++v)
2922 if (
mask[v] ==
false)
2925 bool has_hn_constraints =
false;
2927 if (is_face ==
false)
2933 [this->first_selected_component])
2934 for (
unsigned int v = 0; v < n_lanes; ++v)
2939 has_hn_constraints = true;
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2945 const bool use_vectorized_path =
2946 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2948 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
2949 std::array<VectorizedArrayType *, n_components> values_dofs;
2950 for (
unsigned int c = 0; c < n_components; ++c)
2951 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2952 c * dofs_per_component;
2956 [is_face ? this->dof_access_index :
2957 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2958 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2959 IndexStorageVariants::interleaved &&
2960 use_vectorized_path)
2962 const unsigned int *dof_indices =
2964 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2968 [this->first_selected_component] *
2971 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2972 if (n_components == 1 || this->n_fe_components == 1)
2973 for (
unsigned int comp = 0; comp < n_components; ++comp)
2975 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2978 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2980 if (n_components == 1 || this->n_fe_components == 1)
2981 for (
unsigned int i = 0; i < dofs_per_component;
2982 ++i, dof_indices += n_lanes)
2983 for (
unsigned int comp = 0; comp < n_components; ++comp)
2984 operation.process_dof_gather(dof_indices,
2988 values_dofs[comp][i],
2991 for (
unsigned int comp = 0; comp < n_components; ++comp)
2992 for (
unsigned int i = 0; i < dofs_per_component;
2993 ++i, dof_indices += n_lanes)
2994 operation.process_dof_gather(dof_indices,
2998 values_dofs[comp][i],
3005 std::array<const unsigned int *, n_lanes> dof_indices;
3006 dof_indices.fill(
nullptr);
3011 bool has_constraints =
false;
3012 const unsigned int n_components_read =
3013 this->n_fe_components > 1 ? n_components : 1;
3017 for (
unsigned int v = 0; v < n_lanes; ++v)
3023 const std::pair<unsigned int, unsigned int> *my_index_start =
3024 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3025 this->first_selected_component];
3030 if (my_index_start[n_components_read].
second !=
3031 my_index_start[0].
second)
3032 has_constraints =
true;
3035 dof_info.
dof_indices.data() + my_index_start[0].first;
3040 for (
unsigned int v = 0; v < n_lanes; ++v)
3045 const std::pair<unsigned int, unsigned int> *my_index_start =
3046 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3047 this->first_selected_component];
3048 if (my_index_start[n_components_read].
second !=
3049 my_index_start[0].
second)
3050 has_constraints =
true;
3057 dof_info.hanging_node_constraint_masks_comp
3058 [this->active_fe_index][this->first_selected_component])
3059 has_hn_constraints = true;
3062 my_index_start[0].
first ||
3065 my_index_start[0].
first,
3068 dof_info.
dof_indices.data() + my_index_start[0].first;
3072 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3073 return i != numbers::invalid_unsigned_int;
3075 for (
unsigned int comp = 0; comp < n_components; ++comp)
3076 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3077 operation.process_empty(values_dofs[comp][i]);
3081 if (!has_constraints && apply_constraints)
3083 if (n_components == 1 || this->n_fe_components == 1)
3085 for (
unsigned int v = 0; v < n_lanes; ++v)
3090 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3091 for (
unsigned int comp = 0; comp < n_components; ++comp)
3092 operation.process_dof(dof_indices[v][i],
3094 values_dofs[comp][i][v]);
3099 for (
unsigned int comp = 0; comp < n_components; ++comp)
3100 for (
unsigned int v = 0; v < n_lanes; ++v)
3105 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3106 operation.process_dof(
3107 dof_indices[v][comp * dofs_per_component + i],
3109 values_dofs[comp][i][v]);
3121 for (
unsigned int v = 0; v < n_lanes; ++v)
3127 const unsigned int cell_dof_index =
3128 cell_index * this->n_fe_components + this->first_selected_component;
3129 const unsigned int n_components_read =
3130 this->n_fe_components > 1 ? n_components : 1;
3131 unsigned int index_indicators =
3133 unsigned int next_index_indicators =
3134 dof_info.
row_starts[cell_dof_index + 1].second;
3138 if (apply_constraints ==
false &&
3139 (dof_info.
row_starts[cell_dof_index].second !=
3140 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3146 dof_info.hanging_node_constraint_masks_comp
3147 [this->active_fe_index][this->first_selected_component])))
3156 [this->first_selected_component] +
3158 next_index_indicators = index_indicators;
3161 if (n_components == 1 || this->n_fe_components == 1)
3163 unsigned int ind_local = 0;
3164 for (; index_indicators != next_index_indicators; ++index_indicators)
3166 const std::pair<unsigned short, unsigned short> indicator =
3169 for (
unsigned int j = 0; j < indicator.first; ++j)
3170 for (
unsigned int comp = 0; comp < n_components; ++comp)
3171 operation.process_dof(dof_indices[v][j],
3173 values_dofs[comp][ind_local + j][v]);
3175 ind_local += indicator.first;
3176 dof_indices[v] += indicator.first;
3180 Number
value[n_components];
3181 for (
unsigned int comp = 0; comp < n_components; ++comp)
3182 operation.pre_constraints(values_dofs[comp][ind_local][v],
3185 const Number *data_val =
3187 const Number *end_pool =
3189 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3190 for (
unsigned int comp = 0; comp < n_components; ++comp)
3191 operation.process_constraint(*dof_indices[v],
3196 for (
unsigned int comp = 0; comp < n_components; ++comp)
3197 operation.post_constraints(
value[comp],
3198 values_dofs[comp][ind_local][v]);
3204 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3205 for (
unsigned int comp = 0; comp < n_components; ++comp)
3206 operation.process_dof(*dof_indices[v],
3208 values_dofs[comp][ind_local][v]);
3217 for (
unsigned int comp = 0; comp < n_components; ++comp)
3219 unsigned int ind_local = 0;
3222 for (; index_indicators != next_index_indicators;
3225 const std::pair<unsigned short, unsigned short> indicator =
3229 for (
unsigned int j = 0; j < indicator.first; ++j)
3230 operation.process_dof(dof_indices[v][j],
3232 values_dofs[comp][ind_local + j][v]);
3233 ind_local += indicator.first;
3234 dof_indices[v] += indicator.first;
3239 operation.pre_constraints(values_dofs[comp][ind_local][v],
3242 const Number *data_val =
3244 const Number *end_pool =
3247 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3248 operation.process_constraint(*dof_indices[v],
3253 operation.post_constraints(
value,
3254 values_dofs[comp][ind_local][v]);
3261 for (; ind_local < dofs_per_component;
3262 ++dof_indices[v], ++ind_local)
3265 operation.process_dof(*dof_indices[v],
3267 values_dofs[comp][ind_local][v]);
3270 if (apply_constraints ==
true && comp + 1 < n_components)
3271 next_index_indicators =
3272 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3284 typename VectorizedArrayType>
3285template <
typename VectorType,
typename VectorOperation>
3290 const std::array<VectorType *, n_components_> &src)
const
3294 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3295 unsigned int index = this->first_selected_component * dofs_per_component;
3296 for (
unsigned int comp = 0; comp < n_components; ++comp)
3298 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3300 operation.process_empty(
3301 this->values_dofs[comp * dofs_per_component + i]);
3302 operation.process_dof_global(
3303 local_dof_indices[this->
data->lexicographic_numbering[index]],
3305 this->values_dofs[comp * dofs_per_component + i][0]);
3316 typename VectorizedArrayType>
3317template <
typename VectorType,
typename VectorOperation>
3322 const std::array<VectorType *, n_components_> &src,
3325 n_components_> &vectors_sm,
3326 const std::bitset<n_lanes> &
mask)
const
3335 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3338 is_face ? this->dof_access_index :
3340 const unsigned int n_active_lanes =
mask.count();
3343 const std::vector<unsigned int> &dof_indices_cont =
3346 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3347 std::array<VectorizedArrayType *, n_components> values_dofs{{
nullptr}};
3348 for (
unsigned int c = 0; c < n_components; ++c)
3349 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3350 c * dofs_per_component;
3354 const bool accesses_exterior_dofs =
3355 this->dof_access_index ==
3357 this->is_interior_face() ==
false;
3363 interleaved_contiguous &&
3364 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3366 const unsigned int dof_index =
3367 dof_indices_cont[this->cell * n_lanes] +
3370 [this->first_selected_component] *
3372 if (n_components == 1 || this->n_fe_components == 1)
3373 for (
unsigned int comp = 0; comp < n_components; ++comp)
3374 operation.process_dofs_vectorized(dofs_per_component,
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3388 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3392 const unsigned int n_filled_lanes =
3395 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3396 n_active_lanes == n_lanes &&
3397 !accesses_exterior_dofs;
3399 if (vectors_sm[0] !=
nullptr)
3401 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3402 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3405 const auto upper_bound =
3406 std::min<unsigned int>(n_filled_lanes, n_lanes);
3407 for (
unsigned int v = 0; v < upper_bound; ++v)
3409 if (
mask[v] ==
false)
3411 vector_ptrs[v] =
nullptr;
3431 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3432 vectors_sm[comp]->operator[](temp.first).
data() + temp.second +
3434 [this->active_fe_index][this->first_selected_component]);
3436 vector_ptrs[v] =
nullptr;
3438 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3439 vector_ptrs[v] =
nullptr;
3444 if (use_vectorized_path)
3446 if (n_components == 1 || this->n_fe_components == 1)
3448 for (
unsigned int comp = 0; comp < n_components; ++comp)
3450 auto vector_ptrs = compute_vector_ptrs(comp);
3451 operation.process_dofs_vectorized_transpose(
3460 auto vector_ptrs = compute_vector_ptrs(0);
3461 operation.process_dofs_vectorized_transpose(dofs_per_component *
3469 for (
unsigned int comp = 0; comp < n_components; ++comp)
3471 auto vector_ptrs = compute_vector_ptrs(
3472 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3474 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3475 operation.process_empty(values_dofs[comp][i]);
3477 if (n_components == 1 || this->n_fe_components == 1)
3479 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3480 if (
mask[v] ==
true)
3481 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3482 operation.process_dof(vector_ptrs[v][i],
3483 values_dofs[comp][i][v]);
3487 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3488 if (
mask[v] ==
true)
3489 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3490 operation.process_dof(
3491 vector_ptrs[v][i + comp * dofs_per_component],
3492 values_dofs[comp][i][v]);
3498 std::array<unsigned int, n_lanes> dof_indices{
3501 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3505 if (
mask[v] ==
true)
3507 dof_indices_cont[cells[v]] +
3510 [this->first_selected_component] *
3516 if (use_vectorized_path)
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (
unsigned int comp = 0; comp < n_components; ++comp)
3524 operation.process_dofs_vectorized_transpose(dofs_per_component,
3530 operation.process_dofs_vectorized_transpose(dofs_per_component *
3539 interleaved_contiguous_strided)
3541 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (
unsigned int comp = 0; comp < n_components; ++comp)
3545 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3546 src[comp]->
begin());
3549 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3551 if (n_components == 1 || this->n_fe_components == 1)
3552 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3554 for (
unsigned int comp = 0; comp < n_components; ++comp)
3555 operation.process_dof_gather(dof_indices.data(),
3558 src_ptrs[comp] + i * n_lanes,
3559 values_dofs[comp][i],
3563 for (
unsigned int comp = 0; comp < n_components; ++comp)
3564 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3566 operation.process_dof_gather(
3569 (comp * dofs_per_component + i) * n_lanes,
3570 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3571 values_dofs[comp][i],
3579 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3581 std::array<typename VectorType::value_type *, n_components> src_ptrs{
3583 if (n_components == 1 || this->n_fe_components == 1)
3584 for (
unsigned int comp = 0; comp < n_components; ++comp)
3585 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3586 src[comp]->
begin());
3589 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3591 const unsigned int *offsets =
3593 if (n_components == 1 || this->n_fe_components == 1)
3594 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3596 for (
unsigned int comp = 0; comp < n_components; ++comp)
3597 operation.process_dof_gather(dof_indices.data(),
3601 values_dofs[comp][i],
3604 for (
unsigned int v = 0; v < n_lanes; ++v)
3605 dof_indices[v] += offsets[v];
3608 for (
unsigned int comp = 0; comp < n_components; ++comp)
3609 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3611 operation.process_dof_gather(dof_indices.data(),
3615 values_dofs[comp][i],
3618 for (
unsigned int v = 0; v < n_lanes; ++v)
3619 dof_indices[v] += offsets[v];
3624 for (
unsigned int comp = 0; comp < n_components; ++comp)
3626 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3627 operation.process_empty(values_dofs[comp][i]);
3628 if (accesses_exterior_dofs)
3630 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3631 if (
mask[v] ==
true)
3634 [ind][cells[v] / VectorizedArrayType::size()] ==
3638 if (n_components == 1 || this->n_fe_components == 1)
3640 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3641 operation.process_dof(dof_indices[v] + i,
3643 values_dofs[comp][i][v]);
3647 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3651 values_dofs[comp][i][v]);
3656 const unsigned int offset =
3659 if (n_components == 1 || this->n_fe_components == 1)
3661 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3662 operation.process_dof(dof_indices[v] + i * offset,
3664 values_dofs[comp][i][v]);
3668 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3669 operation.process_dof(
3671 (i + comp * dofs_per_component) * offset,
3673 values_dofs[comp][i][v]);
3684 if (n_components == 1 || this->n_fe_components == 1)
3686 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3687 if (
mask[v] ==
true)
3688 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3689 operation.process_dof(dof_indices[v] + i,
3691 values_dofs[comp][i][v]);
3695 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3696 if (
mask[v] ==
true)
3697 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3698 operation.process_dof(dof_indices[v] + i +
3699 comp * dofs_per_component,
3701 values_dofs[comp][i][v]);
3706 const unsigned int *offsets =
3708 [ind][VectorizedArrayType::size() * this->cell];
3709 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3711 if (n_components == 1 || this->n_fe_components == 1)
3712 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3714 if (
mask[v] ==
true)
3715 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3716 operation.process_dof(dof_indices[v] + i * offsets[v],
3718 values_dofs[comp][i][v]);
3722 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3723 if (
mask[v] ==
true)
3724 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3725 operation.process_dof(
3727 (i + comp * dofs_per_component) * offsets[v],
3729 values_dofs[comp][i][v]);
3741 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3742 decltype(std::declval<VectorType>().begin())
3743 get_beginning(VectorType &vec)
3751 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3752 typename VectorType::value_type *
3753 get_beginning(VectorType &)
3759 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3761 const std::vector<ArrayView<const typename VectorType::value_type>> *
3762 get_shared_vector_data(VectorType *vec,
3763 const bool is_valid_mode_for_sm,
3764 const unsigned int active_fe_index,
3768 if (is_valid_mode_for_sm &&
3771 active_fe_index == 0)
3772 return &vec->shared_vector_data();
3778 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3780 const std::vector<ArrayView<const typename VectorType::value_type>> *
3781 get_shared_vector_data(VectorType *,
3789 template <
int n_components,
typename VectorType>
3791 std::array<
typename internal::BlockVectorSelector<
3796 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3800 get_vector_data(VectorType &src,
3801 const unsigned int first_index,
3802 const bool is_valid_mode_for_sm,
3803 const unsigned int active_fe_index,
3809 std::array<
typename internal::BlockVectorSelector<
3815 ArrayView<
const typename internal::BlockVectorSelector<
3821 for (
unsigned int d = 0;
d < n_components; ++
d)
3822 src_data.first[d] = internal::BlockVectorSelector<
3828 for (
unsigned int d = 0;
d < n_components; ++
d)
3829 src_data.second[d] = get_shared_vector_data(
3830 const_cast<typename internal::BlockVectorSelector<
3831 std::remove_const_t<VectorType>,
3833 *>(src_data.first[d]),
3834 is_valid_mode_for_sm,
3848 typename VectorizedArrayType>
3853 if (this->dof_info ==
nullptr ||
3855 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3856 this->dof_info->hanging_node_constraint_masks_comp
3857 [this->active_fe_index][this->first_selected_component] ==
false)
3860 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3861 constraint_mask{{internal::MatrixFreeFunctions::
3862 unconstrained_compressed_constraint_kind}};
3864 bool hn_available =
false;
3866 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3868 for (
unsigned int v = 0; v < n_lanes; ++v)
3880 constraint_mask[v] =
mask;
3886 if (hn_available ==
false)
3891 this->
data->data.front().fe_degree,
3892 this->get_shape_info(),
3904 typename VectorizedArrayType>
3905template <
typename VectorType>
3909 const unsigned int first_index,
3910 const std::bitset<n_lanes> &
mask)
3912 const auto src_data = internal::get_vector_data<n_components_>(
3915 this->dof_info !=
nullptr &&
3916 this->dof_access_index ==
3918 this->active_fe_index,
3922 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3924 apply_hanging_node_constraints(
false);
3928 this->dof_values_initialized =
true;
3938 typename VectorizedArrayType>
3939template <
typename VectorType>
3943 const unsigned int first_index,
3944 const std::bitset<n_lanes> &
mask)
3946 const auto src_data = internal::get_vector_data<n_components_>(
3949 this->dof_access_index ==
3951 this->active_fe_index,
3955 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3959 this->dof_values_initialized =
true;
3969 typename VectorizedArrayType>
3970template <
typename VectorType>
3974 const unsigned int first_index,
3975 const std::bitset<n_lanes> &
mask)
const
3979 Assert(this->dof_values_initialized ==
true,
3983 apply_hanging_node_constraints(
true);
3985 const auto dst_data = internal::get_vector_data<n_components_>(
3988 this->dof_access_index ==
3990 this->active_fe_index,
3995 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
4004 typename VectorizedArrayType>
4005template <
typename VectorType>
4009 const unsigned int first_index,
4010 const std::bitset<n_lanes> &
mask)
const
4014 Assert(this->dof_values_initialized ==
true,
4018 const auto dst_data = internal::get_vector_data<n_components_>(
4021 this->dof_access_index ==
4023 this->active_fe_index,
4027 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4036 typename VectorizedArrayType>
4037template <
typename VectorType>
4041 const unsigned int first_index,
4042 const std::bitset<n_lanes> &
mask)
const
4046 Assert(this->dof_values_initialized ==
true,
4050 const auto dst_data = internal::get_vector_data<n_components_>(
4053 this->dof_access_index ==
4055 this->active_fe_index,
4059 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4072 typename VectorizedArrayType>
4078 VectorizedArrayType>::value_type
4083 if constexpr (n_components == 1)
4084 return this->values_dofs[dof];
4087 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4088 Tensor<1, n_components_, VectorizedArrayType> return_value;
4089 for (
unsigned int comp = 0; comp < n_components; ++comp)
4090 return_value[comp] = this->values_dofs[comp * dofs + dof];
4091 return return_value;
4101 typename VectorizedArrayType>
4107 VectorizedArrayType>::value_type
4109 get_value(
const unsigned int q_point)
const
4113 Assert(this->values_quad_initialized ==
true,
4118 if constexpr (n_components == 1)
4119 return this->values_quad[q_point];
4122 if (n_components == dim &&
4123 this->
data->element_type ==
4129 Assert(this->values_quad_initialized ==
true,
4134 Assert(this->J_value !=
nullptr,
4137 const std::size_t nqp = this->n_quadrature_points;
4145 const VectorizedArrayType inv_det =
4146 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4147 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4148 this->jacobian[0][2][2];
4151 for (
unsigned int comp = 0; comp < n_components; ++comp)
4152 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4153 jac[comp][comp] * inv_det;
4160 this->jacobian[q_point] :
4169 const VectorizedArrayType inv_det =
4170 (is_face && dim == 2 && this->get_face_no() < 2) ?
4174 for (
unsigned int comp = 0; comp < n_components; ++comp)
4176 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4177 for (
unsigned int e = 1;
e < dim; ++
e)
4179 this->values_quad[e * nqp + q_point] * jac[comp][e];
4180 value_out[comp] *= inv_det;
4187 const std::size_t nqp = this->n_quadrature_points;
4189 for (
unsigned int comp = 0; comp < n_components; ++comp)
4190 return_value[comp] = this->values_quad[comp * nqp + q_point];
4191 return return_value;
4202 typename VectorizedArrayType>
4208 VectorizedArrayType>::gradient_type
4214 Assert(this->gradients_quad_initialized ==
true,
4219 Assert(this->jacobian !=
nullptr,
4221 "update_gradients"));
4222 const std::size_t nqp = this->n_quadrature_points;
4224 if constexpr (n_components == dim && dim > 1)
4226 if (this->
data->element_type ==
4232 Assert(this->gradients_quad_initialized ==
true,
4237 Assert(this->jacobian !=
nullptr,
4239 "update_gradients"));
4240 const std::size_t nqp = this->n_quadrature_points;
4241 const std::size_t nqp_d = nqp * dim;
4244 this->gradients_quad + q_point * dim;
4255 const VectorizedArrayType inv_det =
4256 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4257 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4258 this->jacobian[0][2][2];
4261 for (
unsigned int d = 0;
d < dim; ++
d)
4262 for (
unsigned int comp = 0; comp < n_components; ++comp)
4263 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4265 (jac[comp][comp] * inv_det);
4277 const VectorizedArrayType inv_det =
4278 (is_face && dim == 2 && this->get_face_no() < 2) ?
4282 VectorizedArrayType tmp[dim][dim];
4284 for (
unsigned int d = 0;
d < dim; ++
d)
4285 for (
unsigned int e = 0;
e < dim; ++
e)
4288 for (
unsigned int f = 1; f < dim; ++f)
4289 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4291 for (
unsigned int comp = 0; comp < n_components; ++comp)
4292 for (
unsigned int d = 0;
d < dim; ++
d)
4294 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4295 for (
unsigned int f = 1; f < dim; ++f)
4296 res += jac[comp][f] * tmp[d][f];
4298 grad_out[comp][
d] = res * inv_det;
4309 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4311 "update_hessians"));
4313 const auto jac_grad =
4314 this->jacobian_gradients_non_inverse[q_point];
4316 this->jacobian[q_point];
4320 const VectorizedArrayType inv_det =
4321 (is_face && dim == 2 && this->get_face_no() < 2) ?
4328 VectorizedArrayType tmp[dim][dim];
4329 for (
unsigned int d = 0;
d < dim; ++
d)
4330 for (
unsigned int e = 0;
e < dim; ++
e)
4333 for (
unsigned int f = 1; f < dim; ++f)
4334 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4339 for (
unsigned int d = 0;
d < dim; ++
d)
4341 for (
unsigned int e = 0;
e < dim; ++
e)
4343 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4344 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4345 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4348 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4350 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4355 for (
unsigned int d = 0;
d < dim; ++
d)
4356 for (
unsigned int e = 0;
e < dim; ++
e)
4358 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4359 for (
unsigned int f = 1; f < dim; ++f)
4360 res += tmp[f][d] * inv_t_jac[e][f];
4361 grad_out[
d][
e] = res;
4367 VectorizedArrayType tmp3[dim], tmp4[dim];
4368 for (
unsigned int d = 0;
d < dim; ++
d)
4370 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4371 for (
unsigned int e = 1;
e < dim; ++
e)
4372 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4374 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4375 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4376 for (
unsigned int d = 0;
d < dim; ++
d)
4378 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4379 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4381 for (
unsigned int d = 0;
d < dim; ++
d)
4383 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4384 for (
unsigned int e = 1;
e < dim; ++
e)
4385 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4388 VectorizedArrayType tmp2[dim];
4389 for (
unsigned int d = 0;
d < dim; ++
d)
4391 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4392 for (
unsigned e = 1;
e < dim; ++
e)
4394 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4397 for (
unsigned int d = 0;
d < dim; ++
d)
4398 for (
unsigned int e = 0;
e < dim; ++
e)
4400 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4404 grad_out[
d][
e] *= inv_det;
4415 for (
unsigned int comp = 0; comp < n_components; ++comp)
4416 for (
unsigned int d = 0;
d < dim; ++
d)
4418 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4419 this->jacobian[0][
d][
d];
4428 for (
unsigned int comp = 0; comp < n_components; ++comp)
4429 for (
unsigned int d = 0;
d < dim; ++
d)
4432 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4433 for (
unsigned int e = 1;
e < dim; ++
e)
4434 grad_out[comp][d] +=
4436 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4439 if constexpr (n_components == 1)
4451 typename VectorizedArrayType>
4457 VectorizedArrayType>::value_type
4464 Assert(this->gradients_quad_initialized ==
true,
4468 Assert(this->normal_x_jacobian !=
nullptr,
4470 "update_gradients"));
4472 const std::size_t nqp = this->n_quadrature_points;
4476 for (
unsigned int comp = 0; comp < n_components; ++comp)
4478 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4479 (this->normal_x_jacobian[0][dim - 1]);
4482 const std::size_t
index =
4484 for (
unsigned int comp = 0; comp < n_components; ++comp)
4486 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4487 this->normal_x_jacobian[
index][0];
4488 for (
unsigned int d = 1;
d < dim; ++
d)
4490 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4491 this->normal_x_jacobian[
index][
d];
4494 if constexpr (n_components == 1)
4506 template <
typename VectorizedArrayType>
4509 const VectorizedArrayType *
const hessians,
4511 VectorizedArrayType (&tmp)[1][1])
4513 tmp[0][0] = jac[0][0] *
hessians[0];
4516 template <
typename VectorizedArrayType>
4519 const VectorizedArrayType *
const hessians,
4520 const unsigned int nqp,
4521 VectorizedArrayType (&tmp)[2][2])
4523 for (
unsigned int d = 0;
d < 2; ++
d)
4531 template <
typename VectorizedArrayType>
4534 const VectorizedArrayType *
const hessians,
4535 const unsigned int nqp,
4536 VectorizedArrayType (&tmp)[3][3])
4538 for (
unsigned int d = 0;
d < 3; ++
d)
4559 typename VectorizedArrayType>
4564 VectorizedArrayType>::hessian_type
4570 Assert(this->hessians_quad_initialized ==
true,
4575 Assert(this->jacobian !=
nullptr,
4585 const std::size_t nqp = this->n_quadrature_points;
4586 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4591 for (
unsigned int comp = 0; comp < n_components; ++comp)
4593 for (
unsigned int d = 0;
d < dim; ++
d)
4594 hessian_out[comp][d][d] =
4595 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4596 (jac[
d][
d] * jac[
d][
d]);
4602 hessian_out[comp][0][1] =
4603 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4604 (jac[0][0] * jac[1][1]);
4607 hessian_out[comp][0][1] =
4608 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4609 (jac[0][0] * jac[1][1]);
4610 hessian_out[comp][0][2] =
4611 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4612 (jac[0][0] * jac[2][2]);
4613 hessian_out[comp][1][2] =
4614 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4615 (jac[1][1] * jac[2][2]);
4620 for (
unsigned int d = 0;
d < dim; ++
d)
4621 for (
unsigned int e = d + 1;
e < dim; ++
e)
4622 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4628 for (
unsigned int comp = 0; comp < n_components; ++comp)
4630 VectorizedArrayType tmp[dim][dim];
4631 internal::hessian_unit_times_jac(
4632 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4635 for (
unsigned int d = 0;
d < dim; ++
d)
4636 for (
unsigned int e = d;
e < dim; ++
e)
4638 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4639 for (
unsigned int f = 1; f < dim; ++f)
4640 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4647 for (
unsigned int d = 0;
d < dim; ++
d)
4648 for (
unsigned int e = d + 1;
e < dim; ++
e)
4649 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4655 const auto &jac_grad = this->jacobian_gradients[q_point];
4656 for (
unsigned int comp = 0; comp < n_components; ++comp)
4658 VectorizedArrayType tmp[dim][dim];
4659 internal::hessian_unit_times_jac(
4660 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4663 for (
unsigned int d = 0;
d < dim; ++
d)
4664 for (
unsigned int e = d;
e < dim; ++
e)
4666 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4667 for (
unsigned int f = 1; f < dim; ++f)
4668 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4672 for (
unsigned int d = 0;
d < dim; ++
d)
4673 for (
unsigned int e = 0;
e < dim; ++
e)
4674 hessian_out[comp][d][d] +=
4676 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4679 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4680 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4681 for (
unsigned int f = 0; f < dim; ++f)
4682 hessian_out[comp][d][e] +=
4683 jac_grad[count][f] *
4684 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4687 for (
unsigned int d = 0;
d < dim; ++
d)
4688 for (
unsigned int e = d + 1;
e < dim; ++
e)
4689 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4692 if constexpr (n_components == 1)
4693 return hessian_out[0];
4704 typename VectorizedArrayType>
4709 VectorizedArrayType>::gradient_type
4716 Assert(this->hessians_quad_initialized ==
true,
4727 const std::size_t nqp = this->n_quadrature_points;
4728 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4734 for (
unsigned int comp = 0; comp < n_components; ++comp)
4735 for (
unsigned int d = 0;
d < dim; ++
d)
4736 hessian_out[comp][d] =
4737 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4738 (jac[
d][
d] * jac[
d][
d]);
4743 for (
unsigned int comp = 0; comp < n_components; ++comp)
4747 VectorizedArrayType tmp[dim][dim];
4748 internal::hessian_unit_times_jac(
4749 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4753 for (
unsigned int d = 0;
d < dim; ++
d)
4755 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4756 for (
unsigned int f = 1; f < dim; ++f)
4757 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4764 const auto &jac_grad = this->jacobian_gradients[q_point];
4765 for (
unsigned int comp = 0; comp < n_components; ++comp)
4769 VectorizedArrayType tmp[dim][dim];
4770 internal::hessian_unit_times_jac(
4771 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4775 for (
unsigned int d = 0;
d < dim; ++
d)
4777 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4778 for (
unsigned int f = 1; f < dim; ++f)
4779 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4782 for (
unsigned int d = 0;
d < dim; ++
d)
4783 for (
unsigned int e = 0;
e < dim; ++
e)
4784 hessian_out[comp][d] +=
4786 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4790 if constexpr (n_components == 1)
4791 return hessian_out[0];
4802 typename VectorizedArrayType>
4807 VectorizedArrayType>::value_type
4814 Assert(this->hessians_quad_initialized ==
true,
4819 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4820 if constexpr (n_components == 1)
4822 VectorizedArrayType
sum = hess_diag[0];
4823 for (
unsigned int d = 1;
d < dim; ++
d)
4824 sum += hess_diag[d];
4830 for (
unsigned int comp = 0; comp < n_components; ++comp)
4832 laplacian_out[comp] = hess_diag[comp][0];
4833 for (
unsigned int d = 1;
d < dim; ++
d)
4834 laplacian_out[comp] += hess_diag[comp][d];
4836 return laplacian_out;
4846 typename VectorizedArrayType>
4851 VectorizedArrayType>::value_type
4857 Assert(this->hessians_quad_initialized ==
true,
4862 Assert(this->normal_x_jacobian !=
nullptr,
4864 "update_hessians"));
4868 const std::size_t nqp = this->n_quadrature_points;
4869 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4873 const auto nxj = this->normal_x_jacobian[0];
4875 for (
unsigned int comp = 0; comp < n_components; ++comp)
4877 for (
unsigned int d = 0;
d < dim; ++
d)
4878 hessian_out[comp] +=
4879 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4880 (nxj[
d]) * (nxj[d]);
4887 hessian_out[comp] +=
4888 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4892 hessian_out[comp] +=
4893 2. * this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4895 hessian_out[comp] +=
4896 2. * this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4898 hessian_out[comp] +=
4899 2. * this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4910 const auto normal = this->normal_vector(q_point);
4911 const auto hessian = get_hessian(q_point);
4913 if constexpr (n_components == 1)
4914 hessian_out[0] =
hessian * normal * normal;
4916 for (
unsigned int comp = 0; comp < n_components; ++comp)
4917 hessian_out[comp] =
hessian[comp] * normal * normal;
4919 if constexpr (n_components == 1)
4920 return hessian_out[0];
4931 typename VectorizedArrayType>
4938 this->dof_values_initialized =
true;
4940 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4942 for (
unsigned int comp = 0; comp < n_components; ++comp)
4943 if constexpr (n_components == 1)
4944 this->values_dofs[comp * dofs + dof] = val_in;
4946 this->values_dofs[comp * dofs + dof] = val_in[comp];
4955 typename VectorizedArrayType>
4958 submit_value(
const value_type val_in,
const unsigned int q_point)
4965 Assert(this->J_value !=
nullptr,
4970 this->values_quad_submitted =
true;
4973 const std::size_t nqp = this->n_quadrature_points;
4974 VectorizedArrayType *
values = this->values_quad + q_point;
4976 const VectorizedArrayType JxW =
4978 this->J_value[0] * this->quadrature_weights[q_point] :
4979 this->J_value[q_point];
4980 if constexpr (n_components == 1)
4981 values[0] = val_in * JxW;
4984 if (n_components == dim &&
4985 this->
data->element_type ==
4990 Assert(this->J_value !=
nullptr,
4996 this->values_quad_submitted =
true;
4999 VectorizedArrayType *
values = this->values_quad + q_point;
5000 const std::size_t nqp = this->n_quadrature_points;
5006 const VectorizedArrayType weight =
5007 this->quadrature_weights[q_point];
5009 for (
unsigned int comp = 0; comp < n_components; ++comp)
5010 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
5017 this->jacobian[q_point] :
5022 const VectorizedArrayType fac =
5024 this->quadrature_weights[q_point] :
5026 this->J_value[q_point] :
5027 this->J_value[0] * this->quadrature_weights[q_point]) *
5028 ((dim == 2 && this->get_face_no() < 2) ?
5037 for (
unsigned int comp = 0; comp < n_components; ++comp)
5039 values[comp * nqp] = val_in[0] * jac[0][comp];
5040 for (
unsigned int e = 1;
e < dim; ++
e)
5041 values[comp * nqp] += val_in[e] * jac[e][comp];
5042 values[comp * nqp] *= fac;
5047 for (
unsigned int comp = 0; comp < n_components; ++comp)
5048 values[comp * nqp] = val_in[comp] * JxW;
5058 typename VectorizedArrayType>
5059template <
int,
typename>
5063 const unsigned int q_point)
5065 static_assert(n_components == 1,
5066 "Do not try to modify the default template parameters used for"
5067 " selectively enabling this function via std::enable_if!");
5068 submit_value(val_in[0], q_point);
5077 typename VectorizedArrayType>
5080 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
5087 Assert(this->J_value !=
nullptr,
5089 "update_gradients"));
5090 Assert(this->jacobian !=
nullptr,
5092 "update_gradients"));
5095 this->gradients_quad_submitted =
true;
5098 if constexpr (dim > 1 && n_components == dim)
5100 if (this->
data->element_type ==
5110 Assert(this->J_value !=
nullptr,
5112 "update_gradients"));
5113 Assert(this->jacobian !=
nullptr,
5115 "update_gradients"));
5118 this->gradients_quad_submitted =
true;
5121 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5122 VectorizedArrayType *
values =
5123 this->values_from_gradients_quad + q_point;
5124 const std::size_t nqp = this->n_quadrature_points;
5125 const std::size_t nqp_d = nqp * dim;
5135 const VectorizedArrayType weight =
5136 this->quadrature_weights[q_point];
5137 for (
unsigned int d = 0;
d < dim; ++
d)
5138 for (
unsigned int comp = 0; comp < n_components; ++comp)
5139 gradients[comp * nqp_d + d] = grad_in[comp][d] *
5141 (jac[comp][comp] * weight);
5153 const VectorizedArrayType fac =
5155 this->quadrature_weights[q_point] :
5156 this->J_value[0] * this->quadrature_weights[q_point] *
5157 ((dim == 2 && this->get_face_no() < 2) ?
5162 VectorizedArrayType tmp[dim][dim];
5163 for (
unsigned int d = 0;
d < dim; ++
d)
5164 for (
unsigned int e = 0;
e < dim; ++
e)
5166 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
5167 for (
unsigned int f = 1; f < dim; ++f)
5168 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
5170 for (
unsigned int comp = 0; comp < n_components; ++comp)
5171 for (
unsigned int d = 0;
d < dim; ++
d)
5173 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
5174 for (
unsigned int f = 1; f < dim; ++f)
5175 res += jac[f][comp] * tmp[d][f];
5184 const auto jac_grad =
5185 this->jacobian_gradients_non_inverse[q_point];
5187 this->jacobian[q_point];
5191 const VectorizedArrayType fac =
5192 (!is_face) ? this->quadrature_weights[q_point] :
5193 this->J_value[q_point] *
5194 ((dim == 2 && this->get_face_no() < 2) ?
5203 VectorizedArrayType tmp3[dim], tmp4[dim];
5204 for (
unsigned int d = 0;
d < dim; ++
d)
5206 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5207 for (
unsigned int e = 1;
e < dim; ++
e)
5208 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5210 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5211 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5212 for (
unsigned int d = 0;
d < dim; ++
d)
5214 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5215 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5217 for (
unsigned int d = 0;
d < dim; ++
d)
5219 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5220 for (
unsigned int e = 1;
e < dim; ++
e)
5221 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5227 VectorizedArrayType tmp[dim][dim];
5230 for (
unsigned int d = 0;
d < dim; ++
d)
5231 for (
unsigned int e = 0;
e < dim; ++
e)
5233 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5234 for (
unsigned int f = 1; f < dim; ++f)
5235 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5238 for (
unsigned int d = 0;
d < dim; ++
d)
5239 for (
unsigned int e = 0;
e < dim; ++
e)
5241 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5242 for (
unsigned int f = 1; f < dim; ++f)
5243 res += t_jac[d][f] * tmp[e][f];
5250 VectorizedArrayType
value[dim];
5251 for (
unsigned int d = 0;
d < dim; ++
d)
5253 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5254 for (
unsigned int e = 1;
e < dim; ++
e)
5255 value[d] += tmp[d][e] * jac_grad[d][e];
5257 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5258 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5259 for (
unsigned int d = 0;
d < dim; ++
d)
5261 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5262 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5267 for (
unsigned int d = 0;
d < dim; ++
d)
5269 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5270 for (
unsigned int e = 1;
e < dim; ++
e)
5271 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5272 for (
unsigned int e = 0;
e < dim; ++
e)
5273 value[e] -= t_jac[e][d] * tmp2;
5276 for (
unsigned int d = 0;
d < dim; ++
d)
5277 values[d * nqp] =
value[d];
5283 const std::size_t nqp_d = this->n_quadrature_points * dim;
5284 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5288 const VectorizedArrayType JxW =
5289 this->J_value[0] * this->quadrature_weights[q_point];
5295 std::array<VectorizedArrayType, dim> jac;
5296 for (
unsigned int d = 0;
d < dim; ++
d)
5297 jac[d] = this->jacobian[0][d][d];
5299 for (
unsigned int d = 0;
d < dim; ++
d)
5301 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5302 if constexpr (n_components == 1)
5303 gradients[d] = grad_in[d] * factor;
5305 for (
unsigned int comp = 0; comp < n_components; ++comp)
5306 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5313 this->jacobian[q_point] :
5315 const VectorizedArrayType JxW =
5317 this->J_value[q_point] :
5318 this->J_value[0] * this->quadrature_weights[q_point];
5319 if constexpr (n_components == 1)
5320 for (
unsigned int d = 0;
d < dim; ++
d)
5322 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5323 for (
unsigned int e = 1;
e < dim; ++
e)
5324 new_val += (jac[e][d] * grad_in[e]);
5328 for (
unsigned int comp = 0; comp < n_components; ++comp)
5329 for (
unsigned int d = 0;
d < dim; ++
d)
5331 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5332 for (
unsigned int e = 1;
e < dim; ++
e)
5333 new_val += (jac[e][d] * grad_in[comp][e]);
5345 typename VectorizedArrayType>
5346template <
int,
typename>
5350 const unsigned int q_point)
5352 static_assert(n_components == 1 && dim == 1,
5353 "Do not try to modify the default template parameters used for"
5354 " selectively enabling this function via std::enable_if!");
5355 submit_gradient(grad_in[0], q_point);
5364 typename VectorizedArrayType>
5370 Assert(this->normal_x_jacobian !=
nullptr,
5372 "update_gradients"));
5375 this->gradients_quad_submitted =
true;
5378 const std::size_t nqp_d = this->n_quadrature_points * dim;
5379 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5383 const VectorizedArrayType JxW_jac = this->J_value[0] *
5384 this->quadrature_weights[q_point] *
5385 this->normal_x_jacobian[0][dim - 1];
5386 for (
unsigned int comp = 0; comp < n_components; ++comp)
5388 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5389 gradients[comp * nqp_d + d] = VectorizedArrayType();
5390 if constexpr (n_components == 1)
5391 gradients[dim - 1] = grad_in * JxW_jac;
5393 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5398 const unsigned int index =
5401 this->normal_x_jacobian[
index];
5402 const VectorizedArrayType JxW =
5404 this->J_value[
index] * this->quadrature_weights[q_point] :
5405 this->J_value[
index];
5406 for (
unsigned int comp = 0; comp < n_components; ++comp)
5407 for (
unsigned int d = 0;
d < dim; ++
d)
5408 if constexpr (n_components == 1)
5411 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5421 typename VectorizedArrayType>
5424 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5431 Assert(this->J_value !=
nullptr,
5433 "update_hessians"));
5434 Assert(this->jacobian !=
nullptr,
5436 "update_hessians"));
5439 this->hessians_quad_submitted =
true;
5443 const std::size_t nqp = this->n_quadrature_points;
5444 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5447 const VectorizedArrayType JxW =
5448 this->J_value[0] * this->quadrature_weights[q_point];
5451 for (
unsigned int d = 0;
d < dim; ++
d)
5453 const auto jac_d = this->jacobian[0][
d][
d];
5454 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5455 for (
unsigned int comp = 0; comp < n_components; ++comp)
5456 if constexpr (n_components == 1)
5457 this->hessians_quad[
d * nqp + q_point] =
5458 hessian_in[
d][
d] * factor;
5460 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5461 hessian_in[comp][d][d] * factor;
5465 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5466 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5468 const auto jac_d = this->jacobian[0][
d][
d];
5469 const auto jac_e = this->jacobian[0][
e][
e];
5470 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5471 for (
unsigned int comp = 0; comp < n_components; ++comp)
5472 if constexpr (n_components == 1)
5473 this->hessians_quad[off_dia * nqp + q_point] =
5474 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5476 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5477 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5484 const VectorizedArrayType JxW =
5485 this->J_value[0] * this->quadrature_weights[q_point];
5486 for (
unsigned int comp = 0; comp < n_components; ++comp)
5489 if constexpr (n_components == 1)
5490 hessian_c = hessian_in;
5492 hessian_c = hessian_in[comp];
5495 VectorizedArrayType tmp[dim][dim];
5496 for (
unsigned int i = 0; i < dim; ++i)
5497 for (
unsigned int j = 0; j < dim; ++j)
5499 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5500 for (
unsigned int k = 1; k < dim; ++k)
5501 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5505 VectorizedArrayType tmp2[dim][dim];
5506 for (
unsigned int i = 0; i < dim; ++i)
5507 for (
unsigned int j = 0; j < dim; ++j)
5509 tmp2[i][j] = jac[0][i] * tmp[0][j];
5510 for (
unsigned int k = 1; k < dim; ++k)
5511 tmp2[i][j] += jac[k][i] * tmp[k][j];
5515 for (
unsigned int d = 0;
d < dim; ++
d)
5516 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5520 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5521 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5522 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5523 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5529 const VectorizedArrayType JxW = this->J_value[q_point];
5530 const auto &jac_grad = this->jacobian_gradients[q_point];
5531 for (
unsigned int comp = 0; comp < n_components; ++comp)
5534 if constexpr (n_components == 1)
5535 hessian_c = hessian_in;
5537 hessian_c = hessian_in[comp];
5540 VectorizedArrayType tmp[dim][dim];
5541 for (
unsigned int i = 0; i < dim; ++i)
5542 for (
unsigned int j = 0; j < dim; ++j)
5544 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5545 for (
unsigned int k = 1; k < dim; ++k)
5546 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5550 VectorizedArrayType tmp2[dim][dim];
5551 for (
unsigned int i = 0; i < dim; ++i)
5552 for (
unsigned int j = 0; j < dim; ++j)
5554 tmp2[i][j] = jac[0][i] * tmp[0][j];
5555 for (
unsigned int k = 1; k < dim; ++k)
5556 tmp2[i][j] += jac[k][i] * tmp[k][j];
5560 for (
unsigned int d = 0;
d < dim; ++
d)
5561 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5565 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5566 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5567 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5568 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5571 for (
unsigned int d = 0;
d < dim; ++
d)
5573 VectorizedArrayType
sum = 0;
5574 for (
unsigned int e = 0;
e < dim; ++
e)
5575 sum += hessian_c[e][e] * jac_grad[e][d];
5576 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5577 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5579 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5580 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5593 typename VectorizedArrayType>
5597 const unsigned int q_point)
5604 Assert(this->J_value !=
nullptr,
5606 "update_hessians"));
5607 Assert(this->jacobian !=
nullptr,
5609 "update_hessians"));
5612 this->hessians_quad_submitted =
true;
5616 const std::size_t nqp = this->n_quadrature_points;
5617 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5620 const VectorizedArrayType JxW =
5621 this->J_value[0] * this->quadrature_weights[q_point];
5623 const auto nxj = this->normal_x_jacobian[0];
5626 for (
unsigned int d = 0;
d < dim; ++
d)
5628 const auto nxj_d = nxj[
d];
5629 const VectorizedArrayType factor = nxj_d * nxj_d * JxW;
5630 for (
unsigned int comp = 0; comp < n_components; ++comp)
5631 if constexpr (n_components == 1)
5632 this->hessians_quad[
d * nqp + q_point] =
5633 normal_hessian_in * factor;
5635 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5636 normal_hessian_in[comp] * factor;
5640 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5641 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5643 const auto jac_d = nxj[
d];
5644 const auto jac_e = nxj[
e];
5645 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5646 for (
unsigned int comp = 0; comp < n_components; ++comp)
5647 if constexpr (n_components == 1)
5648 this->hessians_quad[off_dia * nqp + q_point] =
5649 2. * normal_hessian_in * factor;
5651 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5652 2. * normal_hessian_in[comp] * factor;
5657 const auto normal = this->normal_vector(q_point);
5658 const auto normal_projector =
outer_product(normal, normal);
5659 if constexpr (n_components == 1)
5660 submit_hessian(normal_hessian_in * normal_projector, q_point);
5664 for (
unsigned int comp = 0; comp < n_components; ++comp)
5665 tmp[comp] = normal_hessian_in[comp] * normal_projector;
5666 submit_hessian(tmp, q_point);
5677 typename VectorizedArrayType>
5682 VectorizedArrayType>::value_type
5689 Assert(this->values_quad_submitted ==
true,
5694 const std::size_t nqp = this->n_quadrature_points;
5695 for (
unsigned int q = 0; q < nqp; ++q)
5696 for (
unsigned int comp = 0; comp < n_components; ++comp)
5697 return_value[comp] += this->values_quad[comp * nqp + q];
5698 if constexpr (n_components == 1)
5699 return return_value[0];
5701 return return_value;
5710 typename VectorizedArrayType>
5711template <
int,
typename>
5716 static_assert(n_components == dim,
5717 "Do not try to modify the default template parameters used for"
5718 " selectively enabling this function via std::enable_if!");
5722 Assert(this->gradients_quad_initialized ==
true,
5726 Assert(this->jacobian !=
nullptr,
5728 "update_gradients"));
5730 VectorizedArrayType divergence;
5731 const std::size_t nqp = this->n_quadrature_points;
5734 this->
data->element_type ==
5737 VectorizedArrayType inv_det =
5740 this->jacobian[0][0][0] *
5741 ((dim == 2) ? this->jacobian[0][1][1] :
5742 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5750 if (is_face && dim == 2 && this->get_face_no() < 2)
5754 divergence = this->gradients_quad[q_point * dim];
5755 for (
unsigned int d = 1;
d < dim; ++
d)
5756 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5757 divergence *= inv_det;
5766 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5767 for (
unsigned int d = 1;
d < dim; ++
d)
5768 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5769 this->jacobian[0][
d][
d];
5776 this->jacobian[q_point] :
5778 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5779 for (
unsigned int e = 1;
e < dim; ++
e)
5780 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5781 for (
unsigned int d = 1;
d < dim; ++
d)
5782 for (
unsigned int e = 0;
e < dim; ++
e)
5784 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5796 typename VectorizedArrayType>
5797template <
int,
typename>
5802 static_assert(n_components == dim,
5803 "Do not try to modify the default template parameters used for"
5804 " selectively enabling this function via std::enable_if!");
5807 const auto grad = get_gradient(q_point);
5808 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5809 VectorizedArrayType half = Number(0.5);
5810 for (
unsigned int d = 0;
d < dim; ++
d)
5811 symmetrized[d] = grad[d][d];
5817 symmetrized[2] = grad[0][1] + grad[1][0];
5818 symmetrized[2] *= half;
5821 symmetrized[3] = grad[0][1] + grad[1][0];
5822 symmetrized[3] *= half;
5823 symmetrized[4] = grad[0][2] + grad[2][0];
5824 symmetrized[4] *= half;
5825 symmetrized[5] = grad[1][2] + grad[2][1];
5826 symmetrized[5] *= half;
5840 typename VectorizedArrayType>
5841template <
int,
typename>
5843 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5845 get_curl(const unsigned
int q_point) const
5847 static_assert(dim > 1 && n_components == dim,
5848 "Do not try to modify the default template parameters used for"
5849 " selectively enabling this function via std::enable_if!");
5853 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5857 curl[0] = grad[1][0] - grad[0][1];
5860 curl[0] = grad[2][1] - grad[1][2];
5861 curl[1] = grad[0][2] - grad[2][0];
5862 curl[2] = grad[1][0] - grad[0][1];
5876 typename VectorizedArrayType>
5877template <
int,
typename>
5881 const unsigned int q_point)
5883 static_assert(n_components == dim,
5884 "Do not try to modify the default template parameters used for"
5885 " selectively enabling this function via std::enable_if!");
5892 Assert(this->J_value !=
nullptr,
5894 "update_gradients"));
5895 Assert(this->jacobian !=
nullptr,
5897 "update_gradients"));
5900 this->gradients_quad_submitted =
true;
5903 const std::size_t nqp_d = this->n_quadrature_points * dim;
5904 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5906 if (this->
data->element_type ==
5913 const VectorizedArrayType fac =
5915 this->quadrature_weights[q_point] * div_in :
5917 this->J_value[q_point] :
5918 this->J_value[0] * this->quadrature_weights[q_point]) *
5921 this->jacobian[this->cell_type >
5925 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5927 for (
unsigned int d = 0;
d < dim; ++
d)
5929 for (
unsigned int e = 0;
e < dim; ++
e)
5930 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5932 this->divergence_is_requested =
true;
5939 const VectorizedArrayType fac =
5940 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5941 for (
unsigned int d = 0;
d < dim; ++
d)
5943 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5944 for (
unsigned int e = 0;
e < dim; ++
e)
5945 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5952 this->jacobian[q_point] :
5954 const VectorizedArrayType fac =
5956 this->J_value[q_point] :
5957 this->J_value[0] * this->quadrature_weights[q_point]) *
5959 for (
unsigned int d = 0;
d < dim; ++
d)
5961 for (
unsigned int e = 0;
e < dim; ++
e)
5962 gradients[d * nqp_d + e] = jac[d][e] * fac;
5974 typename VectorizedArrayType>
5975template <
int,
typename>
5980 const unsigned int q_point)
5982 static_assert(n_components == dim,
5983 "Do not try to modify the default template parameters used for"
5984 " selectively enabling this function via std::enable_if!");
5987 this->
data->element_type !=
5999 Assert(this->J_value !=
nullptr,
6001 "update_gradients"));
6002 Assert(this->jacobian !=
nullptr,
6004 "update_gradients"));
6007 this->gradients_quad_submitted =
true;
6010 const std::size_t nqp_d = this->n_quadrature_points * dim;
6011 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
6014 const VectorizedArrayType JxW =
6015 this->J_value[0] * this->quadrature_weights[q_point];
6017 for (
unsigned int d = 0;
d < dim; ++
d)
6018 gradients[d * nqp_d + d] =
6020 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6021 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
6023 const VectorizedArrayType
value =
6032 const VectorizedArrayType JxW =
6034 this->J_value[q_point] :
6035 this->J_value[0] * this->quadrature_weights[q_point];
6038 this->jacobian[q_point] :
6040 VectorizedArrayType weighted[dim][dim];
6041 for (
unsigned int i = 0; i < dim; ++i)
6043 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6044 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6046 const VectorizedArrayType
value =
6048 weighted[i][j] =
value;
6049 weighted[j][i] =
value;
6051 for (
unsigned int comp = 0; comp < dim; ++comp)
6052 for (
unsigned int d = 0;
d < dim; ++
d)
6054 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6055 for (
unsigned int e = 1;
e < dim; ++
e)
6056 new_val += jac[e][d] * weighted[comp][e];
6068 typename VectorizedArrayType>
6069template <
int,
typename>
6073 const unsigned int q_point)
6075 static_assert(n_components == dim,
6076 "Do not try to modify the default template parameters used for"
6077 " selectively enabling this function via std::enable_if!");
6083 grad[1][0] = curl[0];
6084 grad[0][1] = -curl[0];
6087 grad[2][1] = curl[0];
6088 grad[1][2] = -curl[0];
6089 grad[0][2] = curl[1];
6090 grad[2][0] = -curl[1];
6091 grad[1][0] = curl[2];
6092 grad[0][1] = -curl[2];
6097 submit_gradient(grad, q_point);
6110 typename VectorizedArrayType>
6116 VectorizedArrayType>::
6118 const unsigned int fe_no,
6119 const unsigned int quad_no,
6120 const unsigned int first_selected_component,
6121 const unsigned int active_fe_index,
6122 const unsigned int active_quad_index)
6123 : BaseClass(matrix_free,
6125 first_selected_component,
6133 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6134 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6135 , n_q_points(this->
data->n_q_points)
6137 check_template_arguments(fe_no, 0);
6147 typename VectorizedArrayType>
6153 VectorizedArrayType>::
6155 const std::pair<unsigned int, unsigned int> &range,
6156 const unsigned int dof_no,
6157 const unsigned int quad_no,
6158 const unsigned int first_selected_component)
6162 first_selected_component,
6163 matrix_free.get_cell_active_fe_index(range))
6173 typename VectorizedArrayType>
6179 VectorizedArrayType>::
6184 const unsigned int first_selected_component)
6185 : BaseClass(mapping,
6189 first_selected_component,
6191 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6192 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6193 , n_q_points(this->
data->n_q_points)
6205 typename VectorizedArrayType>
6211 VectorizedArrayType>::
6215 const unsigned int first_selected_component)
6220 first_selected_component,
6222 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6223 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6224 , n_q_points(this->
data->n_q_points)
6236 typename VectorizedArrayType>
6242 VectorizedArrayType>::
6245 const unsigned int first_selected_component)
6246 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6248 other.mapped_geometry->get_quadrature(),
6249 other.mapped_geometry->get_fe_values().get_update_flags(),
6250 first_selected_component,
6252 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6253 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6254 , n_q_points(this->
data->n_q_points)
6266 typename VectorizedArrayType>
6275 , dofs_per_component(this->
data->dofs_per_component_on_cell)
6276 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
6277 , n_q_points(this->
data->n_q_points)
6289 typename VectorizedArrayType>
6295 VectorizedArrayType> &
6301 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6303 BaseClass::operator=(other);
6315 typename VectorizedArrayType>
6322 VectorizedArrayType>::
6323 check_template_arguments(
const unsigned int dof_no,
6324 const unsigned int first_selected_component)
6327 (void)first_selected_component;
6330 this->
data->dofs_per_component_on_cell > 0,
6332 "There is nothing useful you can do with an FEEvaluation object with "
6333 "FE_Nothing, i.e., without DoFs! If you have passed to "
6334 "MatrixFree::reinit() a collection of finite elements also containing "
6335 "FE_Nothing, please check - before creating FEEvaluation - the category "
6336 "of the current range by calling either "
6337 "MatrixFree::get_cell_range_category(range) or "
6338 "MatrixFree::get_face_range_category(range). The returned category "
6339 "is the index of the active FE, which you can use to exclude "
6346 if ((
static_cast<unsigned int>(fe_degree) !=
6348 static_cast<unsigned int>(fe_degree) !=
6349 this->
data->data.front().fe_degree) ||
6350 n_q_points != this->n_quadrature_points)
6352 std::string message =
6353 "-------------------------------------------------------\n";
6355 "Illegal arguments in constructor/wrong template arguments!\n";
6356 message +=
" Called --> FEEvaluation<dim,";
6360 message +=
",Number>(data";
6376 if (
static_cast<unsigned int>(fe_degree) ==
6377 this->
data->data.front().fe_degree)
6379 proposed_dof_comp = dof_no;
6380 proposed_fe_comp = first_selected_component;
6383 for (
unsigned int no = 0;
6386 for (
unsigned int nf = 0;
6389 if (this->matrix_free
6390 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6392 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6394 proposed_dof_comp = no;
6395 proposed_fe_comp = nf;
6399 this->mapping_data->descriptor[this->active_quad_index]
6401 proposed_quad_comp = this->quad_no;
6403 for (
unsigned int no = 0;
6409 .descriptor[this->active_quad_index]
6410 .n_q_points == n_q_points)
6412 proposed_quad_comp = no;
6419 if (proposed_dof_comp != first_selected_component)
6420 message +=
"Wrong vector component selection:\n";
6422 message +=
"Wrong quadrature formula selection:\n";
6423 message +=
" Did you mean FEEvaluation<dim,";
6427 message +=
",Number>(data";
6437 std::string correct_pos;
6438 if (proposed_dof_comp != dof_no)
6439 correct_pos =
" ^ ";
6442 if (proposed_quad_comp != this->quad_no)
6443 correct_pos +=
" ^ ";
6446 if (proposed_fe_comp != first_selected_component)
6447 correct_pos +=
" ^\n";
6449 correct_pos +=
" \n";
6456 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6457 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6458 message +=
"Wrong template arguments:\n";
6459 message +=
" Did you mean FEEvaluation<dim,";
6464 message +=
",Number>(data";
6473 std::string correct_pos;
6474 if (this->
data->data.front().fe_degree !=
6475 static_cast<unsigned int>(fe_degree))
6479 if (proposed_n_q_points_1d != n_q_points_1d)
6480 correct_pos +=
" ^\n";
6482 correct_pos +=
" \n";
6483 message +=
" " + correct_pos;
6485 Assert(
static_cast<unsigned int>(fe_degree) ==
6486 this->
data->data.front().fe_degree &&
6487 n_q_points == this->n_quadrature_points,
6493 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6504 typename VectorizedArrayType>
6513 Assert(this->matrix_free !=
nullptr,
6514 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6515 " Integer indexing is not possible."));
6523 const unsigned int offsets =
6524 this->mapping_data->data_index_offsets[
cell_index];
6525 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6526 this->J_value = &this->mapping_data->JxW_values[offsets];
6527 if (!this->mapping_data->jacobian_gradients[0].empty())
6529 this->jacobian_gradients =
6530 this->mapping_data->jacobian_gradients[0].data() + offsets;
6531 this->jacobian_gradients_non_inverse =
6532 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6538 for (
unsigned int i = 0; i < n_lanes; ++i)
6539 this->cell_ids[i] =
cell_index * n_lanes + i;
6546 this->cell_ids[i] =
cell_index * n_lanes + i;
6547 for (; i < n_lanes; ++i)
6551 if (this->mapping_data->quadrature_points.empty() ==
false)
6553 &this->mapping_data->quadrature_points
6554 [this->mapping_data->quadrature_point_offsets[this->cell]];
6558 this->is_reinitialized =
true;
6559 this->dof_values_initialized =
false;
6560 this->values_quad_initialized =
false;
6561 this->gradients_quad_initialized =
false;
6562 this->hessians_quad_initialized =
false;
6573 typename VectorizedArrayType>
6580 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6587 this->cell_ids = cell_ids;
6592 for (
unsigned int v = 0; v < n_lanes; ++v)
6606 if (this->mapped_geometry ==
nullptr)
6607 this->mapped_geometry =
6608 std::make_shared<internal::MatrixFreeFunctions::
6609 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6611 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6613 auto &this_jacobian_data = mapping_storage.jacobians[0];
6614 auto &this_J_value_data = mapping_storage.JxW_values;
6615 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6616 auto &this_jacobian_gradients_non_inverse_data =
6617 mapping_storage.jacobian_gradients_non_inverse[0];
6618 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6622 if (this_jacobian_data.size() != 2)
6623 this_jacobian_data.resize_fast(2);
6625 if (this_J_value_data.size() != 1)
6626 this_J_value_data.resize_fast(1);
6628 const auto &update_flags_cells =
6632 this_jacobian_gradients_data.size() != 1)
6634 this_jacobian_gradients_data.resize_fast(1);
6635 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6639 this_quadrature_points_data.size() != 1)
6640 this_quadrature_points_data.resize_fast(1);
6644 if (this_jacobian_data.size() != this->n_quadrature_points)
6645 this_jacobian_data.resize_fast(this->n_quadrature_points);
6647 if (this_J_value_data.size() != this->n_quadrature_points)
6648 this_J_value_data.resize_fast(this->n_quadrature_points);
6650 const auto &update_flags_cells =
6654 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6656 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6657 this_jacobian_gradients_non_inverse_data.resize_fast(
6658 this->n_quadrature_points);
6662 this_quadrature_points_data.size() != this->n_quadrature_points)
6663 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6667 this->jacobian = this_jacobian_data.data();
6668 this->J_value = this_J_value_data.data();
6669 this->jacobian_gradients = this_jacobian_gradients_data.data();
6670 this->jacobian_gradients_non_inverse =
6671 this_jacobian_gradients_non_inverse_data.data();
6675 for (
unsigned int v = 0; v < n_lanes; ++v)
6682 const unsigned int cell_batch_index =
cell_index / n_lanes;
6683 const unsigned int offsets =
6684 this->mapping_data->data_index_offsets[cell_batch_index];
6685 const unsigned int lane =
cell_index % n_lanes;
6687 if (this->cell_type <=
6691 for (
unsigned int q = 0; q < 2; ++q)
6692 for (
unsigned int i = 0; i < dim; ++i)
6693 for (
unsigned int j = 0; j < dim; ++j)
6694 this_jacobian_data[q][i][j][v] =
6695 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6697 const unsigned int q = 0;
6699 this_J_value_data[q][v] =
6700 this->mapping_data->JxW_values[offsets + q][lane];
6702 const auto &update_flags_cells =
6707 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6708 for (
unsigned int j = 0; j < dim; ++j)
6709 this_jacobian_gradients_data[q][i][j][v] =
6711 ->jacobian_gradients[0][offsets + q][i][j][lane];
6713 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6714 for (
unsigned int j = 0; j < dim; ++j)
6715 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6717 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6722 for (
unsigned int i = 0; i < dim; ++i)
6723 this_quadrature_points_data[q][i][v] =
6724 this->mapping_data->quadrature_points
6726 ->quadrature_point_offsets[cell_batch_index] +
6732 const auto cell_type =
6736 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6738 const unsigned int q_src =
6744 this_J_value_data[q][v] =
6745 this->mapping_data->JxW_values[offsets + q_src][lane];
6747 for (
unsigned int i = 0; i < dim; ++i)
6748 for (
unsigned int j = 0; j < dim; ++j)
6749 this_jacobian_data[q][i][j][v] =
6751 ->jacobians[0][offsets + q_src][i][j][lane];
6753 const auto &update_flags_cells =
6758 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6759 for (
unsigned int j = 0; j < dim; ++j)
6760 this_jacobian_gradients_data[q][i][j][v] =
6762 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6764 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6765 for (
unsigned int j = 0; j < dim; ++j)
6766 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6768 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6781 this->mapping_data->quadrature_points
6783 ->quadrature_point_offsets[cell_batch_index] +
6787 this->mapping_data->jacobians[0][offsets + 1];
6789 for (
unsigned int d = 0;
d < dim; ++
d)
6792 static_cast<Number
>(
6793 this->descriptor->quadrature.point(q)[
d]);
6795 for (
unsigned int d = 0;
d < dim; ++
d)
6796 for (
unsigned int e = 0;
e < dim; ++
e)
6799 static_cast<Number
>(
6800 this->descriptor->quadrature.point(q)[
e]);
6802 for (
unsigned int i = 0; i < dim; ++i)
6803 this_quadrature_points_data[q][i][v] = point[i][lane];
6808 for (
unsigned int i = 0; i < dim; ++i)
6809 this_quadrature_points_data[q][i][v] =
6810 this->mapping_data->quadrature_points
6812 ->quadrature_point_offsets[cell_batch_index] +
6822 this->is_reinitialized =
true;
6823 this->dof_values_initialized =
false;
6824 this->values_quad_initialized =
false;
6825 this->gradients_quad_initialized =
false;
6826 this->hessians_quad_initialized =
false;
6837 typename VectorizedArrayType>
6838template <
bool level_dof_access>
6845 VectorizedArrayType>
::
6848 Assert(this->matrix_free ==
nullptr,
6849 ExcMessage(
"Cannot use initialization from cell iterator if "
6850 "initialized from MatrixFree object. Use variant for "
6851 "on the fly computation with arguments as for FEValues "
6854 this->mapped_geometry->reinit(
6856 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6857 if (level_dof_access)
6858 cell->get_mg_dof_indices(this->local_dof_indices);
6860 cell->get_dof_indices(this->local_dof_indices);
6864 this->is_reinitialized =
true;
6875 typename VectorizedArrayType>
6882 VectorizedArrayType>
::
6885 Assert(this->matrix_free == 0,
6886 ExcMessage(
"Cannot use initialization from cell iterator if "
6887 "initialized from MatrixFree object. Use variant for "
6888 "on the fly computation with arguments as for FEValues "
6891 this->mapped_geometry->reinit(cell);
6895 this->is_reinitialized =
true;
6906 typename VectorizedArrayType>
6913 VectorizedArrayType>::
6918 Assert(this->dof_values_initialized ==
true,
6921 evaluate(this->values_dofs, evaluation_flags);
6931 typename VectorizedArrayType>
6938 VectorizedArrayType>::
6939 evaluate(
const VectorizedArrayType *values_array,
6942 const bool hessians_on_general_cells =
6946 if (hessians_on_general_cells)
6949 if (this->
data->element_type ==
6955 if constexpr (fe_degree > -1)
6958 template run<fe_degree, n_q_points_1d>(n_components,
6959 evaluation_flag_actual,
6967 evaluation_flag_actual,
6968 const_cast<VectorizedArrayType *
>(values_array),
6974 this->values_quad_initialized =
6976 this->gradients_quad_initialized =
6978 this->hessians_quad_initialized =
6989 template <
typename Number,
6990 typename VectorizedArrayType,
6992 typename EvaluatorType,
6993 std::enable_if_t<internal::has_begin<VectorType> &&
6996 VectorizedArrayType *
6997 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
7003 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7004 const auto &dof_info = fe_eval.get_dof_info();
7011 if (std::is_same_v<typename VectorType::value_type, Number> &&
7015 interleaved_contiguous &&
7016 reinterpret_cast<
std::size_t>(
7018 dof_info.dof_indices_contiguous
7019 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
7020 [cell * VectorizedArrayType::
size()]) %
7021 sizeof(VectorizedArrayType) ==
7024 return reinterpret_cast<VectorizedArrayType *
>(
7028 [cell * VectorizedArrayType::size()] +
7030 [fe_eval.get_active_fe_index()]
7031 [fe_eval.get_first_selected_component()] *
7032 VectorizedArrayType::size());
7041 template <
typename Number,
7042 typename VectorizedArrayType,
7044 typename EvaluatorType,
7045 std::enable_if_t<!internal::has_begin<VectorType> ||
7048 VectorizedArrayType *
7049 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7062 typename VectorizedArrayType>
7063template <
typename VectorType>
7070 VectorizedArrayType>::
7071 gather_evaluate(
const VectorType &input_vector,
7074 const VectorizedArrayType *src_ptr =
7075 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7076 *
this, input_vector);
7077 if (src_ptr !=
nullptr)
7078 evaluate(src_ptr, evaluation_flag);
7081 this->read_dof_values(input_vector);
7082 evaluate(this->begin_dof_values(), evaluation_flag);
7093 typename VectorizedArrayType>
7100 VectorizedArrayType>::
7103 integrate(integration_flag, this->values_dofs);
7107 this->dof_values_initialized =
true;
7118 typename VectorizedArrayType>
7125 VectorizedArrayType>::
7127 VectorizedArrayType *values_array,
7128 const bool sum_into_values_array)
7133 Assert(this->values_quad_submitted ==
true,
7136 Assert(this->gradients_quad_submitted ==
true,
7139 Assert(this->hessians_quad_submitted ==
true,
7142 Assert(this->matrix_free !=
nullptr ||
7143 this->mapped_geometry->is_initialized(),
7149 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7150 "EvaluationFlags::hessians are supported."));
7156 unsigned int size = n_components * dim * n_q_points;
7159 for (
unsigned int i = 0; i <
size; ++i)
7160 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7164 for (
unsigned int i = 0; i <
size; ++i)
7165 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7170 if (n_components == dim &&
7171 this->
data->element_type ==
7175 this->divergence_is_requested ==
false)
7177 unsigned int size = n_components * n_q_points;
7180 for (
unsigned int i = 0; i <
size; ++i)
7181 this->values_quad[i] += this->values_from_gradients_quad[i];
7185 for (
unsigned int i = 0; i <
size; ++i)
7186 this->values_quad[i] = this->values_from_gradients_quad[i];
7191 if constexpr (fe_degree > -1)
7194 template run<fe_degree, n_q_points_1d>(n_components,
7195 integration_flag_actual,
7198 sum_into_values_array);
7204 integration_flag_actual,
7207 sum_into_values_array);
7212 this->dof_values_initialized =
true;
7223 typename VectorizedArrayType>
7224template <
typename VectorType>
7231 VectorizedArrayType>::
7233 VectorType &destination)
7235 VectorizedArrayType *dst_ptr =
7236 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7237 *
this, destination);
7238 if (dst_ptr !=
nullptr)
7239 integrate(integration_flag, dst_ptr,
true);
7242 integrate(integration_flag, this->begin_dof_values());
7243 this->distribute_local_to_global(destination);
7254 typename VectorizedArrayType>
7261 VectorizedArrayType>::dof_indices()
const
7278 typename VectorizedArrayType>
7284 VectorizedArrayType>::
7287 const bool is_interior_face,
7288 const unsigned int dof_no,
7289 const unsigned int quad_no,
7290 const unsigned int first_selected_component,
7291 const unsigned int active_fe_index,
7292 const unsigned int active_quad_index,
7293 const unsigned int face_type)
7294 : BaseClass(matrix_free,
7296 first_selected_component,
7304 , dofs_per_component(this->
data->dofs_per_component_on_cell)
7305 , dofs_per_cell(this->
data->dofs_per_component_on_cell * n_components_)
7306 , n_q_points(this->n_quadrature_points)
7316 typename VectorizedArrayType>
7322 VectorizedArrayType>::
7325 const std::pair<unsigned int, unsigned int> &range,
7326 const bool is_interior_face,
7327 const unsigned int dof_no,
7328 const unsigned int quad_no,
7329 const unsigned int first_selected_component)
7334 first_selected_component,
7335 matrix_free.get_face_active_fe_index(range,
7338 matrix_free.get_face_info(range.
first).face_type)
7348 typename VectorizedArrayType>
7355 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7357 Assert(this->mapped_geometry ==
nullptr,
7358 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7359 " Integer indexing is not possible"));
7360 if (this->mapped_geometry !=
nullptr)
7363 this->cell = face_index;
7364 this->dof_access_index =
7365 this->is_interior_face() ?
7373 this->matrix_free->get_task_info().boundary_partition_data.back())
7374 Assert(this->is_interior_face(),
7376 "Boundary faces do not have a neighbor. When looping over "
7377 "boundary faces use FEFaceEvaluation with the parameter "
7378 "is_interior_face set to true. "));
7380 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7385 this->face_ids[i] = face_index * n_lanes + i;
7386 for (; i < n_lanes; ++i)
7389 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7390 const unsigned int offsets =
7391 this->mapping_data->data_index_offsets[face_index];
7392 this->J_value = &this->mapping_data->JxW_values[offsets];
7393 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7395 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7396 this->normal_x_jacobian =
7398 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7399 this->jacobian_gradients =
7400 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7402 this->jacobian_gradients_non_inverse =
7404 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7408 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7411 this->mapping_data->quadrature_point_offsets.size());
7413 this->mapping_data->quadrature_points.data() +
7414 this->mapping_data->quadrature_point_offsets[this->cell];
7419 this->is_reinitialized =
true;
7420 this->dof_values_initialized =
false;
7421 this->values_quad_initialized =
false;
7422 this->gradients_quad_initialized =
false;
7423 this->hessians_quad_initialized =
false;
7434 typename VectorizedArrayType>
7442 const unsigned int face_number)
7448 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7452 Assert(this->mapped_geometry ==
nullptr,
7453 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7454 " Integer indexing is not possible"));
7455 if (this->mapped_geometry !=
nullptr)
7460 .faces_by_cells_type[
cell_index][face_number];
7463 this->dof_access_index =
7466 if (this->is_interior_face() ==
false)
7471 for (
unsigned int i = 0; i < n_lanes; ++i)
7474 const unsigned int cell_this =
cell_index * n_lanes + i;
7476 unsigned int face_index =
7481 this->face_ids[i] = face_index;
7486 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7487 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7494 auto cell_m = faces.cells_interior[face_index % n_lanes];
7495 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7497 const bool face_identifies_as_interior = cell_m != cell_this;
7499 Assert(cell_m == cell_this || cell_p == cell_this,
7503 if (face_identifies_as_interior)
7505 this->cell_ids[i] = cell_m;
7506 this->face_numbers[i] = faces.interior_face_no;
7510 this->cell_ids[i] = cell_p;
7511 this->face_numbers[i] = faces.exterior_face_no;
7514 const bool orientation_interior_face = faces.face_orientation >= 8;
7515 unsigned int face_orientation = faces.face_orientation % 8;
7516 if (face_identifies_as_interior != orientation_interior_face)
7518 constexpr std::array<std::uint8_t, 8> table{
7519 {0, 1, 6, 3, 4, 5, 2, 7}};
7520 face_orientation = table[face_orientation];
7522 this->face_orientations[i] = face_orientation;
7527 this->face_orientations[0] = 0;
7528 this->face_numbers[0] = face_number;
7533 for (
unsigned int i = 0; i < n_lanes; ++i)
7534 this->cell_ids[i] =
cell_index * n_lanes + i;
7542 this->cell_ids[i] =
cell_index * n_lanes + i;
7543 for (; i < n_lanes; ++i)
7546 for (
unsigned int i = 0; i < n_lanes; ++i)
7553 const unsigned int offsets =
7555 .face_data_by_cells[this->quad_no]
7560 .face_data_by_cells[this->quad_no]
7561 .JxW_values.size());
7563 .face_data_by_cells[this->quad_no]
7564 .JxW_values[offsets];
7566 .face_data_by_cells[this->quad_no]
7567 .normal_vectors[offsets];
7569 .face_data_by_cells[this->quad_no]
7570 .jacobians[!this->is_interior_face()][offsets];
7571 this->normal_x_jacobian =
7573 .face_data_by_cells[this->quad_no]
7574 .normals_times_jacobians[!this->is_interior_face()][offsets];
7575 this->jacobian_gradients =
7576 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7578 this->jacobian_gradients_non_inverse =
7580 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7585 .face_data_by_cells[this->quad_no]
7586 .quadrature_point_offsets.empty() ==
false)
7588 const unsigned int index =
7592 .face_data_by_cells[this->quad_no]
7593 .quadrature_point_offsets.size());
7595 .face_data_by_cells[this->quad_no]
7596 .quadrature_points.data() +
7598 .face_data_by_cells[this->quad_no]
7599 .quadrature_point_offsets[
index];
7604 this->is_reinitialized =
true;
7605 this->dof_values_initialized =
false;
7606 this->values_quad_initialized =
false;
7607 this->gradients_quad_initialized =
false;
7608 this->hessians_quad_initialized =
false;
7619 typename VectorizedArrayType>
7626 VectorizedArrayType>::
7634 evaluate(this->values_dofs, evaluation_flag);
7644 typename VectorizedArrayType>
7651 VectorizedArrayType>::
7652 evaluate(
const VectorizedArrayType *values_array,
7655 Assert((evaluation_flag &
7658 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7659 "and EvaluationFlags::hessians are supported."));
7661 const bool hessians_on_general_cells =
7665 if (hessians_on_general_cells)
7668 if (this->
data->element_type ==
7674 if constexpr (fe_degree > -1)
7676 template
run<fe_degree, n_q_points_1d>(n_components,
7677 evaluation_flag_actual,
7682 n_components, evaluation_flag_actual, values_array, *
this);
7686 this->values_quad_initialized =
7688 this->gradients_quad_initialized =
7690 this->hessians_quad_initialized =
7702 typename VectorizedArrayType>
7709 VectorizedArrayType>::
7717 project_to_face(this->values_dofs, evaluation_flag);
7727 typename VectorizedArrayType>
7734 VectorizedArrayType>::
7735 project_to_face(
const VectorizedArrayType *values_array,
7738 Assert((evaluation_flag &
7741 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7742 "and EvaluationFlags::hessians are supported."));
7744 const bool hessians_on_general_cells =
7748 if (hessians_on_general_cells)
7751 if (this->
data->element_type ==
7757 if constexpr (fe_degree > -1)
7760 VectorizedArrayType>::template run<fe_degree>(n_components,
7761 evaluation_flag_actual,
7767 evaluation_flag_actual,
7781 typename VectorizedArrayType>
7788 VectorizedArrayType>::
7791 Assert((evaluation_flag &
7794 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7795 "and EvaluationFlags::hessians are supported."));
7797 const bool hessians_on_general_cells =
7801 if (hessians_on_general_cells)
7804 if (this->
data->element_type ==
7810 if constexpr (fe_degree > -1)
7813 VectorizedArrayType>::template run<fe_degree>(n_components,
7814 evaluation_flag_actual,
7822 this->values_quad_initialized =
7824 this->gradients_quad_initialized =
7826 this->hessians_quad_initialized =
7838 typename VectorizedArrayType>
7845 VectorizedArrayType>::
7847 const bool sum_into_values)
7849 integrate(integration_flag, this->values_dofs, sum_into_values);
7853 this->dof_values_initialized =
true;
7864 typename VectorizedArrayType>
7871 VectorizedArrayType>::
7873 VectorizedArrayType *values_array,
7874 const bool sum_into_values)
7876 Assert((integration_flag &
7879 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7880 "and EvaluationFlags::hessians are supported."));
7886 unsigned int size = n_components * dim * n_q_points;
7889 for (
unsigned int i = 0; i <
size; ++i)
7890 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7894 for (
unsigned int i = 0; i <
size; ++i)
7895 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7900 if (this->
data->element_type ==
7904 this->divergence_is_requested ==
false)
7906 unsigned int size = n_components * n_q_points;
7909 for (
unsigned int i = 0; i <
size; ++i)
7910 this->values_quad[i] += this->values_from_gradients_quad[i];
7914 for (
unsigned int i = 0; i <
size; ++i)
7915 this->values_quad[i] = this->values_from_gradients_quad[i];
7920 if constexpr (fe_degree > -1)
7922 template
run<fe_degree, n_q_points_1d>(n_components,
7923 integration_flag_actual,
7930 integration_flag_actual,
7943 typename VectorizedArrayType>
7950 VectorizedArrayType>::
7953 Assert((integration_flag &
7956 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7957 "and EvaluationFlags::hessians are supported."));
7963 unsigned int size = n_components * dim * n_q_points;
7966 for (
unsigned int i = 0; i <
size; ++i)
7967 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7971 for (
unsigned int i = 0; i <
size; ++i)
7972 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7977 if (this->
data->element_type ==
7981 this->divergence_is_requested ==
false)
7983 unsigned int size = n_components * n_q_points;
7986 for (
unsigned int i = 0; i <
size; ++i)
7987 this->values_quad[i] += this->values_from_gradients_quad[i];
7991 for (
unsigned int i = 0; i <
size; ++i)
7992 this->values_quad[i] = this->values_from_gradients_quad[i];
7997 if constexpr (fe_degree > -1)
8000 VectorizedArrayType>::template run<fe_degree>(n_components,
8001 integration_flag_actual,
8017 typename VectorizedArrayType>
8024 VectorizedArrayType>::
8026 const bool sum_into_values)
8028 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
8032 this->dof_values_initialized =
true;
8043 typename VectorizedArrayType>
8050 VectorizedArrayType>::
8052 VectorizedArrayType *values_array,
8053 const bool sum_into_values)
8055 Assert((integration_flag &
8058 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8059 "and EvaluationFlags::hessians are supported."));
8066 if (this->
data->element_type ==
8070 this->divergence_is_requested ==
false)
8073 if constexpr (fe_degree > -1)
8076 VectorizedArrayType>::template run<fe_degree>(n_components,
8077 integration_flag_actual,
8084 integration_flag_actual,
8097 typename VectorizedArrayType>
8098template <
typename VectorType>
8105 VectorizedArrayType>::
8106 gather_evaluate(
const VectorType &input_vector,
8109 Assert((evaluation_flag &
8112 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8113 "and EvaluationFlags::hessians are supported."));
8115 const auto shared_vector_data = internal::get_shared_vector_data(
8117 this->dof_access_index ==
8119 this->active_fe_index,
8122 if (this->
data->data.front().fe_degree > 0 &&
8123 fast_evaluation_supported(this->
data->data.front().fe_degree,
8124 this->data->data.front().n_q_points_1d) &&
8127 typename VectorType::value_type,
8128 VectorizedArrayType>::
8129 supports(evaluation_flag,
8133 this->dof_info->index_storage_variants[this->dof_access_index]
8136 if constexpr (fe_degree > -1)
8140 typename VectorType::value_type,
8141 VectorizedArrayType>::template
run<fe_degree,
8145 internal::get_beginning<typename VectorType::value_type>(
8154 typename VectorType::value_type,
8155 VectorizedArrayType>::evaluate(n_components,
8157 internal::get_beginning<
8158 typename VectorType::value_type>(
8166 this->read_dof_values(input_vector);
8167 this->evaluate(evaluation_flag);
8173 this->gradients_quad_initialized =
8175 this->hessians_quad_initialized =
8187 typename VectorizedArrayType>
8188template <
typename VectorType>
8196 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8197 const bool integrate_gradients,
8198 VectorType &destination)
8205 integrate_scatter(flag, destination);
8215 typename VectorizedArrayType>
8216template <
typename VectorType>
8223 VectorizedArrayType>::
8225 VectorType &destination)
8227 Assert((this->dof_access_index ==
8229 this->is_interior_face() ==
false) ==
false,
8232 const auto shared_vector_data = internal::get_shared_vector_data(
8234 this->dof_access_index ==
8236 this->active_fe_index,
8239 if (this->
data->data.front().fe_degree > 0 &&
8240 fast_evaluation_supported(this->
data->data.front().fe_degree,
8241 this->data->data.front().n_q_points_1d) &&
8244 typename VectorType::value_type,
8245 VectorizedArrayType>::
8246 supports(integration_flag,
8250 this->dof_info->index_storage_variants[this->dof_access_index]
8253 if constexpr (fe_degree > -1)
8257 typename VectorType::value_type,
8258 VectorizedArrayType>::template
run<fe_degree,
8262 internal::get_beginning<typename VectorType::value_type>(
8271 typename VectorType::value_type,
8272 VectorizedArrayType>::integrate(n_components,
8274 internal::get_beginning<
8275 typename VectorType::value_type>(
8283 integrate(integration_flag);
8284 this->distribute_local_to_global(destination);
8295 typename VectorizedArrayType>
8302 VectorizedArrayType>::dof_indices()
const
8315 typename VectorizedArrayType>
8322 VectorizedArrayType>::
8323 fast_evaluation_supported(
const unsigned int given_degree,
8324 const unsigned int given_n_q_points_1d)
8326 return fe_degree == -1 ?
8339 typename VectorizedArrayType>
8346 VectorizedArrayType>::
8347 fast_evaluation_supported(
const unsigned int given_degree,
8348 const unsigned int given_n_q_points_1d)
8350 return fe_degree == -1 ?
8363 typename VectorizedArrayType>
8370 VectorizedArrayType>::at_boundary()
const
8372 Assert(this->dof_access_index !=
8376 if (this->is_interior_face() ==
false)
8381 this->matrix_free->n_boundary_face_batches()))
8394 typename VectorizedArrayType>
8401 VectorizedArrayType>::boundary_id()
const
8403 Assert(this->dof_access_index !=
8420 typename VectorizedArrayType>
8428 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8430 return this->
data->dofs_per_component_on_face;
8440 typename VectorizedArrayType>
8447 VectorizedArrayType>::get_dofs_projected_to_face()
8449 return this->
data->dofs_per_component_on_face * n_components_;
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > val_in, const unsigned int q_point)
void submit_normal_hessian(const value_type normal_hessian_in, const unsigned int q_point)
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
std::conditional_t< n_components_==1, Tensor< 2, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 3, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > > > hessian_type
std::conditional_t< n_components_==1, Tensor< 1, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 2, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > > > gradient_type
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
VectorizedArrayType get_divergence(const unsigned int q_point) const
hessian_type get_hessian(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
static constexpr unsigned int n_components
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
std::conditional_t< n_components_==1, VectorizedArrayType, Tensor< 1, n_components_, VectorizedArrayType > > value_type
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
static constexpr unsigned int n_lanes
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
value_type get_normal_hessian(const unsigned int q_point) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
const unsigned int quad_no
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const unsigned int active_fe_index
std::array< unsigned int, n_lanes > cell_ids
const unsigned int active_quad_index
bool is_interior_face() const
FEEvaluationData & operator=(const FEEvaluationData &other)
const unsigned int first_selected_component
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
void reinit(const std::array< unsigned int, n_lanes > &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
static constexpr unsigned int n_lanes
FEEvaluation(const FEEvaluation &other)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
void project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
unsigned int get_dofs_per_component_projected_to_face()
void project_to_face(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
static constexpr unsigned int static_dofs_per_component
static constexpr unsigned int n_lanes
static constexpr unsigned int static_n_q_points
unsigned int get_dofs_projected_to_face()
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
types::boundary_id boundary_id() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
void evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::ShapeInfo< Number > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_jacobian_grads
Gradient of volume element.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static constexpr unsigned int value
@ dof_access_face_exterior
@ dof_access_face_interior
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
std::vector< std::pair< unsigned int, unsigned int > > row_starts
std::vector< std::vector< unsigned int > > component_dof_indices_offset
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< unsigned int > dof_indices
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
std::vector< unsigned int > row_starts_plain_indices
std::vector< unsigned int > component_to_base_index
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
std::vector< unsigned int > plain_dof_indices
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
std::vector< unsigned int > dof_indices_interleaved
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
std::vector< QuadratureDescriptor > descriptor
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)