16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
92 typename VectorizedArrayType>
118 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
156 template <
typename VectorType>
159 const VectorType &src,
161 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().
flip());
191 template <
typename VectorType>
194 const VectorType &src,
196 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().
flip());
229 template <
typename VectorType>
234 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().
flip())
const;
274 template <
typename VectorType>
279 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().
flip())
const;
284 template <
typename VectorType>
289 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().
flip())
const;
365 typename = std::enable_if_t<n_components == n_components_local>>
419 template <
int dim_ = dim,
420 typename = std::enable_if_t<dim_ == 1 && n_components == dim_>>
527 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
546 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
559 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
578 template <
int dim_ = dim,
typename = std::enable_if_t<n_components_ == dim_>>
592 template <
int dim_ = dim,
593 typename = std::enable_if_t<n_components_ == dim_ && dim_ != 1>>
594 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
610 template <
int dim_ = dim,
611 typename = std::enable_if_t<n_components_ == dim_ && dim != 1>>
655 const unsigned int dof_no,
658 const unsigned int fe_degree,
659 const unsigned int n_q_points,
663 const unsigned int face_type);
737 template <
typename VectorType,
typename VectorOperation>
741 const std::array<VectorType *, n_components_> &vectors,
745 const std::bitset<n_lanes> &mask,
755 template <
typename VectorType,
typename VectorOperation>
759 const std::array<VectorType *, n_components_> &vectors,
763 const std::bitset<n_lanes> &mask)
const;
772 template <
typename VectorType,
typename VectorOperation>
776 const std::array<VectorType *, n_components_> &vectors)
const;
1380 typename VectorizedArrayType>
1385 VectorizedArrayType>
1388 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1389 "Type of Number and of VectorizedArrayType do not match.");
1431 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1516 const unsigned int dof_no = 0,
1517 const unsigned int quad_no = 0,
1530 const std::pair<unsigned int, unsigned int> &
range,
1531 const unsigned int dof_no = 0,
1532 const unsigned int quad_no = 0,
1642 template <
bool level_dof_access>
1707 template <
typename VectorType>
1753 template <
typename VectorType>
1840 int n_q_points_1d = fe_degree + 1,
1842 typename Number = double,
1848 VectorizedArrayType>
1851 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1852 "Type of Number and of VectorizedArrayType do not match.");
1894 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1996 const unsigned int dof_no = 0,
1997 const unsigned int quad_no = 0,
2012 const std::pair<unsigned int, unsigned int> &
range,
2014 const unsigned int dof_no = 0,
2015 const unsigned int quad_no = 0,
2110 template <
typename VectorType>
2179 template <
typename VectorType>
2187 template <
typename VectorType>
2190 const bool integrate_gradients,
2267 namespace MatrixFreeFunctions
2271 template <
int dim,
int degree>
2280 template <
int degree>
2283 static constexpr unsigned int value = degree + 1;
2299 template <
bool is_face,
2302 typename VectorizedArrayType>
2307 const unsigned int dof_no,
2308 const unsigned int first_selected_component,
2309 const unsigned int quad_no,
2310 const unsigned int fe_degree,
2311 const unsigned int n_q_points,
2314 const unsigned int face_type)
2321 &internal::MatrixFreeFunctions::
2322 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2323 matrix_free.get_mapping_info(), quad_no);
2327 init_data.dof_info->fe_index_from_degree(first_selected_component,
2337 std::min<unsigned int>(
2339 init_data.mapping_data->descriptor.size() /
2340 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2342 init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2344 init_data.shape_info = &matrix_free.get_shape_info(
2347 init_data.dof_info->component_to_base_index[first_selected_component],
2353 (
init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2369 typename VectorizedArrayType>
2374 VectorizedArrayType>
::
2377 const unsigned int dof_no,
2378 const unsigned int first_selected_component,
2379 const unsigned int quad_no,
2380 const unsigned int fe_degree,
2381 const unsigned int n_q_points,
2382 const bool is_interior_face,
2383 const unsigned int active_fe_index,
2384 const unsigned int active_quad_index,
2385 const unsigned int face_type)
2389 first_selected_component,
2398 first_selected_component)
2399 , scratch_data_array(matrix_free.acquire_scratch_data())
2400 , matrix_free(&matrix_free)
2404 this->dof_info->start_components.back() == 1 ||
2407 this->dof_info->start_components
2408 [
this->dof_info->component_to_base_index[first_selected_component] +
2410 first_selected_component,
2412 "You tried to construct a vector-valued evaluator with " +
2413 std::to_string(n_components) +
2414 " components. However, "
2415 "the current base element has only " +
2417 this->dof_info->start_components
2418 [
this->dof_info->component_to_base_index[first_selected_component] +
2420 first_selected_component) +
2421 " components left when starting from local element index " +
2423 first_selected_component -
2424 this->dof_info->start_components
2425 [
this->dof_info->component_to_base_index[first_selected_component]]) +
2426 " (global index " + std::to_string(first_selected_component) +
")"));
2438 typename VectorizedArrayType>
2443 VectorizedArrayType>
::
2449 const unsigned int first_selected_component,
2453 other->mapped_geometry->get_quadrature() == quadrature ?
2454 other->mapped_geometry :
2456 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2461 first_selected_component)
2466 fe.component_to_base_index(first_selected_component).first;
2469 first_selected_component >=
2471 ExcMessage(
"The underlying element must at least contain as many "
2472 "components as requested by this class"));
2477 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2479 fe.component_to_base_index(first_selected_component).
first);
2490 typename VectorizedArrayType>
2495 VectorizedArrayType>
::
2500 VectorizedArrayType> &
other)
2504 other.matrix_free->acquire_scratch_data())
2505 , matrix_free(
other.matrix_free)
2507 if (
other.matrix_free ==
nullptr)
2514 this->mapped_geometry =
2515 std::make_shared<internal::MatrixFreeFunctions::
2516 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2517 other.mapped_geometry->get_fe_values().get_mapping(),
2518 other.mapped_geometry->get_quadrature(),
2519 other.mapped_geometry->get_fe_values().get_update_flags());
2521 if constexpr (is_face ==
false)
2522 this->mapping_data = &this->mapped_geometry->get_data_storage();
2526 "face evaluators is not currently "
2532 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2534 this->mapped_geometry->get_data_storage().JxW_values.begin();
2535 this->jacobian_gradients =
2536 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2537 this->jacobian_gradients_non_inverse =
2538 this->mapped_geometry->get_data_storage()
2539 .jacobian_gradients_non_inverse[0]
2542 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2554 typename VectorizedArrayType>
2559 VectorizedArrayType> &
2565 VectorizedArrayType> &
other)
2568 if (matrix_free ==
nullptr)
2571 delete scratch_data_array;
2575 matrix_free->release_scratch_data(scratch_data_array);
2580 matrix_free =
other.matrix_free;
2582 if (
other.matrix_free ==
nullptr)
2590 this->mapped_geometry =
2591 std::make_shared<internal::MatrixFreeFunctions::
2592 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2593 other.mapped_geometry->get_fe_values().get_mapping(),
2594 other.mapped_geometry->get_quadrature(),
2595 other.mapped_geometry->get_fe_values().get_update_flags());
2597 if constexpr (is_face ==
false)
2598 this->mapping_data = &this->mapped_geometry->get_data_storage();
2602 "face evaluators is not currently "
2607 this->mapped_geometry->get_data_storage().jacobians[0].begin();
2609 this->mapped_geometry->get_data_storage().JxW_values.begin();
2610 this->jacobian_gradients =
2611 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2612 this->jacobian_gradients_non_inverse =
2613 this->mapped_geometry->get_data_storage()
2614 .jacobian_gradients_non_inverse[0]
2617 this->mapped_geometry->get_data_storage().quadrature_points.begin();
2621 scratch_data_array = matrix_free->acquire_scratch_data();
2635 typename VectorizedArrayType>
2640 VectorizedArrayType>::~FEEvaluationBase()
2642 if (matrix_free !=
nullptr)
2646 matrix_free->release_scratch_data(scratch_data_array);
2653 delete scratch_data_array;
2664 typename VectorizedArrayType>
2669 Assert(matrix_free !=
nullptr,
2671 "FEEvaluation was not initialized with a MatrixFree object!"));
2672 return *matrix_free;
2681 template <
typename VectorType,
bool>
2684 template <
typename VectorType>
2690 template <
typename VectorType>
2699 template <
typename VectorType,
bool>
2702 template <
typename VectorType>
2713 return &
vec.block(component);
2717 template <
typename VectorType>
2742 template <
typename VectorType>
2749 const unsigned int component)
2752 return &
vec[component];
2756 template <
typename VectorType>
2763 const unsigned int component)
2766 return &
vec[component];
2770 template <
typename VectorType>
2777 const unsigned int component)
2780 return vec[component];
2784 template <
typename VectorType>
2791 const unsigned int component)
2794 return vec[component];
2798 template <
typename VectorType, std::
size_t N>
2805 const unsigned int component)
2808 return vec[component];
2819 typename VectorizedArrayType>
2820template <
typename VectorType,
typename VectorOperation>
2825 const std::array<VectorType *, n_components_> &src,
2829 const std::bitset<n_lanes> &
mask,
2835 if (this->matrix_free ==
nullptr)
2837 read_write_operation_global(
operation, src);
2844 if (this->n_fe_components == 1)
2845 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
2848 ExcMessage(
"The finite element underlying this FEEvaluation "
2849 "object is scalar, but you requested " +
2850 std::to_string(n_components) +
2851 " components via the template argument in "
2852 "FEEvaluation. In that case, you must pass an "
2853 "std::vector<VectorType> or a BlockVector to " +
2854 "read_dof_values and distribute_local_to_global."));
2867 this->dof_access_index ==
2869 this->is_interior_face() ==
false;
2879 bool is_contiguous =
true;
2883 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
2890 if (
mask[v] ==
true &&
2893 [cells[v] / n_lanes] <
2896 is_contiguous =
false;
2900 this->dof_access_index :
2901 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2902 [
this->cell] <
internal::MatrixFreeFunctions::DoFInfo::
2905 is_contiguous =
false;
2917 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
2921 for (
unsigned int v = 0; v < n_lanes; ++v)
2922 if (
mask[v] ==
false)
2927 if (is_face ==
false)
2933 [
this->first_selected_component])
2934 for (
unsigned int v = 0; v < n_lanes; ++v)
2942 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
2948 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
2949 std::array<VectorizedArrayType *, n_components> values_dofs;
2950 for (
unsigned int c = 0; c < n_components; ++c)
2951 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
2952 c * dofs_per_component;
2956 [is_face ?
this->dof_access_index :
2957 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
2958 [
this->cell] ==
internal::MatrixFreeFunctions::DoFInfo::
2959 IndexStorageVariants::interleaved &&
2962 const unsigned int *dof_indices =
2964 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
2968 [this->first_selected_component] *
2971 std::array<typename VectorType::value_type *, n_components>
src_ptrs;
2972 if (n_components == 1 || this->n_fe_components == 1)
2973 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
2978 const_cast<typename VectorType::value_type *
>(src[0]->begin());
2980 if (n_components == 1 || this->n_fe_components == 1)
2981 for (
unsigned int i = 0; i < dofs_per_component;
2982 ++i, dof_indices += n_lanes)
2983 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
2984 operation.process_dof_gather(dof_indices,
2988 values_dofs[
comp][i],
2991 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
2992 for (
unsigned int i = 0; i < dofs_per_component;
2993 ++i, dof_indices += n_lanes)
2994 operation.process_dof_gather(dof_indices,
2998 values_dofs[
comp][i],
3005 std::array<const unsigned int *, n_lanes> dof_indices;
3006 dof_indices.fill(
nullptr);
3013 this->n_fe_components > 1 ? n_components : 1;
3017 for (
unsigned int v = 0; v < n_lanes; ++v)
3024 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3025 this->first_selected_component];
3040 for (
unsigned int v = 0; v < n_lanes; ++v)
3046 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3047 this->first_selected_component];
3057 dof_info.hanging_node_constraint_masks_comp
3058 [
this->active_fe_index][
this->first_selected_component])
3072 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3073 return i != numbers::invalid_unsigned_int;
3075 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3076 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3083 if (n_components == 1 || this->n_fe_components == 1)
3085 for (
unsigned int v = 0; v < n_lanes; ++v)
3090 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3091 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3092 operation.process_dof(dof_indices[v][i],
3094 values_dofs[
comp][i][v]);
3099 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3100 for (
unsigned int v = 0; v < n_lanes; ++v)
3105 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3107 dof_indices[v][
comp * dofs_per_component + i],
3109 values_dofs[
comp][i][v]);
3121 for (
unsigned int v = 0; v < n_lanes; ++v)
3128 cell_index * this->n_fe_components + this->first_selected_component;
3130 this->n_fe_components > 1 ? n_components : 1;
3146 dof_info.hanging_node_constraint_masks_comp
3147 [
this->active_fe_index][
this->first_selected_component])))
3156 [this->first_selected_component] +
3161 if (n_components == 1 || this->n_fe_components == 1)
3166 const std::pair<unsigned short, unsigned short>
indicator =
3170 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3180 Number
value[n_components];
3181 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3186 this->matrix_free->constraint_pool_begin(
indicator.second);
3188 this->matrix_free->constraint_pool_end(
indicator.second);
3190 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3191 operation.process_constraint(*dof_indices[v],
3196 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3205 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3217 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3225 const std::pair<unsigned short, unsigned short>
indicator =
3243 this->matrix_free->constraint_pool_begin(
indicator.second);
3245 this->matrix_free->constraint_pool_end(
indicator.second);
3248 operation.process_constraint(*dof_indices[v],
3284 typename VectorizedArrayType>
3285template <
typename VectorType,
typename VectorOperation>
3290 const std::array<VectorType *, n_components_> &src)
const
3294 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3295 unsigned int index = this->first_selected_component * dofs_per_component;
3296 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3298 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3301 this->values_dofs[
comp * dofs_per_component + i]);
3303 local_dof_indices[this->
data->lexicographic_numbering[index]],
3305 this->values_dofs[
comp * dofs_per_component + i][0]);
3316 typename VectorizedArrayType>
3317template <
typename VectorType,
typename VectorOperation>
3322 const std::array<VectorType *, n_components_> &src,
3326 const std::bitset<n_lanes> &
mask)
const
3335 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3338 is_face ? this->dof_access_index :
3346 const std::size_t dofs_per_component = this->
data->dofs_per_component_on_cell;
3347 std::array<VectorizedArrayType *, n_components> values_dofs{{
nullptr}};
3348 for (
unsigned int c = 0; c < n_components; ++c)
3349 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3350 c * dofs_per_component;
3355 this->dof_access_index ==
3357 this->is_interior_face() ==
false;
3363 interleaved_contiguous &&
3366 const unsigned int dof_index =
3370 [this->first_selected_component] *
3372 if (n_components == 1 || this->n_fe_components == 1)
3373 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3374 operation.process_dofs_vectorized(dofs_per_component,
3380 operation.process_dofs_vectorized(dofs_per_component * n_components,
3388 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3402 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs{
3405 const auto upper_bound =
3407 for (
unsigned int v = 0; v < upper_bound; ++v)
3409 if (
mask[v] ==
false)
3411 vector_ptrs[v] =
nullptr;
3431 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3434 [this->active_fe_index][this->first_selected_component]);
3436 vector_ptrs[v] =
nullptr;
3439 vector_ptrs[v] =
nullptr;
3446 if (n_components == 1 || this->n_fe_components == 1)
3448 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3451 operation.process_dofs_vectorized_transpose(
3461 operation.process_dofs_vectorized_transpose(dofs_per_component *
3469 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3472 (n_components == 1 || this->n_fe_components == 1) ?
comp : 0);
3474 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3477 if (n_components == 1 || this->n_fe_components == 1)
3480 if (
mask[v] ==
true)
3481 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3482 operation.process_dof(vector_ptrs[v][i],
3483 values_dofs[
comp][i][v]);
3488 if (
mask[v] ==
true)
3489 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3491 vector_ptrs[v][i +
comp * dofs_per_component],
3492 values_dofs[
comp][i][v]);
3498 std::array<unsigned int, n_lanes> dof_indices{
3505 if (
mask[v] ==
true)
3510 [this->first_selected_component] *
3522 if (n_components == 1 || this->n_fe_components == 1)
3523 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3524 operation.process_dofs_vectorized_transpose(dofs_per_component,
3530 operation.process_dofs_vectorized_transpose(dofs_per_component *
3539 interleaved_contiguous_strided)
3541 std::array<typename VectorType::value_type *, n_components>
src_ptrs{
3543 if (n_components == 1 || this->n_fe_components == 1)
3544 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3549 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3551 if (n_components == 1 || this->n_fe_components == 1)
3552 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3554 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3555 operation.process_dof_gather(dof_indices.data(),
3559 values_dofs[
comp][i],
3563 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3564 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3569 (
comp * dofs_per_component + i) * n_lanes,
3570 src_ptrs[0] + (
comp * dofs_per_component + i) * n_lanes,
3571 values_dofs[
comp][i],
3579 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3581 std::array<typename VectorType::value_type *, n_components>
src_ptrs{
3583 if (n_components == 1 || this->n_fe_components == 1)
3584 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3589 const_cast<typename VectorType::value_type *
>(src[0]->begin());
3593 if (n_components == 1 || this->n_fe_components == 1)
3594 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3596 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3597 operation.process_dof_gather(dof_indices.data(),
3601 values_dofs[
comp][i],
3604 for (
unsigned int v = 0; v < n_lanes; ++v)
3608 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3609 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3611 operation.process_dof_gather(dof_indices.data(),
3615 values_dofs[
comp][i],
3618 for (
unsigned int v = 0; v < n_lanes; ++v)
3624 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
3626 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3631 if (
mask[v] ==
true)
3634 [
ind][cells[v] / VectorizedArrayType::size()] ==
3638 if (n_components == 1 || this->n_fe_components == 1)
3640 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3641 operation.process_dof(dof_indices[v] + i,
3643 values_dofs[
comp][i][v]);
3647 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3648 operation.process_dof(dof_indices[v] + i +
3649 comp * dofs_per_component,
3651 values_dofs[
comp][i][v]);
3656 const unsigned int offset =
3659 if (n_components == 1 || this->n_fe_components == 1)
3661 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3662 operation.process_dof(dof_indices[v] + i * offset,
3664 values_dofs[
comp][i][v]);
3668 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3671 (i +
comp * dofs_per_component) * offset,
3673 values_dofs[
comp][i][v]);
3684 if (n_components == 1 || this->n_fe_components == 1)
3687 if (
mask[v] ==
true)
3688 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3689 operation.process_dof(dof_indices[v] + i,
3691 values_dofs[
comp][i][v]);
3696 if (
mask[v] ==
true)
3697 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3698 operation.process_dof(dof_indices[v] + i +
3699 comp * dofs_per_component,
3701 values_dofs[
comp][i][v]);
3708 [
ind][VectorizedArrayType::size() * this->cell];
3711 if (n_components == 1 || this->n_fe_components == 1)
3714 if (
mask[v] ==
true)
3715 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3718 values_dofs[
comp][i][v]);
3723 if (
mask[v] ==
true)
3724 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3729 values_dofs[
comp][i][v]);
3741 std::enable_if_t<!IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3742 decltype(std::declval<VectorType>().begin())
3751 std::enable_if_t<IsBlockVector<VectorType>::value,
VectorType> * =
nullptr>
3752 typename VectorType::value_type *
3759 std::enable_if_t<has_shared_vector_data<VectorType>,
VectorType> * =
3761 const std::vector<ArrayView<const typename VectorType::value_type>> *
3764 const unsigned int active_fe_index,
3771 active_fe_index == 0)
3772 return &
vec->shared_vector_data();
3778 std::enable_if_t<!has_shared_vector_data<VectorType>,
VectorType>
3780 const std::vector<ArrayView<const typename VectorType::value_type>> *
3789 template <
int n_components,
typename VectorType>
3791 std::array<
typename internal::BlockVectorSelector<
3796 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
3803 const unsigned int active_fe_index,
3809 std::array<
typename internal::BlockVectorSelector<
3815 ArrayView<
const typename internal::BlockVectorSelector<
3821 for (
unsigned int d = 0;
d < n_components; ++
d)
3822 src_data.first[d] = internal::BlockVectorSelector<
3828 for (
unsigned int d = 0;
d < n_components; ++
d)
3830 const_cast<typename internal::BlockVectorSelector<
3831 std::remove_const_t<VectorType>,
3848 typename VectorizedArrayType>
3853 if (this->dof_info ==
nullptr ||
3855 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
3856 this->dof_info->hanging_node_constraint_masks_comp
3857 [
this->active_fe_index][
this->first_selected_component] ==
false)
3860 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
3861 constraint_mask{{internal::MatrixFreeFunctions::
3862 unconstrained_compressed_constraint_kind}};
3866 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3868 for (
unsigned int v = 0; v < n_lanes; ++v)
3880 constraint_mask[v] =
mask;
3891 this->
data->data.front().fe_degree,
3892 this->get_shape_info(),
3904 typename VectorizedArrayType>
3905template <
typename VectorType>
3910 const std::bitset<n_lanes> &
mask)
3912 const auto src_data = internal::get_vector_data<n_components_>(
3915 this->dof_info !=
nullptr &&
3916 this->dof_access_index ==
3918 this->active_fe_index,
3924 apply_hanging_node_constraints(
false);
3927 this->dof_values_initialized =
true;
3937 typename VectorizedArrayType>
3938template <
typename VectorType>
3943 const std::bitset<n_lanes> &
mask)
3945 const auto src_data = internal::get_vector_data<n_components_>(
3948 this->dof_access_index ==
3950 this->active_fe_index,
3957 this->dof_values_initialized =
true;
3967 typename VectorizedArrayType>
3968template <
typename VectorType>
3973 const std::bitset<n_lanes> &
mask)
const
3976 Assert(this->dof_values_initialized ==
true,
3980 apply_hanging_node_constraints(
true);
3982 const auto dst_data = internal::get_vector_data<n_components_>(
3985 this->dof_access_index ==
3987 this->active_fe_index,
4001 typename VectorizedArrayType>
4002template <
typename VectorType>
4007 const std::bitset<n_lanes> &
mask)
const
4010 Assert(this->dof_values_initialized ==
true,
4014 const auto dst_data = internal::get_vector_data<n_components_>(
4017 this->dof_access_index ==
4019 this->active_fe_index,
4032 typename VectorizedArrayType>
4033template <
typename VectorType>
4038 const std::bitset<n_lanes> &
mask)
const
4041 Assert(this->dof_values_initialized ==
true,
4045 const auto dst_data = internal::get_vector_data<n_components_>(
4048 this->dof_access_index ==
4050 this->active_fe_index,
4067 typename VectorizedArrayType>
4073 VectorizedArrayType>::value_type
4078 if constexpr (n_components == 1)
4079 return this->values_dofs[dof];
4082 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4083 Tensor<1, n_components_, VectorizedArrayType> return_value;
4084 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4085 return_value[
comp] = this->values_dofs[
comp * dofs + dof];
4086 return return_value;
4096 typename VectorizedArrayType>
4102 VectorizedArrayType>::value_type
4107 Assert(this->values_quad_initialized ==
true,
4112 if constexpr (n_components == 1)
4113 return this->values_quad[
q_point];
4116 if (n_components == dim &&
4117 this->
data->element_type ==
4122 Assert(this->values_quad_initialized ==
true,
4127 Assert(this->J_value !=
nullptr,
4130 const std::size_t
nqp = this->n_quadrature_points;
4138 const VectorizedArrayType
inv_det =
4139 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4140 this->jacobian[0][0][0] *
this->jacobian[0][1][1] *
4141 this->jacobian[0][2][2];
4144 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4162 const VectorizedArrayType
inv_det =
4163 (is_face && dim == 2 && this->get_face_no() < 2) ?
4167 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4170 for (
unsigned int e = 1;
e < dim; ++
e)
4180 const std::size_t
nqp = this->n_quadrature_points;
4182 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4184 return return_value;
4195 typename VectorizedArrayType>
4201 VectorizedArrayType>::gradient_type
4206 Assert(this->gradients_quad_initialized ==
true,
4211 Assert(this->jacobian !=
nullptr,
4213 "update_gradients"));
4214 const std::size_t
nqp = this->n_quadrature_points;
4216 if constexpr (n_components == dim && dim > 1)
4218 if (this->
data->element_type ==
4223 Assert(this->gradients_quad_initialized ==
true,
4228 Assert(this->jacobian !=
nullptr,
4230 "update_gradients"));
4231 const std::size_t
nqp = this->n_quadrature_points;
4232 const std::size_t
nqp_d =
nqp * dim;
4235 this->gradients_quad +
q_point * dim;
4246 const VectorizedArrayType
inv_det =
4247 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
4248 this->jacobian[0][0][0] *
this->jacobian[0][1][1] *
4249 this->jacobian[0][2][2];
4252 for (
unsigned int d = 0;
d < dim; ++
d)
4253 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4268 const VectorizedArrayType
inv_det =
4269 (is_face && dim == 2 && this->get_face_no() < 2) ?
4273 VectorizedArrayType tmp[dim][dim];
4275 for (
unsigned int d = 0;
d < dim; ++
d)
4276 for (
unsigned int e = 0;
e < dim; ++
e)
4279 for (
unsigned int f = 1; f < dim; ++f)
4282 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4283 for (
unsigned int d = 0;
d < dim; ++
d)
4285 VectorizedArrayType
res =
jac[
comp][0] * tmp[
d][0];
4286 for (
unsigned int f = 1; f < dim; ++f)
4300 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
4302 "update_hessians"));
4305 this->jacobian_gradients_non_inverse[
q_point];
4311 const VectorizedArrayType
inv_det =
4312 (is_face && dim == 2 && this->get_face_no() < 2) ?
4319 VectorizedArrayType tmp[dim][dim];
4320 for (
unsigned int d = 0;
d < dim; ++
d)
4321 for (
unsigned int e = 0;
e < dim; ++
e)
4324 for (
unsigned int f = 1; f < dim; ++f)
4325 tmp[e][d] +=
t_jac[f][d] * gradients[f *
nqp_d + e];
4330 for (
unsigned int d = 0;
d < dim; ++
d)
4332 for (
unsigned int e = 0;
e < dim; ++
e)
4335 for (
unsigned int f = 0,
r = dim; f < dim; ++f)
4336 for (
unsigned int k = f + 1;
k < dim; ++
k, ++
r)
4346 for (
unsigned int d = 0;
d < dim; ++
d)
4347 for (
unsigned int e = 0;
e < dim; ++
e)
4350 for (
unsigned int f = 1; f < dim; ++f)
4358 VectorizedArrayType
tmp3[dim],
tmp4[dim];
4359 for (
unsigned int d = 0;
d < dim; ++
d)
4362 for (
unsigned int e = 1;
e < dim; ++
e)
4365 for (
unsigned int e = 0,
k = dim;
e < dim; ++
e)
4366 for (
unsigned int f = e + 1; f < dim; ++
k, ++f)
4367 for (
unsigned int d = 0;
d < dim; ++
d)
4372 for (
unsigned int d = 0;
d < dim; ++
d)
4375 for (
unsigned int e = 1;
e < dim; ++
e)
4379 VectorizedArrayType
tmp2[dim];
4380 for (
unsigned int d = 0;
d < dim; ++
d)
4383 for (
unsigned e = 1;
e < dim; ++
e)
4388 for (
unsigned int d = 0;
d < dim; ++
d)
4389 for (
unsigned int e = 0;
e < dim; ++
e)
4406 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4407 for (
unsigned int d = 0;
d < dim; ++
d)
4410 this->jacobian[0][
d][
d];
4419 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4420 for (
unsigned int d = 0;
d < dim; ++
d)
4424 for (
unsigned int e = 1;
e < dim; ++
e)
4430 if constexpr (n_components == 1)
4442 typename VectorizedArrayType>
4448 VectorizedArrayType>::value_type
4454 Assert(this->gradients_quad_initialized ==
true,
4458 Assert(this->normal_x_jacobian !=
nullptr,
4460 "update_gradients"));
4462 const std::size_t
nqp = this->n_quadrature_points;
4466 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4469 (this->normal_x_jacobian[0][dim - 1]);
4472 const std::size_t
index =
4474 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4477 this->normal_x_jacobian[
index][0];
4478 for (
unsigned int d = 1;
d < dim; ++
d)
4481 this->normal_x_jacobian[
index][
d];
4484 if constexpr (n_components == 1)
4496 template <
typename VectorizedArrayType>
4499 const VectorizedArrayType *
const hessians,
4501 VectorizedArrayType (&tmp)[1][1])
4506 template <
typename VectorizedArrayType>
4509 const VectorizedArrayType *
const hessians,
4510 const unsigned int nqp,
4511 VectorizedArrayType (&tmp)[2][2])
4513 for (
unsigned int d = 0;
d < 2; ++
d)
4521 template <
typename VectorizedArrayType>
4524 const VectorizedArrayType *
const hessians,
4525 const unsigned int nqp,
4526 VectorizedArrayType (&tmp)[3][3])
4528 for (
unsigned int d = 0;
d < 3; ++
d)
4549 typename VectorizedArrayType>
4554 VectorizedArrayType>::hessian_type
4559 Assert(this->hessians_quad_initialized ==
true,
4564 Assert(this->jacobian !=
nullptr,
4574 const std::size_t
nqp = this->n_quadrature_points;
4575 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4580 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4582 for (
unsigned int d = 0;
d < dim; ++
d)
4609 for (
unsigned int d = 0;
d < dim; ++
d)
4610 for (
unsigned int e = d + 1;
e < dim; ++
e)
4617 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4619 VectorizedArrayType tmp[dim][dim];
4620 internal::hessian_unit_times_jac(
4624 for (
unsigned int d = 0;
d < dim; ++
d)
4625 for (
unsigned int e = d;
e < dim; ++
e)
4628 for (
unsigned int f = 1; f < dim; ++f)
4636 for (
unsigned int d = 0;
d < dim; ++
d)
4637 for (
unsigned int e = d + 1;
e < dim; ++
e)
4645 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4647 VectorizedArrayType tmp[dim][dim];
4648 internal::hessian_unit_times_jac(
4652 for (
unsigned int d = 0;
d < dim; ++
d)
4653 for (
unsigned int e = d;
e < dim; ++
e)
4656 for (
unsigned int f = 1; f < dim; ++f)
4661 for (
unsigned int d = 0;
d < dim; ++
d)
4662 for (
unsigned int e = 0;
e < dim; ++
e)
4668 for (
unsigned int d = 0,
count = dim;
d < dim; ++
d)
4669 for (
unsigned int e = d + 1;
e < dim; ++
e, ++
count)
4670 for (
unsigned int f = 0; f < dim; ++f)
4676 for (
unsigned int d = 0;
d < dim; ++
d)
4677 for (
unsigned int e = d + 1;
e < dim; ++
e)
4681 if constexpr (n_components == 1)
4693 typename VectorizedArrayType>
4698 VectorizedArrayType>::gradient_type
4704 Assert(this->hessians_quad_initialized ==
true,
4715 const std::size_t
nqp = this->n_quadrature_points;
4716 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4722 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4723 for (
unsigned int d = 0;
d < dim; ++
d)
4731 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4735 VectorizedArrayType tmp[dim][dim];
4736 internal::hessian_unit_times_jac(
4741 for (
unsigned int d = 0;
d < dim; ++
d)
4744 for (
unsigned int f = 1; f < dim; ++f)
4753 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4757 VectorizedArrayType tmp[dim][dim];
4758 internal::hessian_unit_times_jac(
4763 for (
unsigned int d = 0;
d < dim; ++
d)
4766 for (
unsigned int f = 1; f < dim; ++f)
4770 for (
unsigned int d = 0;
d < dim; ++
d)
4771 for (
unsigned int e = 0;
e < dim; ++
e)
4778 if constexpr (n_components == 1)
4790 typename VectorizedArrayType>
4795 VectorizedArrayType>::value_type
4801 Assert(this->hessians_quad_initialized ==
true,
4807 if constexpr (n_components == 1)
4810 for (
unsigned int d = 1;
d < dim; ++
d)
4817 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4820 for (
unsigned int d = 1;
d < dim; ++
d)
4833 typename VectorizedArrayType>
4838 VectorizedArrayType>::value_type
4843 Assert(this->hessians_quad_initialized ==
true,
4848 Assert(this->normal_x_jacobian !=
nullptr,
4850 "update_hessians"));
4854 const std::size_t
nqp = this->n_quadrature_points;
4855 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4859 const auto nxj = this->normal_x_jacobian[0];
4861 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4863 for (
unsigned int d = 0;
d < dim; ++
d)
4896 const auto normal = this->normal_vector(
q_point);
4899 if constexpr (n_components == 1)
4902 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4905 if constexpr (n_components == 1)
4917 typename VectorizedArrayType>
4923 this->dof_values_initialized =
true;
4925 const std::size_t dofs = this->
data->dofs_per_component_on_cell;
4927 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
4928 if constexpr (n_components == 1)
4929 this->values_dofs[
comp * dofs + dof] =
val_in;
4940 typename VectorizedArrayType>
4949 Assert(this->J_value !=
nullptr,
4953 this->values_quad_submitted =
true;
4956 const std::size_t
nqp = this->n_quadrature_points;
4959 const VectorizedArrayType JxW =
4961 this->J_value[0] * this->quadrature_weights[
q_point] :
4963 if constexpr (n_components == 1)
4964 values[0] =
val_in * JxW;
4967 if (n_components == dim &&
4968 this->
data->element_type ==
4973 Assert(this->J_value !=
nullptr,
4978 this->values_quad_submitted =
true;
4982 const std::size_t
nqp = this->n_quadrature_points;
4988 const VectorizedArrayType weight =
4989 this->quadrature_weights[
q_point];
4991 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5004 const VectorizedArrayType
fac =
5006 this->quadrature_weights[
q_point] :
5010 ((dim == 2 &&
this->get_face_no() < 2) ?
5019 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5022 for (
unsigned int e = 1;
e < dim; ++
e)
5029 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5040 typename VectorizedArrayType>
5041template <
int,
typename>
5047 static_assert(n_components == 1,
5048 "Do not try to modify the default template parameters used for"
5049 " selectively enabling this function via std::enable_if!");
5059 typename VectorizedArrayType>
5068 Assert(this->J_value !=
nullptr,
5070 "update_gradients"));
5071 Assert(this->jacobian !=
nullptr,
5073 "update_gradients"));
5075 this->gradients_quad_submitted =
true;
5078 if constexpr (dim > 1 && n_components == dim)
5080 if (this->
data->element_type ==
5089 Assert(this->J_value !=
nullptr,
5091 "update_gradients"));
5092 Assert(this->jacobian !=
nullptr,
5094 "update_gradients"));
5096 this->gradients_quad_submitted =
true;
5100 VectorizedArrayType *
values =
5101 this->values_from_gradients_quad +
q_point;
5102 const std::size_t
nqp = this->n_quadrature_points;
5103 const std::size_t
nqp_d =
nqp * dim;
5113 const VectorizedArrayType weight =
5114 this->quadrature_weights[
q_point];
5115 for (
unsigned int d = 0;
d < dim; ++
d)
5116 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5131 const VectorizedArrayType
fac =
5133 this->quadrature_weights[
q_point] :
5135 ((dim == 2 &&
this->get_face_no() < 2) ?
5140 VectorizedArrayType tmp[dim][dim];
5141 for (
unsigned int d = 0;
d < dim; ++
d)
5142 for (
unsigned int e = 0;
e < dim; ++
e)
5145 for (
unsigned int f = 1; f < dim; ++f)
5148 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5149 for (
unsigned int d = 0;
d < dim; ++
d)
5151 VectorizedArrayType
res =
jac[0][
comp] * tmp[
d][0];
5152 for (
unsigned int f = 1; f < dim; ++f)
5163 this->jacobian_gradients_non_inverse[
q_point];
5169 const VectorizedArrayType
fac =
5170 (!is_face) ? this->quadrature_weights[
q_point] :
5172 ((dim == 2 &&
this->get_face_no() < 2) ?
5181 VectorizedArrayType
tmp3[dim],
tmp4[dim];
5182 for (
unsigned int d = 0;
d < dim; ++
d)
5185 for (
unsigned int e = 1;
e < dim; ++
e)
5188 for (
unsigned int e = 0,
k = dim;
e < dim; ++
e)
5189 for (
unsigned int f = e + 1; f < dim; ++
k, ++f)
5190 for (
unsigned int d = 0;
d < dim; ++
d)
5195 for (
unsigned int d = 0;
d < dim; ++
d)
5198 for (
unsigned int e = 1;
e < dim; ++
e)
5205 VectorizedArrayType tmp[dim][dim];
5208 for (
unsigned int d = 0;
d < dim; ++
d)
5209 for (
unsigned int e = 0;
e < dim; ++
e)
5212 for (
unsigned int f = 1; f < dim; ++f)
5216 for (
unsigned int d = 0;
d < dim; ++
d)
5217 for (
unsigned int e = 0;
e < dim; ++
e)
5219 VectorizedArrayType
res =
t_jac[
d][0] * tmp[
e][0];
5220 for (
unsigned int f = 1; f < dim; ++f)
5228 VectorizedArrayType
value[dim];
5229 for (
unsigned int d = 0;
d < dim; ++
d)
5232 for (
unsigned int e = 1;
e < dim; ++
e)
5235 for (
unsigned int e = 0,
k = dim;
e < dim; ++
e)
5236 for (
unsigned int f = e + 1; f < dim; ++
k, ++f)
5237 for (
unsigned int d = 0;
d < dim; ++
d)
5245 for (
unsigned int d = 0;
d < dim; ++
d)
5248 for (
unsigned int e = 1;
e < dim; ++
e)
5250 for (
unsigned int e = 0;
e < dim; ++
e)
5254 for (
unsigned int d = 0;
d < dim; ++
d)
5261 const std::size_t
nqp_d = this->n_quadrature_points * dim;
5266 const VectorizedArrayType JxW =
5267 this->J_value[0] * this->quadrature_weights[
q_point];
5273 std::array<VectorizedArrayType, dim>
jac;
5274 for (
unsigned int d = 0;
d < dim; ++
d)
5275 jac[d] = this->jacobian[0][d][d];
5277 for (
unsigned int d = 0;
d < dim; ++
d)
5279 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
5280 if constexpr (n_components == 1)
5281 gradients[d] =
grad_in[d] * factor;
5283 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5293 const VectorizedArrayType JxW =
5296 this->J_value[0] * this->quadrature_weights[
q_point];
5297 if constexpr (n_components == 1)
5298 for (
unsigned int d = 0;
d < dim; ++
d)
5301 for (
unsigned int e = 1;
e < dim; ++
e)
5306 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5307 for (
unsigned int d = 0;
d < dim; ++
d)
5310 for (
unsigned int e = 1;
e < dim; ++
e)
5323 typename VectorizedArrayType>
5324template <
int,
typename>
5330 static_assert(n_components == 1 && dim == 1,
5331 "Do not try to modify the default template parameters used for"
5332 " selectively enabling this function via std::enable_if!");
5342 typename VectorizedArrayType>
5348 Assert(this->normal_x_jacobian !=
nullptr,
5350 "update_gradients"));
5352 this->gradients_quad_submitted =
true;
5355 const std::size_t
nqp_d = this->n_quadrature_points * dim;
5360 const VectorizedArrayType
JxW_jac = this->J_value[0] *
5361 this->quadrature_weights[
q_point] *
5362 this->normal_x_jacobian[0][dim - 1];
5363 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5365 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5366 gradients[
comp *
nqp_d + d] = VectorizedArrayType();
5367 if constexpr (n_components == 1)
5375 const unsigned int index =
5378 this->normal_x_jacobian[
index];
5379 const VectorizedArrayType JxW =
5381 this->J_value[
index] * this->quadrature_weights[
q_point] :
5383 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5384 for (
unsigned int d = 0;
d < dim; ++
d)
5385 if constexpr (n_components == 1)
5398 typename VectorizedArrayType>
5407 Assert(this->J_value !=
nullptr,
5409 "update_hessians"));
5410 Assert(this->jacobian !=
nullptr,
5412 "update_hessians"));
5414 this->hessians_quad_submitted =
true;
5418 const std::size_t
nqp = this->n_quadrature_points;
5419 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5422 const VectorizedArrayType JxW =
5423 this->J_value[0] * this->quadrature_weights[
q_point];
5426 for (
unsigned int d = 0;
d < dim; ++
d)
5428 const auto jac_d = this->jacobian[0][
d][
d];
5429 const VectorizedArrayType factor =
jac_d *
jac_d * JxW;
5430 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5431 if constexpr (n_components == 1)
5440 for (
unsigned int d = 1,
off_dia = dim;
d < dim; ++
d)
5441 for (
unsigned int e = 0;
e <
d; ++
e, ++
off_dia)
5443 const auto jac_d = this->jacobian[0][
d][
d];
5444 const auto jac_e = this->jacobian[0][
e][
e];
5445 const VectorizedArrayType factor =
jac_d *
jac_e * JxW;
5446 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5447 if constexpr (n_components == 1)
5459 const VectorizedArrayType JxW =
5460 this->J_value[0] * this->quadrature_weights[
q_point];
5461 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5464 if constexpr (n_components == 1)
5470 VectorizedArrayType tmp[dim][dim];
5471 for (
unsigned int i = 0; i < dim; ++i)
5472 for (
unsigned int j = 0;
j < dim; ++
j)
5475 for (
unsigned int k = 1;
k < dim; ++
k)
5480 VectorizedArrayType
tmp2[dim][dim];
5481 for (
unsigned int i = 0; i < dim; ++i)
5482 for (
unsigned int j = 0;
j < dim; ++
j)
5485 for (
unsigned int k = 1;
k < dim; ++
k)
5490 for (
unsigned int d = 0;
d < dim; ++
d)
5495 for (
unsigned int d = 0,
off_diag = dim;
d < dim; ++
d)
5496 for (
unsigned int e = d + 1;
e < dim; ++
e, ++
off_diag)
5504 const VectorizedArrayType JxW = this->J_value[
q_point];
5506 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5509 if constexpr (n_components == 1)
5515 VectorizedArrayType tmp[dim][dim];
5516 for (
unsigned int i = 0; i < dim; ++i)
5517 for (
unsigned int j = 0;
j < dim; ++
j)
5520 for (
unsigned int k = 1;
k < dim; ++
k)
5525 VectorizedArrayType
tmp2[dim][dim];
5526 for (
unsigned int i = 0; i < dim; ++i)
5527 for (
unsigned int j = 0;
j < dim; ++
j)
5530 for (
unsigned int k = 1;
k < dim; ++
k)
5535 for (
unsigned int d = 0;
d < dim; ++
d)
5540 for (
unsigned int d = 0,
off_diag = dim;
d < dim; ++
d)
5541 for (
unsigned int e = d + 1;
e < dim; ++
e, ++
off_diag)
5546 for (
unsigned int d = 0;
d < dim; ++
d)
5548 VectorizedArrayType
sum = 0;
5549 for (
unsigned int e = 0;
e < dim; ++
e)
5551 for (
unsigned int e = 0,
count = dim;
e < dim; ++
e)
5552 for (
unsigned int f = e + 1; f < dim; ++f, ++
count)
5555 this->gradients_from_hessians_quad[(
comp *
nqp +
q_point) * dim +
5568 typename VectorizedArrayType>
5578 Assert(this->J_value !=
nullptr,
5580 "update_hessians"));
5581 Assert(this->jacobian !=
nullptr,
5583 "update_hessians"));
5585 this->hessians_quad_submitted =
true;
5589 const std::size_t
nqp = this->n_quadrature_points;
5590 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5593 const VectorizedArrayType JxW =
5594 this->J_value[0] * this->quadrature_weights[
q_point];
5596 const auto nxj = this->normal_x_jacobian[0];
5599 for (
unsigned int d = 0;
d < dim; ++
d)
5602 const VectorizedArrayType factor =
nxj_d *
nxj_d * JxW;
5603 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5604 if constexpr (n_components == 1)
5613 for (
unsigned int d = 1,
off_dia = dim;
d < dim; ++
d)
5614 for (
unsigned int e = 0;
e <
d; ++
e, ++
off_dia)
5618 const VectorizedArrayType factor =
jac_d *
jac_e * JxW;
5619 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5620 if constexpr (n_components == 1)
5630 const auto normal = this->normal_vector(
q_point);
5632 if constexpr (n_components == 1)
5637 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5650 typename VectorizedArrayType>
5655 VectorizedArrayType>::value_type
5661 Assert(this->values_quad_submitted ==
true,
5666 const std::size_t
nqp = this->n_quadrature_points;
5667 for (
unsigned int q = 0;
q <
nqp; ++
q)
5668 for (
unsigned int comp = 0;
comp < n_components; ++
comp)
5669 return_value[
comp] += this->values_quad[
comp *
nqp +
q];
5670 if constexpr (n_components == 1)
5671 return return_value[0];
5673 return return_value;
5682 typename VectorizedArrayType>
5683template <
int,
typename>
5688 static_assert(n_components == dim,
5689 "Do not try to modify the default template parameters used for"
5690 " selectively enabling this function via std::enable_if!");
5693 Assert(this->gradients_quad_initialized ==
true,
5697 Assert(this->jacobian !=
nullptr,
5699 "update_gradients"));
5701 VectorizedArrayType divergence;
5702 const std::size_t
nqp = this->n_quadrature_points;
5705 this->
data->element_type ==
5711 this->jacobian[0][0][0] *
5712 ((dim == 2) ? this->jacobian[0][1][1] :
5713 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
5721 if (is_face && dim == 2 && this->get_face_no() < 2)
5725 divergence = this->gradients_quad[
q_point * dim];
5726 for (
unsigned int d = 1;
d < dim; ++
d)
5727 divergence += this->gradients_quad[(d *
nqp +
q_point) * dim +
d];
5737 this->gradients_quad[
q_point * dim] * this->jacobian[0][0][0];
5738 for (
unsigned int d = 1;
d < dim; ++
d)
5739 divergence += this->gradients_quad[(d *
nqp +
q_point) * dim +
d] *
5740 this->jacobian[0][
d][
d];
5749 divergence =
jac[0][0] * this->gradients_quad[
q_point * dim];
5750 for (
unsigned int e = 1;
e < dim; ++
e)
5751 divergence +=
jac[0][e] * this->gradients_quad[
q_point * dim + e];
5752 for (
unsigned int d = 1;
d < dim; ++
d)
5753 for (
unsigned int e = 0;
e < dim; ++
e)
5767 typename VectorizedArrayType>
5768template <
int,
typename>
5773 static_assert(n_components == dim,
5774 "Do not try to modify the default template parameters used for"
5775 " selectively enabling this function via std::enable_if!");
5779 VectorizedArrayType
symmetrized[(dim * dim + dim) / 2];
5780 VectorizedArrayType
half = Number(0.5);
5781 for (
unsigned int d = 0;
d < dim; ++
d)
5811 typename VectorizedArrayType>
5812template <
int,
typename>
5814 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
5818 static_assert(dim > 1 && n_components == dim,
5819 "Do not try to modify the default template parameters used for"
5820 " selectively enabling this function via std::enable_if!");
5824 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
5847 typename VectorizedArrayType>
5848template <
int,
typename>
5854 static_assert(n_components == dim,
5855 "Do not try to modify the default template parameters used for"
5856 " selectively enabling this function via std::enable_if!");
5862 Assert(this->J_value !=
nullptr,
5864 "update_gradients"));
5865 Assert(this->jacobian !=
nullptr,
5867 "update_gradients"));
5869 this->gradients_quad_submitted =
true;
5872 const std::size_t
nqp_d = this->n_quadrature_points * dim;
5875 if (this->
data->element_type ==
5882 const VectorizedArrayType
fac =
5894 Number((dim == 2 &&
this->get_face_no() < 2) ? -1 : 1);
5896 for (
unsigned int d = 0;
d < dim; ++
d)
5898 for (
unsigned int e = 0;
e < dim; ++
e)
5899 gradients[d *
nqp_d + e] = (d == e) ?
fac : 0.;
5901 this->divergence_is_requested =
true;
5908 const VectorizedArrayType
fac =
5909 this->J_value[0] * this->quadrature_weights[
q_point] *
div_in;
5910 for (
unsigned int d = 0;
d < dim; ++
d)
5912 const VectorizedArrayType
jac_dd = this->jacobian[0][
d][
d];
5913 for (
unsigned int e = 0;
e < dim; ++
e)
5923 const VectorizedArrayType
fac =
5926 this->J_value[0] * this->quadrature_weights[
q_point]) *
5928 for (
unsigned int d = 0;
d < dim; ++
d)
5930 for (
unsigned int e = 0;
e < dim; ++
e)
5943 typename VectorizedArrayType>
5944template <
int,
typename>
5951 static_assert(n_components == dim,
5952 "Do not try to modify the default template parameters used for"
5953 " selectively enabling this function via std::enable_if!");
5956 this->
data->element_type !=
5967 Assert(this->J_value !=
nullptr,
5969 "update_gradients"));
5970 Assert(this->jacobian !=
nullptr,
5972 "update_gradients"));
5974 this->gradients_quad_submitted =
true;
5977 const std::size_t
nqp_d = this->n_quadrature_points * dim;
5981 const VectorizedArrayType JxW =
5982 this->J_value[0] * this->quadrature_weights[
q_point];
5984 for (
unsigned int d = 0;
d < dim; ++
d)
5985 gradients[d *
nqp_d + d] =
5987 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
5988 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
5990 const VectorizedArrayType
value =
5991 sym_grad.access_raw_entry(counter) * JxW;
5999 const VectorizedArrayType JxW =
6002 this->J_value[0] * this->quadrature_weights[
q_point];
6007 VectorizedArrayType
weighted[dim][dim];
6008 for (
unsigned int i = 0; i < dim; ++i)
6010 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6011 for (
unsigned int j = i + 1;
j < dim; ++
j, ++counter)
6013 const VectorizedArrayType
value =
6014 sym_grad.access_raw_entry(counter) * JxW;
6019 for (
unsigned int d = 0;
d < dim; ++
d)
6022 for (
unsigned int e = 1;
e < dim; ++
e)
6035 typename VectorizedArrayType>
6036template <
int,
typename>
6042 static_assert(n_components == dim,
6043 "Do not try to modify the default template parameters used for"
6044 " selectively enabling this function via std::enable_if!");
6050 grad[1][0] = curl[0];
6051 grad[0][1] = -curl[0];
6054 grad[2][1] = curl[0];
6055 grad[1][2] = -curl[0];
6056 grad[0][2] = curl[1];
6057 grad[2][0] = -curl[1];
6058 grad[1][0] = curl[2];
6059 grad[0][1] = -curl[2];
6077 typename VectorizedArrayType>
6083 VectorizedArrayType>
::
6085 const unsigned int fe_no,
6086 const unsigned int quad_no,
6087 const unsigned int first_selected_component,
6088 const unsigned int active_fe_index,
6089 const unsigned int active_quad_index)
6090 : BaseClass(matrix_free,
6092 first_selected_component,
6100 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
6102 , n_q_points(
this->
data->n_q_points)
6104 check_template_arguments(
fe_no, 0);
6114 typename VectorizedArrayType>
6120 VectorizedArrayType>
::
6122 const std::pair<unsigned int, unsigned int> &
range,
6123 const unsigned int dof_no,
6124 const unsigned int quad_no,
6125 const unsigned int first_selected_component)
6129 first_selected_component,
6130 matrix_free.get_cell_active_fe_index(
range))
6140 typename VectorizedArrayType>
6146 VectorizedArrayType>
::
6151 const unsigned int first_selected_component)
6152 : BaseClass(mapping,
6156 first_selected_component,
6158 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
6160 , n_q_points(
this->
data->n_q_points)
6172 typename VectorizedArrayType>
6178 VectorizedArrayType>
::
6182 const unsigned int first_selected_component)
6187 first_selected_component,
6189 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
6191 , n_q_points(
this->
data->n_q_points)
6203 typename VectorizedArrayType>
6209 VectorizedArrayType>
::
6212 const unsigned int first_selected_component)
6213 : BaseClass(
other.mapped_geometry->get_fe_values().get_mapping(),
6215 other.mapped_geometry->get_quadrature(),
6216 other.mapped_geometry->get_fe_values().get_update_flags(),
6217 first_selected_component,
6219 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
6221 , n_q_points(
this->
data->n_q_points)
6233 typename VectorizedArrayType>
6242 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
6244 , n_q_points(
this->
data->n_q_points)
6256 typename VectorizedArrayType>
6262 VectorizedArrayType> &
6270 BaseClass::operator=(
other);
6282 typename VectorizedArrayType>
6289 VectorizedArrayType>
::
6291 const unsigned int first_selected_component)
6294 (
void)first_selected_component;
6297 this->
data->dofs_per_component_on_cell > 0,
6299 "There is nothing useful you can do with an FEEvaluation object with "
6300 "FE_Nothing, i.e., without DoFs! If you have passed to "
6301 "MatrixFree::reinit() a collection of finite elements also containing "
6302 "FE_Nothing, please check - before creating FEEvaluation - the category "
6303 "of the current range by calling either "
6304 "MatrixFree::get_cell_range_category(range) or "
6305 "MatrixFree::get_face_range_category(range). The returned category "
6306 "is the index of the active FE, which you can use to exclude "
6313 static_cast<unsigned int>(fe_degree) !=
6314 this->
data->data.front().fe_degree) ||
6315 n_q_points !=
this->n_quadrature_points)
6318 "-------------------------------------------------------\n";
6319 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
6320 message +=
" Called --> FEEvaluation<dim,";
6340 if (
static_cast<unsigned int>(fe_degree) ==
6341 this->
data->data.front().fe_degree)
6343 proposed_dof_comp =
dof_no;
6347 for (
unsigned int no = 0;
no < this->matrix_free->n_components();
6349 for (
unsigned int nf = 0;
6350 nf < this->matrix_free->n_base_elements(
no);
6352 if (this->matrix_free
6353 ->get_shape_info(
no, 0,
nf, this->active_fe_index, 0)
6355 .fe_degree ==
static_cast<unsigned int>(fe_degree))
6362 this->mapping_data->descriptor[
this->active_quad_index]
6366 for (
unsigned int no = 0;
6367 no < this->matrix_free->get_mapping_info().cell_data.size();
6369 if (this->matrix_free->get_mapping_info()
6371 .descriptor[this->active_quad_index]
6372 .n_q_points == n_q_points)
6382 message +=
"Wrong vector component selection:\n";
6384 message +=
"Wrong quadrature formula selection:\n";
6385 message +=
" Did you mean FEEvaluation<dim,";
6417 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6418 message +=
"Wrong template arguments:\n";
6419 message +=
" Did you mean FEEvaluation<dim,";
6433 if (this->
data->data.front().fe_degree !=
6434 static_cast<unsigned int>(fe_degree))
6444 Assert(
static_cast<unsigned int>(fe_degree) ==
6445 this->
data->data.front().fe_degree &&
6446 n_q_points ==
this->n_quadrature_points,
6452 this->mapping_data->descriptor[
this->active_quad_index].n_q_points);
6463 typename VectorizedArrayType>
6472 Assert(this->matrix_free !=
nullptr,
6473 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
6474 " Integer indexing is not possible."));
6480 this->matrix_free->get_mapping_info().get_cell_type(
cell_index);
6483 this->mapping_data->data_index_offsets[
cell_index];
6484 this->jacobian = &this->mapping_data->jacobians[0][
offsets];
6485 this->J_value = &this->mapping_data->JxW_values[
offsets];
6486 if (!this->mapping_data->jacobian_gradients[0].empty())
6488 this->jacobian_gradients =
6489 this->mapping_data->jacobian_gradients[0].data() +
offsets;
6490 this->jacobian_gradients_non_inverse =
6491 this->mapping_data->jacobian_gradients_non_inverse[0].data() +
offsets;
6494 if (this->matrix_free->n_active_entries_per_cell_batch(
this->cell) == n_lanes)
6497 for (
unsigned int i = 0; i < n_lanes; ++i)
6498 this->cell_ids[i] =
cell_index * n_lanes + i;
6503 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
6505 this->cell_ids[i] =
cell_index * n_lanes + i;
6506 for (; i < n_lanes; ++i)
6510 if (this->mapping_data->quadrature_points.empty() ==
false)
6512 &this->mapping_data->quadrature_points
6513 [this->mapping_data->quadrature_point_offsets[this->cell]];
6516 this->is_reinitialized =
true;
6517 this->dof_values_initialized =
false;
6518 this->values_quad_initialized =
false;
6519 this->gradients_quad_initialized =
false;
6520 this->hessians_quad_initialized =
false;
6531 typename VectorizedArrayType>
6538 VectorizedArrayType>
::reinit(
const std::array<
unsigned int,
6545 this->cell_ids = cell_ids;
6550 for (
unsigned int v = 0; v < n_lanes; ++v)
6559 this->matrix_free->get_mapping_info().get_cell_type(
6564 if (this->mapped_geometry ==
nullptr)
6565 this->mapped_geometry =
6566 std::make_shared<internal::MatrixFreeFunctions::
6567 MappingDataOnTheFly<dim, VectorizedArrayType>>();
6586 const auto &update_flags_cells =
6587 this->matrix_free->get_mapping_info().update_flags_cells;
6608 const auto &update_flags_cells =
6609 this->matrix_free->get_mapping_info().update_flags_cells;
6616 this->n_quadrature_points);
6628 this->jacobian_gradients_non_inverse =
6633 for (
unsigned int v = 0; v < n_lanes; ++v)
6643 const unsigned int lane =
cell_index % n_lanes;
6645 if (this->cell_type <=
6649 for (
unsigned int q = 0;
q < 2; ++
q)
6650 for (
unsigned int i = 0; i < dim; ++i)
6651 for (
unsigned int j = 0;
j < dim; ++
j)
6653 this->mapping_data->jacobians[0][
offsets +
q][i][
j][lane];
6655 const unsigned int q = 0;
6658 this->mapping_data->JxW_values[
offsets +
q][lane];
6660 const auto &update_flags_cells =
6661 this->matrix_free->get_mapping_info().update_flags_cells;
6665 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6666 for (
unsigned int j = 0;
j < dim; ++
j)
6669 ->jacobian_gradients[0][
offsets +
q][i][
j][lane];
6671 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6672 for (
unsigned int j = 0;
j < dim; ++
j)
6675 ->jacobian_gradients_non_inverse[0][
offsets +
q][i][
j]
6680 for (
unsigned int i = 0; i < dim; ++i)
6682 this->mapping_data->quadrature_points
6690 const auto cell_type =
6691 this->matrix_free->get_mapping_info().get_cell_type(
6694 for (
unsigned int q = 0;
q < this->n_quadrature_points; ++
q)
6696 const unsigned int q_src =
6705 for (
unsigned int i = 0; i < dim; ++i)
6706 for (
unsigned int j = 0;
j < dim; ++
j)
6711 const auto &update_flags_cells =
6712 this->matrix_free->get_mapping_info().update_flags_cells;
6716 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6717 for (
unsigned int j = 0;
j < dim; ++
j)
6722 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
6723 for (
unsigned int j = 0;
j < dim; ++
j)
6739 this->mapping_data->quadrature_points
6745 this->mapping_data->jacobians[0][
offsets + 1];
6747 for (
unsigned int d = 0;
d < dim; ++
d)
6750 static_cast<Number
>(
6751 this->descriptor->quadrature.point(
q)[
d]);
6753 for (
unsigned int d = 0;
d < dim; ++
d)
6754 for (
unsigned int e = 0;
e < dim; ++
e)
6757 static_cast<Number
>(
6758 this->descriptor->quadrature.point(
q)[
e]);
6760 for (
unsigned int i = 0; i < dim; ++i)
6766 for (
unsigned int i = 0; i < dim; ++i)
6768 this->mapping_data->quadrature_points
6779 this->is_reinitialized =
true;
6780 this->dof_values_initialized =
false;
6781 this->values_quad_initialized =
false;
6782 this->gradients_quad_initialized =
false;
6783 this->hessians_quad_initialized =
false;
6794 typename VectorizedArrayType>
6795template <
bool level_dof_access>
6802 VectorizedArrayType>
::
6805 Assert(this->matrix_free ==
nullptr,
6806 ExcMessage(
"Cannot use initialization from cell iterator if "
6807 "initialized from MatrixFree object. Use variant for "
6808 "on the fly computation with arguments as for FEValues "
6811 this->mapped_geometry->reinit(
6813 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6815 cell->get_mg_dof_indices(this->local_dof_indices);
6817 cell->get_dof_indices(this->local_dof_indices);
6820 this->is_reinitialized =
true;
6831 typename VectorizedArrayType>
6838 VectorizedArrayType>
::
6841 Assert(this->matrix_free == 0,
6842 ExcMessage(
"Cannot use initialization from cell iterator if "
6843 "initialized from MatrixFree object. Use variant for "
6844 "on the fly computation with arguments as for FEValues "
6847 this->mapped_geometry->reinit(cell);
6850 this->is_reinitialized =
true;
6861 typename VectorizedArrayType>
6868 VectorizedArrayType>
::
6872 Assert(this->dof_values_initialized ==
true,
6875 evaluate(this->values_dofs, evaluation_flags);
6885 typename VectorizedArrayType>
6892 VectorizedArrayType>
::
6903 if (this->
data->element_type ==
6909 if constexpr (fe_degree > -1)
6927 this->values_quad_initialized =
6929 this->gradients_quad_initialized =
6931 this->hessians_quad_initialized =
6942 template <
typename Number,
6943 typename VectorizedArrayType,
6946 std::enable_if_t<internal::has_begin<VectorType> &&
6949 VectorizedArrayType *
6956 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
6957 const auto &dof_info = fe_eval.get_dof_info();
6964 if (std::is_same_v<typename VectorType::value_type, Number> &&
6968 interleaved_contiguous &&
6971 dof_info.dof_indices_contiguous
6972 [
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
6973 [cell * VectorizedArrayType::
size()]) %
6974 sizeof(VectorizedArrayType) ==
6977 return reinterpret_cast<VectorizedArrayType *
>(
6981 [cell * VectorizedArrayType::size()] +
6983 [fe_eval.get_active_fe_index()]
6984 [fe_eval.get_first_selected_component()] *
6985 VectorizedArrayType::size());
6994 template <
typename Number,
6995 typename VectorizedArrayType,
6998 std::enable_if_t<!internal::has_begin<VectorType> ||
7001 VectorizedArrayType *
7015 typename VectorizedArrayType>
7016template <
typename VectorType>
7023 VectorizedArrayType>
::
7027 const VectorizedArrayType *
src_ptr =
7028 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7046 typename VectorizedArrayType>
7053 VectorizedArrayType>
::
7059 this->dof_values_initialized =
true;
7070 typename VectorizedArrayType>
7077 VectorizedArrayType>
::
7084 Assert(this->values_quad_submitted ==
true,
7087 Assert(this->gradients_quad_submitted ==
true,
7090 Assert(this->hessians_quad_submitted ==
true,
7093 Assert(this->matrix_free !=
nullptr ||
7094 this->mapped_geometry->is_initialized(),
7100 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7101 "EvaluationFlags::hessians are supported."));
7107 unsigned int size = n_components * dim * n_q_points;
7110 for (
unsigned int i = 0; i <
size; ++i)
7111 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7115 for (
unsigned int i = 0; i <
size; ++i)
7116 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7121 if (n_components == dim &&
7122 this->
data->element_type ==
7126 this->divergence_is_requested ==
false)
7128 unsigned int size = n_components * n_q_points;
7131 for (
unsigned int i = 0; i <
size; ++i)
7132 this->values_quad[i] += this->values_from_gradients_quad[i];
7136 for (
unsigned int i = 0; i <
size; ++i)
7137 this->values_quad[i] = this->values_from_gradients_quad[i];
7142 if constexpr (fe_degree > -1)
7162 this->dof_values_initialized =
true;
7173 typename VectorizedArrayType>
7174template <
typename VectorType>
7181 VectorizedArrayType>
::
7185 VectorizedArrayType *
dst_ptr =
7186 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7204 typename VectorizedArrayType>
7211 VectorizedArrayType>::dof_indices()
const
7213 return {0
U, dofs_per_cell};
7227 typename VectorizedArrayType>
7233 VectorizedArrayType>
::
7236 const bool is_interior_face,
7237 const unsigned int dof_no,
7238 const unsigned int quad_no,
7239 const unsigned int first_selected_component,
7240 const unsigned int active_fe_index,
7241 const unsigned int active_quad_index,
7242 const unsigned int face_type)
7243 : BaseClass(matrix_free,
7245 first_selected_component,
7253 , dofs_per_component(
this->
data->dofs_per_component_on_cell)
7255 , n_q_points(
this->n_quadrature_points)
7265 typename VectorizedArrayType>
7271 VectorizedArrayType>
::
7274 const std::pair<unsigned int, unsigned int> &
range,
7275 const bool is_interior_face,
7276 const unsigned int dof_no,
7277 const unsigned int quad_no,
7278 const unsigned int first_selected_component)
7283 first_selected_component,
7284 matrix_free.get_face_active_fe_index(
range,
7287 matrix_free.get_face_info(
range.
first).face_type)
7297 typename VectorizedArrayType>
7304 VectorizedArrayType>
::reinit(
const unsigned int face_index)
7306 Assert(this->mapped_geometry ==
nullptr,
7307 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7308 " Integer indexing is not possible"));
7309 if (this->mapped_geometry !=
nullptr)
7312 this->cell = face_index;
7313 this->dof_access_index =
7314 this->is_interior_face() ?
7320 this->matrix_free->get_task_info().face_partition_data.back() &&
7322 this->matrix_free->get_task_info().boundary_partition_data.back())
7323 Assert(this->is_interior_face(),
7325 "Boundary faces do not have a neighbor. When looping over "
7326 "boundary faces use FEFaceEvaluation with the parameter "
7327 "is_interior_face set to true. "));
7329 this->reinit_face(this->matrix_free->get_face_info(face_index));
7332 for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
7334 this->face_ids[i] = face_index * n_lanes + i;
7335 for (; i < n_lanes; ++i)
7338 this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
7340 this->mapping_data->data_index_offsets[face_index];
7341 this->J_value = &this->mapping_data->JxW_values[
offsets];
7342 this->normal_vectors = &this->mapping_data->normal_vectors[
offsets];
7344 &this->mapping_data->jacobians[!this->is_interior_face()][
offsets];
7345 this->normal_x_jacobian =
7347 ->normals_times_jacobians[!this->is_interior_face()][
offsets];
7348 this->jacobian_gradients =
7349 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7351 this->jacobian_gradients_non_inverse =
7353 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7357 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
7360 this->mapping_data->quadrature_point_offsets.size());
7362 this->mapping_data->quadrature_points.data() +
7363 this->mapping_data->quadrature_point_offsets[this->cell];
7367 this->is_reinitialized =
true;
7368 this->dof_values_initialized =
false;
7369 this->values_quad_initialized =
false;
7370 this->gradients_quad_initialized =
false;
7371 this->hessians_quad_initialized =
false;
7382 typename VectorizedArrayType>
7390 const unsigned int face_number)
7394 this->matrix_free->get_mapping_info().face_data_by_cells.size(),
7396 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7399 this->matrix_free->get_mapping_info().cell_type.size());
7400 Assert(this->mapped_geometry ==
nullptr,
7401 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7402 " Integer indexing is not possible"));
7403 if (this->mapped_geometry !=
nullptr)
7407 this->cell_type = this->matrix_free->get_mapping_info()
7408 .faces_by_cells_type[
cell_index][face_number];
7411 this->dof_access_index =
7414 if (this->is_interior_face() ==
false)
7419 for (
unsigned int i = 0; i < n_lanes; ++i)
7424 unsigned int face_index =
7425 this->matrix_free->get_cell_and_face_to_plain_faces()(
cell_index,
7429 this->face_ids[i] = face_index;
7434 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
7435 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
7440 this->matrix_free->get_face_info(face_index / n_lanes);
7442 auto cell_m = faces.cells_interior[face_index % n_lanes];
7443 auto cell_p = faces.cells_exterior[face_index % n_lanes];
7453 this->cell_ids[i] =
cell_m;
7454 this->face_numbers[i] = faces.interior_face_no;
7458 this->cell_ids[i] =
cell_p;
7459 this->face_numbers[i] = faces.exterior_face_no;
7463 unsigned int face_orientation = faces.face_orientation % 8;
7466 constexpr std::array<std::uint8_t, 8> table{
7467 {0, 1, 6, 3, 4, 5, 2, 7}};
7468 face_orientation = table[face_orientation];
7470 this->face_orientations[i] = face_orientation;
7475 this->face_orientations[0] = 0;
7476 this->face_numbers[0] = face_number;
7477 if (this->matrix_free->n_active_entries_per_cell_batch(
this->cell) ==
7481 for (
unsigned int i = 0; i < n_lanes; ++i)
7482 this->cell_ids[i] =
cell_index * n_lanes + i;
7488 this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7490 this->cell_ids[i] =
cell_index * n_lanes + i;
7491 for (; i < n_lanes; ++i)
7494 for (
unsigned int i = 0; i < n_lanes; ++i)
7496 this->matrix_free->get_cell_and_face_to_plain_faces()(
cell_index,
7502 this->matrix_free->get_mapping_info()
7503 .face_data_by_cells[this->quad_no]
7507 this->matrix_free->get_mapping_info()
7508 .face_data_by_cells[
this->quad_no]
7509 .JxW_values.size());
7510 this->J_value = &this->matrix_free->get_mapping_info()
7511 .face_data_by_cells[this->quad_no]
7513 this->normal_vectors = &this->matrix_free->get_mapping_info()
7514 .face_data_by_cells[this->quad_no]
7516 this->jacobian = &this->matrix_free->get_mapping_info()
7517 .face_data_by_cells[this->quad_no]
7518 .jacobians[!this->is_interior_face()][
offsets];
7519 this->normal_x_jacobian =
7520 &this->matrix_free->get_mapping_info()
7521 .face_data_by_cells[this->quad_no]
7522 .normals_times_jacobians[!this->is_interior_face()][
offsets];
7523 this->jacobian_gradients =
7524 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
7526 this->jacobian_gradients_non_inverse =
7528 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
7532 if (this->matrix_free->get_mapping_info()
7533 .face_data_by_cells[
this->quad_no]
7534 .quadrature_point_offsets.empty() ==
false)
7536 const unsigned int index =
7539 this->matrix_free->get_mapping_info()
7540 .face_data_by_cells[
this->quad_no]
7541 .quadrature_point_offsets.size());
7543 .face_data_by_cells[this->quad_no]
7544 .quadrature_points.data() +
7545 this->matrix_free->get_mapping_info()
7546 .face_data_by_cells[this->quad_no]
7547 .quadrature_point_offsets[
index];
7551 this->is_reinitialized =
true;
7552 this->dof_values_initialized =
false;
7553 this->values_quad_initialized =
false;
7554 this->gradients_quad_initialized =
false;
7555 this->hessians_quad_initialized =
false;
7566 typename VectorizedArrayType>
7573 VectorizedArrayType>
::
7590 typename VectorizedArrayType>
7597 VectorizedArrayType>
::
7604 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7605 "and EvaluationFlags::hessians are supported."));
7614 if (this->
data->element_type ==
7620 if constexpr (fe_degree > -1)
7622 template run<fe_degree, n_q_points_1d>(n_components,
7631 this->values_quad_initialized =
7633 this->gradients_quad_initialized =
7635 this->hessians_quad_initialized =
7647 typename VectorizedArrayType>
7654 VectorizedArrayType>
::
7671 typename VectorizedArrayType>
7678 VectorizedArrayType>
::
7685 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7686 "and EvaluationFlags::hessians are supported."));
7695 if (this->
data->element_type ==
7701 if constexpr (fe_degree > -1)
7725 typename VectorizedArrayType>
7732 VectorizedArrayType>
::
7738 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7739 "and EvaluationFlags::hessians are supported."));
7748 if (this->
data->element_type ==
7754 if constexpr (fe_degree > -1)
7765 this->values_quad_initialized =
7767 this->gradients_quad_initialized =
7769 this->hessians_quad_initialized =
7781 typename VectorizedArrayType>
7788 VectorizedArrayType>
::
7795 this->dof_values_initialized =
true;
7806 typename VectorizedArrayType>
7813 VectorizedArrayType>
::
7821 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7822 "and EvaluationFlags::hessians are supported."));
7828 unsigned int size = n_components * dim * n_q_points;
7831 for (
unsigned int i = 0; i <
size; ++i)
7832 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7836 for (
unsigned int i = 0; i <
size; ++i)
7837 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7842 if (this->
data->element_type ==
7846 this->divergence_is_requested ==
false)
7848 unsigned int size = n_components * n_q_points;
7851 for (
unsigned int i = 0; i <
size; ++i)
7852 this->values_quad[i] += this->values_from_gradients_quad[i];
7856 for (
unsigned int i = 0; i <
size; ++i)
7857 this->values_quad[i] = this->values_from_gradients_quad[i];
7862 if constexpr (fe_degree > -1)
7864 template run<fe_degree, n_q_points_1d>(n_components,
7885 typename VectorizedArrayType>
7892 VectorizedArrayType>
::
7898 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
7899 "and EvaluationFlags::hessians are supported."));
7905 unsigned int size = n_components * dim * n_q_points;
7908 for (
unsigned int i = 0; i <
size; ++i)
7909 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7913 for (
unsigned int i = 0; i <
size; ++i)
7914 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7919 if (this->
data->element_type ==
7923 this->divergence_is_requested ==
false)
7925 unsigned int size = n_components * n_q_points;
7928 for (
unsigned int i = 0; i <
size; ++i)
7929 this->values_quad[i] += this->values_from_gradients_quad[i];
7933 for (
unsigned int i = 0; i <
size; ++i)
7934 this->values_quad[i] = this->values_from_gradients_quad[i];
7939 if constexpr (fe_degree > -1)
7959 typename VectorizedArrayType>
7966 VectorizedArrayType>
::
7973 this->dof_values_initialized =
true;
7984 typename VectorizedArrayType>
7991 VectorizedArrayType>
::
7999 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8000 "and EvaluationFlags::hessians are supported."));
8007 if (this->
data->element_type ==
8011 this->divergence_is_requested ==
false)
8014 if constexpr (fe_degree > -1)
8038 typename VectorizedArrayType>
8039template <
typename VectorType>
8046 VectorizedArrayType>
::
8053 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8054 "and EvaluationFlags::hessians are supported."));
8056 const auto shared_vector_data = internal::get_shared_vector_data(
8058 this->dof_access_index ==
8060 this->active_fe_index,
8063 if (this->
data->data.front().fe_degree > 0 &&
8064 fast_evaluation_supported(this->
data->data.front().fe_degree,
8065 this->data->data.front().n_q_points_1d) &&
8068 typename VectorType::value_type,
8069 VectorizedArrayType>::
8074 this->dof_info->index_storage_variants[
this->dof_access_index]
8077 if constexpr (fe_degree > -1)
8081 typename VectorType::value_type,
8086 internal::get_beginning<typename VectorType::value_type>(
8095 typename VectorType::value_type,
8096 VectorizedArrayType>::evaluate(n_components,
8098 internal::get_beginning<
8099 typename VectorType::value_type>(
8113 this->gradients_quad_initialized =
8126 typename VectorizedArrayType>
8127template <
typename VectorType>
8135 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8136 const bool integrate_gradients,
8154 typename VectorizedArrayType>
8155template <
typename VectorType>
8162 VectorizedArrayType>
::
8166 Assert((this->dof_access_index ==
8168 this->is_interior_face() ==
false) ==
false,
8171 const auto shared_vector_data = internal::get_shared_vector_data(
8173 this->dof_access_index ==
8175 this->active_fe_index,
8178 if (this->
data->data.front().fe_degree > 0 &&
8179 fast_evaluation_supported(this->
data->data.front().fe_degree,
8180 this->data->data.front().n_q_points_1d) &&
8183 typename VectorType::value_type,
8184 VectorizedArrayType>::
8189 this->dof_info->index_storage_variants[
this->dof_access_index]
8192 if constexpr (fe_degree > -1)
8196 typename VectorType::value_type,
8201 internal::get_beginning<typename VectorType::value_type>(
8210 typename VectorType::value_type,
8211 VectorizedArrayType>::integrate(n_components,
8213 internal::get_beginning<
8214 typename VectorType::value_type>(
8234 typename VectorizedArrayType>
8241 VectorizedArrayType>::dof_indices()
const
8243 return {0
U, dofs_per_cell};
8253 typename VectorizedArrayType>
8260 VectorizedArrayType>
::
8264 return fe_degree == -1 ?
8277 typename VectorizedArrayType>
8284 VectorizedArrayType>
::
8288 return fe_degree == -1 ?
8301 typename VectorizedArrayType>
8308 VectorizedArrayType>::at_boundary()
const
8310 Assert(this->dof_access_index !=
8314 if (this->is_interior_face() ==
false)
8316 else if (this->
cell < this->matrix_free->n_inner_face_batches())
8318 else if (this->cell < (this->matrix_free->n_inner_face_batches() +
8319 this->matrix_free->n_boundary_face_batches()))
8332 typename VectorizedArrayType>
8339 VectorizedArrayType>::boundary_id()
const
8341 Assert(this->dof_access_index !=
8346 return this->matrix_free->get_boundary_id(this->cell);
8358 typename VectorizedArrayType>
8366 VectorizedArrayType>::get_dofs_per_component_projected_to_face()
8368 return this->
data->dofs_per_component_on_face;
8378 typename VectorizedArrayType>
8385 VectorizedArrayType>::get_dofs_projected_to_face()
8387 return this->
data->dofs_per_component_on_face * n_components_;
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > val_in, const unsigned int q_point)
void submit_normal_hessian(const value_type normal_hessian_in, const unsigned int q_point)
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
std::conditional_t< n_components_==1, Tensor< 2, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 3, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > > > hessian_type
std::conditional_t< n_components_==1, Tensor< 1, dim, VectorizedArrayType >, std::conditional_t< n_components_==dim, Tensor< 2, dim, VectorizedArrayType >, Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > > > gradient_type
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
VectorizedArrayType get_divergence(const unsigned int q_point) const
hessian_type get_hessian(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
static constexpr unsigned int n_components
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
std::conditional_t< n_components_==1, VectorizedArrayType, Tensor< 1, n_components_, VectorizedArrayType > > value_type
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< n_lanes > &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
static constexpr unsigned int n_lanes
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
value_type get_normal_hessian(const unsigned int q_point) const
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip())
const unsigned int quad_no
const 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::conditional_t< rank_==1, Number, Tensor< rank_ - 1, dim, Number > > value_type
#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()
std::vector< index_type > data
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< Number > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate_in_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, FEEvaluationData< dim, Number, true > &fe_eval)
static void collect_from_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval, const bool sum_into_values)
static void project_to_face(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static constexpr unsigned int value
@ dof_access_face_exterior
@ dof_access_face_interior
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
std::vector< std::pair< unsigned int, unsigned int > > row_starts
std::vector< std::vector< unsigned int > > component_dof_indices_offset
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
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::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
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)