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 " +
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 " +
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;
5752 const VectorizedArrayType inv_det =
5753 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5754 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5755 this->jacobian[0][2][2];
5758 for (
unsigned int comp = 0; comp < n_components; ++comp)
5759 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5760 jac[comp][comp] * inv_det;
5767 this->jacobian[q_point] :
5775 const VectorizedArrayType inv_det =
5776 (is_face && dim == 2 && this->get_face_no() < 2) ?
5780 for (
unsigned int comp = 0; comp < n_components; ++comp)
5783 this->values_quad[q_point] * jac[comp][0] * inv_det;
5784 for (
unsigned int e = 1;
e < dim; ++
e)
5786 this->values_quad[
e * nqp + q_point] * jac[comp][
e] * inv_det;
5794 return BaseClass::get_value(q_point);
5798template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5803 if (this->data->element_type ==
5808 Assert(this->gradients_quad_initialized ==
true,
5813 Assert(this->jacobian !=
nullptr,
5815 "update_gradients"));
5816 const std::size_t nqp = this->n_quadrature_points;
5826 const VectorizedArrayType inv_det =
5827 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5828 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5829 this->jacobian[0][2][2];
5832 for (
unsigned int d = 0;
d < dim; ++
d)
5833 for (
unsigned int comp = 0; comp < n_components; ++comp)
5835 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
5836 inv_t_jac[
d][
d] * jac[comp][comp] * inv_det;
5846 const VectorizedArrayType inv_det =
5847 (is_face && dim == 2 && this->get_face_no() < 2) ?
5851 VectorizedArrayType tmp;
5853 for (
unsigned int comp = 0; comp < n_components; ++comp)
5854 for (
unsigned int d = 0;
d < dim; ++
d)
5857 for (
unsigned int f = 0; f < dim; ++f)
5858 for (
unsigned int e = 0;
e < dim; ++
e)
5859 tmp += jac[comp][f] * inv_t_jac[
d][
e] * inv_det *
5860 this->gradients_quad[(f * dim +
e) * nqp + q_point];
5862 grad_out[comp][
d] = tmp;
5876 return BaseClass::get_gradient(q_point);
5882template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5888 Assert(this->gradients_quad_initialized ==
true,
5892 Assert(this->jacobian !=
nullptr,
5894 "update_gradients"));
5896 VectorizedArrayType divergence;
5897 const std::size_t nqp = this->n_quadrature_points;
5899 if (this->data->element_type ==
5906 const VectorizedArrayType inv_det =
5907 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5908 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5909 this->jacobian[0][2][2];
5912 divergence = this->gradients_quad[q_point] * inv_det;
5913 for (
unsigned int d = 1;
d < dim; ++
d)
5915 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
5921 const VectorizedArrayType inv_det =
5922 (is_face && dim == 2 && this->get_face_no() < 2) ?
5927 divergence = this->gradients_quad[q_point] * inv_det;
5928 for (
unsigned int d = 1;
d < dim; ++
d)
5930 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
5944 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
5945 for (
unsigned int d = 1;
d < dim; ++
d)
5946 divergence += this->gradients_quad[(dim *
d +
d) * nqp + q_point] *
5947 this->jacobian[0][
d][
d];
5954 this->jacobian[q_point] :
5956 divergence = jac[0][0] * this->gradients_quad[q_point];
5957 for (
unsigned int e = 1;
e < dim; ++
e)
5958 divergence += jac[0][
e] * this->gradients_quad[
e * nqp + q_point];
5959 for (
unsigned int d = 1;
d < dim; ++
d)
5960 for (
unsigned int e = 0;
e < dim; ++
e)
5962 jac[
d][
e] * this->gradients_quad[(
d * dim +
e) * nqp + q_point];
5970template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5976 const auto grad = get_gradient(q_point);
5977 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5978 VectorizedArrayType half = Number(0.5);
5979 for (
unsigned int d = 0;
d < dim; ++
d)
5980 symmetrized[
d] = grad[
d][
d];
5986 symmetrized[2] = grad[0][1] + grad[1][0];
5987 symmetrized[2] *= half;
5990 symmetrized[3] = grad[0][1] + grad[1][0];
5991 symmetrized[3] *= half;
5992 symmetrized[4] = grad[0][2] + grad[2][0];
5993 symmetrized[4] *= half;
5994 symmetrized[5] = grad[1][2] + grad[2][1];
5995 symmetrized[5] *= half;
6005template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6007 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6009 const unsigned int q_point)
const
6013 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6019 "Computing the curl in 1d is not a useful operation"));
6022 curl[0] = grad[1][0] - grad[0][1];
6025 curl[0] = grad[2][1] - grad[1][2];
6026 curl[1] = grad[0][2] - grad[2][0];
6027 curl[2] = grad[1][0] - grad[0][1];
6037template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6042 return BaseClass::get_hessian_diagonal(q_point);
6047template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6050 const unsigned int q_point)
const
6053 Assert(this->hessians_quad_initialized ==
true,
6057 return BaseClass::get_hessian(q_point);
6061template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6065 const unsigned int q_point)
6067 if (this->data->element_type ==
6072 Assert(this->J_value !=
nullptr,
6077 this->values_quad_submitted =
true;
6080 const std::size_t nqp = this->n_quadrature_points;
6086 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6088 for (
unsigned int comp = 0; comp < n_components; ++comp)
6089 this->values_quad[comp * nqp + q_point] =
6090 val_in[comp] * weight * jac[comp][comp];
6097 this->jacobian[q_point] :
6106 const VectorizedArrayType fac =
6108 this->quadrature_weights[q_point] :
6110 this->J_value[q_point] :
6111 this->J_value[0] * this->quadrature_weights[q_point]) *
6112 ((dim == 2 && this->get_face_no() < 2) ?
6117 for (
unsigned int comp = 0; comp < n_components; ++comp)
6119 this->values_quad[comp * nqp + q_point] =
6120 val_in[0] * jac[0][comp] * fac;
6121 for (
unsigned int e = 1;
e < dim; ++
e)
6122 this->values_quad[comp * nqp + q_point] +=
6123 val_in[
e] * jac[
e][comp] * fac;
6130 BaseClass::submit_value(val_in, q_point);
6134template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6138 const unsigned int q_point)
6140 if (this->data->element_type ==
6149 Assert(this->J_value !=
nullptr,
6151 "update_gradients"));
6152 Assert(this->jacobian !=
nullptr,
6154 "update_gradients"));
6156 this->gradients_quad_submitted =
true;
6159 const std::size_t nqp = this->n_quadrature_points;
6167 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6168 for (
unsigned int d = 0;
d < dim; ++
d)
6169 for (
unsigned int comp = 0; comp < n_components; ++comp)
6170 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6171 grad_in[comp][
d] * inv_t_jac[
d][
d] * jac[comp][comp] * weight;
6182 const VectorizedArrayType fac =
6183 (!is_face) ? this->quadrature_weights[q_point] :
6184 this->J_value[0] * this->quadrature_weights[q_point] *
6185 ((dim == 2 && this->get_face_no() < 2) ?
6190 for (
unsigned int comp = 0; comp < n_components; ++comp)
6191 for (
unsigned int d = 0;
d < dim; ++
d)
6193 VectorizedArrayType tmp = 0;
6194 for (
unsigned int f = 0; f < dim; ++f)
6195 for (
unsigned int e = 0;
e < dim; ++
e)
6196 tmp += jac[f][comp] * inv_t_jac[
e][
d] * grad_in[f][
e] * fac;
6198 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = tmp;
6209 BaseClass::submit_gradient(grad_in, q_point);
6215template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6220 const unsigned int q_point)
6222 if (this->data->element_type ==
6232 BaseClass::submit_gradient(grad_in, q_point);
6238template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6242 const unsigned int q_point)
6248 Assert(this->J_value !=
nullptr,
6250 "update_gradients"));
6251 Assert(this->jacobian !=
nullptr,
6253 "update_gradients"));
6255 this->gradients_quad_submitted =
true;
6258 const std::size_t nqp = this->n_quadrature_points;
6259 if (this->data->element_type ==
6268 const VectorizedArrayType fac =
6271 this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
6274 this->quadrature_weights[q_point] * div_in;
6276 for (
unsigned int d = 0;
d < dim; ++
d)
6278 this->gradients_quad[(dim *
d +
d) * nqp + q_point] = fac;
6279 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6281 this->gradients_quad[(dim *
d +
e) * nqp + q_point] =
6282 VectorizedArrayType();
6283 this->gradients_quad[(dim *
e +
d) * nqp + q_point] =
6284 VectorizedArrayType();
6299 const VectorizedArrayType fac =
6300 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6301 for (
unsigned int d = 0;
d < dim; ++
d)
6303 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6304 (fac * this->jacobian[0][
d][
d]);
6305 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6307 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6308 VectorizedArrayType();
6309 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6310 VectorizedArrayType();
6318 this->jacobian[q_point] :
6320 const VectorizedArrayType fac =
6322 this->J_value[q_point] :
6323 this->J_value[0] * this->quadrature_weights[q_point]) *
6325 for (
unsigned int d = 0;
d < dim; ++
d)
6327 for (
unsigned int e = 0;
e < dim; ++
e)
6328 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6337template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6342 const unsigned int q_point)
6345 this->data->element_type !=
6356 Assert(this->J_value !=
nullptr,
6358 "update_gradients"));
6359 Assert(this->jacobian !=
nullptr,
6361 "update_gradients"));
6363 this->gradients_quad_submitted =
true;
6366 const std::size_t nqp = this->n_quadrature_points;
6369 const VectorizedArrayType JxW =
6370 this->J_value[0] * this->quadrature_weights[q_point];
6371 for (
unsigned int d = 0;
d < dim; ++
d)
6372 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6374 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6375 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6377 const VectorizedArrayType
value =
6379 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6380 value * this->jacobian[0][
d][
d];
6381 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6382 value * this->jacobian[0][
e][
e];
6388 const VectorizedArrayType JxW =
6390 this->J_value[q_point] :
6391 this->J_value[0] * this->quadrature_weights[q_point];
6394 this->jacobian[q_point] :
6396 VectorizedArrayType weighted[dim][dim];
6397 for (
unsigned int i = 0; i < dim; ++i)
6399 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6400 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6402 const VectorizedArrayType
value =
6404 weighted[i][j] =
value;
6405 weighted[j][i] =
value;
6407 for (
unsigned int comp = 0; comp < dim; ++comp)
6408 for (
unsigned int d = 0;
d < dim; ++
d)
6410 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6411 for (
unsigned int e = 1;
e < dim; ++
e)
6412 new_val += jac[
e][
d] * weighted[comp][
e];
6413 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
6420template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6424 const unsigned int q_point)
6432 "Testing by the curl in 1d is not a useful operation"));
6435 grad[1][0] = curl[0];
6436 grad[0][1] = -curl[0];
6439 grad[2][1] = curl[0];
6440 grad[1][2] = -curl[0];
6441 grad[0][2] = curl[1];
6442 grad[2][0] = -curl[1];
6443 grad[1][0] = curl[2];
6444 grad[0][1] = -curl[2];
6449 submit_gradient(grad, q_point);
6456template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6460 const unsigned int dof_no,
6461 const unsigned int first_selected_component,
6462 const unsigned int quad_no,
6463 const unsigned int fe_degree,
6464 const unsigned int n_q_points,
6465 const bool is_interior_face,
6466 const unsigned int active_fe_index,
6467 const unsigned int active_quad_index,
6468 const unsigned int face_type)
6472 first_selected_component,
6484template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6491 const unsigned int first_selected_component,
6498 first_selected_component,
6504template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6513template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6525template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6528 const unsigned int dof)
const
6531 return this->values_dofs[dof];
6536template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6539 const unsigned int q_point)
const
6542 Assert(this->values_quad_initialized ==
true,
6546 return this->values_quad[q_point];
6551template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6554 const unsigned int q_point)
const
6560 Assert(this->gradients_quad_initialized ==
true,
6567 this->jacobian[q_point] :
6571 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6578template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6581 const unsigned int q_point)
const
6583 return get_gradient(q_point)[0];
6588template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6593 return BaseClass::get_normal_derivative(q_point)[0];
6598template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6601 const unsigned int q_point)
const
6603 return BaseClass::get_hessian(q_point)[0];
6608template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6613 return BaseClass::get_hessian_diagonal(q_point)[0];
6618template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6621 const unsigned int q_point)
const
6623 return BaseClass::get_laplacian(q_point)[0];
6628template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6634 this->dof_values_initialized =
true;
6637 this->values_dofs[dof] = val_in;
6642template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6645 const VectorizedArrayType val_in,
6646 const unsigned int q_point)
6653 this->values_quad_submitted =
true;
6658 const VectorizedArrayType JxW = this->J_value[q_point];
6659 this->values_quad[q_point] = val_in * JxW;
6663 const VectorizedArrayType JxW =
6664 this->J_value[0] * this->quadrature_weights[q_point];
6665 this->values_quad[q_point] = val_in * JxW;
6671template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6675 const unsigned int q_point)
6677 submit_value(val_in[0], q_point);
6682template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6686 const unsigned int q_point)
6688 submit_gradient(grad_in[0], q_point);
6693template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6696 const VectorizedArrayType grad_in,
6697 const unsigned int q_point)
6704 this->gradients_quad_submitted =
true;
6709 this->jacobian[q_point] :
6711 const VectorizedArrayType JxW =
6713 this->J_value[q_point] :
6714 this->J_value[0] * this->quadrature_weights[q_point];
6716 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6721template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6725 const unsigned int q_point)
6727 submit_gradient(grad_in[0][0], q_point);
6732template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6736 const unsigned int q_point)
6740 BaseClass::submit_normal_derivative(grad, q_point);
6745template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6749 const unsigned int q_point)
6751 BaseClass::submit_normal_derivative(grad_in, q_point);
6755template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6759 const unsigned int q_point)
6763 BaseClass::submit_hessian(hessian, q_point);
6767template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6768inline VectorizedArrayType
6772 return BaseClass::integrate_value()[0];
6785 typename VectorizedArrayType>
6791 VectorizedArrayType>::
6793 const unsigned int fe_no,
6794 const unsigned int quad_no,
6795 const unsigned int first_selected_component,
6796 const unsigned int active_fe_index,
6797 const unsigned int active_quad_index)
6798 : BaseClass(matrix_free,
6800 first_selected_component,
6807 , dofs_per_component(this->data->dofs_per_component_on_cell)
6808 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6809 , n_q_points(this->data->n_q_points)
6811 check_template_arguments(fe_no, 0);
6821 typename VectorizedArrayType>
6827 VectorizedArrayType>::
6829 const std::pair<unsigned int, unsigned int> & range,
6830 const unsigned int dof_no,
6831 const unsigned int quad_no,
6832 const unsigned int first_selected_component)
6836 first_selected_component,
6837 matrix_free.get_cell_active_fe_index(range))
6847 typename VectorizedArrayType>
6853 VectorizedArrayType>::
6858 const unsigned int first_selected_component)
6859 : BaseClass(mapping,
6863 first_selected_component,
6865 , dofs_per_component(this->data->dofs_per_component_on_cell)
6866 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6867 , n_q_points(this->data->n_q_points)
6879 typename VectorizedArrayType>
6885 VectorizedArrayType>::
6889 const unsigned int first_selected_component)
6894 first_selected_component,
6896 , dofs_per_component(this->data->dofs_per_component_on_cell)
6897 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6898 , n_q_points(this->data->n_q_points)
6910 typename VectorizedArrayType>
6916 VectorizedArrayType>::
6919 const unsigned int first_selected_component)
6920 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6922 other.mapped_geometry->get_quadrature(),
6923 other.mapped_geometry->get_fe_values().get_update_flags(),
6924 first_selected_component,
6926 , dofs_per_component(this->data->dofs_per_component_on_cell)
6927 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6928 , n_q_points(this->data->n_q_points)
6940 typename VectorizedArrayType>
6949 , dofs_per_component(this->data->dofs_per_component_on_cell)
6950 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6951 , n_q_points(this->data->n_q_points)
6963 typename VectorizedArrayType>
6969 VectorizedArrayType> &
6975 VectorizedArrayType>::operator=(
const FEEvaluation &other)
6977 BaseClass::operator=(other);
6989 typename VectorizedArrayType>
6996 VectorizedArrayType>::
6997 check_template_arguments(
const unsigned int dof_no,
6998 const unsigned int first_selected_component)
7001 (void)first_selected_component;
7004 this->data->dofs_per_component_on_cell > 0,
7006 "There is nothing useful you can do with an FEEvaluation object with "
7007 "FE_Nothing, i.e., without DoFs! If you have passed to "
7008 "MatrixFree::reinit() a collection of finite elements also containing "
7009 "FE_Nothing, please check - before creating FEEvaluation - the category "
7010 "of the current range by calling either "
7011 "MatrixFree::get_cell_range_category(range) or "
7012 "MatrixFree::get_face_range_category(range). The returned category "
7013 "is the index of the active FE, which you can use to exclude "
7020 static_cast<unsigned int>(fe_degree) !=
7021 this->data->data.front().fe_degree) ||
7022 n_q_points != this->n_quadrature_points)
7024 std::string message =
7025 "-------------------------------------------------------\n";
7026 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7027 message +=
" Called --> FEEvaluation<dim,";
7031 message +=
",Number>(data";
7047 if (
static_cast<unsigned int>(fe_degree) ==
7048 this->data->data.front().fe_degree)
7050 proposed_dof_comp = dof_no;
7051 proposed_fe_comp = first_selected_component;
7054 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
7056 for (
unsigned int nf = 0;
7059 if (this->matrix_free
7060 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7062 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7064 proposed_dof_comp = no;
7065 proposed_fe_comp = nf;
7069 this->mapping_data->descriptor[this->active_quad_index]
7071 proposed_quad_comp = this->quad_no;
7073 for (
unsigned int no = 0;
7078 .descriptor[this->active_quad_index]
7079 .n_q_points == n_q_points)
7081 proposed_quad_comp = no;
7088 if (proposed_dof_comp != first_selected_component)
7089 message +=
"Wrong vector component selection:\n";
7091 message +=
"Wrong quadrature formula selection:\n";
7092 message +=
" Did you mean FEEvaluation<dim,";
7096 message +=
",Number>(data";
7105 std::string correct_pos;
7106 if (proposed_dof_comp != dof_no)
7107 correct_pos =
" ^ ";
7110 if (proposed_quad_comp != this->quad_no)
7111 correct_pos +=
" ^ ";
7114 if (proposed_fe_comp != first_selected_component)
7115 correct_pos +=
" ^\n";
7117 correct_pos +=
" \n";
7123 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7124 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7125 message +=
"Wrong template arguments:\n";
7126 message +=
" Did you mean FEEvaluation<dim,";
7131 message +=
",Number>(data";
7139 std::string correct_pos;
7140 if (this->data->data.front().fe_degree !=
7141 static_cast<unsigned int>(fe_degree))
7145 if (proposed_n_q_points_1d != n_q_points_1d)
7146 correct_pos +=
" ^\n";
7148 correct_pos +=
" \n";
7149 message +=
" " + correct_pos;
7151 Assert(
static_cast<unsigned int>(fe_degree) ==
7152 this->data->data.front().fe_degree &&
7153 n_q_points == this->n_quadrature_points,
7159 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7170 typename VectorizedArrayType>
7179 Assert(this->mapped_geometry ==
nullptr,
7180 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7181 " Integer indexing is not possible"));
7182 if (this->mapped_geometry !=
nullptr)
7191 const unsigned int offsets =
7192 this->mapping_data->data_index_offsets[
cell_index];
7193 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7194 this->J_value = &this->mapping_data->JxW_values[offsets];
7195 this->jacobian_gradients =
7196 this->mapping_data->jacobian_gradients[0].data() + offsets;
7201 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7202 for (; i < VectorizedArrayType::size(); ++i)
7205 if (this->mapping_data->quadrature_points.empty() ==
false)
7207 &this->mapping_data->quadrature_points
7208 [this->mapping_data->quadrature_point_offsets[this->cell]];
7211 this->is_reinitialized =
true;
7212 this->dof_values_initialized =
false;
7213 this->values_quad_initialized =
false;
7214 this->gradients_quad_initialized =
false;
7215 this->hessians_quad_initialized =
false;
7226 typename VectorizedArrayType>
7233 VectorizedArrayType>
::
7234 reinit(
const std::array<
unsigned int, VectorizedArrayType::size()> &cell_ids)
7240 this->cell_ids = cell_ids;
7245 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7259 if (this->mapped_geometry ==
nullptr)
7260 this->mapped_geometry =
7261 std::make_shared<internal::MatrixFreeFunctions::
7262 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7264 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7266 auto &this_jacobian_data = mapping_storage.jacobians[0];
7267 auto &this_J_value_data = mapping_storage.JxW_values;
7268 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7269 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7273 if (this->mapping_data->jacobians[0].size() > 0)
7274 this_jacobian_data.resize_fast(2);
7276 if (this->mapping_data->JxW_values.size() > 0)
7277 this_J_value_data.resize_fast(1);
7279 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7280 this_jacobian_gradients_data.resize_fast(1);
7282 if (this->mapping_data->quadrature_points.size() > 0)
7283 this_quadrature_points_data.resize_fast(1);
7287 if (this->mapping_data->jacobians[0].size() > 0)
7288 this_jacobian_data.resize_fast(this->n_quadrature_points);
7290 if (this->mapping_data->JxW_values.size() > 0)
7291 this_J_value_data.resize_fast(this->n_quadrature_points);
7293 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7294 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7296 if (this->mapping_data->quadrature_points.size() > 0)
7297 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7301 this->jacobian = this_jacobian_data.data();
7302 this->J_value = this_J_value_data.data();
7303 this->jacobian_gradients = this_jacobian_gradients_data.data();
7307 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7314 const unsigned int cell_batch_index =
7316 const unsigned int offsets =
7317 this->mapping_data->data_index_offsets[cell_batch_index];
7318 const unsigned int lane =
cell_index % VectorizedArrayType::size();
7320 if (this->cell_type <=
7324 const unsigned int q = 0;
7326 if (this->mapping_data->JxW_values.size() > 0)
7327 this_J_value_data[q][v] =
7328 this->mapping_data->JxW_values[offsets + q][lane];
7330 if (this->mapping_data->jacobians[0].size() > 0)
7331 for (
unsigned int q = 0; q < 2; ++q)
7332 for (
unsigned int i = 0; i < dim; ++i)
7333 for (
unsigned int j = 0; j < dim; ++j)
7334 this_jacobian_data[q][i][j][v] =
7335 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7337 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7338 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7339 for (
unsigned int j = 0; j < dim; ++j)
7340 this_jacobian_gradients_data[q][i][j][v] =
7342 ->jacobian_gradients[0][offsets + q][i][j][lane];
7344 if (this->mapping_data->quadrature_points.size() > 0)
7345 for (
unsigned int i = 0; i < dim; ++i)
7346 this_quadrature_points_data[q][i][v] =
7347 this->mapping_data->quadrature_points
7349 ->quadrature_point_offsets[cell_batch_index] +
7355 const auto cell_type =
7359 for (
unsigned int q = 0; q < this->n_quadrature_points; ++q)
7361 const unsigned int q_src =
7367 if (this->mapping_data->JxW_values.size() > 0)
7368 this_J_value_data[q][v] =
7369 this->mapping_data->JxW_values[offsets + q_src][lane];
7371 if (this->mapping_data->jacobians[0].size() > 0)
7372 for (
unsigned int i = 0; i < dim; ++i)
7373 for (
unsigned int j = 0; j < dim; ++j)
7374 this_jacobian_data[q][i][j][v] =
7376 ->jacobians[0][offsets + q_src][i][j][lane];
7378 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7379 for (
unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7380 for (
unsigned int j = 0; j < dim; ++j)
7381 this_jacobian_gradients_data[q][i][j][v] =
7383 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7385 if (this->mapping_data->quadrature_points.size() > 0)
7394 this->mapping_data->quadrature_points
7396 ->quadrature_point_offsets[cell_batch_index] +
7400 this->mapping_data->jacobians[0][offsets + 1];
7402 for (
unsigned int d = 0;
d < dim; ++
d)
7405 static_cast<Number
>(
7406 this->descriptor->quadrature.point(q)[
d]);
7408 for (
unsigned int d = 0;
d < dim; ++
d)
7409 for (
unsigned int e = 0;
e < dim; ++
e)
7412 static_cast<Number
>(
7413 this->descriptor->quadrature.point(q)[
e]);
7415 for (
unsigned int i = 0; i < dim; ++i)
7416 this_quadrature_points_data[q][i][v] =
point[i][lane];
7421 for (
unsigned int i = 0; i < dim; ++i)
7422 this_quadrature_points_data[q][i][v] =
7423 this->mapping_data->quadrature_points
7425 ->quadrature_point_offsets[cell_batch_index] +
7434 this->is_reinitialized =
true;
7435 this->dof_values_initialized =
false;
7436 this->values_quad_initialized =
false;
7437 this->gradients_quad_initialized =
false;
7438 this->hessians_quad_initialized =
false;
7449 typename VectorizedArrayType>
7450template <
bool level_dof_access>
7457 VectorizedArrayType>
::
7460 Assert(this->matrix_free ==
nullptr,
7461 ExcMessage(
"Cannot use initialization from cell iterator if "
7462 "initialized from MatrixFree object. Use variant for "
7463 "on the fly computation with arguments as for FEValues "
7466 this->mapped_geometry->reinit(
7468 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7469 if (level_dof_access)
7470 cell->get_mg_dof_indices(this->local_dof_indices);
7472 cell->get_dof_indices(this->local_dof_indices);
7475 this->is_reinitialized =
true;
7486 typename VectorizedArrayType>
7493 VectorizedArrayType>
::
7496 Assert(this->matrix_free == 0,
7497 ExcMessage(
"Cannot use initialization from cell iterator if "
7498 "initialized from MatrixFree object. Use variant for "
7499 "on the fly computation with arguments as for FEValues "
7502 this->mapped_geometry->reinit(cell);
7505 this->is_reinitialized =
true;
7516 typename VectorizedArrayType>
7523 VectorizedArrayType>::evaluate(
const bool evaluate_values,
7524 const bool evaluate_gradients,
7525 const bool evaluate_hessians)
7528 Assert(this->dof_values_initialized ==
true,
7531 evaluate(this->values_dofs,
7543 typename VectorizedArrayType>
7550 VectorizedArrayType>::
7554 Assert(this->dof_values_initialized ==
true,
7557 evaluate(this->values_dofs, evaluation_flags);
7567 typename VectorizedArrayType>
7574 VectorizedArrayType>::evaluate(
const VectorizedArrayType
7576 const bool evaluate_values,
7577 const bool evaluate_gradients,
7578 const bool evaluate_hessians)
7587 evaluate(values_array, flag);
7597 typename VectorizedArrayType>
7604 VectorizedArrayType>::
7605 evaluate(
const VectorizedArrayType * values_array,
7608 const bool hessians_on_general_cells =
7612 if (hessians_on_general_cells)
7618 evaluate(n_components, evaluation_flag_actual, values_array, *
this);
7624 evaluation_flag_actual,
7625 const_cast<VectorizedArrayType *
>(values_array),
7631 this->values_quad_initialized =
true;
7633 this->gradients_quad_initialized =
true;
7635 this->hessians_quad_initialized =
true;
7646 typename VectorizedArrayType>
7647template <
typename VectorType>
7655 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
7656 const bool evaluate_values,
7657 const bool evaluate_gradients,
7658 const bool evaluate_hessians)
7667 gather_evaluate(input_vector, flag);
7676 template <
typename Number,
7677 typename VectorizedArrayType,
7678 typename VectorType,
7679 typename EvaluatorType,
7680 typename std::enable_if<internal::has_begin<VectorType> &&
7682 VectorType>::type * =
nullptr>
7683 VectorizedArrayType *
7684 check_vector_access_inplace(
const EvaluatorType &fe_eval, VectorType &vector)
7690 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7691 const auto & dof_info = fe_eval.get_dof_info();
7698 if (std::is_same<typename VectorType::value_type, Number>::value &&
7702 interleaved_contiguous &&
7703 reinterpret_cast<std::size_t
>(
7707 [cell * VectorizedArrayType::size()]) %
7708 sizeof(VectorizedArrayType) ==
7711 return reinterpret_cast<VectorizedArrayType *
>(
7715 [cell * VectorizedArrayType::size()] +
7717 [fe_eval.get_active_fe_index()]
7718 [fe_eval.get_first_selected_component()] *
7719 VectorizedArrayType::size());
7728 template <
typename Number,
7729 typename VectorizedArrayType,
7730 typename VectorType,
7731 typename EvaluatorType,
7732 typename std::enable_if<!internal::has_begin<VectorType> ||
7734 VectorType>::type * =
nullptr>
7735 VectorizedArrayType *
7736 check_vector_access_inplace(
const EvaluatorType &, VectorType &)
7749 typename VectorizedArrayType>
7750template <
typename VectorType>
7757 VectorizedArrayType>::
7758 gather_evaluate(
const VectorType & input_vector,
7761 const VectorizedArrayType *src_ptr =
7762 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7763 *
this, input_vector);
7764 if (src_ptr !=
nullptr)
7765 evaluate(src_ptr, evaluation_flag);
7768 this->read_dof_values(input_vector);
7769 evaluate(this->begin_dof_values(), evaluation_flag);
7780 typename VectorizedArrayType>
7787 VectorizedArrayType>::integrate(
const bool integrate_values,
7788 const bool integrate_gradients)
7790 integrate(integrate_values, integrate_gradients, this->values_dofs);
7793 this->dof_values_initialized =
true;
7804 typename VectorizedArrayType>
7811 VectorizedArrayType>::
7814 integrate(integration_flag, this->values_dofs);
7817 this->dof_values_initialized =
true;
7828 typename VectorizedArrayType>
7835 VectorizedArrayType>::integrate(
const bool integrate_values,
7836 const bool integrate_gradients,
7837 VectorizedArrayType *values_array)
7843 integrate(flag, values_array);
7853 typename VectorizedArrayType>
7860 VectorizedArrayType>::
7862 VectorizedArrayType * values_array,
7863 const bool sum_into_values_array)
7867 Assert(this->values_quad_submitted ==
true,
7870 Assert(this->gradients_quad_submitted ==
true,
7873 Assert(this->hessians_quad_submitted ==
true,
7876 Assert(this->matrix_free !=
nullptr ||
7877 this->mapped_geometry->is_initialized(),
7883 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7884 "EvaluationFlags::hessians are supported."));
7890 unsigned int size = n_components * dim * n_q_points;
7893 for (
unsigned int i = 0; i < size; ++i)
7894 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7898 for (
unsigned int i = 0; i < size; ++i)
7899 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7908 integration_flag_actual,
7911 sum_into_values_array);
7917 integration_flag_actual,
7920 sum_into_values_array);
7924 this->dof_values_initialized =
true;
7935 typename VectorizedArrayType>
7936template <
typename VectorType>
7944 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
7945 const bool integrate_gradients,
7946 VectorType &destination)
7953 integrate_scatter(flag, destination);
7963 typename VectorizedArrayType>
7964template <
typename VectorType>
7971 VectorizedArrayType>::
7973 VectorType & destination)
7975 VectorizedArrayType *dst_ptr =
7976 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7977 *
this, destination);
7978 if (dst_ptr !=
nullptr)
7979 integrate(integration_flag, dst_ptr,
true);
7982 integrate(integration_flag, this->begin_dof_values());
7983 this->distribute_local_to_global(destination);
7994 typename VectorizedArrayType>
8001 VectorizedArrayType>::dof_indices()
const
8003 return {0
U, dofs_per_cell};
8017 typename VectorizedArrayType>
8023 VectorizedArrayType>::
8026 const bool is_interior_face,
8027 const unsigned int dof_no,
8028 const unsigned int quad_no,
8029 const unsigned int first_selected_component,
8030 const unsigned int active_fe_index,
8031 const unsigned int active_quad_index,
8032 const unsigned int face_type)
8033 : BaseClass(matrix_free,
8035 first_selected_component,
8043 , dofs_per_component(this->data->dofs_per_component_on_cell)
8044 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
8045 , n_q_points(this->n_quadrature_points)
8055 typename VectorizedArrayType>
8061 VectorizedArrayType>::
8064 const std::pair<unsigned int, unsigned int> & range,
8065 const bool is_interior_face,
8066 const unsigned int dof_no,
8067 const unsigned int quad_no,
8068 const unsigned int first_selected_component)
8073 first_selected_component,
8074 matrix_free.get_face_active_fe_index(range,
8077 matrix_free.get_face_info(range.
first).face_type)
8087 typename VectorizedArrayType>
8094 VectorizedArrayType>
::reinit(
const unsigned int face_index)
8096 Assert(this->mapped_geometry ==
nullptr,
8097 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8098 " Integer indexing is not possible"));
8099 if (this->mapped_geometry !=
nullptr)
8102 this->cell = face_index;
8103 this->dof_access_index =
8104 this->is_interior_face() ?
8112 this->matrix_free->get_task_info().boundary_partition_data.back())
8113 Assert(this->is_interior_face(),
8115 "Boundary faces do not have a neighbor. When looping over "
8116 "boundary faces use FEFaceEvaluation with the parameter "
8117 "is_interior_face set to true. "));
8119 this->reinit_face(this->matrix_free->
get_face_info(face_index));
8124 this->face_ids[i] = face_index * VectorizedArrayType::size() + i;
8125 for (; i < VectorizedArrayType::size(); ++i)
8128 this->cell_type = this->matrix_free->
get_mapping_info().face_type[face_index];
8129 const unsigned int offsets =
8130 this->mapping_data->data_index_offsets[face_index];
8131 this->J_value = &this->mapping_data->JxW_values[offsets];
8132 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8134 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
8135 this->normal_x_jacobian =
8137 ->normals_times_jacobians[!this->is_interior_face()][offsets];
8138 this->jacobian_gradients =
8139 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8142 if (this->mapping_data->quadrature_point_offsets.empty() ==
false)
8145 this->mapping_data->quadrature_point_offsets.size());
8147 this->mapping_data->quadrature_points.data() +
8148 this->mapping_data->quadrature_point_offsets[this->cell];
8152 this->is_reinitialized =
true;
8153 this->dof_values_initialized =
false;
8154 this->values_quad_initialized =
false;
8155 this->gradients_quad_initialized =
false;
8156 this->hessians_quad_initialized =
false;
8167 typename VectorizedArrayType>
8175 const unsigned int face_number)
8181 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8185 Assert(this->mapped_geometry ==
nullptr,
8186 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
8187 " Integer indexing is not possible"));
8188 if (this->mapped_geometry !=
nullptr)
8193 .faces_by_cells_type[
cell_index][face_number];
8196 this->dof_access_index =
8199 constexpr unsigned int n_lanes = VectorizedArrayType::size();
8201 if (this->is_interior_face() ==
false)
8206 for (
unsigned int i = 0; i < n_lanes; ++i)
8209 const unsigned int cell_this =
cell_index * n_lanes + i;
8211 unsigned int face_index =
8216 this->face_ids[i] = face_index;
8221 this->face_numbers[i] =
static_cast<std::uint8_t
>(-1);
8222 this->face_orientations[i] =
static_cast<std::uint8_t
>(-1);
8229 auto cell_m = faces.cells_interior[face_index % n_lanes];
8230 auto cell_p = faces.cells_exterior[face_index % n_lanes];
8232 const bool face_identifies_as_interior = cell_m != cell_this;
8234 Assert(cell_m == cell_this || cell_p == cell_this,
8238 if (face_identifies_as_interior)
8240 this->cell_ids[i] = cell_m;
8241 this->face_numbers[i] = faces.interior_face_no;
8245 this->cell_ids[i] = cell_p;
8246 this->face_numbers[i] = faces.exterior_face_no;
8249 const bool orientation_interior_face = faces.face_orientation >= 8;
8250 unsigned int face_orientation = faces.face_orientation % 8;
8251 if (face_identifies_as_interior != orientation_interior_face)
8253 constexpr std::array<std::uint8_t, 8> table{
8254 {0, 1, 2, 3, 6, 5, 4, 7}};
8255 face_orientation = table[face_orientation];
8257 this->face_orientations[i] = face_orientation;
8262 this->face_orientations[0] = 0;
8263 this->face_numbers[0] = face_number;
8264 for (
unsigned int i = 0; i < n_lanes; ++i)
8265 this->cell_ids[i] =
cell_index * n_lanes + i;
8266 for (
unsigned int i = 0; i < n_lanes; ++i)
8273 const unsigned int offsets =
8275 .face_data_by_cells[this->quad_no]
8280 .face_data_by_cells[this->quad_no]
8281 .JxW_values.size());
8283 .face_data_by_cells[this->quad_no]
8284 .JxW_values[offsets];
8286 .face_data_by_cells[this->quad_no]
8287 .normal_vectors[offsets];
8289 .face_data_by_cells[this->quad_no]
8290 .jacobians[!this->is_interior_face()][offsets];
8291 this->normal_x_jacobian =
8293 .face_data_by_cells[this->quad_no]
8294 .normals_times_jacobians[!this->is_interior_face()][offsets];
8295 this->jacobian_gradients =
8296 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8300 .face_data_by_cells[this->quad_no]
8301 .quadrature_point_offsets.empty() ==
false)
8303 const unsigned int index =
8307 .face_data_by_cells[this->quad_no]
8308 .quadrature_point_offsets.size());
8310 .face_data_by_cells[this->quad_no]
8311 .quadrature_points.data() +
8313 .face_data_by_cells[this->quad_no]
8314 .quadrature_point_offsets[
index];
8318 this->is_reinitialized =
true;
8319 this->dof_values_initialized =
false;
8320 this->values_quad_initialized =
false;
8321 this->gradients_quad_initialized =
false;
8322 this->hessians_quad_initialized =
false;
8333 typename VectorizedArrayType>
8340 VectorizedArrayType>::evaluate(
const bool evaluate_values,
8341 const bool evaluate_gradients)
8347 evaluate(this->values_dofs, evaluate_values, evaluate_gradients);
8357 typename VectorizedArrayType>
8364 VectorizedArrayType>::
8371 evaluate(this->values_dofs, evaluation_flag);
8381 typename VectorizedArrayType>
8388 VectorizedArrayType>::evaluate(
const VectorizedArrayType
8390 const bool evaluate_values,
8391 const bool evaluate_gradients)
8398 evaluate(values_array, flag);
8408 typename VectorizedArrayType>
8415 VectorizedArrayType>::
8416 evaluate(
const VectorizedArrayType * values_array,
8419 Assert((evaluation_flag &
8422 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8423 "and EvaluationFlags::hessians are supported."));
8425 const bool hessians_on_general_cells =
8429 if (hessians_on_general_cells)
8434 template run<fe_degree, n_q_points_1d>(n_components,
8435 evaluation_flag_actual,
8440 n_components, evaluation_flag_actual, values_array, *
this);
8444 this->values_quad_initialized =
true;
8446 this->gradients_quad_initialized =
true;
8448 this->hessians_quad_initialized =
true;
8459 typename VectorizedArrayType>
8466 VectorizedArrayType>::
8469 integrate(integration_flag, this->values_dofs);
8472 this->dof_values_initialized =
true;
8483 typename VectorizedArrayType>
8490 VectorizedArrayType>::integrate(
const bool integrate_values,
8491 const bool integrate_gradients)
8493 integrate(integrate_values, integrate_gradients, this->values_dofs);
8496 this->dof_values_initialized =
true;
8507 typename VectorizedArrayType>
8514 VectorizedArrayType>::integrate(
const bool integrate_values,
8515 const bool integrate_gradients,
8524 integrate(flag, values_array);
8534 typename VectorizedArrayType>
8541 VectorizedArrayType>::
8543 VectorizedArrayType * values_array)
8545 Assert((integration_flag &
8548 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8549 "and EvaluationFlags::hessians are supported."));
8555 unsigned int size = n_components * dim * n_q_points;
8558 for (
unsigned int i = 0; i < size; ++i)
8559 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
8563 for (
unsigned int i = 0; i < size; ++i)
8564 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
8571 template run<fe_degree, n_q_points_1d>(n_components,
8572 integration_flag_actual,
8577 n_components, integration_flag_actual, values_array, *
this);
8587 typename VectorizedArrayType>
8588template <
typename VectorType>
8596 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
8597 const bool evaluate_values,
8598 const bool evaluate_gradients)
8605 gather_evaluate(input_vector, flag);
8615 typename VectorizedArrayType>
8616template <
typename VectorType>
8623 VectorizedArrayType>::
8624 gather_evaluate(
const VectorType & input_vector,
8627 Assert((evaluation_flag &
8630 ExcMessage(
"Only EvaluationFlags::values, EvaluationFlags::gradients, "
8631 "and EvaluationFlags::hessians are supported."));
8633 const auto shared_vector_data = internal::get_shared_vector_data(
8635 this->dof_access_index ==
8637 this->active_fe_index,
8640 if (this->data->data.front().fe_degree > 0 &&
8641 fast_evaluation_supported(this->data->data.front().fe_degree,
8642 this->data->data.front().n_q_points_1d) &&
8645 typename VectorType::value_type,
8646 VectorizedArrayType>::
8647 supports(evaluation_flag,
8649 internal::get_beginning<typename VectorType::value_type>(
8658 typename VectorType::value_type,
8659 VectorizedArrayType>::template
run<fe_degree,
8663 internal::get_beginning<typename VectorType::value_type>(
8672 typename VectorType::value_type,
8673 VectorizedArrayType>::evaluate(n_components,
8675 internal::get_beginning<
8676 typename VectorType::value_type>(
8684 this->read_dof_values(input_vector);
8685 this->evaluate(evaluation_flag);
8690 this->values_quad_initialized =
true;
8692 this->gradients_quad_initialized =
true;
8694 this->hessians_quad_initialized =
true;
8705 typename VectorizedArrayType>
8706template <
typename VectorType>
8714 VectorizedArrayType>::integrate_scatter(
const bool integrate_values,
8715 const bool integrate_gradients,
8716 VectorType &destination)
8723 integrate_scatter(flag, destination);
8733 typename VectorizedArrayType>
8734template <
typename VectorType>
8741 VectorizedArrayType>::
8743 VectorType & destination)
8745 Assert((this->dof_access_index ==
8747 this->is_interior_face() ==
false) ==
false,
8750 const auto shared_vector_data = internal::get_shared_vector_data(
8752 this->dof_access_index ==
8754 this->active_fe_index,
8757 if (this->data->data.front().fe_degree > 0 &&
8758 fast_evaluation_supported(this->data->data.front().fe_degree,
8759 this->data->data.front().n_q_points_1d) &&
8762 typename VectorType::value_type,
8763 VectorizedArrayType>::
8764 supports(integration_flag,
8766 internal::get_beginning<typename VectorType::value_type>(
8775 typename VectorType::value_type,
8776 VectorizedArrayType>::template
run<fe_degree,
8780 internal::get_beginning<typename VectorType::value_type>(
8789 typename VectorType::value_type,
8790 VectorizedArrayType>::integrate(n_components,
8792 internal::get_beginning<
8793 typename VectorType::value_type>(
8801 integrate(integration_flag);
8802 this->distribute_local_to_global(destination);
8813 typename VectorizedArrayType>
8820 VectorizedArrayType>::dof_indices()
const
8822 return {0
U, dofs_per_cell};
8832 typename VectorizedArrayType>
8839 VectorizedArrayType>::
8840 fast_evaluation_supported(
const unsigned int given_degree,
8841 const unsigned int give_n_q_points_1d)
8843 return fe_degree == -1 ?
8856 typename VectorizedArrayType>
8863 VectorizedArrayType>::
8864 fast_evaluation_supported(
const unsigned int given_degree,
8865 const unsigned int give_n_q_points_1d)
8867 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)
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)
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)
std::string to_string(const T &t)
#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 > &)