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)