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;
350 template <
int dim_ = dim,
351 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
354 const unsigned int q_point);
396 template <
int dim_ = dim,
397 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
400 const unsigned int q_point);
418 const unsigned int q_point);
466 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
478 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
488 template <
int dim_ = dim,
489 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
490 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
509 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
512 const unsigned int q_point);
530 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
534 const unsigned int q_point);
548 template <
int dim_ = dim,
549 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
552 const unsigned int q_point);
593 const unsigned int dof_no,
596 const unsigned int fe_degree,
597 const unsigned int n_q_points,
601 const unsigned int face_type);
675 template <
typename VectorType,
typename VectorOperation>
679 const std::array<VectorType *, n_components_> &vectors,
682 n_components_> &vectors_sm,
683 const std::bitset<n_lanes> &mask,
684 const bool apply_constraints =
true)
const;
693 template <
typename VectorType,
typename VectorOperation>
697 const std::array<VectorType *, n_components_> &vectors,
700 n_components_> &vectors_sm,
701 const std::bitset<n_lanes> &mask)
const;
710 template <
typename VectorType,
typename VectorOperation>
714 const std::array<VectorType *, n_components_> &vectors)
const;
1318 typename VectorizedArrayType>
1323 VectorizedArrayType>
1326 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1327 "Type of Number and of VectorizedArrayType do not match.");
1369 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1454 const unsigned int dof_no = 0,
1455 const unsigned int quad_no = 0,
1468 const std::pair<unsigned int, unsigned int> &range,
1469 const unsigned int dof_no = 0,
1470 const unsigned int quad_no = 0,
1580 template <
bool level_dof_access>
1602 const unsigned int given_n_q_points_1d);
1645 template <
typename VectorType>
1676 VectorizedArrayType *values_array,
1677 const bool sum_into_values =
false);
1692 template <
typename VectorType>
1695 VectorType &output_vector);
1779 int n_q_points_1d = fe_degree + 1,
1780 int n_components_ = 1,
1781 typename Number = double,
1787 VectorizedArrayType>
1790 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1791 "Type of Number and of VectorizedArrayType do not match.");
1833 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1935 const unsigned int dof_no = 0,
1936 const unsigned int quad_no = 0,
1951 const std::pair<unsigned int, unsigned int> &range,
1953 const unsigned int dof_no = 0,
1954 const unsigned int quad_no = 0,
1968 reinit(
const unsigned int face_batch_number);
1978 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
1985 const unsigned int given_n_q_points_1d);
2049 template <
typename VectorType>
2065 const bool sum_into_values =
false);
2077 VectorizedArrayType *values_array,
2078 const bool sum_into_values =
false);
2095 const bool sum_into_values =
false);
2103 VectorizedArrayType *values_array,
2104 const bool sum_into_values =
false);
2117 template <
typename VectorType>
2120 VectorType &output_vector);
2125 template <
typename VectorType>
2128 const bool integrate_gradients,
2129 VectorType &output_vector);
2205 namespace MatrixFreeFunctions
2209 template <
int dim,
int degree>
2218 template <
int degree>
2221 static constexpr unsigned int value = degree + 1;
2237 template <
bool is_face,
2240 typename VectorizedArrayType>
2243 extract_initialization_data(
2245 const unsigned int dof_no,
2246 const unsigned int first_selected_component,
2247 const unsigned int quad_no,
2248 const unsigned int fe_degree,
2249 const unsigned int n_q_points,
2250 const unsigned int active_fe_index_given,
2251 const unsigned int active_quad_index_given,
2252 const unsigned int face_type)
2255 InitializationData init_data;
2259 &internal::MatrixFreeFunctions::
2260 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2268 active_fe_index_given :
2274 active_quad_index_given :
2275 std::min<unsigned int>(
2278 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2307 typename VectorizedArrayType>
2312 VectorizedArrayType>::
2315 const unsigned int dof_no,
2316 const unsigned int first_selected_component,
2317 const unsigned int quad_no,
2318 const unsigned int fe_degree,
2319 const unsigned int n_q_points,
2320 const bool is_interior_face,
2321 const unsigned int active_fe_index,
2322 const unsigned int active_quad_index,
2323 const unsigned int face_type)
2325 internal::extract_initialization_data<is_face>(matrix_free,
2327 first_selected_component,
2336 first_selected_component)
2337 , scratch_data_array(matrix_free.acquire_scratch_data())
2338 , matrix_free(&matrix_free)
2340 this->set_data_pointers(scratch_data_array, n_components_);
2342 this->dof_info->start_components.back() == 1 ||
2343 static_cast<int>(n_components_) <=
2345 this->dof_info->start_components
2346 [this->dof_info->component_to_base_index[first_selected_component] +
2348 first_selected_component,
2350 "You tried to construct a vector-valued evaluator with " +
2351 std::to_string(n_components) +
2352 " components. However, "
2353 "the current base element has only " +
2355 this->dof_info->start_components
2356 [this->dof_info->component_to_base_index[first_selected_component] +
2358 first_selected_component) +
2359 " components left when starting from local element index " +
2361 first_selected_component -
2362 this->dof_info->start_components
2363 [this->dof_info->component_to_base_index[first_selected_component]]) +
2364 " (global index " + std::to_string(first_selected_component) +
")"));
2376 typename VectorizedArrayType>
2381 VectorizedArrayType>::
2387 const unsigned int first_selected_component,
2391 other->mapped_geometry->get_quadrature() == quadrature ?
2392 other->mapped_geometry :
2394 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2399 first_selected_component)
2400 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2401 , matrix_free(nullptr)
2403 const unsigned int base_element_number =
2407 first_selected_component >=
2409 ExcMessage(
"The underlying element must at least contain as many "
2410 "components as requested by this class"));
2411 (void)base_element_number;
2416 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2418 fe.component_to_base_index(first_selected_component).
first);
2420 this->set_data_pointers(scratch_data_array, n_components_);
2429 typename VectorizedArrayType>
2434 VectorizedArrayType>::
2439 VectorizedArrayType> &other)
2441 , scratch_data_array(other.matrix_free == nullptr ?
2443 other.matrix_free->acquire_scratch_data())
2444 , matrix_free(other.matrix_free)
2446 if (other.matrix_free ==
nullptr)
2454 this->mapped_geometry =
2455 std::make_shared<internal::MatrixFreeFunctions::
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2460 this->mapping_data = &this->mapped_geometry->get_data_storage();
2464 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2466 this->mapped_geometry->get_data_storage().JxW_values.begin();
2467 this->jacobian_gradients =
2468 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2469 this->jacobian_gradients_non_inverse =
2470 this->mapped_geometry->get_data_storage()
2471 .jacobian_gradients_non_inverse[0]
2474 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2477 this->set_data_pointers(scratch_data_array, n_components_);
2486 typename VectorizedArrayType>
2491 VectorizedArrayType> &
2497 VectorizedArrayType> &other)
2500 if (matrix_free ==
nullptr)
2503 delete scratch_data_array;
2512 matrix_free = other.matrix_free;
2514 if (other.matrix_free ==
nullptr)
2523 this->mapped_geometry =
2524 std::make_shared<internal::MatrixFreeFunctions::
2525 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2530 this->mapping_data = &this->mapped_geometry->get_data_storage();
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();
2549 this->set_data_pointers(scratch_data_array, n_components_);
2560 typename VectorizedArrayType>
2565 VectorizedArrayType>::~FEEvaluationBase()
2567 if (matrix_free !=
nullptr)
2578 delete scratch_data_array;
2589 typename VectorizedArrayType>
2594 Assert(matrix_free !=
nullptr,
2596 "FEEvaluation was not initialized with a MatrixFree object!"));
2597 return *matrix_free;
2606 template <
typename VectorType,
bool>
2607 struct ConstBlockVectorSelector;
2609 template <
typename VectorType>
2610 struct ConstBlockVectorSelector<VectorType, true>
2612 using BaseVectorType =
const typename VectorType::BlockType;
2615 template <
typename VectorType>
2616 struct ConstBlockVectorSelector<VectorType, false>
2618 using BaseVectorType =
typename VectorType::BlockType;
2624 template <
typename VectorType,
bool>
2625 struct BlockVectorSelector;
2627 template <
typename VectorType>
2628 struct BlockVectorSelector<VectorType, true>
2630 using BaseVectorType =
typename ConstBlockVectorSelector<
2632 std::is_const_v<VectorType>>::BaseVectorType;
2634 static BaseVectorType *
2635 get_vector_component(VectorType &vec,
const unsigned int component)
2638 return &vec.block(component);
2642 template <
typename VectorType>
2643 struct BlockVectorSelector<VectorType, false>
2645 using BaseVectorType = VectorType;
2647 static BaseVectorType *
2648 get_vector_component(VectorType &vec,
const unsigned int component)
2667 template <
typename VectorType>
2668 struct BlockVectorSelector<
std::vector<VectorType>, false>
2670 using BaseVectorType = VectorType;
2672 static BaseVectorType *
2673 get_vector_component(std::vector<VectorType> &vec,
2674 const unsigned int component)
2677 return &vec[component];
2681 template <
typename VectorType>
2682 struct BlockVectorSelector<const
std::vector<VectorType>, false>
2684 using BaseVectorType =
const VectorType;
2686 static const BaseVectorType *
2687 get_vector_component(
const std::vector<VectorType> &vec,
2688 const unsigned int component)
2691 return &vec[component];
2695 template <
typename VectorType>
2696 struct BlockVectorSelector<
std::vector<VectorType *>, false>
2698 using BaseVectorType = VectorType;
2700 static BaseVectorType *
2701 get_vector_component(std::vector<VectorType *> &vec,
2702 const unsigned int component)
2705 return vec[component];
2709 template <
typename VectorType>
2710 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
2712 using BaseVectorType =
const VectorType;
2714 static const BaseVectorType *
2715 get_vector_component(
const std::vector<VectorType *> &vec,
2716 const unsigned int component)
2719 return vec[component];
2730 typename VectorizedArrayType>
2731template <
typename VectorType,
typename VectorOperation>
2736 const std::array<VectorType *, n_components_> &src,
2739 n_components_> &src_sm,
2740 const std::bitset<n_lanes> &
mask,
2741 const bool apply_constraints)
const
2746 if (this->matrix_free ==
nullptr)
2748 read_write_operation_global(operation, src);
2755 if (this->n_fe_components == 1)
2756 for (
unsigned int comp = 0; comp < n_components; ++comp)
2758 Assert(src[comp] !=
nullptr,
2759 ExcMessage(
"The finite element underlying this FEEvaluation "
2760 "object is scalar, but you requested " +
2761 std::to_string(n_components) +
2762 " components via the template argument in "
2763 "FEEvaluation. In that case, you must pass an "
2764 "std::vector<VectorType> or a BlockVector to " +
2765 "read_dof_values and distribute_local_to_global."));
2777 const bool accesses_exterior_dofs =
2778 this->dof_access_index ==
2780 this->is_interior_face() ==
false;
2790 bool is_contiguous =
true;
2792 if (accesses_exterior_dofs)
2794 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2795 const unsigned int n_filled_lanes =
2800 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
2801 if (
mask[v] ==
true &&
2804 [cells[v] / n_lanes] <
2807 is_contiguous = false;
2811 this->dof_access_index :
2812 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2813 [this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2814 IndexStorageVariants::contiguous)
2816 is_contiguous =
false;
2821 read_write_operation_contiguous(operation, src, src_sm,
mask);
2828 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2830 const bool masking_is_active =
mask.count() < n_lanes;
2831 if (masking_is_active)
2832 for (
unsigned int v = 0; v < n_lanes; ++v)
2833 if (
mask[v] ==
false)
2836 bool has_hn_constraints =
false;
2838 if (is_face ==
false)
2844 [this->first_selected_component])
2845 for (
unsigned int v = 0; v < n_lanes; ++v)
2850 has_hn_constraints = true;
2853 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2856 const bool use_vectorized_path =
2857 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
2859 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
2860 std::array<VectorizedArrayType *, n_components> values_dofs;
2861 for (
unsigned int c = 0; c < n_components; ++c)
2862 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2863 c * dofs_per_component;
2867 [is_face ? this->dof_access_index :
2868 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2869 [this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2870 IndexStorageVariants::interleaved &&
2871 use_vectorized_path)
2873 const unsigned int *dof_indices =
2875 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2879 [this->first_selected_component] *
2882 std::array<typename VectorType::value_type *, n_components> src_ptrs;
2883 if (n_components == 1 || this->n_fe_components == 1)
2884 for (
unsigned int comp = 0; comp < n_components; ++comp)
2886 const_cast<typename VectorType::value_type *
>(src[comp]->
begin());
2889 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2891 if (n_components == 1 || this->n_fe_components == 1)
2892 for (
unsigned int i = 0; i < dofs_per_component;
2893 ++i, dof_indices += n_lanes)
2894 for (
unsigned int comp = 0; comp < n_components; ++comp)
2895 operation.process_dof_gather(dof_indices,
2899 values_dofs[comp][i],
2902 for (
unsigned int comp = 0; comp < n_components; ++comp)
2903 for (
unsigned int i = 0; i < dofs_per_component;
2904 ++i, dof_indices += n_lanes)
2905 operation.process_dof_gather(dof_indices,
2909 values_dofs[comp][i],
2916 std::array<const unsigned int *, n_lanes> dof_indices;
2917 dof_indices.fill(
nullptr);
2922 bool has_constraints =
false;
2923 const unsigned int n_components_read =
2924 this->n_fe_components > 1 ? n_components : 1;
2928 for (
unsigned int v = 0; v < n_lanes; ++v)
2934 const std::pair<unsigned int, unsigned int> *my_index_start =
2935 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2936 this->first_selected_component];
2941 if (my_index_start[n_components_read].
second !=
2942 my_index_start[0].
second)
2943 has_constraints =
true;
2946 dof_info.
dof_indices.data() + my_index_start[0].first;
2951 for (
unsigned int v = 0; v < n_lanes; ++v)
2956 const std::pair<unsigned int, unsigned int> *my_index_start =
2957 &dof_info.
row_starts[cells[v] * this->n_fe_components +
2958 this->first_selected_component];
2959 if (my_index_start[n_components_read].
second !=
2960 my_index_start[0].
second)
2961 has_constraints =
true;
2968 dof_info.hanging_node_constraint_masks_comp
2969 [this->active_fe_index][this->first_selected_component])
2970 has_hn_constraints = true;
2973 my_index_start[0].
first ||
2976 my_index_start[0].
first,
2979 dof_info.
dof_indices.data() + my_index_start[0].first;
2983 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
2984 return i != numbers::invalid_unsigned_int;
2986 for (
unsigned int comp = 0; comp < n_components; ++comp)
2987 for (
unsigned int i = 0; i < dofs_per_component; ++i)
2988 operation.process_empty(values_dofs[comp][i]);
2992 if (!has_constraints && apply_constraints)
2994 if (n_components == 1 || this->n_fe_components == 1)
2996 for (
unsigned int v = 0; v < n_lanes; ++v)
3001 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3002 for (
unsigned int comp = 0; comp < n_components; ++comp)
3003 operation.process_dof(dof_indices[v][i],
3005 values_dofs[comp][i][v]);
3010 for (
unsigned int comp = 0; comp < n_components; ++comp)
3011 for (
unsigned int v = 0; v < n_lanes; ++v)
3016 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3017 operation.process_dof(
3018 dof_indices[v][comp * dofs_per_component + i],
3020 values_dofs[comp][i][v]);
3032 for (
unsigned int v = 0; v < n_lanes; ++v)
3038 const unsigned int cell_dof_index =
3039 cell_index * this->n_fe_components + this->first_selected_component;
3040 const unsigned int n_components_read =
3041 this->n_fe_components > 1 ? n_components : 1;
3042 unsigned int index_indicators =
3044 unsigned int next_index_indicators =
3045 dof_info.
row_starts[cell_dof_index + 1].second;
3049 if (apply_constraints ==
false &&
3050 (dof_info.
row_starts[cell_dof_index].second !=
3051 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3057 dof_info.hanging_node_constraint_masks_comp
3058 [this->active_fe_index][this->first_selected_component])))
3067 [this->first_selected_component] +
3069 next_index_indicators = index_indicators;
3072 if (n_components == 1 || this->n_fe_components == 1)
3074 unsigned int ind_local = 0;
3075 for (; index_indicators != next_index_indicators; ++index_indicators)
3077 const std::pair<unsigned short, unsigned short> indicator =
3080 for (
unsigned int j = 0; j < indicator.first; ++j)
3081 for (
unsigned int comp = 0; comp < n_components; ++comp)
3082 operation.process_dof(dof_indices[v][j],
3084 values_dofs[comp][ind_local + j][v]);
3086 ind_local += indicator.first;
3087 dof_indices[v] += indicator.first;
3091 Number
value[n_components];
3092 for (
unsigned int comp = 0; comp < n_components; ++comp)
3093 operation.pre_constraints(values_dofs[comp][ind_local][v],
3096 const Number *data_val =
3098 const Number *end_pool =
3100 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3101 for (
unsigned int comp = 0; comp < n_components; ++comp)
3102 operation.process_constraint(*dof_indices[v],
3107 for (
unsigned int comp = 0; comp < n_components; ++comp)
3108 operation.post_constraints(
value[comp],
3109 values_dofs[comp][ind_local][v]);
3115 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3116 for (
unsigned int comp = 0; comp < n_components; ++comp)
3117 operation.process_dof(*dof_indices[v],
3119 values_dofs[comp][ind_local][v]);
3128 for (
unsigned int comp = 0; comp < n_components; ++comp)
3130 unsigned int ind_local = 0;
3133 for (; index_indicators != next_index_indicators;
3136 const std::pair<unsigned short, unsigned short> indicator =
3140 for (
unsigned int j = 0; j < indicator.first; ++j)
3141 operation.process_dof(dof_indices[v][j],
3143 values_dofs[comp][ind_local + j][v]);
3144 ind_local += indicator.first;
3145 dof_indices[v] += indicator.first;
3150 operation.pre_constraints(values_dofs[comp][ind_local][v],
3153 const Number *data_val =
3155 const Number *end_pool =
3158 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3159 operation.process_constraint(*dof_indices[v],
3164 operation.post_constraints(
value,
3165 values_dofs[comp][ind_local][v]);
3172 for (; ind_local < dofs_per_component;
3173 ++dof_indices[v], ++ind_local)
3176 operation.process_dof(*dof_indices[v],
3178 values_dofs[comp][ind_local][v]);
3181 if (apply_constraints ==
true && comp + 1 < n_components)
3182 next_index_indicators =
3183 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3195 typename VectorizedArrayType>
3196template <
typename VectorType,
typename VectorOperation>
3201 const std::array<VectorType *, n_components_> &src)
const
3205 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3206 unsigned int index = this->first_selected_component * dofs_per_component;
3207 for (
unsigned int comp = 0; comp < n_components; ++comp)
3209 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3211 operation.process_empty(
3212 this->values_dofs[comp * dofs_per_component + i]);
3213 operation.process_dof_global(
3214 local_dof_indices[this->data->lexicographic_numbering[
index]],
3216 this->values_dofs[comp * dofs_per_component + i][0]);
3227 typename VectorizedArrayType>
3228template <
typename VectorType,
typename VectorOperation>
3233 const std::array<VectorType *, n_components_> &src,
3236 n_components_> &vectors_sm,
3237 const std::bitset<n_lanes> &
mask)
const
3246 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3249 is_face ? this->dof_access_index :
3251 const unsigned int n_active_lanes =
mask.count();
3254 const std::vector<unsigned int> &dof_indices_cont =
3257 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3258 std::array<VectorizedArrayType *, n_components> values_dofs;
3259 for (
unsigned int c = 0; c < n_components; ++c)
3260 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3261 c * dofs_per_component;
3265 const bool accesses_exterior_dofs =
3266 this->dof_access_index ==
3268 this->is_interior_face() ==
false;
3274 interleaved_contiguous &&
3275 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3277 const unsigned int dof_index =
3278 dof_indices_cont[this->cell * n_lanes] +
3281 [this->first_selected_component] *
3283 if (n_components == 1 || this->n_fe_components == 1)
3284 for (
unsigned int comp = 0; comp < n_components; ++comp)
3285 operation.process_dofs_vectorized(dofs_per_component,
3291 operation.process_dofs_vectorized(dofs_per_component * n_components,
3299 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3303 const unsigned int n_filled_lanes =
3306 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3307 n_active_lanes == n_lanes &&
3308 !accesses_exterior_dofs;
3310 if (vectors_sm[0] !=
nullptr)
3312 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3313 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3315 const auto upper_bound =
3316 std::min<unsigned int>(n_filled_lanes, n_lanes);
3317 for (
unsigned int v = 0; v < upper_bound; ++v)
3319 if (
mask[v] ==
false)
3321 vector_ptrs[v] =
nullptr;
3341 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3342 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3344 [this->active_fe_index][this->first_selected_component]);
3346 vector_ptrs[v] =
nullptr;
3348 for (
unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3349 vector_ptrs[v] =
nullptr;
3354 if (use_vectorized_path)
3356 if (n_components == 1 || this->n_fe_components == 1)
3358 for (
unsigned int comp = 0; comp < n_components; ++comp)
3360 auto vector_ptrs = compute_vector_ptrs(comp);
3361 operation.process_dofs_vectorized_transpose(
3370 auto vector_ptrs = compute_vector_ptrs(0);
3371 operation.process_dofs_vectorized_transpose(dofs_per_component *
3379 for (
unsigned int comp = 0; comp < n_components; ++comp)
3381 auto vector_ptrs = compute_vector_ptrs(
3382 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3384 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3385 operation.process_empty(values_dofs[comp][i]);
3387 if (n_components == 1 || this->n_fe_components == 1)
3389 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3390 if (
mask[v] ==
true)
3391 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3392 operation.process_dof(vector_ptrs[v][i],
3393 values_dofs[comp][i][v]);
3397 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3398 if (
mask[v] ==
true)
3399 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3400 operation.process_dof(
3401 vector_ptrs[v][i + comp * dofs_per_component],
3402 values_dofs[comp][i][v]);
3408 std::array<unsigned int, n_lanes> dof_indices;
3409 std::fill(dof_indices.begin(),
3414 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3418 if (
mask[v] ==
true)
3420 dof_indices_cont[cells[v]] +
3423 [this->first_selected_component] *
3429 if (use_vectorized_path)
3435 if (n_components == 1 || this->n_fe_components == 1)
3436 for (
unsigned int comp = 0; comp < n_components; ++comp)
3437 operation.process_dofs_vectorized_transpose(dofs_per_component,
3443 operation.process_dofs_vectorized_transpose(dofs_per_component *
3452 interleaved_contiguous_strided)
3454 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3455 if (n_components == 1 || this->n_fe_components == 1)
3456 for (
unsigned int comp = 0; comp < n_components; ++comp)
3457 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3458 src[comp]->
begin());
3461 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3463 if (n_components == 1 || this->n_fe_components == 1)
3464 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3466 for (
unsigned int comp = 0; comp < n_components; ++comp)
3467 operation.process_dof_gather(dof_indices.data(),
3470 src_ptrs[comp] + i * n_lanes,
3471 values_dofs[comp][i],
3475 for (
unsigned int comp = 0; comp < n_components; ++comp)
3476 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3478 operation.process_dof_gather(
3481 (comp * dofs_per_component + i) * n_lanes,
3482 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
3483 values_dofs[comp][i],
3491 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3493 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3494 if (n_components == 1 || this->n_fe_components == 1)
3495 for (
unsigned int comp = 0; comp < n_components; ++comp)
3496 src_ptrs[comp] =
const_cast<typename VectorType::value_type *
>(
3497 src[comp]->
begin());
3500 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3502 const unsigned int *offsets =
3504 if (n_components == 1 || this->n_fe_components == 1)
3505 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3507 for (
unsigned int comp = 0; comp < n_components; ++comp)
3508 operation.process_dof_gather(dof_indices.data(),
3512 values_dofs[comp][i],
3515 for (
unsigned int v = 0; v < n_lanes; ++v)
3516 dof_indices[v] += offsets[v];
3519 for (
unsigned int comp = 0; comp < n_components; ++comp)
3520 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3522 operation.process_dof_gather(dof_indices.data(),
3526 values_dofs[comp][i],
3529 for (
unsigned int v = 0; v < n_lanes; ++v)
3530 dof_indices[v] += offsets[v];
3535 for (
unsigned int comp = 0; comp < n_components; ++comp)
3537 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3538 operation.process_empty(values_dofs[comp][i]);
3539 if (accesses_exterior_dofs)
3541 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3542 if (
mask[v] ==
true)
3545 [ind][cells[v] / VectorizedArrayType::size()] ==
3549 if (n_components == 1 || this->n_fe_components == 1)
3551 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3552 operation.process_dof(dof_indices[v] + i,
3554 values_dofs[comp][i][v]);
3558 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3559 operation.process_dof(dof_indices[v] + i +
3560 comp * dofs_per_component,
3562 values_dofs[comp][i][v]);
3567 const unsigned int offset =
3570 if (n_components == 1 || this->n_fe_components == 1)
3572 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3573 operation.process_dof(dof_indices[v] + i * offset,
3575 values_dofs[comp][i][v]);
3579 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3580 operation.process_dof(
3582 (i + comp * dofs_per_component) * offset,
3584 values_dofs[comp][i][v]);
3595 if (n_components == 1 || this->n_fe_components == 1)
3597 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3598 if (
mask[v] ==
true)
3599 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3600 operation.process_dof(dof_indices[v] + i,
3602 values_dofs[comp][i][v]);
3606 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3607 if (
mask[v] ==
true)
3608 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3609 operation.process_dof(dof_indices[v] + i +
3610 comp * dofs_per_component,
3612 values_dofs[comp][i][v]);
3617 const unsigned int *offsets =
3619 [ind][VectorizedArrayType::size() * this->cell];
3620 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3622 if (n_components == 1 || this->n_fe_components == 1)
3623 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3625 if (
mask[v] ==
true)
3626 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3627 operation.process_dof(dof_indices[v] + i * offsets[v],
3629 values_dofs[comp][i][v]);
3633 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3634 if (
mask[v] ==
true)
3635 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3636 operation.process_dof(
3638 (i + comp * dofs_per_component) * offsets[v],
3640 values_dofs[comp][i][v]);
3651 typename VectorType,
3652 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
3653 decltype(std::declval<VectorType>().begin())
3654 get_beginning(VectorType &vec)
3661 typename VectorType,
3662 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
3663 typename VectorType::value_type *
3664 get_beginning(VectorType &)
3669 template <
typename VectorType,
3670 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
3672 const std::vector<ArrayView<const typename VectorType::value_type>> *
3673 get_shared_vector_data(VectorType *vec,
3674 const bool is_valid_mode_for_sm,
3675 const unsigned int active_fe_index,
3679 if (is_valid_mode_for_sm &&
3682 active_fe_index == 0)
3683 return &vec->shared_vector_data();
3688 template <
typename VectorType,
3689 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
3691 const std::vector<ArrayView<const typename VectorType::value_type>> *
3692 get_shared_vector_data(VectorType *,
3700 template <
int n_components,
typename VectorType>
3702 std::array<
typename internal::BlockVectorSelector<
3707 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3711 get_vector_data(VectorType &src,
3712 const unsigned int first_index,
3713 const bool is_valid_mode_for_sm,
3714 const unsigned int active_fe_index,
3720 std::array<
typename internal::BlockVectorSelector<
3726 ArrayView<
const typename internal::BlockVectorSelector<
3732 for (
unsigned int d = 0;
d < n_components; ++
d)
3733 src_data.first[d] = internal::BlockVectorSelector<
3739 for (
unsigned int d = 0;
d < n_components; ++
d)
3740 src_data.second[d] = get_shared_vector_data(
3741 const_cast<typename internal::BlockVectorSelector<
3742 std::remove_const_t<VectorType>,
3744 *>(src_data.first[d]),
3745 is_valid_mode_for_sm,
3759 typename VectorizedArrayType>
3764 if (this->dof_info ==
nullptr ||
3766 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3767 this->dof_info->hanging_node_constraint_masks_comp
3768 [this->active_fe_index][this->first_selected_component] ==
false)
3771 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3774 bool hn_available =
false;
3776 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3778 for (
unsigned int v = 0; v < n_lanes; ++v)
3790 constraint_mask[v] =
mask;
3796 if (hn_available ==
false)
3801 this->data->data.front().fe_degree,
3802 this->get_shape_info(),
3814 typename VectorizedArrayType>
3815template <
typename VectorType>
3819 const unsigned int first_index,
3820 const std::bitset<n_lanes> &
mask)
3822 const auto src_data = internal::get_vector_data<n_components_>(
3825 this->dof_access_index ==
3827 this->active_fe_index,
3831 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
3833 apply_hanging_node_constraints(
false);
3836 this->dof_values_initialized =
true;
3846 typename VectorizedArrayType>
3847template <
typename VectorType>
3851 const unsigned int first_index,
3852 const std::bitset<n_lanes> &
mask)
3854 const auto src_data = internal::get_vector_data<n_components_>(
3857 this->dof_access_index ==
3859 this->active_fe_index,
3863 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
3866 this->dof_values_initialized =
true;
3876 typename VectorizedArrayType>
3877template <
typename VectorType>
3881 const unsigned int first_index,
3882 const std::bitset<n_lanes> &
mask)
const
3885 Assert(this->dof_values_initialized ==
true,
3889 apply_hanging_node_constraints(
true);
3891 const auto dst_data = internal::get_vector_data<n_components_>(
3894 this->dof_access_index ==
3896 this->active_fe_index,
3901 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
3910 typename VectorizedArrayType>
3911template <
typename VectorType>
3915 const unsigned int first_index,
3916 const std::bitset<n_lanes> &
mask)
const
3919 Assert(this->dof_values_initialized ==
true,
3923 const auto dst_data = internal::get_vector_data<n_components_>(
3926 this->dof_access_index ==
3928 this->active_fe_index,
3932 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
3941 typename VectorizedArrayType>
3942template <
typename VectorType>
3946 const unsigned int first_index,
3947 const std::bitset<n_lanes> &
mask)
const
3950 Assert(this->dof_values_initialized ==
true,
3954 const auto dst_data = internal::get_vector_data<n_components_>(
3957 this->dof_access_index ==
3959 this->active_fe_index,
3963 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
3976 typename VectorizedArrayType>
3982 VectorizedArrayType>::value_type
3987 if constexpr (n_components == 1)
3988 return this->values_dofs[dof];
3991 const std::size_t dofs = this->data->dofs_per_component_on_cell;
3993 for (
unsigned int comp = 0; comp < n_components; ++comp)
3994 return_value[comp] = this->values_dofs[comp * dofs + dof];
3995 return return_value;
4005 typename VectorizedArrayType>
4011 VectorizedArrayType>::value_type
4013 get_value(
const unsigned int q_point)
const
4016 Assert(this->values_quad_initialized ==
true,
4021 if constexpr (n_components == 1)
4022 return this->values_quad[q_point];
4025 if (n_components == dim &&
4026 this->data->element_type ==
4031 Assert(this->values_quad_initialized ==
true,
4036 Assert(this->J_value !=
nullptr,
4039 const std::size_t nqp = this->n_quadrature_points;
4047 const VectorizedArrayType inv_det =
4048 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4049 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4050 this->jacobian[0][2][2];
4053 for (
unsigned int comp = 0; comp < n_components; ++comp)
4054 value_out[comp] = this->values_quad[comp * nqp + q_point] *
4055 jac[comp][comp] * inv_det;
4062 this->jacobian[q_point] :
4071 const VectorizedArrayType inv_det =
4072 (is_face && dim == 2 && this->get_face_no() < 2) ?
4076 for (
unsigned int comp = 0; comp < n_components; ++comp)
4078 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
4079 for (
unsigned int e = 1;
e < dim; ++
e)
4081 this->values_quad[e * nqp + q_point] * jac[comp][e];
4082 value_out[comp] *= inv_det;
4089 const std::size_t nqp = this->n_quadrature_points;
4091 for (
unsigned int comp = 0; comp < n_components; ++comp)
4092 return_value[comp] = this->values_quad[comp * nqp + q_point];
4093 return return_value;
4104 typename VectorizedArrayType>
4110 VectorizedArrayType>::gradient_type
4115 Assert(this->gradients_quad_initialized ==
true,
4120 Assert(this->jacobian !=
nullptr,
4122 "update_gradients"));
4123 const std::size_t nqp = this->n_quadrature_points;
4125 if constexpr (n_components == dim && dim > 1)
4127 if (this->data->element_type ==
4132 Assert(this->gradients_quad_initialized ==
true,
4137 Assert(this->jacobian !=
nullptr,
4139 "update_gradients"));
4140 const std::size_t nqp = this->n_quadrature_points;
4141 const std::size_t nqp_d = nqp * dim;
4144 this->gradients_quad + q_point * dim;
4155 const VectorizedArrayType inv_det =
4156 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4157 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
4158 this->jacobian[0][2][2];
4161 for (
unsigned int d = 0;
d < dim; ++
d)
4162 for (
unsigned int comp = 0; comp < n_components; ++comp)
4163 grad_out[comp][d] = gradients[comp * nqp_d + d] *
4165 (jac[comp][comp] * inv_det);
4177 const VectorizedArrayType inv_det =
4178 (is_face && dim == 2 && this->get_face_no() < 2) ?
4182 VectorizedArrayType tmp[dim][dim];
4184 for (
unsigned int d = 0;
d < dim; ++
d)
4185 for (
unsigned int e = 0;
e < dim; ++
e)
4188 for (
unsigned int f = 1; f < dim; ++f)
4189 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
4191 for (
unsigned int comp = 0; comp < n_components; ++comp)
4192 for (
unsigned int d = 0;
d < dim; ++
d)
4194 VectorizedArrayType res = jac[comp][0] * tmp[
d][0];
4195 for (
unsigned int f = 1; f < dim; ++f)
4196 res += jac[comp][f] * tmp[d][f];
4198 grad_out[comp][
d] = res * inv_det;
4209 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4211 "update_hessians"));
4213 const auto jac_grad =
4214 this->jacobian_gradients_non_inverse[q_point];
4216 this->jacobian[q_point];
4220 const VectorizedArrayType inv_det =
4221 (is_face && dim == 2 && this->get_face_no() < 2) ?
4228 VectorizedArrayType tmp[dim][dim];
4229 for (
unsigned int d = 0;
d < dim; ++
d)
4230 for (
unsigned int e = 0;
e < dim; ++
e)
4233 for (
unsigned int f = 1; f < dim; ++f)
4234 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
4239 for (
unsigned int d = 0;
d < dim; ++
d)
4241 for (
unsigned int e = 0;
e < dim; ++
e)
4243 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
4244 for (
unsigned int f = 0, r = dim; f < dim; ++f)
4245 for (
unsigned int k = f + 1; k < dim; ++k, ++r)
4248 jac_grad[r][
d] * this->values_quad[f * nqp + q_point];
4250 jac_grad[r][
d] * this->values_quad[k * nqp + q_point];
4255 for (
unsigned int d = 0;
d < dim; ++
d)
4256 for (
unsigned int e = 0;
e < dim; ++
e)
4258 VectorizedArrayType res = tmp[0][
d] * inv_t_jac[
e][0];
4259 for (
unsigned int f = 1; f < dim; ++f)
4260 res += tmp[f][d] * inv_t_jac[e][f];
4261 grad_out[
d][
e] = res;
4267 VectorizedArrayType tmp3[dim], tmp4[dim];
4268 for (
unsigned int d = 0;
d < dim; ++
d)
4270 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
4271 for (
unsigned int e = 1;
e < dim; ++
e)
4272 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
4274 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
4275 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
4276 for (
unsigned int d = 0;
d < dim; ++
d)
4278 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
4279 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
4281 for (
unsigned int d = 0;
d < dim; ++
d)
4283 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
4284 for (
unsigned int e = 1;
e < dim; ++
e)
4285 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
4288 VectorizedArrayType tmp2[dim];
4289 for (
unsigned int d = 0;
d < dim; ++
d)
4291 tmp2[
d] = t_jac[0][
d] * this->values_quad[q_point];
4292 for (
unsigned e = 1;
e < dim; ++
e)
4294 t_jac[e][d] * this->values_quad[e * nqp + q_point];
4297 for (
unsigned int d = 0;
d < dim; ++
d)
4298 for (
unsigned int e = 0;
e < dim; ++
e)
4300 grad_out[
d][
e] -= tmp4[
e] * tmp2[
d];
4304 grad_out[
d][
e] *= inv_det;
4315 for (
unsigned int comp = 0; comp < n_components; ++comp)
4316 for (
unsigned int d = 0;
d < dim; ++
d)
4318 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4319 this->jacobian[0][
d][
d];
4328 for (
unsigned int comp = 0; comp < n_components; ++comp)
4329 for (
unsigned int d = 0;
d < dim; ++
d)
4332 jac[
d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4333 for (
unsigned int e = 1;
e < dim; ++
e)
4334 grad_out[comp][d] +=
4336 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4339 if constexpr (n_components == 1)
4351 typename VectorizedArrayType>
4357 VectorizedArrayType>::value_type
4363 Assert(this->gradients_quad_initialized ==
true,
4367 Assert(this->normal_x_jacobian !=
nullptr,
4369 "update_gradients"));
4371 const std::size_t nqp = this->n_quadrature_points;
4375 for (
unsigned int comp = 0; comp < n_components; ++comp)
4377 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4378 (this->normal_x_jacobian[0][dim - 1]);
4381 const std::size_t
index =
4383 for (
unsigned int comp = 0; comp < n_components; ++comp)
4385 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4386 this->normal_x_jacobian[
index][0];
4387 for (
unsigned int d = 1;
d < dim; ++
d)
4389 this->gradients_quad[(comp * nqp + q_point) * dim +
d] *
4390 this->normal_x_jacobian[
index][
d];
4393 if constexpr (n_components == 1)
4405 template <
typename VectorizedArrayType>
4408 const VectorizedArrayType *
const hessians,
4410 VectorizedArrayType (&tmp)[1][1])
4412 tmp[0][0] = jac[0][0] *
hessians[0];
4415 template <
typename VectorizedArrayType>
4418 const VectorizedArrayType *
const hessians,
4419 const unsigned int nqp,
4420 VectorizedArrayType (&tmp)[2][2])
4422 for (
unsigned int d = 0;
d < 2; ++
d)
4430 template <
typename VectorizedArrayType>
4433 const VectorizedArrayType *
const hessians,
4434 const unsigned int nqp,
4435 VectorizedArrayType (&tmp)[3][3])
4437 for (
unsigned int d = 0;
d < 3; ++
d)
4458 typename VectorizedArrayType>
4463 VectorizedArrayType>::hessian_type
4468 Assert(this->hessians_quad_initialized ==
true,
4473 Assert(this->jacobian !=
nullptr,
4483 const std::size_t nqp = this->n_quadrature_points;
4484 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4489 for (
unsigned int comp = 0; comp < n_components; ++comp)
4491 for (
unsigned int d = 0;
d < dim; ++
d)
4492 hessian_out[comp][d][d] =
4493 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4494 (jac[
d][
d] * jac[
d][
d]);
4500 hessian_out[comp][0][1] =
4501 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4502 (jac[0][0] * jac[1][1]);
4505 hessian_out[comp][0][1] =
4506 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4507 (jac[0][0] * jac[1][1]);
4508 hessian_out[comp][0][2] =
4509 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4510 (jac[0][0] * jac[2][2]);
4511 hessian_out[comp][1][2] =
4512 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4513 (jac[1][1] * jac[2][2]);
4518 for (
unsigned int d = 0;
d < dim; ++
d)
4519 for (
unsigned int e = d + 1;
e < dim; ++
e)
4520 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4526 for (
unsigned int comp = 0; comp < n_components; ++comp)
4528 VectorizedArrayType tmp[dim][dim];
4529 internal::hessian_unit_times_jac(
4530 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4533 for (
unsigned int d = 0;
d < dim; ++
d)
4534 for (
unsigned int e = d;
e < dim; ++
e)
4536 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4537 for (
unsigned int f = 1; f < dim; ++f)
4538 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4545 for (
unsigned int d = 0;
d < dim; ++
d)
4546 for (
unsigned int e = d + 1;
e < dim; ++
e)
4547 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4553 const auto &jac_grad = this->jacobian_gradients[q_point];
4554 for (
unsigned int comp = 0; comp < n_components; ++comp)
4556 VectorizedArrayType tmp[dim][dim];
4557 internal::hessian_unit_times_jac(
4558 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4561 for (
unsigned int d = 0;
d < dim; ++
d)
4562 for (
unsigned int e = d;
e < dim; ++
e)
4564 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4565 for (
unsigned int f = 1; f < dim; ++f)
4566 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4570 for (
unsigned int d = 0;
d < dim; ++
d)
4571 for (
unsigned int e = 0;
e < dim; ++
e)
4572 hessian_out[comp][d][d] +=
4574 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4577 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4578 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4579 for (
unsigned int f = 0; f < dim; ++f)
4580 hessian_out[comp][d][e] +=
4581 jac_grad[count][f] *
4582 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4585 for (
unsigned int d = 0;
d < dim; ++
d)
4586 for (
unsigned int e = d + 1;
e < dim; ++
e)
4587 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4590 if constexpr (n_components == 1)
4591 return hessian_out[0];
4602 typename VectorizedArrayType>
4607 VectorizedArrayType>::gradient_type
4613 Assert(this->hessians_quad_initialized ==
true,
4624 const std::size_t nqp = this->n_quadrature_points;
4625 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4631 for (
unsigned int comp = 0; comp < n_components; ++comp)
4632 for (
unsigned int d = 0;
d < dim; ++
d)
4633 hessian_out[comp][d] =
4634 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4635 (jac[
d][
d] * jac[
d][
d]);
4640 for (
unsigned int comp = 0; comp < n_components; ++comp)
4644 VectorizedArrayType tmp[dim][dim];
4645 internal::hessian_unit_times_jac(
4646 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4650 for (
unsigned int d = 0;
d < dim; ++
d)
4652 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4653 for (
unsigned int f = 1; f < dim; ++f)
4654 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4661 const auto &jac_grad = this->jacobian_gradients[q_point];
4662 for (
unsigned int comp = 0; comp < n_components; ++comp)
4666 VectorizedArrayType tmp[dim][dim];
4667 internal::hessian_unit_times_jac(
4668 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4672 for (
unsigned int d = 0;
d < dim; ++
d)
4674 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4675 for (
unsigned int f = 1; f < dim; ++f)
4676 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4679 for (
unsigned int d = 0;
d < dim; ++
d)
4680 for (
unsigned int e = 0;
e < dim; ++
e)
4681 hessian_out[comp][d] +=
4683 this->gradients_quad[(comp * nqp + q_point) * dim +
e];
4687 if constexpr (n_components == 1)
4688 return hessian_out[0];
4699 typename VectorizedArrayType>
4704 VectorizedArrayType>::value_type
4710 Assert(this->hessians_quad_initialized ==
true,
4715 const gradient_type hess_diag = get_hessian_diagonal(q_point);
4716 if constexpr (n_components == 1)
4718 VectorizedArrayType
sum = hess_diag[0];
4719 for (
unsigned int d = 1;
d < dim; ++
d)
4720 sum += hess_diag[d];
4726 for (
unsigned int comp = 0; comp < n_components; ++comp)
4728 laplacian_out[comp] = hess_diag[comp][0];
4729 for (
unsigned int d = 1;
d < dim; ++
d)
4730 laplacian_out[comp] += hess_diag[comp][d];
4732 return laplacian_out;
4742 typename VectorizedArrayType>
4748 this->dof_values_initialized =
true;
4750 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4752 for (
unsigned int comp = 0; comp < n_components; ++comp)
4753 if constexpr (n_components == 1)
4754 this->values_dofs[comp * dofs + dof] = val_in;
4756 this->values_dofs[comp * dofs + dof] = val_in[comp];
4765 typename VectorizedArrayType>
4768 submit_value(
const value_type val_in,
const unsigned int q_point)
4774 Assert(this->J_value !=
nullptr,
4778 this->values_quad_submitted =
true;
4781 const std::size_t nqp = this->n_quadrature_points;
4782 VectorizedArrayType *
values = this->values_quad + q_point;
4784 const VectorizedArrayType JxW =
4786 this->J_value[0] * this->quadrature_weights[q_point] :
4787 this->J_value[q_point];
4788 if constexpr (n_components == 1)
4789 values[0] = val_in * JxW;
4792 if (n_components == dim &&
4793 this->data->element_type ==
4798 Assert(this->J_value !=
nullptr,
4803 this->values_quad_submitted =
true;
4806 VectorizedArrayType *
values = this->values_quad + q_point;
4807 const std::size_t nqp = this->n_quadrature_points;
4813 const VectorizedArrayType weight =
4814 this->quadrature_weights[q_point];
4816 for (
unsigned int comp = 0; comp < n_components; ++comp)
4817 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
4824 this->jacobian[q_point] :
4829 const VectorizedArrayType fac =
4831 this->quadrature_weights[q_point] :
4833 this->J_value[q_point] :
4834 this->J_value[0] * this->quadrature_weights[q_point]) *
4835 ((dim == 2 && this->get_face_no() < 2) ?
4844 for (
unsigned int comp = 0; comp < n_components; ++comp)
4846 values[comp * nqp] = val_in[0] * jac[0][comp];
4847 for (
unsigned int e = 1;
e < dim; ++
e)
4848 values[comp * nqp] += val_in[e] * jac[e][comp];
4849 values[comp * nqp] *= fac;
4854 for (
unsigned int comp = 0; comp < n_components; ++comp)
4855 values[comp * nqp] = val_in[comp] * JxW;
4865 typename VectorizedArrayType>
4866template <
int,
typename>
4870 const unsigned int q_point)
4872 submit_value(val_in[0], q_point);
4881 typename VectorizedArrayType>
4884 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point)
4890 Assert(this->J_value !=
nullptr,
4892 "update_gradients"));
4893 Assert(this->jacobian !=
nullptr,
4895 "update_gradients"));
4897 this->gradients_quad_submitted =
true;
4900 if constexpr (dim > 1 && n_components == dim)
4902 if (this->data->element_type ==
4911 Assert(this->J_value !=
nullptr,
4913 "update_gradients"));
4914 Assert(this->jacobian !=
nullptr,
4916 "update_gradients"));
4918 this->gradients_quad_submitted =
true;
4921 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
4922 VectorizedArrayType *
values =
4923 this->values_from_gradients_quad + q_point;
4924 const std::size_t nqp = this->n_quadrature_points;
4925 const std::size_t nqp_d = nqp * dim;
4935 const VectorizedArrayType weight =
4936 this->quadrature_weights[q_point];
4937 for (
unsigned int d = 0;
d < dim; ++
d)
4938 for (
unsigned int comp = 0; comp < n_components; ++comp)
4939 gradients[comp * nqp_d + d] = grad_in[comp][d] *
4941 (jac[comp][comp] * weight);
4953 const VectorizedArrayType fac =
4955 this->quadrature_weights[q_point] :
4956 this->J_value[0] * this->quadrature_weights[q_point] *
4957 ((dim == 2 && this->get_face_no() < 2) ?
4962 VectorizedArrayType tmp[dim][dim];
4963 for (
unsigned int d = 0;
d < dim; ++
d)
4964 for (
unsigned int e = 0;
e < dim; ++
e)
4966 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in[
e][0];
4967 for (
unsigned int f = 1; f < dim; ++f)
4968 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
4970 for (
unsigned int comp = 0; comp < n_components; ++comp)
4971 for (
unsigned int d = 0;
d < dim; ++
d)
4973 VectorizedArrayType res = jac[0][comp] * tmp[
d][0];
4974 for (
unsigned int f = 1; f < dim; ++f)
4975 res += jac[f][comp] * tmp[d][f];
4984 const auto jac_grad =
4985 this->jacobian_gradients_non_inverse[q_point];
4987 this->jacobian[q_point];
4991 const VectorizedArrayType fac =
4992 (!is_face) ? this->quadrature_weights[q_point] :
4993 this->J_value[q_point] *
4994 ((dim == 2 && this->get_face_no() < 2) ?
5003 VectorizedArrayType tmp3[dim], tmp4[dim];
5004 for (
unsigned int d = 0;
d < dim; ++
d)
5006 tmp3[
d] = inv_t_jac[0][
d] * jac_grad[
d][0];
5007 for (
unsigned int e = 1;
e < dim; ++
e)
5008 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
5010 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5011 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5012 for (
unsigned int d = 0;
d < dim; ++
d)
5014 tmp3[f] += inv_t_jac[
d][
e] * jac_grad[k][
d];
5015 tmp3[
e] += inv_t_jac[
d][f] * jac_grad[k][
d];
5017 for (
unsigned int d = 0;
d < dim; ++
d)
5019 tmp4[
d] = tmp3[0] * inv_t_jac[
d][0];
5020 for (
unsigned int e = 1;
e < dim; ++
e)
5021 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
5027 VectorizedArrayType tmp[dim][dim];
5030 for (
unsigned int d = 0;
d < dim; ++
d)
5031 for (
unsigned int e = 0;
e < dim; ++
e)
5033 tmp[
d][
e] = inv_t_jac[0][
d] * grad_in_scaled[
e][0];
5034 for (
unsigned int f = 1; f < dim; ++f)
5035 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
5038 for (
unsigned int d = 0;
d < dim; ++
d)
5039 for (
unsigned int e = 0;
e < dim; ++
e)
5041 VectorizedArrayType res = t_jac[
d][0] * tmp[
e][0];
5042 for (
unsigned int f = 1; f < dim; ++f)
5043 res += t_jac[d][f] * tmp[e][f];
5050 VectorizedArrayType
value[dim];
5051 for (
unsigned int d = 0;
d < dim; ++
d)
5053 value[
d] = tmp[
d][0] * jac_grad[
d][0];
5054 for (
unsigned int e = 1;
e < dim; ++
e)
5055 value[d] += tmp[d][e] * jac_grad[d][e];
5057 for (
unsigned int e = 0, k = dim;
e < dim; ++
e)
5058 for (
unsigned int f = e + 1; f < dim; ++k, ++f)
5059 for (
unsigned int d = 0;
d < dim; ++
d)
5061 value[
e] += tmp[f][
d] * jac_grad[k][
d];
5062 value[f] += tmp[
e][
d] * jac_grad[k][
d];
5067 for (
unsigned int d = 0;
d < dim; ++
d)
5069 VectorizedArrayType tmp2 = grad_in_scaled[
d][0] * tmp4[0];
5070 for (
unsigned int e = 1;
e < dim; ++
e)
5071 tmp2 += grad_in_scaled[d][e] * tmp4[e];
5072 for (
unsigned int e = 0;
e < dim; ++
e)
5073 value[e] -= t_jac[e][d] * tmp2;
5076 for (
unsigned int d = 0;
d < dim; ++
d)
5077 values[d * nqp] =
value[d];
5083 const std::size_t nqp_d = this->n_quadrature_points * dim;
5084 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5088 const VectorizedArrayType JxW =
5089 this->J_value[0] * this->quadrature_weights[q_point];
5095 std::array<VectorizedArrayType, dim> jac;
5096 for (
unsigned int d = 0;
d < dim; ++
d)
5097 jac[d] = this->jacobian[0][d][d];
5099 for (
unsigned int d = 0;
d < dim; ++
d)
5101 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5102 if constexpr (n_components == 1)
5103 gradients[d] = grad_in[d] * factor;
5105 for (
unsigned int comp = 0; comp < n_components; ++comp)
5106 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5113 this->jacobian[q_point] :
5115 const VectorizedArrayType JxW =
5117 this->J_value[q_point] :
5118 this->J_value[0] * this->quadrature_weights[q_point];
5119 if constexpr (n_components == 1)
5120 for (
unsigned int d = 0;
d < dim; ++
d)
5122 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5123 for (
unsigned int e = 1;
e < dim; ++
e)
5124 new_val += (jac[e][d] * grad_in[e]);
5128 for (
unsigned int comp = 0; comp < n_components; ++comp)
5129 for (
unsigned int d = 0;
d < dim; ++
d)
5131 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5132 for (
unsigned int e = 1;
e < dim; ++
e)
5133 new_val += (jac[e][d] * grad_in[comp][e]);
5145 typename VectorizedArrayType>
5146template <
int,
typename>
5150 const unsigned int q_point)
5152 submit_gradient(grad_in[0], q_point);
5161 typename VectorizedArrayType>
5167 Assert(this->normal_x_jacobian !=
nullptr,
5169 "update_gradients"));
5171 this->gradients_quad_submitted =
true;
5174 const std::size_t nqp_d = this->n_quadrature_points * dim;
5175 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5179 const VectorizedArrayType JxW_jac = this->J_value[0] *
5180 this->quadrature_weights[q_point] *
5181 this->normal_x_jacobian[0][dim - 1];
5182 for (
unsigned int comp = 0; comp < n_components; ++comp)
5184 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5185 gradients[comp * nqp_d + d] = VectorizedArrayType();
5186 if constexpr (n_components == 1)
5187 gradients[dim - 1] = grad_in * JxW_jac;
5189 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5194 const unsigned int index =
5197 this->normal_x_jacobian[
index];
5198 const VectorizedArrayType JxW =
5200 this->J_value[
index] * this->quadrature_weights[q_point] :
5201 this->J_value[
index];
5202 for (
unsigned int comp = 0; comp < n_components; ++comp)
5203 for (
unsigned int d = 0;
d < dim; ++
d)
5204 if constexpr (n_components == 1)
5207 gradients[comp * nqp_d +
d] = (grad_in[comp] * JxW) * jac[d];
5217 typename VectorizedArrayType>
5220 submit_hessian(
const hessian_type hessian_in,
const unsigned int q_point)
5226 Assert(this->J_value !=
nullptr,
5228 "update_hessians"));
5229 Assert(this->jacobian !=
nullptr,
5231 "update_hessians"));
5233 this->hessians_quad_submitted =
true;
5237 const std::size_t nqp = this->n_quadrature_points;
5238 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5241 const VectorizedArrayType JxW =
5242 this->J_value[0] * this->quadrature_weights[q_point];
5245 for (
unsigned int d = 0;
d < dim; ++
d)
5247 const auto jac_d = this->jacobian[0][
d][
d];
5248 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5249 for (
unsigned int comp = 0; comp < n_components; ++comp)
5250 if constexpr (n_components == 1)
5251 this->hessians_quad[
d * nqp + q_point] =
5252 hessian_in[
d][
d] * factor;
5254 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5255 hessian_in[comp][d][d] * factor;
5259 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5260 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5262 const auto jac_d = this->jacobian[0][
d][
d];
5263 const auto jac_e = this->jacobian[0][
e][
e];
5264 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5265 for (
unsigned int comp = 0; comp < n_components; ++comp)
5266 if constexpr (n_components == 1)
5267 this->hessians_quad[off_dia * nqp + q_point] =
5268 (hessian_in[
d][
e] + hessian_in[
e][
d]) * factor;
5270 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5271 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5278 const VectorizedArrayType JxW =
5279 this->J_value[0] * this->quadrature_weights[q_point];
5280 for (
unsigned int comp = 0; comp < n_components; ++comp)
5283 if constexpr (n_components == 1)
5284 hessian_c = hessian_in;
5286 hessian_c = hessian_in[comp];
5289 VectorizedArrayType tmp[dim][dim];
5290 for (
unsigned int i = 0; i < dim; ++i)
5291 for (
unsigned int j = 0; j < dim; ++j)
5293 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5294 for (
unsigned int k = 1; k < dim; ++k)
5295 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5299 VectorizedArrayType tmp2[dim][dim];
5300 for (
unsigned int i = 0; i < dim; ++i)
5301 for (
unsigned int j = 0; j < dim; ++j)
5303 tmp2[i][j] = jac[0][i] * tmp[0][j];
5304 for (
unsigned int k = 1; k < dim; ++k)
5305 tmp2[i][j] += jac[k][i] * tmp[k][j];
5309 for (
unsigned int d = 0;
d < dim; ++
d)
5310 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5314 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5315 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5316 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5317 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5323 const VectorizedArrayType JxW = this->J_value[q_point];
5324 const auto &jac_grad = this->jacobian_gradients[q_point];
5325 for (
unsigned int comp = 0; comp < n_components; ++comp)
5328 if constexpr (n_components == 1)
5329 hessian_c = hessian_in;
5331 hessian_c = hessian_in[comp];
5334 VectorizedArrayType tmp[dim][dim];
5335 for (
unsigned int i = 0; i < dim; ++i)
5336 for (
unsigned int j = 0; j < dim; ++j)
5338 tmp[i][j] = hessian_c[i][0] * jac[0][j];
5339 for (
unsigned int k = 1; k < dim; ++k)
5340 tmp[i][j] += hessian_c[i][k] * jac[k][j];
5344 VectorizedArrayType tmp2[dim][dim];
5345 for (
unsigned int i = 0; i < dim; ++i)
5346 for (
unsigned int j = 0; j < dim; ++j)
5348 tmp2[i][j] = jac[0][i] * tmp[0][j];
5349 for (
unsigned int k = 1; k < dim; ++k)
5350 tmp2[i][j] += jac[k][i] * tmp[k][j];
5354 for (
unsigned int d = 0;
d < dim; ++
d)
5355 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5359 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5360 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5361 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5362 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5365 for (
unsigned int d = 0;
d < dim; ++
d)
5367 VectorizedArrayType
sum = 0;
5368 for (
unsigned int e = 0;
e < dim; ++
e)
5369 sum += hessian_c[e][e] * jac_grad[e][d];
5370 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5371 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5373 (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][
d];
5374 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5387 typename VectorizedArrayType>
5392 VectorizedArrayType>::value_type
5398 Assert(this->values_quad_submitted ==
true,
5403 const std::size_t nqp = this->n_quadrature_points;
5404 for (
unsigned int q = 0; q < nqp; ++q)
5405 for (
unsigned int comp = 0; comp < n_components; ++comp)
5406 return_value[comp] += this->values_quad[comp * nqp + q];
5407 if constexpr (n_components == 1)
5408 return return_value[0];
5410 return return_value;
5419 typename VectorizedArrayType>
5420template <
int,
typename>
5426 Assert(this->gradients_quad_initialized ==
true,
5430 Assert(this->jacobian !=
nullptr,
5432 "update_gradients"));
5434 VectorizedArrayType divergence;
5435 const std::size_t nqp = this->n_quadrature_points;
5438 this->data->element_type ==
5441 VectorizedArrayType inv_det =
5444 this->jacobian[0][0][0] *
5445 ((dim == 2) ? this->jacobian[0][1][1] :
5446 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5454 if (is_face && dim == 2 && this->get_face_no() < 2)
5458 divergence = this->gradients_quad[q_point * dim];
5459 for (
unsigned int d = 1;
d < dim; ++
d)
5460 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d];
5461 divergence *= inv_det;
5470 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
5471 for (
unsigned int d = 1;
d < dim; ++
d)
5472 divergence += this->gradients_quad[(d * nqp + q_point) * dim +
d] *
5473 this->jacobian[0][
d][
d];
5480 this->jacobian[q_point] :
5482 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
5483 for (
unsigned int e = 1;
e < dim; ++
e)
5484 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
5485 for (
unsigned int d = 1;
d < dim; ++
d)
5486 for (
unsigned int e = 0;
e < dim; ++
e)
5488 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim +
e];
5500 typename VectorizedArrayType>
5501template <
int,
typename>
5507 const auto grad = get_gradient(q_point);
5508 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5509 VectorizedArrayType half = Number(0.5);
5510 for (
unsigned int d = 0;
d < dim; ++
d)
5511 symmetrized[d] = grad[d][d];
5517 symmetrized[2] = grad[0][1] + grad[1][0];
5518 symmetrized[2] *= half;
5521 symmetrized[3] = grad[0][1] + grad[1][0];
5522 symmetrized[3] *= half;
5523 symmetrized[4] = grad[0][2] + grad[2][0];
5524 symmetrized[4] *= half;
5525 symmetrized[5] = grad[1][2] + grad[2][1];
5526 symmetrized[5] *= half;
5540 typename VectorizedArrayType>
5541template <
int,
typename>
5543 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5545 get_curl(const unsigned
int q_point) const
5549 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5553 curl[0] = grad[1][0] - grad[0][1];
5556 curl[0] = grad[2][1] - grad[1][2];
5557 curl[1] = grad[0][2] - grad[2][0];
5558 curl[2] = grad[1][0] - grad[0][1];
5572 typename VectorizedArrayType>
5573template <
int,
typename>
5577 const unsigned int q_point)
5583 Assert(this->J_value !=
nullptr,
5585 "update_gradients"));
5586 Assert(this->jacobian !=
nullptr,
5588 "update_gradients"));
5590 this->gradients_quad_submitted =
true;
5593 const std::size_t nqp_d = this->n_quadrature_points * dim;
5594 VectorizedArrayType *
gradients = this->gradients_quad + q_point * dim;
5596 if (this->data->element_type ==
5603 const VectorizedArrayType fac =
5605 this->quadrature_weights[q_point] * div_in :
5607 this->J_value[q_point] :
5608 this->J_value[0] * this->quadrature_weights[q_point]) *
5611 this->jacobian[this->cell_type >
5615 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
5617 for (
unsigned int d = 0;
d < dim; ++
d)
5619 for (
unsigned int e = 0;
e < dim; ++
e)
5620 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
5622 this->divergence_is_requested =
true;
5629 const VectorizedArrayType fac =
5630 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
5631 for (
unsigned int d = 0;
d < dim; ++
d)
5633 const VectorizedArrayType jac_dd = this->jacobian[0][
d][
d];
5634 for (
unsigned int e = 0;
e < dim; ++
e)
5635 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
5642 this->jacobian[q_point] :
5644 const VectorizedArrayType fac =
5646 this->J_value[q_point] :
5647 this->J_value[0] * this->quadrature_weights[q_point]) *
5649 for (
unsigned int d = 0;
d < dim; ++
d)
5651 for (
unsigned int e = 0;
e < dim; ++
e)
5652 gradients[d * nqp_d + e] = jac[d][e] * fac;
5664 typename VectorizedArrayType>
5665template <
int,
typename>
5670 const unsigned int q_point)
5673 this->data->element_type !=
5684 Assert(this->J_value !=
nullptr,
5686 "update_gradients"));
5687 Assert(this->jacobian !=
nullptr,
5689 "update_gradients"));
5691 this->gradients_quad_submitted =
true;
5694 const std::size_t nqp_d = this->n_quadrature_points * dim;
5695 VectorizedArrayType *
gradients = this->gradients_quad + dim * q_point;
5698 const VectorizedArrayType JxW =
5699 this->J_value[0] * this->quadrature_weights[q_point];
5701 for (
unsigned int d = 0;
d < dim; ++
d)
5702 gradients[d * nqp_d + d] =
5704 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5705 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5707 const VectorizedArrayType
value =
5716 const VectorizedArrayType JxW =
5718 this->J_value[q_point] :
5719 this->J_value[0] * this->quadrature_weights[q_point];
5722 this->jacobian[q_point] :
5724 VectorizedArrayType weighted[dim][dim];
5725 for (
unsigned int i = 0; i < dim; ++i)
5727 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
5728 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
5730 const VectorizedArrayType
value =
5732 weighted[i][j] =
value;
5733 weighted[j][i] =
value;
5735 for (
unsigned int comp = 0; comp < dim; ++comp)
5736 for (
unsigned int d = 0;
d < dim; ++
d)
5738 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
5739 for (
unsigned int e = 1;
e < dim; ++
e)
5740 new_val += jac[e][d] * weighted[comp][e];
5752 typename VectorizedArrayType>
5753template <
int,
typename>
5757 const unsigned int q_point)
5763 grad[1][0] = curl[0];
5764 grad[0][1] = -curl[0];
5767 grad[2][1] = curl[0];
5768 grad[1][2] = -curl[0];
5769 grad[0][2] = curl[1];
5770 grad[2][0] = -curl[1];
5771 grad[1][0] = curl[2];
5772 grad[0][1] = -curl[2];
5777 submit_gradient(grad, q_point);
5790 typename VectorizedArrayType>
5796 VectorizedArrayType>::
5798 const unsigned int fe_no,
5799 const unsigned int quad_no,
5800 const unsigned int first_selected_component,
5801 const unsigned int active_fe_index,
5802 const unsigned int active_quad_index)
5803 : BaseClass(matrix_free,
5805 first_selected_component,
5812 numbers::invalid_unsigned_int )
5813 , dofs_per_component(this->data->dofs_per_component_on_cell)
5814 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5815 , n_q_points(this->data->n_q_points)
5817 check_template_arguments(fe_no, 0);
5827 typename VectorizedArrayType>
5833 VectorizedArrayType>::
5835 const std::pair<unsigned int, unsigned int> &range,
5836 const unsigned int dof_no,
5837 const unsigned int quad_no,
5838 const unsigned int first_selected_component)
5842 first_selected_component,
5843 matrix_free.get_cell_active_fe_index(range))
5853 typename VectorizedArrayType>
5859 VectorizedArrayType>::
5864 const unsigned int first_selected_component)
5869 first_selected_component,
5871 , dofs_per_component(this->data->dofs_per_component_on_cell)
5872 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5873 , n_q_points(this->data->n_q_points)
5885 typename VectorizedArrayType>
5891 VectorizedArrayType>::
5895 const unsigned int first_selected_component)
5900 first_selected_component,
5902 , dofs_per_component(this->data->dofs_per_component_on_cell)
5903 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5904 , n_q_points(this->data->n_q_points)
5916 typename VectorizedArrayType>
5922 VectorizedArrayType>::
5925 const unsigned int first_selected_component)
5926 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
5928 other.mapped_geometry->get_quadrature(),
5929 other.mapped_geometry->get_fe_values().get_update_flags(),
5930 first_selected_component,
5932 , dofs_per_component(this->data->dofs_per_component_on_cell)
5933 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5934 , n_q_points(this->data->n_q_points)
5946 typename VectorizedArrayType>
5955 , dofs_per_component(this->data->dofs_per_component_on_cell)
5956 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
5957 , n_q_points(this->data->n_q_points)
5969 typename VectorizedArrayType>
5975 VectorizedArrayType> &
5981 VectorizedArrayType>::operator=(
const FEEvaluation &other)
5983 BaseClass::operator=(other);
5995 typename VectorizedArrayType>
6002 VectorizedArrayType>::
6003 check_template_arguments(
const unsigned int dof_no,
6004 const unsigned int first_selected_component)
6007 (void)first_selected_component;
6010 this->data->dofs_per_component_on_cell > 0,
6012 "There is nothing useful you can do with an FEEvaluation object with "
6013 "FE_Nothing, i.e., without DoFs! If you have passed to "
6014 "MatrixFree::reinit() a collection of finite elements also containing "
6015 "FE_Nothing, please check - before creating FEEvaluation - the category "
6016 "of the current range by calling either "
6017 "MatrixFree::get_cell_range_category(range) or "
6018 "MatrixFree::get_face_range_category(range). The returned category "
6019 "is the index of the active FE, which you can use to exclude "
6026 static_cast<unsigned int>(fe_degree) !=
6027 this->data->data.front().fe_degree) ||
6028 n_q_points != this->n_quadrature_points)
6030 std::string message =
6031 "-------------------------------------------------------\n";
6032 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6033 message +=
" Called --> FEEvaluation<dim,";
6037 message +=
",Number>(data";
6053 if (
static_cast<unsigned int>(fe_degree) ==
6054 this->data->data.front().fe_degree)
6056 proposed_dof_comp = dof_no;
6057 proposed_fe_comp = first_selected_component;
6060 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
6062 for (
unsigned int nf = 0;
6065 if (this->matrix_free
6066 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6068 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6070 proposed_dof_comp = no;
6071 proposed_fe_comp = nf;
6075 this->mapping_data->descriptor[this->active_quad_index]
6077 proposed_quad_comp = this->quad_no;
6079 for (
unsigned int no = 0;
6084 .descriptor[this->active_quad_index]
6085 .n_q_points == n_q_points)
6087 proposed_quad_comp = no;
6094 if (proposed_dof_comp != first_selected_component)
6095 message +=
"Wrong vector component selection:\n";
6097 message +=
"Wrong quadrature formula selection:\n";
6098 message +=
" Did you mean FEEvaluation<dim,";
6102 message +=
",Number>(data";
6111 std::string correct_pos;
6112 if (proposed_dof_comp != dof_no)
6113 correct_pos =
" ^ ";
6116 if (proposed_quad_comp != this->quad_no)
6117 correct_pos +=
" ^ ";
6120 if (proposed_fe_comp != first_selected_component)
6121 correct_pos +=
" ^\n";
6123 correct_pos +=
" \n";
6129 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
6130 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6131 message +=
"Wrong template arguments:\n";
6132 message +=
" Did you mean FEEvaluation<dim,";
6137 message +=
",Number>(data";
6145 std::string correct_pos;
6146 if (this->data->data.front().fe_degree !=
6147 static_cast<unsigned int>(fe_degree))
6151 if (proposed_n_q_points_1d != n_q_points_1d)
6152 correct_pos +=
" ^\n";
6154 correct_pos +=
" \n";
6155 message +=
" " + correct_pos;
6157 Assert(
static_cast<unsigned int>(fe_degree) ==
6158 this->data->data.front().fe_degree &&
6159 n_q_points == this->n_quadrature_points,
6165 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6176 typename VectorizedArrayType>
6185 Assert(this->matrix_free !=
nullptr,
6186 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6187 " Integer indexing is not possible."));
6195 const unsigned int offsets =
6196 this->mapping_data->data_index_offsets[
cell_index];
6197 this->jacobian = &this->mapping_data->jacobians[0][offsets];
6198 this->J_value = &this->mapping_data->JxW_values[offsets];
6199 if (!this->mapping_data->jacobian_gradients[0].empty())
6201 this->jacobian_gradients =
6202 this->mapping_data->jacobian_gradients[0].data() + offsets;
6203 this->jacobian_gradients_non_inverse =
6204 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
6210 for (
unsigned int i = 0; i < n_lanes; ++i)
6211 this->cell_ids[i] =
cell_index * n_lanes + i;
6218 this->cell_ids[i] =
cell_index * n_lanes + i;
6219 for (; i < n_lanes; ++i)
6223 if (this->mapping_data->quadrature_points.empty() ==
false)
6225 &this->mapping_data->quadrature_points
6226 [this->mapping_data->quadrature_point_offsets[this->cell]];
6229 this->is_reinitialized =
true;
6230 this->dof_values_initialized =
false;
6231 this->values_quad_initialized =
false;
6232 this->gradients_quad_initialized =
false;
6233 this->hessians_quad_initialized =
false;
6244 typename VectorizedArrayType>
6251 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6258 this->cell_ids = cell_ids;
6263 for (
unsigned int v = 0; v < n_lanes; ++v)
6277 if (this->mapped_geometry ==
nullptr)
6278 this->mapped_geometry =
6279 std::make_shared<internal::MatrixFreeFunctions::
6280 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6282 auto &mapping_storage = this->mapped_geometry->get_data_storage();
6284 auto &this_jacobian_data = mapping_storage.jacobians[0];
6285 auto &this_J_value_data = mapping_storage.JxW_values;
6286 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
6287 auto &this_jacobian_gradients_non_inverse_data =
6288 mapping_storage.jacobian_gradients_non_inverse[0];
6289 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
6293 if (this_jacobian_data.size() != 2)
6294 this_jacobian_data.resize_fast(2);
6296 if (this_J_value_data.size() != 1)
6297 this_J_value_data.resize_fast(1);
6299 const auto &update_flags_cells =
6303 this_jacobian_gradients_data.size() != 1)
6305 this_jacobian_gradients_data.resize_fast(1);
6306 this_jacobian_gradients_non_inverse_data.resize_fast(1);
6310 this_quadrature_points_data.size() != 1)
6311 this_quadrature_points_data.resize_fast(1);
6315 if (this_jacobian_data.size() != this->n_quadrature_points)
6316 this_jacobian_data.resize_fast(this->n_quadrature_points);
6318 if (this_J_value_data.size() != this->n_quadrature_points)
6319 this_J_value_data.resize_fast(this->n_quadrature_points);
6321 const auto &update_flags_cells =
6325 this_jacobian_gradients_data.size() != this->n_quadrature_points)
6327 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
6328 this_jacobian_gradients_non_inverse_data.resize_fast(
6329 this->n_quadrature_points);
6333 this_quadrature_points_data.size() != this->n_quadrature_points)
6334 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
6338 this->jacobian = this_jacobian_data.data();
6339 this->J_value = this_J_value_data.data();
6340 this->jacobian_gradients = this_jacobian_gradients_data.data();
6341 this->jacobian_gradients_non_inverse =
6342 this_jacobian_gradients_non_inverse_data.data();
6346 for (
unsigned int v = 0; v < n_lanes; ++v)
6353 const unsigned int cell_batch_index =
cell_index / n_lanes;
6354 const unsigned int offsets =
6355 this->mapping_data->data_index_offsets[cell_batch_index];
6356 const unsigned int lane =
cell_index % n_lanes;
6358 if (this->cell_type <=
6362 for (
unsigned int q = 0; q < 2; ++q)
6363 for (
unsigned int i = 0; i < dim; ++i)
6364 for (
unsigned int j = 0; j < dim; ++j)
6365 this_jacobian_data[q][i][j][v] =
6366 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
6368 const unsigned int q = 0;
6370 this_J_value_data[q][v] =
6371 this->mapping_data->JxW_values[offsets + q][lane];
6373 const auto &update_flags_cells =
6378 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6379 for (
unsigned int j = 0; j < dim; ++j)
6380 this_jacobian_gradients_data[q][i][j][v] =
6382 ->jacobian_gradients[0][offsets + q][i][j][lane];
6384 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6385 for (
unsigned int j = 0; j < dim; ++j)
6386 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6388 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
6393 for (
unsigned int i = 0; i < dim; ++i)
6394 this_quadrature_points_data[q][i][v] =
6395 this->mapping_data->quadrature_points
6397 ->quadrature_point_offsets[cell_batch_index] +
6403 const auto cell_type =
6407 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
6409 const unsigned int q_src =
6415 this_J_value_data[q][v] =
6416 this->mapping_data->JxW_values[offsets + q_src][lane];
6418 for (
unsigned int i = 0; i < dim; ++i)
6419 for (
unsigned int j = 0; j < dim; ++j)
6420 this_jacobian_data[q][i][j][v] =
6422 ->jacobians[0][offsets + q_src][i][j][lane];
6424 const auto &update_flags_cells =
6429 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6430 for (
unsigned int j = 0; j < dim; ++j)
6431 this_jacobian_gradients_data[q][i][j][v] =
6433 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
6435 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6436 for (
unsigned int j = 0; j < dim; ++j)
6437 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
6439 ->jacobian_gradients_non_inverse[0][offsets + q_src]
6452 this->mapping_data->quadrature_points
6454 ->quadrature_point_offsets[cell_batch_index] +
6458 this->mapping_data->jacobians[0][offsets + 1];
6460 for (
unsigned int d = 0;
d < dim; ++
d)
6463 static_cast<Number
>(
6464 this->descriptor->quadrature.point(q)[
d]);
6466 for (
unsigned int d = 0;
d < dim; ++
d)
6467 for (
unsigned int e = 0;
e < dim; ++
e)
6470 static_cast<Number
>(
6471 this->descriptor->quadrature.point(q)[
e]);
6473 for (
unsigned int i = 0; i < dim; ++i)
6474 this_quadrature_points_data[q][i][v] = point[i][lane];
6479 for (
unsigned int i = 0; i < dim; ++i)
6480 this_quadrature_points_data[q][i][v] =
6481 this->mapping_data->quadrature_points
6483 ->quadrature_point_offsets[cell_batch_index] +
6492 this->is_reinitialized =
true;
6493 this->dof_values_initialized =
false;
6494 this->values_quad_initialized =
false;
6495 this->gradients_quad_initialized =
false;
6496 this->hessians_quad_initialized =
false;
6507 typename VectorizedArrayType>
6508template <
bool level_dof_access>
6515 VectorizedArrayType>
::
6518 Assert(this->matrix_free ==
nullptr,
6519 ExcMessage(
"Cannot use initialization from cell iterator if "
6520 "initialized from MatrixFree object. Use variant for "
6521 "on the fly computation with arguments as for FEValues "
6524 this->mapped_geometry->reinit(
6526 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6527 if (level_dof_access)
6528 cell->get_mg_dof_indices(this->local_dof_indices);
6530 cell->get_dof_indices(this->local_dof_indices);
6533 this->is_reinitialized =
true;
6544 typename VectorizedArrayType>
6551 VectorizedArrayType>
::
6554 Assert(this->matrix_free == 0,
6555 ExcMessage(
"Cannot use initialization from cell iterator if "
6556 "initialized from MatrixFree object. Use variant for "
6557 "on the fly computation with arguments as for FEValues "
6560 this->mapped_geometry->reinit(cell);
6563 this->is_reinitialized =
true;
6574 typename VectorizedArrayType>
6581 VectorizedArrayType>::
6585 Assert(this->dof_values_initialized ==
true,
6588 evaluate(this->values_dofs, evaluation_flags);
6598 typename VectorizedArrayType>
6605 VectorizedArrayType>::
6606 evaluate(
const VectorizedArrayType *values_array,
6609 const bool hessians_on_general_cells =
6613 if (hessians_on_general_cells)
6616 if (this->data->element_type ==
6622 if constexpr (fe_degree > -1)
6625 template run<fe_degree, n_q_points_1d>(n_components,
6626 evaluation_flag_actual,
6634 evaluation_flag_actual,
6635 const_cast<VectorizedArrayType *
>(values_array),
6640 this->values_quad_initialized =
6642 this->gradients_quad_initialized =
6644 this->hessians_quad_initialized =
6655 template <
typename Number,
6656 typename VectorizedArrayType,
6657 typename VectorType,
6658 typename EvaluatorType,
6659 std::enable_if_t<internal::has_begin<VectorType> &&
6661 VectorType> * =
nullptr>
6662 VectorizedArrayType *
6663 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
6669 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6670 const auto &dof_info = fe_eval.get_dof_info();
6677 if (std::is_same_v<typename VectorType::value_type, Number> &&
6681 interleaved_contiguous &&
6682 reinterpret_cast<
std::size_t>(
6684 dof_info.dof_indices_contiguous
6685 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6686 [cell * VectorizedArrayType::size()]) %
6687 sizeof(VectorizedArrayType) ==
6690 return reinterpret_cast<VectorizedArrayType *
>(
6694 [cell * VectorizedArrayType::size()] +
6696 [fe_eval.get_active_fe_index()]
6697 [fe_eval.get_first_selected_component()] *
6698 VectorizedArrayType::size());
6707 template <
typename Number,
6708 typename VectorizedArrayType,
6709 typename VectorType,
6710 typename EvaluatorType,
6711 std::enable_if_t<!internal::has_begin<VectorType> ||
6713 VectorType> * =
nullptr>
6714 VectorizedArrayType *
6715 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
6728 typename VectorizedArrayType>
6729template <
typename VectorType>
6736 VectorizedArrayType>::
6737 gather_evaluate(
const VectorType &input_vector,
6740 const VectorizedArrayType *src_ptr =
6741 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
6742 *
this, input_vector);
6743 if (src_ptr !=
nullptr)
6744 evaluate(src_ptr, evaluation_flag);
6747 this->read_dof_values(input_vector);
6748 evaluate(this->begin_dof_values(), evaluation_flag);
6759 typename VectorizedArrayType>
6766 VectorizedArrayType>::
6769 integrate(integration_flag, this->values_dofs);
6772 this->dof_values_initialized =
true;
6783 typename VectorizedArrayType>
6790 VectorizedArrayType>::
6792 VectorizedArrayType *values_array,
6793 const bool sum_into_values_array)
6797 Assert(this->values_quad_submitted ==
true,
6800 Assert(this->gradients_quad_submitted ==
true,
6803 Assert(this->hessians_quad_submitted ==
true,
6806 Assert(this->matrix_free !=
nullptr ||
6807 this->mapped_geometry->is_initialized(),
6813 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
6814 "EvaluationFlags::hessians are supported."));
6820 unsigned int size = n_components * dim * n_q_points;
6823 for (
unsigned int i = 0; i < size; ++i)
6824 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
6828 for (
unsigned int i = 0; i < size; ++i)
6829 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
6834 if (n_components == dim &&
6835 this->data->element_type ==
6839 this->divergence_is_requested ==
false)
6841 unsigned int size = n_components * n_q_points;
6844 for (
unsigned int i = 0; i < size; ++i)
6845 this->values_quad[i] += this->values_from_gradients_quad[i];
6849 for (
unsigned int i = 0; i < size; ++i)
6850 this->values_quad[i] = this->values_from_gradients_quad[i];
6855 if constexpr (fe_degree > -1)
6858 template run<fe_degree, n_q_points_1d>(n_components,
6859 integration_flag_actual,
6862 sum_into_values_array);
6868 integration_flag_actual,
6871 sum_into_values_array);
6875 this->dof_values_initialized =
true;
6886 typename VectorizedArrayType>
6887template <
typename VectorType>
6894 VectorizedArrayType>::
6896 VectorType &destination)
6898 VectorizedArrayType *dst_ptr =
6899 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
6900 *
this, destination);
6901 if (dst_ptr !=
nullptr)
6902 integrate(integration_flag, dst_ptr,
true);
6905 integrate(integration_flag, this->begin_dof_values());
6906 this->distribute_local_to_global(destination);
6917 typename VectorizedArrayType>
6924 VectorizedArrayType>::dof_indices()
const
6926 return {0U, dofs_per_cell};
6940 typename VectorizedArrayType>
6946 VectorizedArrayType>::
6949 const bool is_interior_face,
6950 const unsigned int dof_no,
6951 const unsigned int quad_no,
6952 const unsigned int first_selected_component,
6953 const unsigned int active_fe_index,
6954 const unsigned int active_quad_index,
6955 const unsigned int face_type)
6956 : BaseClass(matrix_free,
6958 first_selected_component,
6966 , dofs_per_component(this->data->dofs_per_component_on_cell)
6967 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6968 , n_q_points(this->n_quadrature_points)
6978 typename VectorizedArrayType>
6984 VectorizedArrayType>::
6987 const std::pair<unsigned int, unsigned int> &range,
6988 const bool is_interior_face,
6989 const unsigned int dof_no,
6990 const unsigned int quad_no,
6991 const unsigned int first_selected_component)
6996 first_selected_component,
6997 matrix_free.get_face_active_fe_index(range,
6999 numbers::invalid_unsigned_int,
7000 matrix_free.get_face_info(range.
first).face_type)
7010 typename VectorizedArrayType>
7017 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7019 Assert(this->mapped_geometry ==
nullptr,
7020 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7021 " Integer indexing is not possible"));
7022 if (this->mapped_geometry !=
nullptr)
7025 this->cell = face_index;
7026 this->dof_access_index =
7027 this->is_interior_face() ?
7035 this->matrix_free->get_task_info().boundary_partition_data.back())
7036 Assert(this->is_interior_face(),
7038 "Boundary faces do not have a neighbor. When looping over "
7039 "boundary faces use FEFaceEvaluation with the parameter "
7040 "is_interior_face set to true. "));
7042 this->reinit_face(this->matrix_free->
get_face_info(face_index));
7047 this->face_ids[i] = face_index * n_lanes + i;
7048 for (; i < n_lanes; ++i)
7051 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
7052 const unsigned int offsets =
7053 this->mapping_data->data_index_offsets[face_index];
7054 this->J_value = &this->mapping_data->JxW_values[offsets];
7055 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7057 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
7058 this->normal_x_jacobian =
7060 ->normals_times_jacobians[!this->is_interior_face()][offsets];
7061 this->jacobian_gradients =
7062 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7064 this->jacobian_gradients_non_inverse =
7066 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7070 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7073 this->mapping_data->quadrature_point_offsets.size());
7075 this->mapping_data->quadrature_points.data() +
7076 this->mapping_data->quadrature_point_offsets[this->cell];
7080 this->is_reinitialized =
true;
7081 this->dof_values_initialized =
false;
7082 this->values_quad_initialized =
false;
7083 this->gradients_quad_initialized =
false;
7084 this->hessians_quad_initialized =
false;
7095 typename VectorizedArrayType>
7103 const unsigned int face_number)
7109 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7113 Assert(this->mapped_geometry ==
nullptr,
7114 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7115 " Integer indexing is not possible"));
7116 if (this->mapped_geometry !=
nullptr)
7121 .faces_by_cells_type[
cell_index][face_number];
7124 this->dof_access_index =
7127 if (this->is_interior_face() ==
false)
7132 for (
unsigned int i = 0; i < n_lanes; ++i)
7135 const unsigned int cell_this =
cell_index * n_lanes + i;
7137 unsigned int face_index =
7142 this->face_ids[i] = face_index;
7147 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7148 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7155 auto cell_m = faces.cells_interior[face_index % n_lanes];
7156 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7158 const bool face_identifies_as_interior = cell_m != cell_this;
7160 Assert(cell_m == cell_this || cell_p == cell_this,
7164 if (face_identifies_as_interior)
7166 this->cell_ids[i] = cell_m;
7167 this->face_numbers[i] = faces.interior_face_no;
7171 this->cell_ids[i] = cell_p;
7172 this->face_numbers[i] = faces.exterior_face_no;
7175 const bool orientation_interior_face = faces.face_orientation >= 8;
7176 unsigned int face_orientation = faces.face_orientation % 8;
7177 if (face_identifies_as_interior != orientation_interior_face)
7179 constexpr std::array<std::uint8_t, 8> table{
7180 {0, 1, 2, 3, 6, 5, 4, 7}};
7181 face_orientation = table[face_orientation];
7183 this->face_orientations[i] = face_orientation;
7188 this->face_orientations[0] = 0;
7189 this->face_numbers[0] = face_number;
7194 for (
unsigned int i = 0; i < n_lanes; ++i)
7195 this->cell_ids[i] =
cell_index * n_lanes + i;
7203 this->cell_ids[i] =
cell_index * n_lanes + i;
7204 for (; i < n_lanes; ++i)
7207 for (
unsigned int i = 0; i < n_lanes; ++i)
7214 const unsigned int offsets =
7216 .face_data_by_cells[this->quad_no]
7221 .face_data_by_cells[this->quad_no]
7222 .JxW_values.size());
7224 .face_data_by_cells[this->quad_no]
7225 .JxW_values[offsets];
7227 .face_data_by_cells[this->quad_no]
7228 .normal_vectors[offsets];
7230 .face_data_by_cells[this->quad_no]
7231 .jacobians[!this->is_interior_face()][offsets];
7232 this->normal_x_jacobian =
7234 .face_data_by_cells[this->quad_no]
7235 .normals_times_jacobians[!this->is_interior_face()][offsets];
7236 this->jacobian_gradients =
7237 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7239 this->jacobian_gradients_non_inverse =
7241 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7246 .face_data_by_cells[this->quad_no]
7247 .quadrature_point_offsets.empty() ==
false)
7249 const unsigned int index =
7253 .face_data_by_cells[this->quad_no]
7254 .quadrature_point_offsets.size());
7256 .face_data_by_cells[this->quad_no]
7257 .quadrature_points.data() +
7259 .face_data_by_cells[this->quad_no]
7260 .quadrature_point_offsets[
index];
7264 this->is_reinitialized =
true;
7265 this->dof_values_initialized =
false;
7266 this->values_quad_initialized =
false;
7267 this->gradients_quad_initialized =
false;
7268 this->hessians_quad_initialized =
false;
7279 typename VectorizedArrayType>
7286 VectorizedArrayType>::
7293 evaluate(this->values_dofs, evaluation_flag);
7303 typename VectorizedArrayType>
7310 VectorizedArrayType>::
7311 evaluate(
const VectorizedArrayType *values_array,
7314 Assert((evaluation_flag &
7317 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7318 "and EvaluationFlags::hessians are supported."));
7320 const bool hessians_on_general_cells =
7324 if (hessians_on_general_cells)
7327 if (this->data->element_type ==
7333 if constexpr (fe_degree > -1)
7335 template
run<fe_degree, n_q_points_1d>(n_components,
7336 evaluation_flag_actual,
7341 n_components, evaluation_flag_actual, values_array, *
this);
7344 this->values_quad_initialized =
7346 this->gradients_quad_initialized =
7348 this->hessians_quad_initialized =
7360 typename VectorizedArrayType>
7367 VectorizedArrayType>::
7374 project_to_face(this->values_dofs, evaluation_flag);
7384 typename VectorizedArrayType>
7391 VectorizedArrayType>::
7392 project_to_face(
const VectorizedArrayType *values_array,
7395 Assert((evaluation_flag &
7398 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7399 "and EvaluationFlags::hessians are supported."));
7401 const bool hessians_on_general_cells =
7405 if (hessians_on_general_cells)
7408 if (this->data->element_type ==
7414 if constexpr (fe_degree > -1)
7417 VectorizedArrayType>::template run<fe_degree>(n_components,
7418 evaluation_flag_actual,
7424 evaluation_flag_actual,
7438 typename VectorizedArrayType>
7445 VectorizedArrayType>::
7448 Assert((evaluation_flag &
7451 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7452 "and EvaluationFlags::hessians are supported."));
7454 const bool hessians_on_general_cells =
7458 if (hessians_on_general_cells)
7461 if (this->data->element_type ==
7467 if constexpr (fe_degree > -1)
7470 VectorizedArrayType>::template run<fe_degree>(n_components,
7471 evaluation_flag_actual,
7478 this->values_quad_initialized =
7480 this->gradients_quad_initialized =
7482 this->hessians_quad_initialized =
7494 typename VectorizedArrayType>
7501 VectorizedArrayType>::
7503 const bool sum_into_values)
7505 integrate(integration_flag, this->values_dofs, sum_into_values);
7508 this->dof_values_initialized =
true;
7519 typename VectorizedArrayType>
7526 VectorizedArrayType>::
7528 VectorizedArrayType *values_array,
7529 const bool sum_into_values)
7531 Assert((integration_flag &
7534 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7535 "and EvaluationFlags::hessians are supported."));
7541 unsigned int size = n_components * dim * n_q_points;
7544 for (
unsigned int i = 0; i < size; ++i)
7545 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7549 for (
unsigned int i = 0; i < size; ++i)
7550 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7555 if (this->data->element_type ==
7559 this->divergence_is_requested ==
false)
7561 unsigned int size = n_components * n_q_points;
7564 for (
unsigned int i = 0; i < size; ++i)
7565 this->values_quad[i] += this->values_from_gradients_quad[i];
7569 for (
unsigned int i = 0; i < size; ++i)
7570 this->values_quad[i] = this->values_from_gradients_quad[i];
7575 if constexpr (fe_degree > -1)
7577 template
run<fe_degree, n_q_points_1d>(n_components,
7578 integration_flag_actual,
7585 integration_flag_actual,
7598 typename VectorizedArrayType>
7605 VectorizedArrayType>::
7608 Assert((integration_flag &
7611 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7612 "and EvaluationFlags::hessians are supported."));
7618 unsigned int size = n_components * dim * n_q_points;
7621 for (
unsigned int i = 0; i < size; ++i)
7622 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7626 for (
unsigned int i = 0; i < size; ++i)
7627 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7632 if (this->data->element_type ==
7636 this->divergence_is_requested ==
false)
7638 unsigned int size = n_components * n_q_points;
7641 for (
unsigned int i = 0; i < size; ++i)
7642 this->values_quad[i] += this->values_from_gradients_quad[i];
7646 for (
unsigned int i = 0; i < size; ++i)
7647 this->values_quad[i] = this->values_from_gradients_quad[i];
7652 if constexpr (fe_degree > -1)
7655 VectorizedArrayType>::template run<fe_degree>(n_components,
7656 integration_flag_actual,
7672 typename VectorizedArrayType>
7679 VectorizedArrayType>::
7681 const bool sum_into_values)
7683 collect_from_face(integration_flag, this->values_dofs, sum_into_values);
7686 this->dof_values_initialized =
true;
7697 typename VectorizedArrayType>
7704 VectorizedArrayType>::
7706 VectorizedArrayType *values_array,
7707 const bool sum_into_values)
7709 Assert((integration_flag &
7712 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7713 "and EvaluationFlags::hessians are supported."));
7720 if (this->data->element_type ==
7724 this->divergence_is_requested ==
false)
7727 if constexpr (fe_degree > -1)
7730 VectorizedArrayType>::template run<fe_degree>(n_components,
7731 integration_flag_actual,
7738 integration_flag_actual,
7751 typename VectorizedArrayType>
7752template <
typename VectorType>
7759 VectorizedArrayType>::
7760 gather_evaluate(
const VectorType &input_vector,
7763 Assert((evaluation_flag &
7766 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7767 "and EvaluationFlags::hessians are supported."));
7769 const auto shared_vector_data = internal::get_shared_vector_data(
7771 this->dof_access_index ==
7773 this->active_fe_index,
7776 if (this->data->data.front().fe_degree > 0 &&
7777 fast_evaluation_supported(this->data->data.front().fe_degree,
7778 this->data->data.front().n_q_points_1d) &&
7781 typename VectorType::value_type,
7782 VectorizedArrayType>::
7783 supports(evaluation_flag,
7785 internal::get_beginning<typename VectorType::value_type>(
7787 this->dof_info->index_storage_variants[this->dof_access_index]
7790 if constexpr (fe_degree > -1)
7794 typename VectorType::value_type,
7795 VectorizedArrayType>::template
run<fe_degree,
7799 internal::get_beginning<typename VectorType::value_type>(
7808 typename VectorType::value_type,
7809 VectorizedArrayType>::evaluate(n_components,
7811 internal::get_beginning<
7812 typename VectorType::value_type>(
7820 this->read_dof_values(input_vector);
7821 this->evaluate(evaluation_flag);
7826 this->gradients_quad_initialized =
7839 typename VectorizedArrayType>
7840template <
typename VectorType>
7848 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
7849 const bool integrate_gradients,
7850 VectorType &destination)
7857 integrate_scatter(flag, destination);
7867 typename VectorizedArrayType>
7868template <
typename VectorType>
7875 VectorizedArrayType>::
7877 VectorType &destination)
7879 Assert((this->dof_access_index ==
7881 this->is_interior_face() ==
false) ==
false,
7884 const auto shared_vector_data = internal::get_shared_vector_data(
7886 this->dof_access_index ==
7888 this->active_fe_index,
7891 if (this->data->data.front().fe_degree > 0 &&
7892 fast_evaluation_supported(this->data->data.front().fe_degree,
7893 this->data->data.front().n_q_points_1d) &&
7896 typename VectorType::value_type,
7897 VectorizedArrayType>::
7898 supports(integration_flag,
7900 internal::get_beginning<typename VectorType::value_type>(
7902 this->dof_info->index_storage_variants[this->dof_access_index]
7905 if constexpr (fe_degree > -1)
7909 typename VectorType::value_type,
7910 VectorizedArrayType>::template
run<fe_degree,
7914 internal::get_beginning<typename VectorType::value_type>(
7923 typename VectorType::value_type,
7924 VectorizedArrayType>::integrate(n_components,
7926 internal::get_beginning<
7927 typename VectorType::value_type>(
7935 integrate(integration_flag);
7936 this->distribute_local_to_global(destination);
7947 typename VectorizedArrayType>
7954 VectorizedArrayType>::dof_indices()
const
7956 return {0U, dofs_per_cell};
7966 typename VectorizedArrayType>
7973 VectorizedArrayType>::
7974 fast_evaluation_supported(
const unsigned int given_degree,
7975 const unsigned int given_n_q_points_1d)
7977 return fe_degree == -1 ?
7990 typename VectorizedArrayType>
7997 VectorizedArrayType>::
7998 fast_evaluation_supported(
const unsigned int given_degree,
7999 const unsigned int given_n_q_points_1d)
8001 return fe_degree == -1 ?
8014 typename VectorizedArrayType>
8021 VectorizedArrayType>::at_boundary()
const
8023 Assert(this->dof_access_index !=
8027 if (this->is_interior_face() ==
false)
8032 this->matrix_free->n_boundary_face_batches()))
8045 typename VectorizedArrayType>
8052 VectorizedArrayType>::boundary_id()
const
8054 Assert(this->dof_access_index !=
8071 typename VectorizedArrayType>
8079 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8081 return this->data->dofs_per_component_on_face;
8091 typename VectorizedArrayType>
8098 VectorizedArrayType>::get_dofs_projected_to_face()
8100 return this->data->dofs_per_component_on_face * n_components_;
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > val_in, const unsigned int q_point)
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
std::conditional_t< n_components_==1, Tensor< 2, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 3, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > > > hessian_type
std::conditional_t< n_components_==1, Tensor< 1, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 2, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > > > gradient_type
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
VectorizedArrayType get_divergence(const unsigned int q_point) const
hessian_type get_hessian(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
static constexpr unsigned int n_components
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
std::conditional_t< n_components_==1, VectorizedArrayType, Tensor< 1, n_components_, VectorizedArrayType > > value_type
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
static constexpr unsigned int n_lanes
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
const unsigned int quad_no
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const unsigned int active_fe_index
std::array< unsigned int, n_lanes > cell_ids
const unsigned int active_quad_index
bool is_interior_face() const
FEEvaluationData & operator=(const FEEvaluationData &other)
const unsigned int first_selected_component
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
void reinit(const std::array< unsigned int, n_lanes > &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
static constexpr unsigned int n_lanes
FEEvaluation(const FEEvaluation &other)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
void project_to_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
unsigned int get_dofs_per_component_projected_to_face()
void project_to_face(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void collect_from_face(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
static constexpr unsigned int static_dofs_per_component
static constexpr unsigned int n_lanes
static constexpr unsigned int static_n_q_points
unsigned int get_dofs_projected_to_face()
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
types::boundary_id boundary_id() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, const bool sum_into_values=false)
void evaluate_in_face(const EvaluationFlags::EvaluationFlags evaluation_flag)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate_in_face(const EvaluationFlags::EvaluationFlags integration_flag)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_DEPRECATED
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_jacobian_grads
Gradient of volume element.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const types::boundary_id internal_face_boundary_id
static const unsigned int invalid_unsigned_int
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static constexpr unsigned int value
@ dof_access_face_exterior
@ dof_access_face_interior
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
std::vector< std::pair< unsigned int, unsigned int > > row_starts
std::vector< std::vector< unsigned int > > component_dof_indices_offset
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< unsigned int > dof_indices
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
std::vector< unsigned int > row_starts_plain_indices
std::vector< unsigned int > component_to_base_index
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
std::vector< unsigned int > plain_dof_indices
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
std::vector< unsigned int > dof_indices_interleaved
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
std::vector< QuadratureDescriptor > descriptor
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)