17#ifndef dealii_matrix_free_fe_evaluation_h
18#define dealii_matrix_free_fe_evaluation_h
93 typename VectorizedArrayType>
143 template <
typename VectorType>
146 const unsigned int first_index = 0,
147 const std::bitset<VectorizedArrayType::size()> &mask =
148 std::bitset<VectorizedArrayType::size()>().flip());
178 template <
typename VectorType>
181 const unsigned int first_index = 0,
182 const std::bitset<VectorizedArrayType::size()> &mask =
183 std::bitset<VectorizedArrayType::size()>().flip());
216 template <
typename VectorType>
220 const unsigned int first_index = 0,
221 const std::bitset<VectorizedArrayType::size()> &mask =
222 std::bitset<VectorizedArrayType::size()>().flip())
const;
262 template <
typename VectorType>
265 const unsigned int first_index = 0,
266 const std::bitset<VectorizedArrayType::size()> &mask =
267 std::bitset<VectorizedArrayType::size()>().flip())
const;
272 template <
typename VectorType>
276 const unsigned int first_index = 0,
277 const std::bitset<VectorizedArrayType::size()> &mask =
278 std::bitset<VectorizedArrayType::size()>().flip())
const;
413 const unsigned int q_point);
501 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
521 const unsigned int q_point);
542 const unsigned int q_point);
558 const unsigned int q_point);
600 const unsigned int dof_no,
603 const unsigned int fe_degree,
604 const unsigned int n_q_points,
608 const unsigned int face_type);
682 template <
typename VectorType,
typename VectorOperation>
686 const std::array<VectorType *, n_components_> &vectors,
689 n_components_> & vectors_sm,
690 const std::bitset<VectorizedArrayType::size()> &mask,
691 const bool apply_constraints =
true)
const;
700 template <
typename VectorType,
typename VectorOperation>
704 const std::array<VectorType *, n_components_> &vectors,
707 n_components_> & vectors_sm,
708 const std::bitset<VectorizedArrayType::size()> &mask)
const;
717 template <
typename VectorType,
typename VectorOperation>
721 const std::array<VectorType *, n_components_> &vectors)
const;
767 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
768 "Type of Number and of VectorizedArrayType do not match.");
790 const unsigned int dof_no,
793 const unsigned int fe_degree,
794 const unsigned int n_q_points,
834template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
839 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
840 "Type of Number and of VectorizedArrayType do not match.");
880 const unsigned int q_point);
905 const unsigned int q_point);
947 const unsigned int dof_no,
950 const unsigned int fe_degree,
951 const unsigned int n_q_points,
992template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
997 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
998 "Type of Number and of VectorizedArrayType do not match.");
1041 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1061 const unsigned int q_point);
1080 const unsigned int q_point);
1092 const unsigned int q_point);
1105 const unsigned int q_point);
1113 const unsigned int q_point);
1125 const unsigned int dof_no,
1128 const unsigned int dofs_per_cell,
1129 const unsigned int n_q_points,
1168template <
typename Number,
bool is_face,
typename VectorizedArrayType>
1173 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1174 "Type of Number and of VectorizedArrayType do not match.");
1250 const unsigned int q_point);
1257 const unsigned int q_point);
1264 const unsigned int q_point);
1306 const unsigned int dof_no,
1309 const unsigned int fe_degree,
1310 const unsigned int n_q_points,
1902 typename VectorizedArrayType>
1907 VectorizedArrayType>
1910 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1911 "Type of Number and of VectorizedArrayType do not match.");
2025 const unsigned int dof_no = 0,
2026 const unsigned int quad_no = 0,
2039 const std::pair<unsigned int, unsigned int> & range,
2040 const unsigned int dof_no = 0,
2041 const unsigned int quad_no = 0,
2151 template <
bool level_dof_access>
2173 const unsigned int give_n_q_points_1d);
2191 DEAL_II_DEPRECATED_EARLY
void
2193 const bool evaluate_gradients,
2194 const bool evaluate_hessians =
false);
2216 DEAL_II_DEPRECATED_EARLY
void
2218 const bool evaluate_values,
2219 const bool evaluate_gradients,
2220 const bool evaluate_hessians =
false);
2235 template <
typename VectorType>
2243 template <
typename VectorType>
2244 DEAL_II_DEPRECATED_EARLY
void
2246 const bool evaluate_values,
2247 const bool evaluate_gradients,
2248 const bool evaluate_hessians =
false);
2266 DEAL_II_DEPRECATED_EARLY
void
2267 integrate(
const bool integrate_values,
const bool integrate_gradients);
2282 VectorizedArrayType * values_array,
2283 const bool sum_into_values =
false);
2288 DEAL_II_DEPRECATED_EARLY
void
2290 const bool integrate_gradients,
2291 VectorizedArrayType *values_array);
2306 template <
typename VectorType>
2309 VectorType & output_vector);
2314 template <
typename VectorType>
2315 DEAL_II_DEPRECATED_EARLY
void
2317 const bool integrate_gradients,
2318 VectorType &output_vector);
2402 int n_q_points_1d = fe_degree + 1,
2403 int n_components_ = 1,
2404 typename Number = double,
2410 VectorizedArrayType>
2413 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2414 "Type of Number and of VectorizedArrayType do not match.");
2543 const unsigned int dof_no = 0,
2544 const unsigned int quad_no = 0,
2559 const std::pair<unsigned int, unsigned int> & range,
2561 const unsigned int dof_no = 0,
2562 const unsigned int quad_no = 0,
2576 reinit(
const unsigned int face_batch_number);
2586 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2593 const unsigned int give_n_q_points_1d);
2611 DEAL_II_DEPRECATED_EARLY
void
2612 evaluate(
const bool evaluate_values,
const bool evaluate_gradients);
2633 DEAL_II_DEPRECATED_EARLY
void
2635 const bool evaluate_values,
2636 const bool evaluate_gradients);
2649 template <
typename VectorType>
2657 template <
typename VectorType>
2658 DEAL_II_DEPRECATED_EARLY
void
2660 const bool evaluate_values,
2661 const bool evaluate_gradients);
2678 DEAL_II_DEPRECATED_EARLY
void
2679 integrate(
const bool integrate_values,
const bool integrate_gradients);
2691 VectorizedArrayType * values_array);
2696 DEAL_II_DEPRECATED_EARLY
void
2698 const bool integrate_gradients,
2699 VectorizedArrayType *values_array);
2712 template <
typename VectorType>
2715 VectorType & output_vector);
2720 template <
typename VectorType>
2723 const bool integrate_gradients,
2724 VectorType &output_vector);
2764 namespace MatrixFreeFunctions
2768 template <
int dim,
int degree>
2777 template <
int degree>
2780 static constexpr unsigned int value = degree + 1;
2796 template <
bool is_face,
2799 typename VectorizedArrayType>
2802 extract_initialization_data(
2804 const unsigned int dof_no,
2805 const unsigned int first_selected_component,
2806 const unsigned int quad_no,
2807 const unsigned int fe_degree,
2808 const unsigned int n_q_points,
2809 const unsigned int active_fe_index_given,
2810 const unsigned int active_quad_index_given,
2811 const unsigned int face_type)
2814 InitializationData init_data;
2818 &internal::MatrixFreeFunctions::
2819 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2827 active_fe_index_given :
2832 active_quad_index_given :
2863 typename VectorizedArrayType>
2868 VectorizedArrayType>::
2871 const unsigned int dof_no,
2872 const unsigned int first_selected_component,
2873 const unsigned int quad_no,
2874 const unsigned int fe_degree,
2875 const unsigned int n_q_points,
2876 const bool is_interior_face,
2877 const unsigned int active_fe_index,
2878 const unsigned int active_quad_index,
2879 const unsigned int face_type)
2881 internal::extract_initialization_data<is_face>(matrix_free,
2883 first_selected_component,
2892 first_selected_component)
2893 , scratch_data_array(matrix_free.acquire_scratch_data())
2894 , matrix_free(&matrix_free)
2896 this->set_data_pointers(scratch_data_array, n_components_);
2898 this->dof_info->start_components.back() == 1 ||
2899 static_cast<int>(n_components_) <=
2901 this->dof_info->start_components
2902 [this->dof_info->component_to_base_index[first_selected_component] +
2904 first_selected_component,
2906 "You tried to construct a vector-valued evaluator with " +
2907 std::to_string(n_components) +
2908 " components. However, "
2909 "the current base element has only " +
2911 this->dof_info->start_components
2912 [this->dof_info->component_to_base_index[first_selected_component] +
2914 first_selected_component) +
2915 " components left when starting from local element index " +
2917 first_selected_component -
2918 this->dof_info->start_components
2919 [this->dof_info->component_to_base_index[first_selected_component]]) +
2920 " (global index " + std::to_string(first_selected_component) +
")"));
2932 typename VectorizedArrayType>
2937 VectorizedArrayType>::
2943 const unsigned int first_selected_component,
2947 other->mapped_geometry->get_quadrature() == quadrature ?
2948 other->mapped_geometry :
2950 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2955 first_selected_component)
2956 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2957 , matrix_free(nullptr)
2959 const unsigned int base_element_number =
2963 first_selected_component >=
2965 ExcMessage(
"The underlying element must at least contain as many "
2966 "components as requested by this class"));
2967 (void)base_element_number;
2972 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2976 this->set_data_pointers(scratch_data_array, n_components_);
2985 typename VectorizedArrayType>
2990 VectorizedArrayType>::
2995 VectorizedArrayType> &other)
2997 , scratch_data_array(other.matrix_free == nullptr ?
2999 other.matrix_free->acquire_scratch_data())
3000 , matrix_free(other.matrix_free)
3002 if (other.matrix_free ==
nullptr)
3010 this->mapped_geometry =
3011 std::make_shared<internal::MatrixFreeFunctions::
3012 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3016 this->mapping_data = &this->mapped_geometry->get_data_storage();
3020 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3022 this->mapped_geometry->get_data_storage().JxW_values.begin();
3023 this->jacobian_gradients =
3024 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3026 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3029 this->set_data_pointers(scratch_data_array, n_components_);
3038 typename VectorizedArrayType>
3043 VectorizedArrayType> &
3049 VectorizedArrayType> &other)
3052 if (matrix_free ==
nullptr)
3055 delete scratch_data_array;
3064 matrix_free = other.matrix_free;
3066 if (other.matrix_free ==
nullptr)
3075 this->mapped_geometry =
3076 std::make_shared<internal::MatrixFreeFunctions::
3077 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3082 this->mapping_data = &this->mapped_geometry->get_data_storage();
3084 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3086 this->mapped_geometry->get_data_storage().JxW_values.begin();
3087 this->jacobian_gradients =
3088 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3090 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3097 this->set_data_pointers(scratch_data_array, n_components_);
3108 typename VectorizedArrayType>
3113 VectorizedArrayType>::~FEEvaluationBase()
3115 if (matrix_free !=
nullptr)
3126 delete scratch_data_array;
3137 typename VectorizedArrayType>
3142 Assert(matrix_free !=
nullptr,
3144 "FEEvaluation was not initialized with a MatrixFree object!"));
3145 return *matrix_free;
3154 template <
typename VectorType,
bool>
3155 struct ConstBlockVectorSelector;
3157 template <
typename VectorType>
3158 struct ConstBlockVectorSelector<VectorType, true>
3160 using BaseVectorType =
const typename VectorType::BlockType;
3163 template <
typename VectorType>
3164 struct ConstBlockVectorSelector<VectorType, false>
3166 using BaseVectorType =
typename VectorType::BlockType;
3172 template <
typename VectorType,
bool>
3173 struct BlockVectorSelector;
3175 template <
typename VectorType>
3176 struct BlockVectorSelector<VectorType, true>
3178 using BaseVectorType =
typename ConstBlockVectorSelector<
3180 std::is_const<VectorType>::value>::BaseVectorType;
3182 static BaseVectorType *
3183 get_vector_component(VectorType &vec,
const unsigned int component)
3186 return &vec.block(component);
3190 template <
typename VectorType>
3191 struct BlockVectorSelector<VectorType, false>
3193 using BaseVectorType = VectorType;
3195 static BaseVectorType *
3196 get_vector_component(VectorType &vec,
const unsigned int component)
3215 template <
typename VectorType>
3216 struct BlockVectorSelector<
std::vector<VectorType>, false>
3218 using BaseVectorType = VectorType;
3220 static BaseVectorType *
3221 get_vector_component(std::vector<VectorType> &vec,
3222 const unsigned int component)
3225 return &vec[component];
3229 template <
typename VectorType>
3230 struct BlockVectorSelector<const
std::vector<VectorType>, false>
3232 using BaseVectorType =
const VectorType;
3234 static const BaseVectorType *
3235 get_vector_component(
const std::vector<VectorType> &vec,
3236 const unsigned int component)
3239 return &vec[component];
3243 template <
typename VectorType>
3244 struct BlockVectorSelector<
std::vector<VectorType *>, false>
3246 using BaseVectorType = VectorType;
3248 static BaseVectorType *
3249 get_vector_component(std::vector<VectorType *> &vec,
3250 const unsigned int component)
3253 return vec[component];
3257 template <
typename VectorType>
3258 struct BlockVectorSelector<const
std::vector<VectorType *>, false>
3260 using BaseVectorType =
const VectorType;
3262 static const BaseVectorType *
3263 get_vector_component(
const std::vector<VectorType *> &vec,
3264 const unsigned int component)
3267 return vec[component];
3278 typename VectorizedArrayType>
3279template <
typename VectorType,
typename VectorOperation>
3284 const std::array<VectorType *, n_components_> &src,
3287 n_components_> & src_sm,
3288 const std::bitset<VectorizedArrayType::size()> &
mask,
3289 const bool apply_constraints)
const
3294 if (this->matrix_free ==
nullptr)
3296 read_write_operation_global(operation, src);
3302 if (this->n_fe_components == 1)
3303 for (
unsigned int comp = 0; comp < n_components; ++comp)
3305 Assert(src[comp] !=
nullptr,
3306 ExcMessage(
"The finite element underlying this FEEvaluation "
3307 "object is scalar, but you requested " +
3308 std::to_string(n_components) +
3309 " components via the template argument in "
3310 "FEEvaluation. In that case, you must pass an "
3311 "std::vector<VectorType> or a BlockVector to " +
3312 "read_dof_values and distribute_local_to_global."));
3330 this->dof_info->index_storage_variants[this->dof_access_index].size());
3331 if (this->dof_info->index_storage_variants
3332 [is_face ? this->dof_access_index :
3335 IndexStorageVariants::contiguous)
3337 read_write_operation_contiguous(operation, src, src_sm,
mask);
3344 constexpr unsigned int n_lanes = VectorizedArrayType::size();
3347 "non-contiguous DoF storage"));
3349 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
3350 this->get_cell_ids();
3352 bool has_hn_constraints =
false;
3354 if (is_face ==
false)
3356 for (
unsigned int v = 0; v < n_lanes; ++v)
3358 this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3359 this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3360 this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3363 this->dof_info->hanging_node_constraint_masks_comp
3364 [this->active_fe_index][this->first_selected_component])
3365 has_hn_constraints =
true;
3368 std::integral_constant<
bool,
3369 internal::is_vectorizable<VectorType, Number>::value>
3372 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3373 std::array<VectorizedArrayType *, n_components> values_dofs;
3374 for (
unsigned int c = 0; c < n_components; ++c)
3375 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3376 c * dofs_per_component;
3379 this->dof_info->index_storage_variants
3380 [is_face ? this->dof_access_index :
3383 IndexStorageVariants::interleaved &&
3384 (has_hn_constraints ==
false))
3386 const unsigned int *dof_indices =
3387 this->dof_info->dof_indices_interleaved.data() +
3388 this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes]
3391 ->component_dof_indices_offset[this->active_fe_index]
3392 [this->first_selected_component] *
3394 if (n_components == 1 || this->n_fe_components == 1)
3395 for (
unsigned int i = 0; i < dofs_per_component;
3396 ++i, dof_indices += n_lanes)
3397 for (
unsigned int comp = 0; comp < n_components; ++comp)
3398 operation.process_dof_gather(dof_indices,
3401 values_dofs[comp][i],
3404 for (
unsigned int comp = 0; comp < n_components; ++comp)
3405 for (
unsigned int i = 0; i < dofs_per_component;
3406 ++i, dof_indices += n_lanes)
3407 operation.process_dof_gather(
3408 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3414 std::array<const unsigned int *, n_lanes> dof_indices;
3415 dof_indices.fill(
nullptr);
3420 bool has_constraints =
false;
3421 const unsigned int n_components_read =
3422 this->n_fe_components > 1 ? n_components : 1;
3426 for (
unsigned int v = 0; v < n_lanes; ++v)
3431 Assert(cells[v] < this->dof_info->row_starts.size() - 1,
3433 const std::pair<unsigned int, unsigned int> *my_index_start =
3434 &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3435 this->first_selected_component];
3440 if (my_index_start[n_components_read].
second !=
3441 my_index_start[0].
second)
3442 has_constraints =
true;
3445 this->dof_info->dof_indices.data() + my_index_start[0].first;
3450 for (
unsigned int v = 0; v < n_lanes; ++v)
3455 const std::pair<unsigned int, unsigned int> *my_index_start =
3456 &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3457 this->first_selected_component];
3458 if (my_index_start[n_components_read].
second !=
3459 my_index_start[0].
second)
3460 has_constraints =
true;
3462 if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3463 this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3464 this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3467 this->dof_info->hanging_node_constraint_masks_comp
3468 [this->active_fe_index][this->first_selected_component])
3469 has_hn_constraints =
true;
3472 my_index_start[0].
first ||
3473 my_index_start[0].first < this->dof_info->dof_indices.size(),
3475 my_index_start[0].
first,
3476 this->dof_info->dof_indices.size()));
3478 this->dof_info->dof_indices.data() + my_index_start[0].first;
3482 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3483 return i != numbers::invalid_unsigned_int;
3485 for (
unsigned int comp = 0; comp < n_components; ++comp)
3486 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3487 operation.process_empty(values_dofs[comp][i]);
3491 if (!has_constraints && apply_constraints)
3493 if (n_components == 1 || this->n_fe_components == 1)
3495 for (
unsigned int v = 0; v < n_lanes; ++v)
3500 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3501 for (
unsigned int comp = 0; comp < n_components; ++comp)
3502 operation.process_dof(dof_indices[v][i],
3504 values_dofs[comp][i][v]);
3509 for (
unsigned int comp = 0; comp < n_components; ++comp)
3510 for (
unsigned int v = 0; v < n_lanes; ++v)
3515 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3516 operation.process_dof(
3517 dof_indices[v][comp * dofs_per_component + i],
3519 values_dofs[comp][i][v]);
3531 for (
unsigned int v = 0; v < n_lanes; ++v)
3537 const unsigned int cell_dof_index =
3538 cell_index * this->n_fe_components + this->first_selected_component;
3539 const unsigned int n_components_read =
3540 this->n_fe_components > 1 ? n_components : 1;
3541 unsigned int index_indicators =
3542 this->dof_info->row_starts[cell_dof_index].second;
3543 unsigned int next_index_indicators =
3544 this->dof_info->row_starts[cell_dof_index + 1].second;
3548 if (apply_constraints ==
false &&
3549 (this->dof_info->row_starts[cell_dof_index].second !=
3550 this->dof_info->row_starts[cell_dof_index + n_components_read]
3552 ((this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3553 this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3554 this->dof_info->hanging_node_constraint_masks[
cell_index] !=
3557 this->dof_info->hanging_node_constraint_masks_comp
3558 [this->active_fe_index][this->first_selected_component])))
3564 this->dof_info->plain_dof_indices.data() +
3566 ->component_dof_indices_offset[this->active_fe_index]
3567 [this->first_selected_component] +
3568 this->dof_info->row_starts_plain_indices[
cell_index];
3569 next_index_indicators = index_indicators;
3572 if (n_components == 1 || this->n_fe_components == 1)
3574 unsigned int ind_local = 0;
3575 for (; index_indicators != next_index_indicators; ++index_indicators)
3577 const std::pair<unsigned short, unsigned short> indicator =
3578 this->dof_info->constraint_indicator[index_indicators];
3580 for (
unsigned int j = 0; j < indicator.first; ++j)
3581 for (
unsigned int comp = 0; comp < n_components; ++comp)
3582 operation.process_dof(dof_indices[v][j],
3584 values_dofs[comp][ind_local + j][v]);
3586 ind_local += indicator.first;
3587 dof_indices[v] += indicator.first;
3591 Number
value[n_components];
3592 for (
unsigned int comp = 0; comp < n_components; ++comp)
3593 operation.pre_constraints(values_dofs[comp][ind_local][v],
3596 const Number *data_val =
3598 const Number *end_pool =
3600 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3601 for (
unsigned int comp = 0; comp < n_components; ++comp)
3602 operation.process_constraint(*dof_indices[v],
3607 for (
unsigned int comp = 0; comp < n_components; ++comp)
3608 operation.post_constraints(
value[comp],
3609 values_dofs[comp][ind_local][v]);
3615 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3616 for (
unsigned int comp = 0; comp < n_components; ++comp)
3617 operation.process_dof(*dof_indices[v],
3619 values_dofs[comp][ind_local][v]);
3628 for (
unsigned int comp = 0; comp < n_components; ++comp)
3630 unsigned int ind_local = 0;
3633 for (; index_indicators != next_index_indicators;
3636 const std::pair<unsigned short, unsigned short> indicator =
3637 this->dof_info->constraint_indicator[index_indicators];
3640 for (
unsigned int j = 0; j < indicator.first; ++j)
3641 operation.process_dof(dof_indices[v][j],
3643 values_dofs[comp][ind_local + j][v]);
3644 ind_local += indicator.first;
3645 dof_indices[v] += indicator.first;
3650 operation.pre_constraints(values_dofs[comp][ind_local][v],
3653 const Number *data_val =
3655 const Number *end_pool =
3658 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3659 operation.process_constraint(*dof_indices[v],
3664 operation.post_constraints(
value,
3665 values_dofs[comp][ind_local][v]);
3672 for (; ind_local < dofs_per_component;
3673 ++dof_indices[v], ++ind_local)
3676 operation.process_dof(*dof_indices[v],
3678 values_dofs[comp][ind_local][v]);
3681 if (apply_constraints ==
true && comp + 1 < n_components)
3682 next_index_indicators =
3683 this->dof_info->row_starts[cell_dof_index + comp + 2].second;
3695 typename VectorizedArrayType>
3696template <
typename VectorType,
typename VectorOperation>
3701 const std::array<VectorType *, n_components_> &src)
const
3705 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3706 unsigned int index = this->first_selected_component * dofs_per_component;
3707 for (
unsigned int comp = 0; comp < n_components; ++comp)
3709 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3711 operation.process_empty(
3712 this->values_dofs[comp * dofs_per_component + i]);
3713 operation.process_dof_global(
3714 local_dof_indices[this->data->lexicographic_numbering[
index]],
3716 this->values_dofs[comp * dofs_per_component + i][0]);
3727 typename VectorizedArrayType>
3728template <
typename VectorType,
typename VectorOperation>
3733 const std::array<VectorType *, n_components_> &src,
3736 n_components_> & vectors_sm,
3737 const std::bitset<VectorizedArrayType::size()> &
mask)
const
3746 std::integral_constant<
bool,
3747 internal::is_vectorizable<VectorType, Number>::value>
3750 is_face ? this->dof_access_index :
3752 const unsigned int n_lanes =
mask.count();
3754 const std::vector<unsigned int> &dof_indices_cont =
3755 this->dof_info->dof_indices_contiguous[ind];
3757 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3758 std::array<VectorizedArrayType *, n_components> values_dofs;
3759 for (
unsigned int c = 0; c < n_components; ++c)
3760 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3761 c * dofs_per_component;
3767 if ((this->dof_info->index_storage_variants[ind][this->cell] ==
3769 interleaved_contiguous &&
3770 n_lanes == VectorizedArrayType::size()) &&
3772 this->dof_access_index ==
3774 this->is_interior_face() ==
false) &&
3775 !(!is_face && !this->is_interior_face()))
3777 const unsigned int dof_index =
3778 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3780 ->component_dof_indices_offset[this->active_fe_index]
3781 [this->first_selected_component] *
3782 VectorizedArrayType::size();
3783 if (n_components == 1 || this->n_fe_components == 1)
3784 for (
unsigned int comp = 0; comp < n_components; ++comp)
3785 operation.process_dofs_vectorized(dofs_per_component,
3791 operation.process_dofs_vectorized(dofs_per_component * n_components,
3799 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
3800 this->get_cell_or_face_ids();
3804 const unsigned int n_filled_lanes =
3805 this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
3808 (this->dof_access_index ==
3810 this->is_interior_face() ==
false) ||
3811 (!is_face && !this->is_interior_face());
3813 if (vectors_sm[0] !=
nullptr)
3815 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3816 std::array<
typename VectorType::value_type *,
3817 VectorizedArrayType::size()>
3820 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3822 if (
mask[v] ==
false)
3824 vector_ptrs[v] =
nullptr;
3830 Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
3832 ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
3834 this->dof_info->dof_indices_contiguous_sm[ind].size(),
3838 this->dof_info->dof_indices_contiguous_sm[ind].size()));
3841 this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
3844 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3845 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3846 this->dof_info->component_dof_indices_offset
3847 [this->active_fe_index][this->first_selected_component]);
3849 vector_ptrs[v] =
nullptr;
3851 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3853 vector_ptrs[v] =
nullptr;
3858 if (n_filled_lanes == VectorizedArrayType::size() &&
3859 n_lanes == VectorizedArrayType::size() && !is_ecl)
3861 if (n_components == 1 || this->n_fe_components == 1)
3863 for (
unsigned int comp = 0; comp < n_components; ++comp)
3865 auto vector_ptrs = compute_vector_ptrs(comp);
3866 operation.process_dofs_vectorized_transpose(
3875 auto vector_ptrs = compute_vector_ptrs(0);
3876 operation.process_dofs_vectorized_transpose(dofs_per_component *
3884 for (
unsigned int comp = 0; comp < n_components; ++comp)
3886 auto vector_ptrs = compute_vector_ptrs(
3887 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3889 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3890 operation.process_empty(values_dofs[comp][i]);
3892 if (n_components == 1 || this->n_fe_components == 1)
3894 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3895 if (
mask[v] ==
true)
3896 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3897 operation.process_dof(vector_ptrs[v][i],
3898 values_dofs[comp][i][v]);
3902 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3903 if (
mask[v] ==
true)
3904 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3905 operation.process_dof(
3906 vector_ptrs[v][i + comp * dofs_per_component],
3907 values_dofs[comp][i][v]);
3913 unsigned int dof_indices[VectorizedArrayType::size()];
3915 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3919 dof_indices_cont[cells[v]] +
3921 ->component_dof_indices_offset[this->active_fe_index]
3922 [this->first_selected_component] *
3923 this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
3926 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3931 if (n_filled_lanes == VectorizedArrayType::size() &&
3932 n_lanes == VectorizedArrayType::size() && !is_ecl)
3934 if (this->dof_info->index_storage_variants[ind][this->cell] ==
3938 if (n_components == 1 || this->n_fe_components == 1)
3939 for (
unsigned int comp = 0; comp < n_components; ++comp)
3940 operation.process_dofs_vectorized_transpose(dofs_per_component,
3946 operation.process_dofs_vectorized_transpose(dofs_per_component *
3953 else if (this->dof_info->index_storage_variants[ind][this->cell] ==
3955 interleaved_contiguous_strided)
3957 if (n_components == 1 || this->n_fe_components == 1)
3958 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3960 for (
unsigned int comp = 0; comp < n_components; ++comp)
3961 operation.process_dof_gather(dof_indices,
3963 i * VectorizedArrayType::size(),
3964 values_dofs[comp][i],
3968 for (
unsigned int comp = 0; comp < n_components; ++comp)
3969 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3971 operation.process_dof_gather(dof_indices,
3973 (comp * dofs_per_component + i) *
3974 VectorizedArrayType::size(),
3975 values_dofs[comp][i],
3981 Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
3983 IndexStorageVariants::interleaved_contiguous_mixed_strides,
3985 const unsigned int *offsets =
3986 &this->dof_info->dof_indices_interleave_strides
3987 [ind][VectorizedArrayType::size() * this->cell];
3988 if (n_components == 1 || this->n_fe_components == 1)
3989 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3991 for (
unsigned int comp = 0; comp < n_components; ++comp)
3992 operation.process_dof_gather(dof_indices,
3995 values_dofs[comp][i],
3998 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3999 dof_indices[v] += offsets[v];
4002 for (
unsigned int comp = 0; comp < n_components; ++comp)
4003 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4005 operation.process_dof_gather(dof_indices,
4008 values_dofs[comp][i],
4011 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4012 dof_indices[v] += offsets[v];
4017 for (
unsigned int comp = 0; comp < n_components; ++comp)
4019 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4020 operation.process_empty(values_dofs[comp][i]);
4021 if (this->dof_info->index_storage_variants[ind][this->cell] ==
4025 if (n_components == 1 || this->n_fe_components == 1)
4027 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4028 if (
mask[v] ==
true)
4029 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4030 operation.process_dof(dof_indices[v] + i,
4032 values_dofs[comp][i][v]);
4036 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4037 if (
mask[v] ==
true)
4038 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4039 operation.process_dof(dof_indices[v] + i +
4040 comp * dofs_per_component,
4042 values_dofs[comp][i][v]);
4047 const unsigned int *offsets =
4048 &this->dof_info->dof_indices_interleave_strides
4049 [ind][VectorizedArrayType::size() * this->cell];
4050 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4052 if (n_components == 1 || this->n_fe_components == 1)
4053 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4055 if (
mask[v] ==
true)
4056 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4057 operation.process_dof(dof_indices[v] + i * offsets[v],
4059 values_dofs[comp][i][v]);
4063 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4064 if (
mask[v] ==
true)
4065 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4066 operation.process_dof(dof_indices[v] +
4067 (i + comp * dofs_per_component) *
4070 values_dofs[comp][i][v]);
4078 template <
typename Number,
4079 typename VectorType,
4080 typename std::enable_if<!IsBlockVector<VectorType>::value,
4081 VectorType>::type * =
nullptr>
4082 decltype(std::declval<VectorType>().begin())
4083 get_beginning(VectorType &vec)
4088 template <
typename Number,
4089 typename VectorType,
4090 typename std::enable_if<IsBlockVector<VectorType>::value,
4091 VectorType>::type * =
nullptr>
4092 typename VectorType::value_type *
4093 get_beginning(VectorType &)
4098 template <
typename VectorType,
4099 typename std::enable_if<has_shared_vector_data<VectorType>,
4100 VectorType>::type * =
nullptr>
4101 const std::vector<ArrayView<const typename VectorType::value_type>> *
4102 get_shared_vector_data(VectorType & vec,
4103 const bool is_valid_mode_for_sm,
4104 const unsigned int active_fe_index,
4108 if (is_valid_mode_for_sm &&
4111 active_fe_index == 0)
4112 return &vec.shared_vector_data();
4117 template <
typename VectorType,
4118 typename std::enable_if<!has_shared_vector_data<VectorType>,
4119 VectorType>::type * =
nullptr>
4120 const std::vector<ArrayView<const typename VectorType::value_type>> *
4121 get_shared_vector_data(VectorType &,
4129 template <
int n_components,
typename VectorType>
4131 std::array<
typename internal::BlockVectorSelector<
4136 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
4140 get_vector_data(VectorType & src,
4141 const unsigned int first_index,
4142 const bool is_valid_mode_for_sm,
4143 const unsigned int active_fe_index,
4149 std::array<
typename internal::BlockVectorSelector<
4155 ArrayView<
const typename internal::BlockVectorSelector<
4161 for (
unsigned int d = 0;
d < n_components; ++
d)
4162 src_data.first[d] = internal::BlockVectorSelector<
4168 for (
unsigned int d = 0;
d < n_components; ++
d)
4169 src_data.second[d] = get_shared_vector_data(
4170 const_cast<typename internal::BlockVectorSelector<
4171 typename std::remove_const<VectorType>::type,
4173 BaseVectorType &>(*src_data.first[d]),
4174 is_valid_mode_for_sm,
4188 typename VectorizedArrayType>
4193 if (this->dof_info ==
nullptr ||
4195 this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4196 this->dof_info->hanging_node_constraint_masks_comp
4197 [this->active_fe_index][this->first_selected_component] ==
false)
4200 constexpr unsigned int n_lanes = VectorizedArrayType::size();
4201 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4204 bool hn_available =
false;
4206 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
4207 this->get_cell_ids();
4209 for (
unsigned int v = 0; v < n_lanes; ++v)
4221 constraint_mask[v] =
mask;
4227 if (hn_available ==
false)
4232 this->data->data.front().fe_degree,
4233 this->get_shape_info(),
4245 typename VectorizedArrayType>
4246template <
typename VectorType>
4250 const unsigned int first_index,
4251 const std::bitset<VectorizedArrayType::size()> &
mask)
4253 const auto src_data = internal::get_vector_data<n_components_>(
4256 this->dof_access_index ==
4258 this->active_fe_index,
4262 read_write_operation(reader, src_data.first, src_data.second,
mask,
true);
4264 apply_hanging_node_constraints(
false);
4267 this->dof_values_initialized =
true;
4277 typename VectorizedArrayType>
4278template <
typename VectorType>
4282 const unsigned int first_index,
4283 const std::bitset<VectorizedArrayType::size()> &
mask)
4285 const auto src_data = internal::get_vector_data<n_components_>(
4288 this->dof_access_index ==
4290 this->active_fe_index,
4294 read_write_operation(reader, src_data.first, src_data.second,
mask,
false);
4297 this->dof_values_initialized =
true;
4307 typename VectorizedArrayType>
4308template <
typename VectorType>
4313 const unsigned int first_index,
4314 const std::bitset<VectorizedArrayType::size()> &
mask)
const
4317 Assert(this->dof_values_initialized ==
true,
4321 apply_hanging_node_constraints(
true);
4323 const auto dst_data = internal::get_vector_data<n_components_>(
4326 this->dof_access_index ==
4328 this->active_fe_index,
4333 read_write_operation(distributor, dst_data.first, dst_data.second,
mask);
4342 typename VectorizedArrayType>
4343template <
typename VectorType>
4347 const unsigned int first_index,
4348 const std::bitset<VectorizedArrayType::size()> &
mask)
const
4351 Assert(this->dof_values_initialized ==
true,
4355 const auto dst_data = internal::get_vector_data<n_components_>(
4358 this->dof_access_index ==
4360 this->active_fe_index,
4364 read_write_operation(setter, dst_data.first, dst_data.second,
mask);
4373 typename VectorizedArrayType>
4374template <
typename VectorType>
4379 const unsigned int first_index,
4380 const std::bitset<VectorizedArrayType::size()> &
mask)
const
4383 Assert(this->dof_values_initialized ==
true,
4387 const auto dst_data = internal::get_vector_data<n_components_>(
4390 this->dof_access_index ==
4392 this->active_fe_index,
4396 read_write_operation(setter, dst_data.first, dst_data.second,
mask,
false);
4409 typename VectorizedArrayType>
4415 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4417 for (
unsigned int comp = 0; comp < n_components; ++comp)
4418 return_value[comp] = this->values_dofs[comp * dofs + dof];
4419 return return_value;
4428 typename VectorizedArrayType>
4431 get_value(
const unsigned int q_point)
const
4434 Assert(this->values_quad_initialized ==
true,
4439 const std::size_t nqp = this->n_quadrature_points;
4441 for (
unsigned int comp = 0; comp < n_components; ++comp)
4442 return_value[comp] = this->values_quad[comp * nqp + q_point];
4443 return return_value;
4452 typename VectorizedArrayType>
4459 Assert(this->gradients_quad_initialized ==
true,
4464 Assert(this->jacobian !=
nullptr,
4466 "update_gradients"));
4467 const std::size_t nqp = this->n_quadrature_points;
4473 for (
unsigned int d = 0;
d < dim; ++
d)
4474 for (
unsigned int comp = 0; comp < n_components; ++comp)
4476 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4477 this->jacobian[0][
d][
d];
4486 for (
unsigned int comp = 0; comp < n_components; ++comp)
4487 for (
unsigned int d = 0;
d < dim; ++
d)
4490 jac[
d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4491 for (
unsigned int e = 1;
e < dim; ++
e)
4492 grad_out[comp][d] +=
4494 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4506 typename VectorizedArrayType>
4513 Assert(this->gradients_quad_initialized ==
true,
4517 Assert(this->normal_x_jacobian !=
nullptr,
4519 "update_gradients"));
4521 const std::size_t nqp = this->n_quadrature_points;
4525 for (
unsigned int comp = 0; comp < n_components; ++comp)
4527 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4528 (this->normal_x_jacobian[0][dim - 1]);
4531 const std::size_t
index =
4533 for (
unsigned int comp = 0; comp < n_components; ++comp)
4535 grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4536 this->normal_x_jacobian[
index][0];
4537 for (
unsigned int d = 1;
d < dim; ++
d)
4539 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4540 this->normal_x_jacobian[
index][
d];
4552 template <
typename VectorizedArrayType>
4555 const VectorizedArrayType *
const hessians,
4557 VectorizedArrayType (&tmp)[1][1])
4559 tmp[0][0] = jac[0][0] *
hessians[0];
4562 template <
typename VectorizedArrayType>
4565 const VectorizedArrayType *
const hessians,
4566 const unsigned int nqp,
4567 VectorizedArrayType (&tmp)[2][2])
4569 for (
unsigned int d = 0;
d < 2; ++
d)
4577 template <
typename VectorizedArrayType>
4580 const VectorizedArrayType *
const hessians,
4581 const unsigned int nqp,
4582 VectorizedArrayType (&tmp)[3][3])
4584 for (
unsigned int d = 0;
d < 3; ++
d)
4605 typename VectorizedArrayType>
4611 Assert(this->hessians_quad_initialized ==
true,
4616 Assert(this->jacobian !=
nullptr,
4626 const std::size_t nqp = this->n_quadrature_points;
4627 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4632 for (
unsigned int comp = 0; comp < n_components; ++comp)
4634 for (
unsigned int d = 0;
d < dim; ++
d)
4635 hessian_out[comp][d][d] =
4636 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4637 (jac[
d][
d] * jac[
d][
d]);
4643 hessian_out[comp][0][1] =
4644 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4645 (jac[0][0] * jac[1][1]);
4648 hessian_out[comp][0][1] =
4649 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4650 (jac[0][0] * jac[1][1]);
4651 hessian_out[comp][0][2] =
4652 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4653 (jac[0][0] * jac[2][2]);
4654 hessian_out[comp][1][2] =
4655 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4656 (jac[1][1] * jac[2][2]);
4661 for (
unsigned int d = 0;
d < dim; ++
d)
4662 for (
unsigned int e = d + 1;
e < dim; ++
e)
4663 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4669 for (
unsigned int comp = 0; comp < n_components; ++comp)
4671 VectorizedArrayType tmp[dim][dim];
4672 internal::hessian_unit_times_jac(
4673 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4676 for (
unsigned int d = 0;
d < dim; ++
d)
4677 for (
unsigned int e = d;
e < dim; ++
e)
4679 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4680 for (
unsigned int f = 1; f < dim; ++f)
4681 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4688 for (
unsigned int d = 0;
d < dim; ++
d)
4689 for (
unsigned int e = d + 1;
e < dim; ++
e)
4690 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4696 const auto &jac_grad = this->jacobian_gradients[q_point];
4697 for (
unsigned int comp = 0; comp < n_components; ++comp)
4699 VectorizedArrayType tmp[dim][dim];
4700 internal::hessian_unit_times_jac(
4701 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4704 for (
unsigned int d = 0;
d < dim; ++
d)
4705 for (
unsigned int e = d;
e < dim; ++
e)
4707 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4708 for (
unsigned int f = 1; f < dim; ++f)
4709 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4713 for (
unsigned int d = 0;
d < dim; ++
d)
4714 for (
unsigned int e = 0;
e < dim; ++
e)
4715 hessian_out[comp][d][d] +=
4717 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4720 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4721 for (
unsigned int e = d + 1;
e < dim; ++
e, ++count)
4722 for (
unsigned int f = 0; f < dim; ++f)
4723 hessian_out[comp][d][e] +=
4724 jac_grad[count][f] *
4725 this->gradients_quad[(comp * dim + f) * nqp + q_point];
4728 for (
unsigned int d = 0;
d < dim; ++
d)
4729 for (
unsigned int e = d + 1;
e < dim; ++
e)
4730 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4742 typename VectorizedArrayType>
4749 Assert(this->hessians_quad_initialized ==
true,
4760 const std::size_t nqp = this->n_quadrature_points;
4761 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4767 for (
unsigned int comp = 0; comp < n_components; ++comp)
4768 for (
unsigned int d = 0;
d < dim; ++
d)
4769 hessian_out[comp][d] =
4770 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4771 (jac[
d][
d] * jac[
d][
d]);
4776 for (
unsigned int comp = 0; comp < n_components; ++comp)
4780 VectorizedArrayType tmp[dim][dim];
4781 internal::hessian_unit_times_jac(
4782 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4786 for (
unsigned int d = 0;
d < dim; ++
d)
4788 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4789 for (
unsigned int f = 1; f < dim; ++f)
4790 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4797 const auto &jac_grad = this->jacobian_gradients[q_point];
4798 for (
unsigned int comp = 0; comp < n_components; ++comp)
4802 VectorizedArrayType tmp[dim][dim];
4803 internal::hessian_unit_times_jac(
4804 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4808 for (
unsigned int d = 0;
d < dim; ++
d)
4810 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4811 for (
unsigned int f = 1; f < dim; ++f)
4812 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4815 for (
unsigned int d = 0;
d < dim; ++
d)
4816 for (
unsigned int e = 0;
e < dim; ++
e)
4817 hessian_out[comp][d] +=
4819 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4831 typename VectorizedArrayType>
4838 Assert(this->hessians_quad_initialized ==
true,
4844 const auto hess_diag = get_hessian_diagonal(q_point);
4845 for (
unsigned int comp = 0; comp < n_components; ++comp)
4847 laplacian_out[comp] = hess_diag[comp][0];
4848 for (
unsigned int d = 1;
d < dim; ++
d)
4849 laplacian_out[comp] += hess_diag[comp][d];
4851 return laplacian_out;
4860 typename VectorizedArrayType>
4864 const unsigned int dof)
4867 this->dof_values_initialized =
true;
4869 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4871 for (
unsigned int comp = 0; comp < n_components; ++comp)
4872 this->values_dofs[comp * dofs + dof] = val_in[comp];
4881 typename VectorizedArrayType>
4885 const unsigned int q_point)
4891 Assert(this->J_value !=
nullptr,
4895 this->values_quad_submitted =
true;
4898 const std::size_t nqp = this->n_quadrature_points;
4901 const VectorizedArrayType JxW =
4902 this->J_value[0] * this->quadrature_weights[q_point];
4903 for (
unsigned int comp = 0; comp < n_components; ++comp)
4904 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4908 const VectorizedArrayType JxW = this->J_value[q_point];
4909 for (
unsigned int comp = 0; comp < n_components; ++comp)
4910 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4920 typename VectorizedArrayType>
4925 const unsigned int q_point)
4931 Assert(this->J_value !=
nullptr,
4933 "update_gradients"));
4934 Assert(this->jacobian !=
nullptr,
4936 "update_gradients"));
4938 this->gradients_quad_submitted =
true;
4941 const std::size_t nqp = this->n_quadrature_points;
4944 const VectorizedArrayType JxW =
4945 this->J_value[0] * this->quadrature_weights[q_point];
4946 for (
unsigned int d = 0;
d < dim; ++
d)
4948 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
4949 for (
unsigned int comp = 0; comp < n_components; ++comp)
4950 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
4951 grad_in[comp][
d] * factor;
4958 this->jacobian[q_point] :
4960 const VectorizedArrayType JxW =
4962 this->J_value[q_point] :
4963 this->J_value[0] * this->quadrature_weights[q_point];
4964 for (
unsigned int comp = 0; comp < n_components; ++comp)
4965 for (
unsigned int d = 0;
d < dim; ++
d)
4967 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
4968 for (
unsigned int e = 1;
e < dim; ++
e)
4969 new_val += (jac[e][d] * grad_in[comp][e]);
4970 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
4982 typename VectorizedArrayType>
4987 const unsigned int q_point)
4990 Assert(this->normal_x_jacobian !=
nullptr,
4992 "update_gradients"));
4994 this->gradients_quad_submitted =
true;
4997 const std::size_t nqp = this->n_quadrature_points;
4999 for (
unsigned int comp = 0; comp < n_components; ++comp)
5001 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5002 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5003 VectorizedArrayType();
5004 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5006 (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5007 this->quadrature_weights[q_point]);
5011 const unsigned int index =
5014 this->normal_x_jacobian[
index];
5015 for (
unsigned int comp = 0; comp < n_components; ++comp)
5017 VectorizedArrayType factor = grad_in[comp] * this->J_value[
index];
5019 factor = factor * this->quadrature_weights[q_point];
5020 for (
unsigned int d = 0;
d < dim; ++
d)
5021 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5033 typename VectorizedArrayType>
5039 const unsigned int q_point)
5045 Assert(this->J_value !=
nullptr,
5047 "update_hessians"));
5048 Assert(this->jacobian !=
nullptr,
5050 "update_hessians"));
5052 this->hessians_quad_submitted =
true;
5056 const std::size_t nqp = this->n_quadrature_points;
5057 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5060 const VectorizedArrayType JxW =
5061 this->J_value[0] * this->quadrature_weights[q_point];
5064 for (
unsigned int d = 0;
d < dim; ++
d)
5066 const auto jac_d = this->jacobian[0][
d][
d];
5067 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5068 for (
unsigned int comp = 0; comp < n_components; ++comp)
5069 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5070 hessian_in[comp][
d][
d] * factor;
5074 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5075 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5077 const auto jac_d = this->jacobian[0][
d][
d];
5078 const auto jac_e = this->jacobian[0][
e][
e];
5079 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5080 for (
unsigned int comp = 0; comp < n_components; ++comp)
5081 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5082 (hessian_in[comp][
d][
e] + hessian_in[comp][
e][
d]) * factor;
5089 const VectorizedArrayType JxW =
5090 this->J_value[0] * this->quadrature_weights[q_point];
5091 for (
unsigned int comp = 0; comp < n_components; ++comp)
5094 VectorizedArrayType tmp[dim][dim];
5095 for (
unsigned int i = 0; i < dim; ++i)
5096 for (
unsigned int j = 0; j < dim; ++j)
5098 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5099 for (
unsigned int k = 1; k < dim; ++k)
5100 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5104 VectorizedArrayType tmp2[dim][dim];
5105 for (
unsigned int i = 0; i < dim; ++i)
5106 for (
unsigned int j = 0; j < dim; ++j)
5108 tmp2[i][j] = jac[0][i] * tmp[0][j];
5109 for (
unsigned int k = 1; k < dim; ++k)
5110 tmp2[i][j] += jac[k][i] * tmp[k][j];
5114 for (
unsigned int d = 0;
d < dim; ++
d)
5115 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5119 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5120 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5121 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5122 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5128 const VectorizedArrayType JxW = this->J_value[q_point];
5129 const auto &jac_grad = this->jacobian_gradients[q_point];
5130 for (
unsigned int comp = 0; comp < n_components; ++comp)
5133 VectorizedArrayType tmp[dim][dim];
5134 for (
unsigned int i = 0; i < dim; ++i)
5135 for (
unsigned int j = 0; j < dim; ++j)
5137 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5138 for (
unsigned int k = 1; k < dim; ++k)
5139 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5143 VectorizedArrayType tmp2[dim][dim];
5144 for (
unsigned int i = 0; i < dim; ++i)
5145 for (
unsigned int j = 0; j < dim; ++j)
5147 tmp2[i][j] = jac[0][i] * tmp[0][j];
5148 for (
unsigned int k = 1; k < dim; ++k)
5149 tmp2[i][j] += jac[k][i] * tmp[k][j];
5153 for (
unsigned int d = 0;
d < dim; ++
d)
5154 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5158 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5159 for (
unsigned int e = d + 1;
e < dim; ++
e, ++off_diag)
5160 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5161 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5164 for (
unsigned int d = 0;
d < dim; ++
d)
5166 VectorizedArrayType
sum = 0;
5167 for (
unsigned int e = 0;
e < dim; ++
e)
5168 sum += hessian_in[comp][e][e] * jac_grad[e][d];
5169 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5170 for (
unsigned int f = e + 1; f < dim; ++f, ++count)
5171 sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5173 this->gradients_from_hessians_quad[(comp * dim +
d) * nqp +
5174 q_point] = sum * JxW;
5186 typename VectorizedArrayType>
5193 Assert(this->values_quad_submitted ==
true,
5198 const std::size_t nqp = this->n_quadrature_points;
5199 for (
unsigned int q = 0; q < nqp; ++q)
5200 for (
unsigned int comp = 0; comp < n_components; ++comp)
5201 return_value[comp] += this->values_quad[comp * nqp + q];
5202 return (return_value);
5214 typename VectorizedArrayType>
5219 VectorizedArrayType>::
5222 const unsigned int dof_no,
5223 const unsigned int first_selected_component,
5224 const unsigned int quad_no,
5225 const unsigned int fe_degree,
5226 const unsigned int n_q_points,
5227 const bool is_interior_face,
5228 const unsigned int active_fe_index,
5229 const unsigned int active_quad_index,
5230 const unsigned int face_type)
5231 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5234 first_selected_component,
5250 typename VectorizedArrayType>
5255 VectorizedArrayType>::
5261 const unsigned int first_selected_component,
5263 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5268 first_selected_component,
5278 typename VectorizedArrayType>
5283 VectorizedArrayType>::
5288 VectorizedArrayType> &other)
5289 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5299 typename VectorizedArrayType>
5304 VectorizedArrayType> &
5310 VectorizedArrayType> &other)
5316 VectorizedArrayType>::operator=(other);
5325template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5329 const unsigned int dof_no,
5330 const unsigned int first_selected_component,
5331 const unsigned int quad_no,
5332 const unsigned int fe_degree,
5333 const unsigned int n_q_points,
5334 const bool is_interior_face,
5335 const unsigned int active_fe_index,
5336 const unsigned int active_quad_index,
5337 const unsigned int face_type)
5341 first_selected_component,
5353template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5360 const unsigned int first_selected_component,
5367 first_selected_component,
5373template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5383template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5389 ->
FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5396template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5399 const unsigned int dof)
const
5402 return this->values_dofs[dof];
5407template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5410 const unsigned int q_point)
const
5413 Assert(this->values_quad_initialized ==
true,
5417 return this->values_quad[q_point];
5422template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5427 return BaseClass::get_normal_derivative(q_point)[0];
5432template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5435 const unsigned int q_point)
const
5441 Assert(this->gradients_quad_initialized ==
true,
5446 Assert(this->jacobian !=
nullptr,
5448 "update_gradients"));
5452 const std::size_t nqp = this->n_quadrature_points;
5455 for (
unsigned int d = 0;
d < dim; ++
d)
5457 this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5466 for (
unsigned int d = 0;
d < dim; ++
d)
5468 grad_out[
d] = jac[
d][0] * this->gradients_quad[q_point];
5469 for (
unsigned int e = 1;
e < dim; ++
e)
5470 grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5478template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5481 const unsigned int q_point)
const
5483 return BaseClass::get_hessian(q_point)[0];
5488template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5493 return BaseClass::get_hessian_diagonal(q_point)[0];
5498template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5499inline VectorizedArrayType
5501 const unsigned int q_point)
const
5503 return BaseClass::get_laplacian(q_point)[0];
5508template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5514 this->dof_values_initialized =
true;
5517 this->values_dofs[dof] = val_in;
5522template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5525 const VectorizedArrayType val_in,
5526 const unsigned int q_point)
5532 Assert(this->J_value !=
nullptr,
5536 this->values_quad_submitted =
true;
5541 const VectorizedArrayType JxW =
5542 this->J_value[0] * this->quadrature_weights[q_point];
5543 this->values_quad[q_point] = val_in * JxW;
5547 this->values_quad[q_point] = val_in * this->J_value[q_point];
5553template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5557 const unsigned int q_point)
5559 submit_value(val_in[0], q_point);
5564template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5568 const unsigned int q_point)
5572 BaseClass::submit_normal_derivative(grad, q_point);
5577template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5581 const unsigned int q_point)
5587 Assert(this->J_value !=
nullptr,
5589 "update_gradients"));
5590 Assert(this->jacobian !=
nullptr,
5592 "update_gradients"));
5594 this->gradients_quad_submitted =
true;
5597 const std::size_t nqp = this->n_quadrature_points;
5600 const VectorizedArrayType JxW =
5601 this->J_value[0] * this->quadrature_weights[q_point];
5602 for (
unsigned int d = 0;
d < dim; ++
d)
5603 this->gradients_quad[d * nqp + q_point] =
5604 (grad_in[d] * this->jacobian[0][d][d] * JxW);
5611 this->jacobian[q_point] :
5613 const VectorizedArrayType JxW =
5615 this->J_value[q_point] :
5616 this->J_value[0] * this->quadrature_weights[q_point];
5617 for (
unsigned int d = 0;
d < dim; ++
d)
5619 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5620 for (
unsigned int e = 1;
e < dim; ++
e)
5621 new_val += jac[e][d] * grad_in[e];
5622 this->gradients_quad[
d * nqp + q_point] = new_val * JxW;
5629template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5633 const unsigned int q_point)
5637 BaseClass::submit_hessian(
hessian, q_point);
5642template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5643inline VectorizedArrayType
5647 return BaseClass::integrate_value()[0];
5655template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5659 const unsigned int dof_no,
5660 const unsigned int first_selected_component,
5661 const unsigned int quad_no,
5662 const unsigned int fe_degree,
5663 const unsigned int n_q_points,
5664 const bool is_interior_face,
5665 const unsigned int active_fe_index,
5666 const unsigned int active_quad_index,
5667 const unsigned int face_type)
5671 first_selected_component,
5683template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5690 const unsigned int first_selected_component,
5697 first_selected_component,
5703template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5713template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5725template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5728 const unsigned int q_point)
const
5730 if (this->data->element_type ==
5735 Assert(this->values_quad_initialized ==
true,
5740 Assert(this->J_value !=
nullptr,
5743 const std::size_t nqp = this->n_quadrature_points;
5751 const VectorizedArrayType inv_det =
5752 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5753 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5754 this->jacobian[0][2][2];
5757 for (
unsigned int comp = 0; comp < n_components; ++comp)
5758 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5759 jac[comp][comp] * inv_det;
5766 this->jacobian[q_point] :
5774 const VectorizedArrayType inv_det =
5775 (is_face && dim == 2 && this->get_face_no() < 2) ?
5779 for (
unsigned int comp = 0; comp < n_components; ++comp)
5782 this->values_quad[q_point] * jac[comp][0] * inv_det;
5783 for (
unsigned int e = 1;
e < dim; ++
e)
5785 this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
5793 return BaseClass::get_value(q_point);
5797template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5802 if (this->data->element_type ==
5807 Assert(this->gradients_quad_initialized ==
true,
5812 Assert(this->jacobian !=
nullptr,
5814 "update_gradients"));
5815 const std::size_t nqp = this->n_quadrature_points;
5825 const VectorizedArrayType inv_det =
5826 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5827 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5828 this->jacobian[0][2][2];
5831 for (
unsigned int d = 0;
d < dim; ++
d)
5832 for (
unsigned int comp = 0; comp < n_components; ++comp)
5834 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
5835 inv_t_jac[
d][
d] * jac[comp][comp] * inv_det;
5845 const VectorizedArrayType inv_det =
5846 (is_face && dim == 2 && this->get_face_no() < 2) ?
5850 VectorizedArrayType tmp;
5852 for (
unsigned int comp = 0; comp < n_components; ++comp)
5853 for (
unsigned int d = 0;
d < dim; ++
d)
5856 for (
unsigned int f = 0; f < dim; ++f)
5857 for (
unsigned int e = 0;
e < dim; ++
e)
5858 tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
5859 this->gradients_quad[(f * dim + e) * nqp + q_point];
5861 grad_out[comp][
d] = tmp;
5875 return BaseClass::get_gradient(q_point);
5881template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5887 Assert(this->gradients_quad_initialized ==
true,
5891 Assert(this->jacobian !=
nullptr,
5893 "update_gradients"));
5895 VectorizedArrayType divergence;
5896 const std::size_t nqp = this->n_quadrature_points;
5898 if (this->data->element_type ==
5905 const VectorizedArrayType inv_det =
5906 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5907 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5908 this->jacobian[0][2][2];
5911 divergence = this->gradients_quad[q_point] * inv_det;
5912 for (
unsigned int d = 1;
d < dim; ++
d)
5914 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5920 const VectorizedArrayType inv_det =
5921 (is_face && dim == 2 && this->get_face_no() < 2) ?
5926 divergence = this->gradients_quad[q_point] * inv_det;
5927 for (
unsigned int d = 1;
d < dim; ++
d)
5929 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5943 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
5944 for (
unsigned int d = 1;
d < dim; ++
d)
5945 divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
5946 this->jacobian[0][
d][
d];
5953 this->jacobian[q_point] :
5955 divergence = jac[0][0] * this->gradients_quad[q_point];
5956 for (
unsigned int e = 1;
e < dim; ++
e)
5957 divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
5958 for (
unsigned int d = 1;
d < dim; ++
d)
5959 for (
unsigned int e = 0;
e < dim; ++
e)
5961 jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
5969template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5975 const auto grad = get_gradient(q_point);
5976 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5977 VectorizedArrayType half = Number(0.5);
5978 for (
unsigned int d = 0;
d < dim; ++
d)
5979 symmetrized[d] = grad[d][d];
5985 symmetrized[2] = grad[0][1] + grad[1][0];
5986 symmetrized[2] *= half;
5989 symmetrized[3] = grad[0][1] + grad[1][0];
5990 symmetrized[3] *= half;
5991 symmetrized[4] = grad[0][2] + grad[2][0];
5992 symmetrized[4] *= half;
5993 symmetrized[5] = grad[1][2] + grad[2][1];
5994 symmetrized[5] *= half;
6004template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6006 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6008 const unsigned int q_point)
const
6012 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6018 "Computing the curl in 1d is not a useful operation"));
6021 curl[0] = grad[1][0] - grad[0][1];
6024 curl[0] = grad[2][1] - grad[1][2];
6025 curl[1] = grad[0][2] - grad[2][0];
6026 curl[2] = grad[1][0] - grad[0][1];
6036template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6041 return BaseClass::get_hessian_diagonal(q_point);
6046template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6049 const unsigned int q_point)
const
6052 Assert(this->hessians_quad_initialized ==
true,
6056 return BaseClass::get_hessian(q_point);
6060template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6064 const unsigned int q_point)
6066 if (this->data->element_type ==
6071 Assert(this->J_value !=
nullptr,
6076 this->values_quad_submitted =
true;
6079 const std::size_t nqp = this->n_quadrature_points;
6084 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6086 for (
unsigned int comp = 0; comp < n_components; ++comp)
6087 this->values_quad[comp * nqp + q_point] =
6088 val_in[comp] * weight * jac[comp][comp];
6095 this->jacobian[q_point] :
6104 const VectorizedArrayType fac =
6106 this->quadrature_weights[q_point] :
6108 this->J_value[q_point] :
6109 this->J_value[0] * this->quadrature_weights[q_point]) *
6110 ((dim == 2 && this->get_face_no() < 2) ?
6115 for (
unsigned int comp = 0; comp < n_components; ++comp)
6117 this->values_quad[comp * nqp + q_point] =
6118 val_in[0] * jac[0][comp] * fac;
6119 for (
unsigned int e = 1;
e < dim; ++
e)
6120 this->values_quad[comp * nqp + q_point] +=
6121 val_in[e] * jac[e][comp] * fac;
6128 BaseClass::submit_value(val_in, q_point);
6132template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6136 const unsigned int q_point)
6138 if (this->data->element_type ==
6147 Assert(this->J_value !=
nullptr,
6149 "update_gradients"));
6150 Assert(this->jacobian !=
nullptr,
6152 "update_gradients"));
6154 this->gradients_quad_submitted =
true;
6157 const std::size_t nqp = this->n_quadrature_points;
6165 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6166 for (
unsigned int d = 0;
d < dim; ++
d)
6167 for (
unsigned int comp = 0; comp < n_components; ++comp)
6168 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6169 grad_in[comp][
d] * inv_t_jac[
d][
d] * jac[comp][comp] * weight;
6180 const VectorizedArrayType fac =
6181 (!is_face) ? this->quadrature_weights[q_point] :
6182 this->J_value[0] * this->quadrature_weights[q_point] *
6183 ((dim == 2 && this->get_face_no() < 2) ?
6188 for (
unsigned int comp = 0; comp < n_components; ++comp)
6189 for (
unsigned int d = 0;
d < dim; ++
d)
6191 VectorizedArrayType tmp = 0;
6192 for (
unsigned int f = 0; f < dim; ++f)
6193 for (
unsigned int e = 0;
e < dim; ++
e)
6194 tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
6196 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = tmp;
6207 BaseClass::submit_gradient(grad_in, q_point);
6213template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6218 const unsigned int q_point)
6220 if (this->data->element_type ==
6230 BaseClass::submit_gradient(grad_in, q_point);
6236template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6240 const unsigned int q_point)
6246 Assert(this->J_value !=
nullptr,
6248 "update_gradients"));
6249 Assert(this->jacobian !=
nullptr,
6251 "update_gradients"));
6253 this->gradients_quad_submitted =
true;
6256 const std::size_t nqp = this->n_quadrature_points;
6257 if (this->data->element_type ==
6266 const VectorizedArrayType fac =
6269 this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
6272 this->quadrature_weights[q_point] * div_in;
6274 for (
unsigned int d = 0;
d < dim; ++
d)
6276 this->gradients_quad[(dim *
d +
d) * nqp + q_point] = fac;
6277 for (
unsigned int e = d + 1;
e < dim; ++
e)
6279 this->gradients_quad[(dim *
d +
e) * nqp + q_point] =
6280 VectorizedArrayType();
6281 this->gradients_quad[(dim *
e +
d) * nqp + q_point] =
6282 VectorizedArrayType();
6297 const VectorizedArrayType fac =
6298 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6299 for (
unsigned int d = 0;
d < dim; ++
d)
6301 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6302 (fac * this->jacobian[0][d][d]);
6303 for (
unsigned int e = d + 1;
e < dim; ++
e)
6305 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6306 VectorizedArrayType();
6307 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6308 VectorizedArrayType();
6316 this->jacobian[q_point] :
6318 const VectorizedArrayType fac =
6320 this->J_value[q_point] :
6321 this->J_value[0] * this->quadrature_weights[q_point]) *
6323 for (
unsigned int d = 0;
d < dim; ++
d)
6325 for (
unsigned int e = 0;
e < dim; ++
e)
6326 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6335template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6340 const unsigned int q_point)
6343 this->data->element_type !=
6354 Assert(this->J_value !=
nullptr,
6356 "update_gradients"));
6357 Assert(this->jacobian !=
nullptr,
6359 "update_gradients"));
6361 this->gradients_quad_submitted =
true;
6364 const std::size_t nqp = this->n_quadrature_points;
6367 const VectorizedArrayType JxW =
6368 this->J_value[0] * this->quadrature_weights[q_point];
6369 for (
unsigned int d = 0;
d < dim; ++
d)
6370 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6372 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6373 for (
unsigned int d = e + 1;
d < dim; ++
d, ++counter)
6375 const VectorizedArrayType
value =
6377 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6378 value * this->jacobian[0][d][d];
6379 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6380 value * this->jacobian[0][e][e];
6386 const VectorizedArrayType JxW =
6388 this->J_value[q_point] :
6389 this->J_value[0] * this->quadrature_weights[q_point];
6392 this->jacobian[q_point] :
6394 VectorizedArrayType weighted[dim][dim];
6395 for (
unsigned int i = 0; i < dim; ++i)
6397 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6398 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6400 const VectorizedArrayType
value =
6402 weighted[i][j] =
value;
6403 weighted[j][i] =
value;
6405 for (
unsigned int comp = 0; comp < dim; ++comp)
6406 for (
unsigned int d = 0;
d < dim; ++
d)
6408 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6409 for (
unsigned int e = 1;
e < dim; ++
e)
6410 new_val += jac[e][d] * weighted[comp][e];
6411 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
6418template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6422 const unsigned int q_point)
6430 "Testing by the curl in 1d is not a useful operation"));
6433 grad[1][0] = curl[0];
6434 grad[0][1] = -curl[0];
6437 grad[2][1] = curl[0];
6438 grad[1][2] = -curl[0];
6439 grad[0][2] = curl[1];
6440 grad[2][0] = -curl[1];
6441 grad[1][0] = curl[2];
6442 grad[0][1] = -curl[2];
6447 submit_gradient(grad, q_point);
6454template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6458 const unsigned int dof_no,
6459 const unsigned int first_selected_component,
6460 const unsigned int quad_no,
6461 const unsigned int fe_degree,
6462 const unsigned int n_q_points,
6463 const bool is_interior_face,
6464 const unsigned int active_fe_index,
6465 const unsigned int active_quad_index,
6466 const unsigned int face_type)
6470 first_selected_component,
6482template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6489 const unsigned int first_selected_component,
6496 first_selected_component,
6502template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6511template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6523template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6526 const unsigned int dof)
const
6529 return this->values_dofs[dof];
6534template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6537 const unsigned int q_point)
const
6540 Assert(this->values_quad_initialized ==
true,
6544 return this->values_quad[q_point];
6549template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6552 const unsigned int q_point)
const
6558 Assert(this->gradients_quad_initialized ==
true,
6565 this->jacobian[q_point] :
6569 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6576template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6579 const unsigned int q_point)
const
6581 return get_gradient(q_point)[0];
6586template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6591 return BaseClass::get_normal_derivative(q_point)[0];
6596template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6599 const unsigned int q_point)
const
6601 return BaseClass::get_hessian(q_point)[0];
6606template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6611 return BaseClass::get_hessian_diagonal(q_point)[0];
6616template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6619 const unsigned int q_point)
const
6621 return BaseClass::get_laplacian(q_point)[0];
6626template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6632 this->dof_values_initialized =
true;
6635 this->values_dofs[dof] = val_in;
6640template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6643 const VectorizedArrayType val_in,
6644 const unsigned int q_point)
6651 this->values_quad_submitted =
true;
6656 const VectorizedArrayType JxW = this->J_value[q_point];
6657 this->values_quad[q_point] = val_in * JxW;
6661 const VectorizedArrayType JxW =
6662 this->J_value[0] * this->quadrature_weights[q_point];
6663 this->values_quad[q_point] = val_in * JxW;
6669template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6673 const unsigned int q_point)
6675 submit_value(val_in[0], q_point);
6680template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6684 const unsigned int q_point)
6686 submit_gradient(grad_in[0], q_point);
6691template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6694 const VectorizedArrayType grad_in,
6695 const unsigned int q_point)
6702 this->gradients_quad_submitted =
true;
6707 this->jacobian[q_point] :
6709 const VectorizedArrayType JxW =
6711 this->J_value[q_point] :
6712 this->J_value[0] * this->quadrature_weights[q_point];
6714 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6719template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6723 const unsigned int q_point)
6725 submit_gradient(grad_in[0][0], q_point);
6730template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6734 const unsigned int q_point)
6738 BaseClass::submit_normal_derivative(grad, q_point);
6743template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6747 const unsigned int q_point)
6749 BaseClass::submit_normal_derivative(grad_in, q_point);
6753template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6757 const unsigned int q_point)
6761 BaseClass::submit_hessian(
hessian, q_point);
6765template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6766inline VectorizedArrayType
6770 return BaseClass::integrate_value()[0];
6783 typename VectorizedArrayType>
6789 VectorizedArrayType>::
6791 const unsigned int fe_no,
6792 const unsigned int quad_no,
6793 const unsigned int first_selected_component,
6794 const unsigned int active_fe_index,
6795 const unsigned int active_quad_index)
6796 : BaseClass(matrix_free,
6798 first_selected_component,
6805 , dofs_per_component(this->data->dofs_per_component_on_cell)
6806 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6807 , n_q_points(this->data->n_q_points)
6809 check_template_arguments(fe_no, 0);
6819 typename VectorizedArrayType>
6825 VectorizedArrayType>::
6827 const std::pair<unsigned int, unsigned int> & range,
6828 const unsigned int dof_no,
6829 const unsigned int quad_no,
6830 const unsigned int first_selected_component)
6834 first_selected_component,
6835 matrix_free.get_cell_active_fe_index(range))
6845 typename VectorizedArrayType>
6851 VectorizedArrayType>::
6856 const unsigned int first_selected_component)
6857 : BaseClass(mapping,
6861 first_selected_component,
6863 , dofs_per_component(this->data->dofs_per_component_on_cell)
6864 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6865 , n_q_points(this->data->n_q_points)
6877 typename VectorizedArrayType>
6883 VectorizedArrayType>::
6887 const unsigned int first_selected_component)
6892 first_selected_component,
6894 , dofs_per_component(this->data->dofs_per_component_on_cell)
6895 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6896 , n_q_points(this->data->n_q_points)
6908 typename VectorizedArrayType>
6914 VectorizedArrayType>::
6917 const unsigned int first_selected_component)
6918 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6920 other.mapped_geometry->get_quadrature(),
6921 other.mapped_geometry->get_fe_values().get_update_flags(),
6922 first_selected_component,
6924 , dofs_per_component(this->data->dofs_per_component_on_cell)
6925 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6926 , n_q_points(this->data->n_q_points)
6938 typename VectorizedArrayType>
6947 , dofs_per_component(this->data->dofs_per_component_on_cell)
6948 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6949 , n_q_points(this->data->n_q_points)
6961 typename VectorizedArrayType>
6967 VectorizedArrayType> &
6973 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6975 BaseClass::operator=(other);
6987 typename VectorizedArrayType>
6994 VectorizedArrayType>::
6995 check_template_arguments(
const unsigned int dof_no,
6996 const unsigned int first_selected_component)
6999 (void)first_selected_component;
7002 this->data->dofs_per_component_on_cell > 0,
7004 "There is nothing useful you can do with an FEEvaluation object with "
7005 "FE_Nothing, i.e., without DoFs! If you have passed to "
7006 "MatrixFree::reinit() a collection of finite elements also containing "
7007 "FE_Nothing, please check - before creating FEEvaluation - the category "
7008 "of the current range by calling either "
7009 "MatrixFree::get_cell_range_category(range) or "
7010 "MatrixFree::get_face_range_category(range). The returned category "
7011 "is the index of the active FE, which you can use to exclude "
7018 static_cast<unsigned int>(fe_degree) !=
7019 this->data->data.front().fe_degree) ||
7020 n_q_points != this->n_quadrature_points)
7022 std::string message =
7023 "-------------------------------------------------------\n";
7024 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7025 message +=
" Called --> FEEvaluation<dim,";
7029 message +=
",Number>(data";
7045 if (
static_cast<unsigned int>(fe_degree) ==
7046 this->data->data.front().fe_degree)
7048 proposed_dof_comp = dof_no;
7049 proposed_fe_comp = first_selected_component;
7052 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
7054 for (
unsigned int nf = 0;
7057 if (this->matrix_free
7058 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7060 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7062 proposed_dof_comp = no;
7063 proposed_fe_comp = nf;
7067 this->mapping_data->descriptor[this->active_quad_index]
7069 proposed_quad_comp = this->quad_no;
7071 for (
unsigned int no = 0;
7076 .descriptor[this->active_quad_index]
7077 .n_q_points == n_q_points)
7079 proposed_quad_comp = no;
7086 if (proposed_dof_comp != first_selected_component)
7087 message +=
"Wrong vector component selection:\n";
7089 message +=
"Wrong quadrature formula selection:\n";
7090 message +=
" Did you mean FEEvaluation<dim,";
7094 message +=
",Number>(data";
7103 std::string correct_pos;
7104 if (proposed_dof_comp != dof_no)
7105 correct_pos =
" ^ ";
7108 if (proposed_quad_comp != this->quad_no)
7109 correct_pos +=
" ^ ";
7112 if (proposed_fe_comp != first_selected_component)
7113 correct_pos +=
" ^\n";
7115 correct_pos +=
" \n";
7121 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7122 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7123 message +=
"Wrong template arguments:\n";
7124 message +=
" Did you mean FEEvaluation<dim,";
7129 message +=
",Number>(data";
7137 std::string correct_pos;
7138 if (this->data->data.front().fe_degree !=
7139 static_cast<unsigned int>(fe_degree))
7143 if (proposed_n_q_points_1d != n_q_points_1d)
7144 correct_pos +=
" ^\n";
7146 correct_pos +=
" \n";
7147 message +=
" " + correct_pos;
7149 Assert(
static_cast<unsigned int>(fe_degree) ==
7150 this->data->data.front().fe_degree &&
7151 n_q_points == this->n_quadrature_points,
7157 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7168 typename VectorizedArrayType>
7177 Assert(this->mapped_geometry ==
nullptr,
7178 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7179 " Integer indexing is not possible"));
7180 if (this->mapped_geometry !=
nullptr)
7189 const unsigned int offsets =
7190 this->mapping_data->data_index_offsets[
cell_index];
7191 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7192 this->J_value = &this->mapping_data->JxW_values[offsets];
7193 this->jacobian_gradients =
7194 this->mapping_data->jacobian_gradients[0].data() + offsets;
7199 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7200 for (; i < VectorizedArrayType::size(); ++i)
7203 if (this->mapping_data->quadrature_points.empty() ==
false)
7204 this->quadrature_points =
7205 &this->mapping_data->quadrature_points
7206 [this->mapping_data->quadrature_point_offsets[this->cell]];
7209 this->is_reinitialized =
true;
7210 this->dof_values_initialized =
false;
7211 this->values_quad_initialized =
false;
7212 this->gradients_quad_initialized =
false;
7213 this->hessians_quad_initialized =
false;
7224 typename VectorizedArrayType>
7231 VectorizedArrayType>
::
7232 reinit(
const std::array<
unsigned int, VectorizedArrayType::size()> &cell_ids)
7238 this->cell_ids = cell_ids;
7243 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7257 if (this->mapped_geometry ==
nullptr)
7258 this->mapped_geometry =
7259 std::make_shared<internal::MatrixFreeFunctions::
7260 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7262 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7264 auto &this_jacobian_data = mapping_storage.jacobians[0];
7265 auto &this_J_value_data = mapping_storage.JxW_values;
7266 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7267 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7271 if (this->mapping_data->jacobians[0].size() > 0)
7272 this_jacobian_data.resize_fast(2);
7274 if (this->mapping_data->JxW_values.size() > 0)
7275 this_J_value_data.resize_fast(1);
7277 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7278 this_jacobian_gradients_data.resize_fast(1);
7280 if (this->mapping_data->quadrature_points.size() > 0)
7281 this_quadrature_points_data.resize_fast(1);
7285 if (this->mapping_data->jacobians[0].size() > 0)
7286 this_jacobian_data.resize_fast(this->n_quadrature_points);
7288 if (this->mapping_data->JxW_values.size() > 0)
7289 this_J_value_data.resize_fast(this->n_quadrature_points);
7291 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7292 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7294 if (this->mapping_data->quadrature_points.size() > 0)
7295 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7299 this->jacobian = this_jacobian_data.data();
7300 this->J_value = this_J_value_data.data();
7301 this->jacobian_gradients = this_jacobian_gradients_data.data();
7305 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7312 const unsigned int cell_batch_index =
7314 const unsigned int offsets =
7315 this->mapping_data->data_index_offsets[cell_batch_index];
7316 const unsigned int lane =
cell_index % VectorizedArrayType::size();
7318 if (this->cell_type <=
7322 const unsigned int q = 0;
7324 if (this->mapping_data->JxW_values.size() > 0)
7325 this_J_value_data[q][v] =
7326 this->mapping_data->JxW_values[offsets + q][lane];
7328 if (this->mapping_data->jacobians[0].size() > 0)
7329 for (
unsigned int q = 0; q < 2; ++q)
7330 for (
unsigned int i = 0; i < dim; ++i)
7331 for (
unsigned int j = 0; j < dim; ++j)
7332 this_jacobian_data[q][i][j][v] =
7333 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7335 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7336 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7337 for (
unsigned int j = 0; j < dim; ++j)
7338 this_jacobian_gradients_data[q][i][j][v] =
7340 ->jacobian_gradients[0][offsets + q][i][j][lane];
7342 if (this->mapping_data->quadrature_points.size() > 0)
7343 for (
unsigned int i = 0; i < dim; ++i)
7344 this_quadrature_points_data[q][i][v] =
7345 this->mapping_data->quadrature_points
7347 ->quadrature_point_offsets[cell_batch_index] +
7353 const auto cell_type =
7357 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
7359 const unsigned int q_src =
7365 if (this->mapping_data->JxW_values.size() > 0)
7366 this_J_value_data[q][v] =
7367 this->mapping_data->JxW_values[offsets + q_src][lane];
7369 if (this->mapping_data->jacobians[0].size() > 0)
7370 for (
unsigned int i = 0; i < dim; ++i)
7371 for (
unsigned int j = 0; j < dim; ++j)
7372 this_jacobian_data[q][i][j][v] =
7374 ->jacobians[0][offsets + q_src][i][j][lane];
7376 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7377 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7378 for (
unsigned int j = 0; j < dim; ++j)
7379 this_jacobian_gradients_data[q][i][j][v] =
7381 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7383 if (this->mapping_data->quadrature_points.size() > 0)
7392 this->mapping_data->quadrature_points
7394 ->quadrature_point_offsets[cell_batch_index] +
7398 this->mapping_data->jacobians[0][offsets + 1];
7400 for (
unsigned int d = 0;
d < dim; ++
d)
7403 static_cast<Number
>(
7404 this->descriptor->quadrature.point(q)[
d]);
7406 for (
unsigned int d = 0;
d < dim; ++
d)
7407 for (
unsigned int e = 0;
e < dim; ++
e)
7410 static_cast<Number
>(
7411 this->descriptor->quadrature.point(q)[
e]);
7413 for (
unsigned int i = 0; i < dim; ++i)
7414 this_quadrature_points_data[q][i][v] = point[i][lane];
7419 for (
unsigned int i = 0; i < dim; ++i)
7420 this_quadrature_points_data[q][i][v] =
7421 this->mapping_data->quadrature_points
7423 ->quadrature_point_offsets[cell_batch_index] +
7432 this->is_reinitialized =
true;
7433 this->dof_values_initialized =
false;
7434 this->values_quad_initialized =
false;
7435 this->gradients_quad_initialized =
false;
7436 this->hessians_quad_initialized =
false;
7447 typename VectorizedArrayType>
7448template <
bool level_dof_access>
7455 VectorizedArrayType>
::
7458 Assert(this->matrix_free ==
nullptr,
7459 ExcMessage(
"Cannot use initialization from cell iterator if "
7460 "initialized from MatrixFree object. Use variant for "
7461 "on the fly computation with arguments as for FEValues "
7464 this->mapped_geometry->reinit(
7466 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7467 if (level_dof_access)
7468 cell->get_mg_dof_indices(this->local_dof_indices);
7470 cell->get_dof_indices(this->local_dof_indices);
7473 this->is_reinitialized =
true;
7484 typename VectorizedArrayType>
7491 VectorizedArrayType>
::
7494 Assert(this->matrix_free == 0,
7495 ExcMessage(
"Cannot use initialization from cell iterator if "
7496 "initialized from MatrixFree object. Use variant for "
7497 "on the fly computation with arguments as for FEValues "
7500 this->mapped_geometry->reinit(cell);
7503 this->is_reinitialized =
true;
7514 typename VectorizedArrayType>
7521 VectorizedArrayType>::evaluate(
const bool evaluate_values,
7522 const bool evaluate_gradients,
7523 const bool evaluate_hessians)
7526 Assert(this->dof_values_initialized ==
true,
7529 evaluate(this->values_dofs,
7541 typename VectorizedArrayType>
7548 VectorizedArrayType>::
7552 Assert(this->dof_values_initialized ==
true,
7555 evaluate(this->values_dofs, evaluation_flags);
7565 typename VectorizedArrayType>
7572 VectorizedArrayType>::evaluate(
const VectorizedArrayType
7574 const bool evaluate_values,
7575 const bool evaluate_gradients,
7576 const bool evaluate_hessians)
7585 evaluate(values_array, flag);
7595 typename VectorizedArrayType>
7602 VectorizedArrayType>::
7603 evaluate(
const VectorizedArrayType * values_array,
7606 const bool hessians_on_general_cells =
7610 if (hessians_on_general_cells)
7616 evaluate(n_components, evaluation_flag_actual, values_array, *
this);
7622 evaluation_flag_actual,
7623 const_cast<VectorizedArrayType *
>(values_array),
7629 this->values_quad_initialized =
true;
7631 this->gradients_quad_initialized =
true;
7633 this->hessians_quad_initialized =
true;
7644 typename VectorizedArrayType>
7645template <
typename VectorType>
7653 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
7654 const bool evaluate_values,
7655 const bool evaluate_gradients,
7656 const bool evaluate_hessians)
7665 gather_evaluate(input_vector, flag);
7674 template <
typename Number,
7675 typename VectorizedArrayType,
7676 typename VectorType,
7677 typename EvaluatorType,
7678 typename std::enable_if<internal::has_begin<VectorType> &&
7680 VectorType>::type * =
nullptr>
7681 VectorizedArrayType *
7682 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
7688 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7689 const auto & dof_info = fe_eval.get_dof_info();
7696 if (std::is_same<typename VectorType::value_type, Number>::value &&
7700 interleaved_contiguous &&
7701 reinterpret_cast<std::size_t
>(
7705 [cell * VectorizedArrayType::size()]) %
7706 sizeof(VectorizedArrayType) ==
7709 return reinterpret_cast<VectorizedArrayType *
>(
7713 [cell * VectorizedArrayType::size()] +
7715 [fe_eval.get_active_fe_index()]
7716 [fe_eval.get_first_selected_component()] *
7717 VectorizedArrayType::size());
7726 template <
typename Number,
7727 typename VectorizedArrayType,
7728 typename VectorType,
7729 typename EvaluatorType,
7730 typename std::enable_if<!internal::has_begin<VectorType> ||
7732 VectorType>::type * =
nullptr>
7733 VectorizedArrayType *
7734 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7747 typename VectorizedArrayType>
7748template <
typename VectorType>
7755 VectorizedArrayType>::
7756 gather_evaluate(
const VectorType & input_vector,
7759 const VectorizedArrayType *src_ptr =
7760 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7761 *
this, input_vector);
7762 if (src_ptr !=
nullptr)
7763 evaluate(src_ptr, evaluation_flag);
7766 this->read_dof_values(input_vector);
7767 evaluate(this->begin_dof_values(), evaluation_flag);
7778 typename VectorizedArrayType>
7785 VectorizedArrayType>::integrate(
const bool integrate_values,
7786 const bool integrate_gradients)
7788 integrate(integrate_values, integrate_gradients, this->values_dofs);
7791 this->dof_values_initialized =
true;
7802 typename VectorizedArrayType>
7809 VectorizedArrayType>::
7812 integrate(integration_flag, this->values_dofs);
7815 this->dof_values_initialized =
true;
7826 typename VectorizedArrayType>
7833 VectorizedArrayType>::integrate(
const bool integrate_values,
7834 const bool integrate_gradients,
7835 VectorizedArrayType *values_array)
7841 integrate(flag, values_array);
7851 typename VectorizedArrayType>
7858 VectorizedArrayType>::
7860 VectorizedArrayType * values_array,
7861 const bool sum_into_values_array)
7865 Assert(this->values_quad_submitted ==
true,
7868 Assert(this->gradients_quad_submitted ==
true,
7871 Assert(this->hessians_quad_submitted ==
true,
7874 Assert(this->matrix_free !=
nullptr ||
7875 this->mapped_geometry->is_initialized(),
7881 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7882 "EvaluationFlags::hessians are supported."));
7888 unsigned int size = n_components * dim * n_q_points;
7891 for (
unsigned int i = 0; i < size; ++i)
7892 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7896 for (
unsigned int i = 0; i < size; ++i)
7897 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7906 integration_flag_actual,
7909 sum_into_values_array);
7915 integration_flag_actual,
7918 sum_into_values_array);
7922 this->dof_values_initialized =
true;
7933 typename VectorizedArrayType>
7934template <
typename VectorType>
7942 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
7943 const bool integrate_gradients,
7944 VectorType &destination)
7951 integrate_scatter(flag, destination);
7961 typename VectorizedArrayType>
7962template <
typename VectorType>
7969 VectorizedArrayType>::
7971 VectorType & destination)
7973 VectorizedArrayType *dst_ptr =
7974 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7975 *
this, destination);
7976 if (dst_ptr !=
nullptr)
7977 integrate(integration_flag, dst_ptr,
true);
7980 integrate(integration_flag, this->begin_dof_values());
7981 this->distribute_local_to_global(destination);
7992 typename VectorizedArrayType>
7999 VectorizedArrayType>::dof_indices()
const
8001 return {0U, dofs_per_cell};
8015 typename VectorizedArrayType>
8021 VectorizedArrayType>::
8024 const bool is_interior_face,
8025 const unsigned int dof_no,
8026 const unsigned int quad_no,
8027 const unsigned int first_selected_component,
8028 const unsigned int active_fe_index,
8029 const unsigned int active_quad_index,
8030 const unsigned int face_type)
8031 : BaseClass(matrix_free,
8033 first_selected_component,
8041 , dofs_per_component(this->data->dofs_per_component_on_cell)
8042 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
8043 , n_q_points(this->n_quadrature_points)
8053 typename VectorizedArrayType>
8059 VectorizedArrayType>::
8062 const std::pair<unsigned int, unsigned int> & range,
8063 const bool is_interior_face,
8064 const unsigned int dof_no,
8065 const unsigned int quad_no,
8066 const unsigned int first_selected_component)
8071 first_selected_component,
8072 matrix_free.get_face_active_fe_index(range,
8074 numbers::invalid_unsigned_int,
8075 matrix_free.get_face_info(range.
first).face_type)
8085 typename VectorizedArrayType>
8092 VectorizedArrayType>
::reinit(
const unsigned int face_index)
8094 Assert(this->mapped_geometry ==
nullptr,
8095 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8096 " Integer indexing is not possible"));
8097 if (this->mapped_geometry !=
nullptr)
8100 this->cell = face_index;
8101 this->dof_access_index =
8102 this->is_interior_face() ?
8110 this->matrix_free->get_task_info().boundary_partition_data.back())
8111 Assert(this->is_interior_face(),
8113 "Boundary faces do not have a neighbor. When looping over "
8114 "boundary faces use FEFaceEvaluation with the parameter "
8115 "is_interior_face set to true. "));
8117 this->reinit_face(this->matrix_free->
get_face_info(face_index));
8122 this->face_ids[i] = face_index * VectorizedArrayType::size() + i;
8123 for (; i < VectorizedArrayType::size(); ++i)
8126 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
8127 const unsigned int offsets =
8128 this->mapping_data->data_index_offsets[face_index];
8129 this->J_value = &this->mapping_data->JxW_values[offsets];
8130 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8132 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
8133 this->normal_x_jacobian =
8135 ->normals_times_jacobians[!this->is_interior_face()][offsets];
8136 this->jacobian_gradients =
8137 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8140 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
8143 this->mapping_data->quadrature_point_offsets.size());
8145 this->mapping_data->quadrature_points.data() +
8146 this->mapping_data->quadrature_point_offsets[this->cell];
8150 this->is_reinitialized =
true;
8151 this->dof_values_initialized =
false;
8152 this->values_quad_initialized =
false;
8153 this->gradients_quad_initialized =
false;
8154 this->hessians_quad_initialized =
false;
8165 typename VectorizedArrayType>
8173 const unsigned int face_number)
8179 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8183 Assert(this->mapped_geometry ==
nullptr,
8184 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8185 " Integer indexing is not possible"));
8186 if (this->mapped_geometry !=
nullptr)
8191 .faces_by_cells_type[
cell_index][face_number];
8194 this->dof_access_index =
8197 constexpr unsigned int n_lanes = VectorizedArrayType::size();
8199 if (this->is_interior_face() ==
false)
8204 for (
unsigned int i = 0; i < n_lanes; ++i)
8207 const unsigned int cell_this =
cell_index * n_lanes + i;
8209 unsigned int face_index =
8214 this->face_ids[i] = face_index;
8219 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
8220 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
8227 auto cell_m = faces.cells_interior[face_index % n_lanes];
8228 auto cell_p = faces.cells_exterior[face_index % n_lanes];
8230 const bool face_identifies_as_interior = cell_m != cell_this;
8232 Assert(cell_m == cell_this || cell_p == cell_this,
8236 if (face_identifies_as_interior)
8238 this->cell_ids[i] = cell_m;
8239 this->face_numbers[i] = faces.interior_face_no;
8243 this->cell_ids[i] = cell_p;
8244 this->face_numbers[i] = faces.exterior_face_no;
8247 const bool orientation_interior_face = faces.face_orientation >= 8;
8248 unsigned int face_orientation = faces.face_orientation % 8;
8249 if (face_identifies_as_interior != orientation_interior_face)
8251 constexpr std::array<std::uint8_t, 8> table{
8252 {0, 1, 2, 3, 6, 5, 4, 7}};
8253 face_orientation = table[face_orientation];
8255 this->face_orientations[i] = face_orientation;
8260 this->face_orientations[0] = 0;
8261 this->face_numbers[0] = face_number;
8262 for (
unsigned int i = 0; i < n_lanes; ++i)
8263 this->cell_ids[i] =
cell_index * n_lanes + i;
8264 for (
unsigned int i = 0; i < n_lanes; ++i)
8271 const unsigned int offsets =
8273 .face_data_by_cells[this->quad_no]
8278 .face_data_by_cells[this->quad_no]
8279 .JxW_values.size());
8281 .face_data_by_cells[this->quad_no]
8282 .JxW_values[offsets];
8284 .face_data_by_cells[this->quad_no]
8285 .normal_vectors[offsets];
8287 .face_data_by_cells[this->quad_no]
8288 .jacobians[!this->is_interior_face()][offsets];
8289 this->normal_x_jacobian =
8291 .face_data_by_cells[this->quad_no]
8292 .normals_times_jacobians[!this->is_interior_face()][offsets];
8293 this->jacobian_gradients =
8294 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8298 .face_data_by_cells[this->quad_no]
8299 .quadrature_point_offsets.empty() ==
false)
8301 const unsigned int index =
8305 .face_data_by_cells[this->quad_no]
8306 .quadrature_point_offsets.size());
8308 .face_data_by_cells[this->quad_no]
8309 .quadrature_points.data() +
8311 .face_data_by_cells[this->quad_no]
8312 .quadrature_point_offsets[
index];
8316 this->is_reinitialized =
true;
8317 this->dof_values_initialized =
false;
8318 this->values_quad_initialized =
false;
8319 this->gradients_quad_initialized =
false;
8320 this->hessians_quad_initialized =
false;
8331 typename VectorizedArrayType>
8338 VectorizedArrayType>::evaluate(
const bool evaluate_values,
8339 const bool evaluate_gradients)
8345 evaluate(this->values_dofs, evaluate_values, evaluate_gradients);
8355 typename VectorizedArrayType>
8362 VectorizedArrayType>::
8369 evaluate(this->values_dofs, evaluation_flag);
8379 typename VectorizedArrayType>
8386 VectorizedArrayType>::evaluate(
const VectorizedArrayType
8388 const bool evaluate_values,
8389 const bool evaluate_gradients)
8396 evaluate(values_array, flag);
8406 typename VectorizedArrayType>
8413 VectorizedArrayType>::
8414 evaluate(
const VectorizedArrayType * values_array,
8417 Assert((evaluation_flag &
8420 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8421 "and EvaluationFlags::hessians are supported."));
8423 const bool hessians_on_general_cells =
8427 if (hessians_on_general_cells)
8432 template run<fe_degree, n_q_points_1d>(n_components,
8433 evaluation_flag_actual,
8438 n_components, evaluation_flag_actual, values_array, *
this);
8442 this->values_quad_initialized =
true;
8444 this->gradients_quad_initialized =
true;
8446 this->hessians_quad_initialized =
true;
8457 typename VectorizedArrayType>
8464 VectorizedArrayType>::
8467 integrate(integration_flag, this->values_dofs);
8470 this->dof_values_initialized =
true;
8481 typename VectorizedArrayType>
8488 VectorizedArrayType>::integrate(
const bool integrate_values,
8489 const bool integrate_gradients)
8491 integrate(integrate_values, integrate_gradients, this->values_dofs);
8494 this->dof_values_initialized =
true;
8505 typename VectorizedArrayType>
8512 VectorizedArrayType>::integrate(
const bool integrate_values,
8513 const bool integrate_gradients,
8522 integrate(flag, values_array);
8532 typename VectorizedArrayType>
8539 VectorizedArrayType>::
8541 VectorizedArrayType * values_array)
8543 Assert((integration_flag &
8546 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8547 "and EvaluationFlags::hessians are supported."));
8553 unsigned int size = n_components * dim * n_q_points;
8556 for (
unsigned int i = 0; i < size; ++i)
8557 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
8561 for (
unsigned int i = 0; i < size; ++i)
8562 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
8569 template run<fe_degree, n_q_points_1d>(n_components,
8570 integration_flag_actual,
8575 n_components, integration_flag_actual, values_array, *
this);
8585 typename VectorizedArrayType>
8586template <
typename VectorType>
8594 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
8595 const bool evaluate_values,
8596 const bool evaluate_gradients)
8603 gather_evaluate(input_vector, flag);
8613 typename VectorizedArrayType>
8614template <
typename VectorType>
8621 VectorizedArrayType>::
8622 gather_evaluate(
const VectorType & input_vector,
8625 Assert((evaluation_flag &
8628 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8629 "and EvaluationFlags::hessians are supported."));
8631 const auto shared_vector_data = internal::get_shared_vector_data(
8633 this->dof_access_index ==
8635 this->active_fe_index,
8638 if (this->data->data.front().fe_degree > 0 &&
8639 fast_evaluation_supported(this->data->data.front().fe_degree,
8640 this->data->data.front().n_q_points_1d) &&
8643 typename VectorType::value_type,
8644 VectorizedArrayType>::
8645 supports(evaluation_flag,
8647 internal::get_beginning<typename VectorType::value_type>(
8656 typename VectorType::value_type,
8657 VectorizedArrayType>::template
run<fe_degree,
8661 internal::get_beginning<typename VectorType::value_type>(
8670 typename VectorType::value_type,
8671 VectorizedArrayType>::evaluate(n_components,
8673 internal::get_beginning<
8674 typename VectorType::value_type>(
8682 this->read_dof_values(input_vector);
8683 this->evaluate(evaluation_flag);
8688 this->values_quad_initialized =
true;
8690 this->gradients_quad_initialized =
true;
8692 this->hessians_quad_initialized =
true;
8703 typename VectorizedArrayType>
8704template <
typename VectorType>
8712 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8713 const bool integrate_gradients,
8714 VectorType &destination)
8721 integrate_scatter(flag, destination);
8731 typename VectorizedArrayType>
8732template <
typename VectorType>
8739 VectorizedArrayType>::
8741 VectorType & destination)
8743 Assert((this->dof_access_index ==
8745 this->is_interior_face() ==
false) ==
false,
8748 const auto shared_vector_data = internal::get_shared_vector_data(
8750 this->dof_access_index ==
8752 this->active_fe_index,
8755 if (this->data->data.front().fe_degree > 0 &&
8756 fast_evaluation_supported(this->data->data.front().fe_degree,
8757 this->data->data.front().n_q_points_1d) &&
8760 typename VectorType::value_type,
8761 VectorizedArrayType>::
8762 supports(integration_flag,
8764 internal::get_beginning<typename VectorType::value_type>(
8773 typename VectorType::value_type,
8774 VectorizedArrayType>::template
run<fe_degree,
8778 internal::get_beginning<typename VectorType::value_type>(
8787 typename VectorType::value_type,
8788 VectorizedArrayType>::integrate(n_components,
8790 internal::get_beginning<
8791 typename VectorType::value_type>(
8799 integrate(integration_flag);
8800 this->distribute_local_to_global(destination);
8811 typename VectorizedArrayType>
8818 VectorizedArrayType>::dof_indices()
const
8820 return {0U, dofs_per_cell};
8830 typename VectorizedArrayType>
8837 VectorizedArrayType>::
8838 fast_evaluation_supported(
const unsigned int given_degree,
8839 const unsigned int give_n_q_points_1d)
8841 return fe_degree == -1 ?
8854 typename VectorizedArrayType>
8861 VectorizedArrayType>::
8862 fast_evaluation_supported(
const unsigned int given_degree,
8863 const unsigned int give_n_q_points_1d)
8865 return fe_degree == -1 ?
FEEvaluationAccess(const FEEvaluationAccess &other)
hessian_type get_hessian(unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
value_type integrate_value() const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
value_type get_laplacian(const unsigned int q_point) const
void submit_value(const gradient_type val_in, const unsigned int q_point)
value_type get_normal_derivative(const unsigned int q_point) const
void submit_normal_derivative(const gradient_type grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
void submit_gradient(const value_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
VectorizedArrayType value_type
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< 1 > &mapping, const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< 1, VectorizedArrayType, is_face > *other)
value_type get_divergence(const unsigned int q_point) const
value_type get_dof_value(const unsigned int dof) const
FEEvaluationAccess(const MatrixFree< 1, 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=true, 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)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
value_type get_dof_value(const unsigned int dof) const
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
FEEvaluationAccess(const FEEvaluationAccess &other)
FEEvaluationAccess(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=true, 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)
value_type get_laplacian(const unsigned int q_point) const
hessian_type get_hessian(unsigned int q_point) const
FEEvaluationAccess(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)
VectorizedArrayType value_type
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
value_type get_normal_derivative(const unsigned int q_point) const
value_type integrate_value() const
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess(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 dofs_per_cell, const unsigned int n_q_points, const bool is_interior_face=true, 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)
FEEvaluationAccess(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)
FEEvaluationAccess(const FEEvaluationAccess &other)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
Tensor< 3, dim, VectorizedArrayType > get_hessian(const unsigned int q_point) const
VectorizedArrayType get_divergence(const unsigned int q_point) const
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 1, dim, Tensor< 1, dim, VectorizedArrayType > > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_value(const Tensor< 1, dim, VectorizedArrayType > val_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
Tensor< 1, n_components_, VectorizedArrayType > value_type
FEEvaluationAccess(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)
static constexpr unsigned int dimension
Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > gradient_type
FEEvaluationAccess(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=true, 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)
FEEvaluationAccess(const FEEvaluationAccess &other)
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
static constexpr unsigned int n_components
value_type get_dof_value(const unsigned int dof) const
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > get_hessian(const unsigned int q_point) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
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 distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
static constexpr unsigned int dimension
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
VectorizedArrayType get_divergence(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)
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
static constexpr unsigned int n_components
gradient_type get_hessian_diagonal(const unsigned int q_point) const
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< VectorizedArrayType::size()> &mask, const bool apply_constraints=true) const
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) 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< VectorizedArrayType::size()> &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
const unsigned int quad_no
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const unsigned int active_fe_index
std::array< unsigned int, n_lanes > cell_ids
const unsigned int active_quad_index
bool is_interior_face() const
FEEvaluationData & operator=(const FEEvaluationData &other)
const unsigned int first_selected_component
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
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
void integrate(const bool integrate_values, const bool integrate_gradients)
void reinit(const std::array< unsigned int, VectorizedArrayType::size()> &cell_ids)
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)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
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)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
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
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
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
FEEvaluation(const FEEvaluation &other)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
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_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
static constexpr unsigned int tensor_dofs_per_cell
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)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
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 integrate(const EvaluationFlags::EvaluationFlags integration_flag)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array)
static constexpr unsigned int static_dofs_per_component
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
static constexpr unsigned int static_n_q_points
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
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate(const bool integrate_values, const bool integrate_gradients)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_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)
EvaluationFlags
The EvaluationFlags enum.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
static const unsigned int invalid_unsigned_int
boost::integer_range< IncrementableType > iota_view
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &eval, const bool sum_into_values_array=false)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, 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::vector< unsigned int > > component_dof_indices_offset
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
std::vector< unsigned int > component_to_base_index
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
std::vector< QuadratureDescriptor > descriptor
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)