17 #ifndef dealii_matrix_free_fe_evaluation_h
18 #define dealii_matrix_free_fe_evaluation_h
45 #include <type_traits>
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,
834 template <
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,
992 template <
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,
1168 template <
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);
2193 const bool evaluate_gradients,
2194 const bool evaluate_hessians =
false);
2218 const bool evaluate_values,
2219 const bool evaluate_gradients,
2220 const bool evaluate_hessians =
false);
2235 template <
typename VectorType>
2243 template <
typename VectorType>
2246 const bool evaluate_values,
2247 const bool evaluate_gradients,
2248 const bool evaluate_hessians =
false);
2267 integrate(
const bool integrate_values,
const bool integrate_gradients);
2282 VectorizedArrayType * values_array,
2283 const bool sum_into_values =
false);
2290 const bool integrate_gradients,
2291 VectorizedArrayType *values_array);
2306 template <
typename VectorType>
2309 VectorType & output_vector);
2314 template <
typename VectorType>
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);
2612 evaluate(
const bool evaluate_values,
const bool evaluate_gradients);
2635 const bool evaluate_values,
2636 const bool evaluate_gradients);
2649 template <
typename VectorType>
2657 template <
typename VectorType>
2660 const bool evaluate_values,
2661 const bool evaluate_gradients);
2679 integrate(
const bool integrate_values,
const bool integrate_gradients);
2691 VectorizedArrayType * values_array);
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);
2786 namespace MatrixFreeFunctions
2790 template <
int dim,
int degree>
2799 template <
int degree>
2802 static constexpr
unsigned int value = degree + 1;
2818 template <
bool is_face,
2821 typename VectorizedArrayType>
2824 extract_initialization_data(
2826 const unsigned int dof_no,
2827 const unsigned int first_selected_component,
2828 const unsigned int quad_no,
2829 const unsigned int fe_degree,
2830 const unsigned int n_q_points,
2831 const unsigned int active_fe_index_given,
2832 const unsigned int active_quad_index_given,
2833 const unsigned int face_type)
2836 InitializationData init_data;
2840 &internal::MatrixFreeFunctions::
2841 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2849 active_fe_index_given :
2854 active_quad_index_given :
2885 typename VectorizedArrayType>
2890 VectorizedArrayType>::
2893 const unsigned int dof_no,
2894 const unsigned int first_selected_component,
2895 const unsigned int quad_no,
2896 const unsigned int fe_degree,
2897 const unsigned int n_q_points,
2898 const bool is_interior_face,
2899 const unsigned int active_fe_index,
2900 const unsigned int active_quad_index,
2901 const unsigned int face_type)
2903 internal::extract_initialization_data<is_face>(matrix_free,
2905 first_selected_component,
2914 first_selected_component)
2915 , scratch_data_array(matrix_free.acquire_scratch_data())
2916 , matrix_free(&matrix_free)
2918 this->set_data_pointers(scratch_data_array, n_components_);
2920 this->dof_info->start_components.back() == 1 ||
2921 static_cast<int>(n_components_) <=
2923 this->dof_info->start_components
2924 [this->dof_info->component_to_base_index[first_selected_component] +
2926 first_selected_component,
2928 "You tried to construct a vector-valued evaluator with " +
2930 " components. However, "
2931 "the current base element has only " +
2933 this->dof_info->start_components
2934 [this->dof_info->component_to_base_index[first_selected_component] +
2936 first_selected_component) +
2937 " components left when starting from local element index " +
2939 first_selected_component -
2940 this->dof_info->start_components
2941 [this->dof_info->component_to_base_index[first_selected_component]]) +
2942 " (global index " +
std::to_string(first_selected_component) +
")"));
2954 typename VectorizedArrayType>
2959 VectorizedArrayType>::
2965 const unsigned int first_selected_component,
2969 other->mapped_geometry->get_quadrature() == quadrature ?
2970 other->mapped_geometry :
2971 std::make_shared<
internal::MatrixFreeFunctions::
2972 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2977 first_selected_component)
2978 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2979 , matrix_free(nullptr)
2981 const unsigned int base_element_number =
2985 first_selected_component >=
2987 ExcMessage(
"The underlying element must at least contain as many "
2988 "components as requested by this class"));
2989 (void)base_element_number;
2994 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2998 this->set_data_pointers(scratch_data_array, n_components_);
3007 typename VectorizedArrayType>
3012 VectorizedArrayType>::
3017 VectorizedArrayType> &other)
3019 , scratch_data_array(other.matrix_free == nullptr ?
3021 other.matrix_free->acquire_scratch_data())
3022 , matrix_free(other.matrix_free)
3024 if (other.matrix_free ==
nullptr)
3032 this->mapped_geometry =
3033 std::make_shared<internal::MatrixFreeFunctions::
3034 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3038 this->mapping_data = &this->mapped_geometry->get_data_storage();
3042 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3044 this->mapped_geometry->get_data_storage().JxW_values.begin();
3045 this->jacobian_gradients =
3046 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3047 this->jacobian_gradients_non_inverse =
3048 this->mapped_geometry->get_data_storage()
3049 .jacobian_gradients_non_inverse[0]
3052 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3055 this->set_data_pointers(scratch_data_array, n_components_);
3064 typename VectorizedArrayType>
3069 VectorizedArrayType> &
3075 VectorizedArrayType> &other)
3078 if (matrix_free ==
nullptr)
3081 delete scratch_data_array;
3090 matrix_free = other.matrix_free;
3092 if (other.matrix_free ==
nullptr)
3101 this->mapped_geometry =
3102 std::make_shared<internal::MatrixFreeFunctions::
3103 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3108 this->mapping_data = &this->mapped_geometry->get_data_storage();
3110 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3112 this->mapped_geometry->get_data_storage().JxW_values.begin();
3113 this->jacobian_gradients =
3114 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3115 this->jacobian_gradients_non_inverse =
3116 this->mapped_geometry->get_data_storage()
3117 .jacobian_gradients_non_inverse[0]
3120 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3127 this->set_data_pointers(scratch_data_array, n_components_);
3138 typename VectorizedArrayType>
3143 VectorizedArrayType>::~FEEvaluationBase()
3145 if (matrix_free !=
nullptr)
3156 delete scratch_data_array;
3167 typename VectorizedArrayType>
3172 Assert(matrix_free !=
nullptr,
3174 "FEEvaluation was not initialized with a MatrixFree object!"));
3175 return *matrix_free;
3184 template <
typename VectorType,
bool>
3185 struct ConstBlockVectorSelector;
3187 template <
typename VectorType>
3188 struct ConstBlockVectorSelector<VectorType, true>
3190 using BaseVectorType =
const typename VectorType::BlockType;
3193 template <
typename VectorType>
3194 struct ConstBlockVectorSelector<VectorType, false>
3196 using BaseVectorType =
typename VectorType::BlockType;
3202 template <
typename VectorType,
bool>
3203 struct BlockVectorSelector;
3205 template <
typename VectorType>
3206 struct BlockVectorSelector<VectorType, true>
3208 using BaseVectorType =
typename ConstBlockVectorSelector<
3210 std::is_const<VectorType>::value>::BaseVectorType;
3212 static BaseVectorType *
3213 get_vector_component(VectorType &vec,
const unsigned int component)
3216 return &vec.block(component);
3220 template <
typename VectorType>
3221 struct BlockVectorSelector<VectorType, false>
3223 using BaseVectorType = VectorType;
3225 static BaseVectorType *
3226 get_vector_component(VectorType &vec,
const unsigned int component)
3245 template <
typename VectorType>
3246 struct BlockVectorSelector<std::vector<VectorType>, false>
3248 using BaseVectorType = VectorType;
3250 static BaseVectorType *
3251 get_vector_component(std::vector<VectorType> &vec,
3252 const unsigned int component)
3255 return &vec[component];
3259 template <
typename VectorType>
3260 struct BlockVectorSelector<const std::vector<VectorType>, false>
3262 using BaseVectorType =
const VectorType;
3264 static const BaseVectorType *
3265 get_vector_component(
const std::vector<VectorType> &vec,
3266 const unsigned int component)
3269 return &vec[component];
3273 template <
typename VectorType>
3274 struct BlockVectorSelector<std::vector<VectorType *>, false>
3276 using BaseVectorType = VectorType;
3278 static BaseVectorType *
3279 get_vector_component(std::vector<VectorType *> &vec,
3280 const unsigned int component)
3283 return vec[component];
3287 template <
typename VectorType>
3288 struct BlockVectorSelector<const std::vector<VectorType *>, false>
3290 using BaseVectorType =
const VectorType;
3292 static const BaseVectorType *
3293 get_vector_component(
const std::vector<VectorType *> &vec,
3294 const unsigned int component)
3297 return vec[component];
3308 typename VectorizedArrayType>
3309 template <
typename VectorType,
typename VectorOperation>
3314 const std::array<VectorType *, n_components_> &src,
3317 n_components_> & src_sm,
3318 const std::bitset<VectorizedArrayType::size()> &mask,
3319 const bool apply_constraints)
const
3324 if (this->matrix_free ==
nullptr)
3326 read_write_operation_global(operation, src);
3333 if (this->n_fe_components == 1)
3334 for (
unsigned int comp = 0; comp < n_components; ++comp)
3336 Assert(src[comp] !=
nullptr,
3337 ExcMessage(
"The finite element underlying this FEEvaluation "
3338 "object is scalar, but you requested " +
3340 " components via the template argument in "
3341 "FEEvaluation. In that case, you must pass an "
3342 "std::vector<VectorType> or a BlockVector to " +
3343 "read_dof_values and distribute_local_to_global."));
3363 [is_face ? this->dof_access_index :
3366 IndexStorageVariants::contiguous)
3368 read_write_operation_contiguous(operation, src, src_sm, mask);
3375 constexpr
unsigned int n_lanes = VectorizedArrayType::size();
3377 std::array<
unsigned int, VectorizedArrayType::size()> cells =
3378 this->get_cell_ids();
3380 const bool masking_is_active =
mask.count() < n_lanes;
3381 if (masking_is_active)
3382 for (
unsigned int v = 0; v < n_lanes; ++v)
3383 if (mask[v] ==
false)
3386 bool has_hn_constraints =
false;
3388 if (is_face ==
false)
3394 [this->first_selected_component])
3395 for (
unsigned int v = 0; v < n_lanes; ++v)
3400 has_hn_constraints =
true;
3403 std::integral_constant<
bool,
3404 internal::is_vectorizable<VectorType, Number>::value>
3407 const bool is_neighbor_cells = !is_face && !this->is_interior_face();
3408 const bool use_vectorized_path =
3409 !(masking_is_active || is_neighbor_cells || has_hn_constraints);
3411 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3412 std::array<VectorizedArrayType *, n_components> values_dofs;
3413 for (
unsigned int c = 0; c < n_components; ++c)
3414 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3415 c * dofs_per_component;
3419 [is_face ? this->dof_access_index :
3422 IndexStorageVariants::interleaved &&
3423 use_vectorized_path)
3425 const unsigned int *dof_indices =
3427 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
3431 [this->first_selected_component] *
3433 if (n_components == 1 || this->n_fe_components == 1)
3434 for (
unsigned int i = 0; i < dofs_per_component;
3435 ++i, dof_indices += n_lanes)
3436 for (
unsigned int comp = 0; comp < n_components; ++comp)
3437 operation.process_dof_gather(dof_indices,
3440 values_dofs[comp][i],
3443 for (
unsigned int comp = 0; comp < n_components; ++comp)
3444 for (
unsigned int i = 0; i < dofs_per_component;
3445 ++i, dof_indices += n_lanes)
3446 operation.process_dof_gather(
3447 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3453 std::array<const unsigned int *, n_lanes> dof_indices;
3454 dof_indices.fill(
nullptr);
3459 bool has_constraints =
false;
3460 const unsigned int n_components_read =
3461 this->n_fe_components > 1 ? n_components : 1;
3465 for (
unsigned int v = 0; v < n_lanes; ++v)
3471 const std::pair<unsigned int, unsigned int> *my_index_start =
3472 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3473 this->first_selected_component];
3478 if (my_index_start[n_components_read].
second !=
3479 my_index_start[0].
second)
3480 has_constraints =
true;
3483 dof_info.
dof_indices.data() + my_index_start[0].first;
3488 for (
unsigned int v = 0; v < n_lanes; ++v)
3493 const std::pair<unsigned int, unsigned int> *my_index_start =
3494 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3495 this->first_selected_component];
3496 if (my_index_start[n_components_read].
second !=
3497 my_index_start[0].
second)
3498 has_constraints =
true;
3506 [this->active_fe_index][this->first_selected_component])
3507 has_hn_constraints =
true;
3510 my_index_start[0].
first ||
3513 my_index_start[0].
first,
3516 dof_info.
dof_indices.data() + my_index_start[0].first;
3520 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3521 return i != numbers::invalid_unsigned_int;
3523 for (
unsigned int comp = 0; comp < n_components; ++comp)
3524 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3525 operation.process_empty(values_dofs[comp][i]);
3529 if (!has_constraints && apply_constraints)
3531 if (n_components == 1 || this->n_fe_components == 1)
3533 for (
unsigned int v = 0; v < n_lanes; ++v)
3538 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3539 for (
unsigned int comp = 0; comp < n_components; ++comp)
3540 operation.process_dof(dof_indices[v][i],
3542 values_dofs[comp][i][v]);
3547 for (
unsigned int comp = 0; comp < n_components; ++comp)
3548 for (
unsigned int v = 0; v < n_lanes; ++v)
3553 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3554 operation.process_dof(
3555 dof_indices[v][comp * dofs_per_component + i],
3557 values_dofs[comp][i][v]);
3569 for (
unsigned int v = 0; v < n_lanes; ++v)
3575 const unsigned int cell_dof_index =
3576 cell_index * this->n_fe_components + this->first_selected_component;
3577 const unsigned int n_components_read =
3578 this->n_fe_components > 1 ? n_components : 1;
3579 unsigned int index_indicators =
3581 unsigned int next_index_indicators =
3582 dof_info.
row_starts[cell_dof_index + 1].second;
3586 if (apply_constraints ==
false &&
3587 (dof_info.
row_starts[cell_dof_index].second !=
3588 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3595 [this->active_fe_index][this->first_selected_component])))
3604 [this->first_selected_component] +
3606 next_index_indicators = index_indicators;
3609 if (n_components == 1 || this->n_fe_components == 1)
3611 unsigned int ind_local = 0;
3612 for (; index_indicators != next_index_indicators; ++index_indicators)
3614 const std::pair<unsigned short, unsigned short> indicator =
3617 for (
unsigned int j = 0; j < indicator.first; ++j)
3618 for (
unsigned int comp = 0; comp < n_components; ++comp)
3619 operation.process_dof(dof_indices[v][j],
3621 values_dofs[comp][ind_local + j][v]);
3623 ind_local += indicator.first;
3624 dof_indices[v] += indicator.first;
3628 Number
value[n_components];
3629 for (
unsigned int comp = 0; comp < n_components; ++comp)
3630 operation.pre_constraints(values_dofs[comp][ind_local][v],
3633 const Number *data_val =
3635 const Number *end_pool =
3637 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3638 for (
unsigned int comp = 0; comp < n_components; ++comp)
3639 operation.process_constraint(*dof_indices[v],
3644 for (
unsigned int comp = 0; comp < n_components; ++comp)
3645 operation.post_constraints(value[comp],
3646 values_dofs[comp][ind_local][v]);
3652 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3653 for (
unsigned int comp = 0; comp < n_components; ++comp)
3654 operation.process_dof(*dof_indices[v],
3656 values_dofs[comp][ind_local][v]);
3665 for (
unsigned int comp = 0; comp < n_components; ++comp)
3667 unsigned int ind_local = 0;
3670 for (; index_indicators != next_index_indicators;
3673 const std::pair<unsigned short, unsigned short> indicator =
3677 for (
unsigned int j = 0; j < indicator.first; ++j)
3678 operation.process_dof(dof_indices[v][j],
3680 values_dofs[comp][ind_local + j][v]);
3681 ind_local += indicator.first;
3682 dof_indices[v] += indicator.first;
3687 operation.pre_constraints(values_dofs[comp][ind_local][v],
3690 const Number *data_val =
3692 const Number *end_pool =
3695 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3696 operation.process_constraint(*dof_indices[v],
3701 operation.post_constraints(value,
3702 values_dofs[comp][ind_local][v]);
3709 for (; ind_local < dofs_per_component;
3710 ++dof_indices[v], ++ind_local)
3713 operation.process_dof(*dof_indices[v],
3715 values_dofs[comp][ind_local][v]);
3718 if (apply_constraints ==
true && comp + 1 < n_components)
3719 next_index_indicators =
3720 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3732 typename VectorizedArrayType>
3733 template <
typename VectorType,
typename VectorOperation>
3738 const std::array<VectorType *, n_components_> &src)
const
3742 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3743 unsigned int index = this->first_selected_component * dofs_per_component;
3744 for (
unsigned int comp = 0; comp < n_components; ++comp)
3746 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3748 operation.process_empty(
3749 this->values_dofs[comp * dofs_per_component + i]);
3750 operation.process_dof_global(
3751 local_dof_indices[this->data->lexicographic_numbering[index]],
3753 this->values_dofs[comp * dofs_per_component + i][0]);
3764 typename VectorizedArrayType>
3765 template <
typename VectorType,
typename VectorOperation>
3770 const std::array<VectorType *, n_components_> &src,
3773 n_components_> & vectors_sm,
3774 const std::bitset<VectorizedArrayType::size()> &mask)
const
3783 std::integral_constant<
bool,
3784 internal::is_vectorizable<VectorType, Number>::value>
3787 is_face ? this->dof_access_index :
3789 const unsigned int n_lanes =
mask.count();
3792 const std::vector<unsigned int> & dof_indices_cont =
3795 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3796 std::array<VectorizedArrayType *, n_components> values_dofs;
3797 for (
unsigned int c = 0; c < n_components; ++c)
3798 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3799 c * dofs_per_component;
3807 interleaved_contiguous &&
3808 n_lanes == VectorizedArrayType::size()) &&
3810 this->dof_access_index ==
3812 this->is_interior_face() ==
false) &&
3813 !(!is_face && !this->is_interior_face()))
3815 const unsigned int dof_index =
3816 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3819 [this->first_selected_component] *
3820 VectorizedArrayType::size();
3821 if (n_components == 1 || this->n_fe_components == 1)
3822 for (
unsigned int comp = 0; comp < n_components; ++comp)
3823 operation.process_dofs_vectorized(dofs_per_component,
3829 operation.process_dofs_vectorized(dofs_per_component * n_components,
3837 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
3838 this->get_cell_or_face_ids();
3842 const unsigned int n_filled_lanes =
3846 (this->dof_access_index ==
3848 this->is_interior_face() ==
false) ||
3849 (!is_face && !this->is_interior_face());
3851 if (vectors_sm[0] !=
nullptr)
3853 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3854 std::array<
typename VectorType::value_type *,
3855 VectorizedArrayType::size()>
3858 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3860 if (mask[v] ==
false)
3862 vector_ptrs[v] =
nullptr;
3882 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3883 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3885 [this->active_fe_index][this->first_selected_component]);
3887 vector_ptrs[v] =
nullptr;
3889 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3891 vector_ptrs[v] =
nullptr;
3896 if (n_filled_lanes == VectorizedArrayType::size() &&
3897 n_lanes == VectorizedArrayType::size() && !is_ecl)
3899 if (n_components == 1 || this->n_fe_components == 1)
3901 for (
unsigned int comp = 0; comp < n_components; ++comp)
3903 auto vector_ptrs = compute_vector_ptrs(comp);
3904 operation.process_dofs_vectorized_transpose(
3913 auto vector_ptrs = compute_vector_ptrs(0);
3914 operation.process_dofs_vectorized_transpose(dofs_per_component *
3922 for (
unsigned int comp = 0; comp < n_components; ++comp)
3924 auto vector_ptrs = compute_vector_ptrs(
3925 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3927 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3928 operation.process_empty(values_dofs[comp][i]);
3930 if (n_components == 1 || this->n_fe_components == 1)
3932 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3933 if (mask[v] ==
true)
3934 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3935 operation.process_dof(vector_ptrs[v][i],
3936 values_dofs[comp][i][v]);
3940 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3941 if (mask[v] ==
true)
3942 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3943 operation.process_dof(
3944 vector_ptrs[v][i + comp * dofs_per_component],
3945 values_dofs[comp][i][v]);
3951 unsigned int dof_indices[VectorizedArrayType::size()];
3953 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3957 if (mask[v] ==
true)
3959 dof_indices_cont[cells[v]] +
3962 [this->first_selected_component] *
3968 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3973 if (n_filled_lanes == VectorizedArrayType::size() &&
3974 n_lanes == VectorizedArrayType::size() && !is_ecl)
3980 if (n_components == 1 || this->n_fe_components == 1)
3981 for (
unsigned int comp = 0; comp < n_components; ++comp)
3982 operation.process_dofs_vectorized_transpose(dofs_per_component,
3988 operation.process_dofs_vectorized_transpose(dofs_per_component *
3997 interleaved_contiguous_strided)
3999 if (n_components == 1 || this->n_fe_components == 1)
4000 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4002 for (
unsigned int comp = 0; comp < n_components; ++comp)
4003 operation.process_dof_gather(dof_indices,
4005 i * VectorizedArrayType::size(),
4006 values_dofs[comp][i],
4010 for (
unsigned int comp = 0; comp < n_components; ++comp)
4011 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4013 operation.process_dof_gather(dof_indices,
4015 (comp * dofs_per_component + i) *
4016 VectorizedArrayType::size(),
4017 values_dofs[comp][i],
4025 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4027 const unsigned int *offsets =
4029 [ind][VectorizedArrayType::size() * this->cell];
4030 if (n_components == 1 || this->n_fe_components == 1)
4031 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4033 for (
unsigned int comp = 0; comp < n_components; ++comp)
4034 operation.process_dof_gather(dof_indices,
4037 values_dofs[comp][i],
4040 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4041 dof_indices[v] += offsets[v];
4044 for (
unsigned int comp = 0; comp < n_components; ++comp)
4045 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4047 operation.process_dof_gather(dof_indices,
4050 values_dofs[comp][i],
4053 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4054 dof_indices[v] += offsets[v];
4059 for (
unsigned int comp = 0; comp < n_components; ++comp)
4061 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4062 operation.process_empty(values_dofs[comp][i]);
4067 if (n_components == 1 || this->n_fe_components == 1)
4069 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4070 if (mask[v] ==
true)
4071 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4072 operation.process_dof(dof_indices[v] + i,
4074 values_dofs[comp][i][v]);
4078 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4079 if (mask[v] ==
true)
4080 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4081 operation.process_dof(dof_indices[v] + i +
4082 comp * dofs_per_component,
4084 values_dofs[comp][i][v]);
4089 const unsigned int *offsets =
4091 [ind][VectorizedArrayType::size() * this->cell];
4092 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4094 if (n_components == 1 || this->n_fe_components == 1)
4095 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4097 if (mask[v] ==
true)
4098 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4099 operation.process_dof(dof_indices[v] + i * offsets[v],
4101 values_dofs[comp][i][v]);
4105 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4106 if (mask[v] ==
true)
4107 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4108 operation.process_dof(dof_indices[v] +
4109 (i + comp * dofs_per_component) *
4112 values_dofs[comp][i][v]);
4122 typename VectorType,
4123 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4124 decltype(std::declval<VectorType>().
begin())
4125 get_beginning(VectorType &vec)
4132 typename VectorType,
4133 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4134 typename VectorType::value_type *
4135 get_beginning(VectorType &)
4140 template <
typename VectorType,
4141 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4143 const std::vector<ArrayView<const typename VectorType::value_type>> *
4144 get_shared_vector_data(VectorType * vec,
4145 const bool is_valid_mode_for_sm,
4146 const unsigned int active_fe_index,
4150 if (is_valid_mode_for_sm &&
4153 active_fe_index == 0)
4154 return &vec->shared_vector_data();
4159 template <
typename VectorType,
4160 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4162 const std::vector<ArrayView<const typename VectorType::value_type>> *
4163 get_shared_vector_data(VectorType *,
4171 template <
int n_components,
typename VectorType>
4173 std::array<
typename internal::BlockVectorSelector<
4178 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
4182 get_vector_data(VectorType & src,
4183 const unsigned int first_index,
4184 const bool is_valid_mode_for_sm,
4185 const unsigned int active_fe_index,
4191 std::array<
typename internal::BlockVectorSelector<
4197 ArrayView<
const typename internal::BlockVectorSelector<
4203 for (
unsigned int d = 0;
d < n_components; ++
d)
4204 src_data.first[
d] = internal::BlockVectorSelector<
4210 for (
unsigned int d = 0;
d < n_components; ++
d)
4211 src_data.second[
d] = get_shared_vector_data(
4212 const_cast<typename internal::BlockVectorSelector<
4213 typename std::remove_const<VectorType>::type,
4215 BaseVectorType *>(src_data.first[
d]),
4216 is_valid_mode_for_sm,
4230 typename VectorizedArrayType>
4235 if (this->dof_info ==
nullptr ||
4237 this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4238 this->dof_info->hanging_node_constraint_masks_comp
4239 [this->active_fe_index][this->first_selected_component] ==
false)
4242 constexpr
unsigned int n_lanes = VectorizedArrayType::size();
4243 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4246 bool hn_available =
false;
4248 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
4249 this->get_cell_ids();
4251 for (
unsigned int v = 0; v < n_lanes; ++v)
4263 constraint_mask[v] =
mask;
4269 if (hn_available ==
false)
4274 this->data->data.front().fe_degree,
4275 this->get_shape_info(),
4287 typename VectorizedArrayType>
4288 template <
typename VectorType>
4292 const unsigned int first_index,
4293 const std::bitset<VectorizedArrayType::size()> &mask)
4295 const auto src_data = internal::get_vector_data<n_components_>(
4298 this->dof_access_index ==
4300 this->active_fe_index,
4304 read_write_operation(reader, src_data.first, src_data.second, mask,
true);
4306 apply_hanging_node_constraints(
false);
4309 this->dof_values_initialized =
true;
4319 typename VectorizedArrayType>
4320 template <
typename VectorType>
4324 const unsigned int first_index,
4325 const std::bitset<VectorizedArrayType::size()> &mask)
4327 const auto src_data = internal::get_vector_data<n_components_>(
4330 this->dof_access_index ==
4332 this->active_fe_index,
4336 read_write_operation(reader, src_data.first, src_data.second, mask,
false);
4339 this->dof_values_initialized =
true;
4349 typename VectorizedArrayType>
4350 template <
typename VectorType>
4355 const unsigned int first_index,
4356 const std::bitset<VectorizedArrayType::size()> &mask)
const
4359 Assert(this->dof_values_initialized ==
true,
4363 apply_hanging_node_constraints(
true);
4365 const auto dst_data = internal::get_vector_data<n_components_>(
4368 this->dof_access_index ==
4370 this->active_fe_index,
4375 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4384 typename VectorizedArrayType>
4385 template <
typename VectorType>
4389 const unsigned int first_index,
4390 const std::bitset<VectorizedArrayType::size()> &mask)
const
4393 Assert(this->dof_values_initialized ==
true,
4397 const auto dst_data = internal::get_vector_data<n_components_>(
4400 this->dof_access_index ==
4402 this->active_fe_index,
4406 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4415 typename VectorizedArrayType>
4416 template <
typename VectorType>
4421 const unsigned int first_index,
4422 const std::bitset<VectorizedArrayType::size()> &mask)
const
4425 Assert(this->dof_values_initialized ==
true,
4429 const auto dst_data = internal::get_vector_data<n_components_>(
4432 this->dof_access_index ==
4434 this->active_fe_index,
4438 read_write_operation(setter, dst_data.first, dst_data.second, mask,
false);
4451 typename VectorizedArrayType>
4457 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4459 for (
unsigned int comp = 0; comp < n_components; ++comp)
4460 return_value[comp] = this->values_dofs[comp * dofs + dof];
4461 return return_value;
4470 typename VectorizedArrayType>
4473 get_value(
const unsigned int q_point)
const
4476 Assert(this->values_quad_initialized ==
true,
4481 const std::size_t nqp = this->n_quadrature_points;
4483 for (
unsigned int comp = 0; comp < n_components; ++comp)
4484 return_value[comp] = this->values_quad[comp * nqp + q_point];
4485 return return_value;
4494 typename VectorizedArrayType>
4501 Assert(this->gradients_quad_initialized ==
true,
4506 Assert(this->jacobian !=
nullptr,
4508 "update_gradients"));
4509 const std::size_t nqp = this->n_quadrature_points;
4515 for (
unsigned int d = 0;
d < dim; ++
d)
4516 for (
unsigned int comp = 0; comp < n_components; ++comp)
4518 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
4519 this->jacobian[0][
d][
d];
4528 for (
unsigned int comp = 0; comp < n_components; ++comp)
4529 for (
unsigned int d = 0;
d < dim; ++
d)
4532 jac[
d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4533 for (
unsigned int e = 1;
e < dim; ++
e)
4534 grad_out[comp][
d] +=
4536 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4548 typename VectorizedArrayType>
4555 Assert(this->gradients_quad_initialized ==
true,
4559 Assert(this->normal_x_jacobian !=
nullptr,
4561 "update_gradients"));
4563 const std::size_t nqp = this->n_quadrature_points;
4567 for (
unsigned int comp = 0; comp < n_components; ++comp)
4569 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4570 (this->normal_x_jacobian[0][dim - 1]);
4573 const std::size_t
index =
4575 for (
unsigned int comp = 0; comp < n_components; ++comp)
4577 grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4578 this->normal_x_jacobian[
index][0];
4579 for (
unsigned int d = 1;
d < dim; ++
d)
4581 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
4582 this->normal_x_jacobian[
index][
d];
4594 template <
typename VectorizedArrayType>
4597 const VectorizedArrayType *
const hessians,
4599 VectorizedArrayType (&tmp)[1][1])
4601 tmp[0][0] = jac[0][0] *
hessians[0];
4604 template <
typename VectorizedArrayType>
4607 const VectorizedArrayType *
const hessians,
4608 const unsigned int nqp,
4609 VectorizedArrayType (&tmp)[2][2])
4611 for (
unsigned int d = 0;
d < 2; ++
d)
4619 template <
typename VectorizedArrayType>
4622 const VectorizedArrayType *
const hessians,
4623 const unsigned int nqp,
4624 VectorizedArrayType (&tmp)[3][3])
4626 for (
unsigned int d = 0;
d < 3; ++
d)
4647 typename VectorizedArrayType>
4653 Assert(this->hessians_quad_initialized ==
true,
4658 Assert(this->jacobian !=
nullptr,
4668 const std::size_t nqp = this->n_quadrature_points;
4669 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4674 for (
unsigned int comp = 0; comp < n_components; ++comp)
4676 for (
unsigned int d = 0;
d < dim; ++
d)
4677 hessian_out[comp][
d][
d] =
4678 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4679 (jac[
d][
d] * jac[
d][
d]);
4685 hessian_out[comp][0][1] =
4686 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4687 (jac[0][0] * jac[1][1]);
4690 hessian_out[comp][0][1] =
4691 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4692 (jac[0][0] * jac[1][1]);
4693 hessian_out[comp][0][2] =
4694 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4695 (jac[0][0] * jac[2][2]);
4696 hessian_out[comp][1][2] =
4697 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4698 (jac[1][1] * jac[2][2]);
4703 for (
unsigned int d = 0;
d < dim; ++
d)
4704 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4705 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4711 for (
unsigned int comp = 0; comp < n_components; ++comp)
4713 VectorizedArrayType tmp[dim][dim];
4714 internal::hessian_unit_times_jac(
4715 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4718 for (
unsigned int d = 0;
d < dim; ++
d)
4719 for (
unsigned int e =
d;
e < dim; ++
e)
4721 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4722 for (
unsigned int f = 1; f < dim; ++f)
4723 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4730 for (
unsigned int d = 0;
d < dim; ++
d)
4731 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4732 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4738 const auto &jac_grad = this->jacobian_gradients[q_point];
4739 for (
unsigned int comp = 0; comp < n_components; ++comp)
4741 VectorizedArrayType tmp[dim][dim];
4742 internal::hessian_unit_times_jac(
4743 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4746 for (
unsigned int d = 0;
d < dim; ++
d)
4747 for (
unsigned int e =
d;
e < dim; ++
e)
4749 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4750 for (
unsigned int f = 1; f < dim; ++f)
4751 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4755 for (
unsigned int d = 0;
d < dim; ++
d)
4756 for (
unsigned int e = 0;
e < dim; ++
e)
4757 hessian_out[comp][
d][
d] +=
4759 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4762 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4763 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++count)
4764 for (
unsigned int f = 0; f < dim; ++f)
4765 hessian_out[comp][
d][
e] +=
4766 jac_grad[count][f] *
4767 this->gradients_quad[(comp * dim + f) * nqp + q_point];
4770 for (
unsigned int d = 0;
d < dim; ++
d)
4771 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4772 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4784 typename VectorizedArrayType>
4791 Assert(this->hessians_quad_initialized ==
true,
4802 const std::size_t nqp = this->n_quadrature_points;
4803 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4809 for (
unsigned int comp = 0; comp < n_components; ++comp)
4810 for (
unsigned int d = 0;
d < dim; ++
d)
4811 hessian_out[comp][
d] =
4812 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4813 (jac[
d][
d] * jac[
d][
d]);
4818 for (
unsigned int comp = 0; comp < n_components; ++comp)
4822 VectorizedArrayType tmp[dim][dim];
4823 internal::hessian_unit_times_jac(
4824 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4828 for (
unsigned int d = 0;
d < dim; ++
d)
4830 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4831 for (
unsigned int f = 1; f < dim; ++f)
4832 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4839 const auto &jac_grad = this->jacobian_gradients[q_point];
4840 for (
unsigned int comp = 0; comp < n_components; ++comp)
4844 VectorizedArrayType tmp[dim][dim];
4845 internal::hessian_unit_times_jac(
4846 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4850 for (
unsigned int d = 0;
d < dim; ++
d)
4852 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4853 for (
unsigned int f = 1; f < dim; ++f)
4854 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4857 for (
unsigned int d = 0;
d < dim; ++
d)
4858 for (
unsigned int e = 0;
e < dim; ++
e)
4859 hessian_out[comp][
d] +=
4861 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4873 typename VectorizedArrayType>
4880 Assert(this->hessians_quad_initialized ==
true,
4886 const auto hess_diag = get_hessian_diagonal(q_point);
4887 for (
unsigned int comp = 0; comp < n_components; ++comp)
4889 laplacian_out[comp] = hess_diag[comp][0];
4890 for (
unsigned int d = 1;
d < dim; ++
d)
4891 laplacian_out[comp] += hess_diag[comp][
d];
4893 return laplacian_out;
4902 typename VectorizedArrayType>
4906 const unsigned int dof)
4909 this->dof_values_initialized =
true;
4911 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4913 for (
unsigned int comp = 0; comp < n_components; ++comp)
4914 this->values_dofs[comp * dofs + dof] = val_in[comp];
4923 typename VectorizedArrayType>
4927 const unsigned int q_point)
4933 Assert(this->J_value !=
nullptr,
4937 this->values_quad_submitted =
true;
4940 const std::size_t nqp = this->n_quadrature_points;
4943 const VectorizedArrayType JxW =
4944 this->J_value[0] * this->quadrature_weights[q_point];
4945 for (
unsigned int comp = 0; comp < n_components; ++comp)
4946 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4950 const VectorizedArrayType JxW = this->J_value[q_point];
4951 for (
unsigned int comp = 0; comp < n_components; ++comp)
4952 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4962 typename VectorizedArrayType>
4967 const unsigned int q_point)
4973 Assert(this->J_value !=
nullptr,
4975 "update_gradients"));
4976 Assert(this->jacobian !=
nullptr,
4978 "update_gradients"));
4980 this->gradients_quad_submitted =
true;
4983 const std::size_t nqp = this->n_quadrature_points;
4986 const VectorizedArrayType JxW =
4987 this->J_value[0] * this->quadrature_weights[q_point];
4988 std::array<VectorizedArrayType, dim> jac;
4989 for (
unsigned int d = 0;
d < dim; ++
d)
4990 jac[
d] = this->jacobian[0][
d][
d];
4991 for (
unsigned int d = 0;
d < dim; ++
d)
4993 const VectorizedArrayType factor = jac[
d] * JxW;
4994 for (
unsigned int comp = 0; comp < n_components; ++comp)
4995 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
4996 grad_in[comp][
d] * factor;
5003 this->jacobian[q_point] :
5005 const VectorizedArrayType JxW =
5007 this->J_value[q_point] :
5008 this->J_value[0] * this->quadrature_weights[q_point];
5009 for (
unsigned int comp = 0; comp < n_components; ++comp)
5010 for (
unsigned int d = 0;
d < dim; ++
d)
5012 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5013 for (
unsigned int e = 1;
e < dim; ++
e)
5014 new_val += (jac[
e][
d] * grad_in[comp][
e]);
5015 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5027 typename VectorizedArrayType>
5032 const unsigned int q_point)
5035 Assert(this->normal_x_jacobian !=
nullptr,
5037 "update_gradients"));
5039 this->gradients_quad_submitted =
true;
5042 const std::size_t nqp = this->n_quadrature_points;
5045 const VectorizedArrayType JxW_jac = this->J_value[0] *
5046 this->quadrature_weights[q_point] *
5047 this->normal_x_jacobian[0][dim - 1];
5048 for (
unsigned int comp = 0; comp < n_components; ++comp)
5050 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5051 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5052 VectorizedArrayType();
5053 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5054 grad_in[comp] * JxW_jac;
5059 const unsigned int index =
5062 this->normal_x_jacobian[
index];
5063 const VectorizedArrayType JxW =
5065 this->J_value[index] * this->quadrature_weights[q_point] :
5066 this->J_value[index];
5067 for (
unsigned int comp = 0; comp < n_components; ++comp)
5069 for (
unsigned int d = 0;
d < dim; ++
d)
5070 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5071 (grad_in[comp] * JxW) * jac[
d];
5082 typename VectorizedArrayType>
5088 const unsigned int q_point)
5094 Assert(this->J_value !=
nullptr,
5096 "update_hessians"));
5097 Assert(this->jacobian !=
nullptr,
5099 "update_hessians"));
5101 this->hessians_quad_submitted =
true;
5105 const std::size_t nqp = this->n_quadrature_points;
5106 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
5109 const VectorizedArrayType JxW =
5110 this->J_value[0] * this->quadrature_weights[q_point];
5113 for (
unsigned int d = 0;
d < dim; ++
d)
5115 const auto jac_d = this->jacobian[0][
d][
d];
5116 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5117 for (
unsigned int comp = 0; comp < n_components; ++comp)
5118 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5119 hessian_in[comp][
d][
d] * factor;
5123 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5124 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5126 const auto jac_d = this->jacobian[0][
d][
d];
5127 const auto jac_e = this->jacobian[0][
e][
e];
5128 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5129 for (
unsigned int comp = 0; comp < n_components; ++comp)
5130 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5131 (hessian_in[comp][
d][
e] + hessian_in[comp][
e][
d]) * factor;
5138 const VectorizedArrayType JxW =
5139 this->J_value[0] * this->quadrature_weights[q_point];
5140 for (
unsigned int comp = 0; comp < n_components; ++comp)
5143 VectorizedArrayType tmp[dim][dim];
5144 for (
unsigned int i = 0; i < dim; ++i)
5145 for (
unsigned int j = 0; j < dim; ++j)
5147 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5148 for (
unsigned int k = 1; k < dim; ++k)
5149 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5153 VectorizedArrayType tmp2[dim][dim];
5154 for (
unsigned int i = 0; i < dim; ++i)
5155 for (
unsigned int j = 0; j < dim; ++j)
5157 tmp2[i][j] = jac[0][i] * tmp[0][j];
5158 for (
unsigned int k = 1; k < dim; ++k)
5159 tmp2[i][j] += jac[k][i] * tmp[k][j];
5163 for (
unsigned int d = 0;
d < dim; ++
d)
5164 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5168 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5169 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5170 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5171 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5177 const VectorizedArrayType JxW = this->J_value[q_point];
5178 const auto &jac_grad = this->jacobian_gradients[q_point];
5179 for (
unsigned int comp = 0; comp < n_components; ++comp)
5182 VectorizedArrayType tmp[dim][dim];
5183 for (
unsigned int i = 0; i < dim; ++i)
5184 for (
unsigned int j = 0; j < dim; ++j)
5186 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5187 for (
unsigned int k = 1; k < dim; ++k)
5188 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5192 VectorizedArrayType tmp2[dim][dim];
5193 for (
unsigned int i = 0; i < dim; ++i)
5194 for (
unsigned int j = 0; j < dim; ++j)
5196 tmp2[i][j] = jac[0][i] * tmp[0][j];
5197 for (
unsigned int k = 1; k < dim; ++k)
5198 tmp2[i][j] += jac[k][i] * tmp[k][j];
5202 for (
unsigned int d = 0;
d < dim; ++
d)
5203 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5207 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5208 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5209 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5210 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5213 for (
unsigned int d = 0;
d < dim; ++
d)
5215 VectorizedArrayType
sum = 0;
5216 for (
unsigned int e = 0;
e < dim; ++
e)
5217 sum += hessian_in[comp][
e][
e] * jac_grad[
e][
d];
5218 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5219 for (
unsigned int f =
e + 1; f < dim; ++f, ++count)
5220 sum += (hessian_in[comp][
e][f] + hessian_in[comp][f][
e]) *
5222 this->gradients_from_hessians_quad[(comp * dim +
d) * nqp +
5223 q_point] =
sum * JxW;
5235 typename VectorizedArrayType>
5242 Assert(this->values_quad_submitted ==
true,
5247 const std::size_t nqp = this->n_quadrature_points;
5248 for (
unsigned int q = 0; q < nqp; ++q)
5249 for (
unsigned int comp = 0; comp < n_components; ++comp)
5250 return_value[comp] += this->values_quad[comp * nqp + q];
5251 return (return_value);
5263 typename VectorizedArrayType>
5268 VectorizedArrayType>::
5271 const unsigned int dof_no,
5272 const unsigned int first_selected_component,
5273 const unsigned int quad_no,
5274 const unsigned int fe_degree,
5275 const unsigned int n_q_points,
5276 const bool is_interior_face,
5277 const unsigned int active_fe_index,
5278 const unsigned int active_quad_index,
5279 const unsigned int face_type)
5280 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5283 first_selected_component,
5299 typename VectorizedArrayType>
5304 VectorizedArrayType>::
5310 const unsigned int first_selected_component,
5312 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5317 first_selected_component,
5327 typename VectorizedArrayType>
5332 VectorizedArrayType>::
5337 VectorizedArrayType> &other)
5338 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5348 typename VectorizedArrayType>
5353 VectorizedArrayType> &
5359 VectorizedArrayType> &other)
5365 VectorizedArrayType>::operator=(other);
5374 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5378 const unsigned int dof_no,
5379 const unsigned int first_selected_component,
5380 const unsigned int quad_no,
5381 const unsigned int fe_degree,
5382 const unsigned int n_q_points,
5383 const bool is_interior_face,
5384 const unsigned int active_fe_index,
5385 const unsigned int active_quad_index,
5386 const unsigned int face_type)
5390 first_selected_component,
5402 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5409 const unsigned int first_selected_component,
5416 first_selected_component,
5422 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5432 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5438 ->
FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5445 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5448 const unsigned int dof)
const
5451 return this->values_dofs[dof];
5456 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5459 const unsigned int q_point)
const
5462 Assert(this->values_quad_initialized ==
true,
5466 return this->values_quad[q_point];
5471 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5476 return BaseClass::get_normal_derivative(q_point)[0];
5481 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5484 const unsigned int q_point)
const
5490 Assert(this->gradients_quad_initialized ==
true,
5495 Assert(this->jacobian !=
nullptr,
5497 "update_gradients"));
5501 const std::size_t nqp = this->n_quadrature_points;
5504 for (
unsigned int d = 0;
d < dim; ++
d)
5506 this->gradients_quad[
d * nqp + q_point] * this->jacobian[0][
d][
d];
5515 for (
unsigned int d = 0;
d < dim; ++
d)
5517 grad_out[
d] = jac[
d][0] * this->gradients_quad[q_point];
5518 for (
unsigned int e = 1;
e < dim; ++
e)
5519 grad_out[
d] += jac[
d][
e] * this->gradients_quad[
e * nqp + q_point];
5527 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5530 const unsigned int q_point)
const
5532 return BaseClass::get_hessian(q_point)[0];
5537 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5542 return BaseClass::get_hessian_diagonal(q_point)[0];
5547 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5548 inline VectorizedArrayType
5550 const unsigned int q_point)
const
5552 return BaseClass::get_laplacian(q_point)[0];
5557 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5563 this->dof_values_initialized =
true;
5566 this->values_dofs[dof] = val_in;
5571 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5574 const VectorizedArrayType val_in,
5575 const unsigned int q_point)
5581 Assert(this->J_value !=
nullptr,
5585 this->values_quad_submitted =
true;
5590 const VectorizedArrayType JxW =
5591 this->J_value[0] * this->quadrature_weights[q_point];
5592 this->values_quad[q_point] = val_in * JxW;
5596 this->values_quad[q_point] = val_in * this->J_value[q_point];
5602 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5606 const unsigned int q_point)
5608 submit_value(val_in[0], q_point);
5613 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5617 const unsigned int q_point)
5621 BaseClass::submit_normal_derivative(grad, q_point);
5626 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5630 const unsigned int q_point)
5636 Assert(this->J_value !=
nullptr,
5638 "update_gradients"));
5639 Assert(this->jacobian !=
nullptr,
5641 "update_gradients"));
5643 this->gradients_quad_submitted =
true;
5646 const std::size_t nqp = this->n_quadrature_points;
5649 const VectorizedArrayType JxW =
5650 this->J_value[0] * this->quadrature_weights[q_point];
5654 std::array<VectorizedArrayType, dim> jac;
5655 for (
unsigned int d = 0;
d < dim; ++
d)
5656 jac[
d] = this->jacobian[0][
d][
d];
5658 for (
unsigned int d = 0;
d < dim; ++
d)
5659 this->gradients_quad[
d * nqp + q_point] = grad_in[
d] * jac[
d] * JxW;
5666 this->jacobian[q_point] :
5668 const VectorizedArrayType JxW =
5670 this->J_value[q_point] :
5671 this->J_value[0] * this->quadrature_weights[q_point];
5672 for (
unsigned int d = 0;
d < dim; ++
d)
5674 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5675 for (
unsigned int e = 1;
e < dim; ++
e)
5676 new_val += jac[
e][
d] * grad_in[
e];
5677 this->gradients_quad[
d * nqp + q_point] = new_val * JxW;
5684 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5688 const unsigned int q_point)
5692 BaseClass::submit_hessian(hessian, q_point);
5697 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5698 inline VectorizedArrayType
5702 return BaseClass::integrate_value()[0];
5710 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5714 const unsigned int dof_no,
5715 const unsigned int first_selected_component,
5716 const unsigned int quad_no,
5717 const unsigned int fe_degree,
5718 const unsigned int n_q_points,
5719 const bool is_interior_face,
5720 const unsigned int active_fe_index,
5721 const unsigned int active_quad_index,
5722 const unsigned int face_type)
5726 first_selected_component,
5738 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5745 const unsigned int first_selected_component,
5752 first_selected_component,
5758 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5768 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5780 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5783 const unsigned int q_point)
const
5785 if (this->data->element_type ==
5790 Assert(this->values_quad_initialized ==
true,
5795 Assert(this->J_value !=
nullptr,
5798 const std::size_t nqp = this->n_quadrature_points;
5806 const VectorizedArrayType inv_det =
5807 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5808 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5809 this->jacobian[0][2][2];
5812 for (
unsigned int comp = 0; comp < n_components; ++comp)
5813 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5814 jac[comp][comp] * inv_det;
5821 this->jacobian[q_point] :
5829 const VectorizedArrayType inv_det =
5830 (is_face && dim == 2 && this->get_face_no() < 2) ?
5834 for (
unsigned int comp = 0; comp < n_components; ++comp)
5837 this->values_quad[q_point] * jac[comp][0] * inv_det;
5838 for (
unsigned int e = 1;
e < dim; ++
e)
5840 this->values_quad[
e * nqp + q_point] * jac[comp][
e] * inv_det;
5848 return BaseClass::get_value(q_point);
5852 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5857 if (this->data->element_type ==
5862 Assert(this->gradients_quad_initialized ==
true,
5867 Assert(this->jacobian !=
nullptr,
5869 "update_gradients"));
5870 const std::size_t nqp = this->n_quadrature_points;
5880 const VectorizedArrayType inv_det =
5881 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5882 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5883 this->jacobian[0][2][2];
5886 for (
unsigned int d = 0;
d < dim; ++
d)
5887 for (
unsigned int comp = 0; comp < n_components; ++comp)
5889 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
5890 inv_t_jac[
d][
d] * jac[comp][comp] * inv_det;
5900 const VectorizedArrayType inv_det =
5901 (is_face && dim == 2 && this->get_face_no() < 2) ?
5905 VectorizedArrayType tmp;
5907 for (
unsigned int comp = 0; comp < n_components; ++comp)
5908 for (
unsigned int d = 0;
d < dim; ++
d)
5911 for (
unsigned int f = 0; f < dim; ++f)
5912 for (
unsigned int e = 0;
e < dim; ++
e)
5913 tmp += jac[comp][f] * inv_t_jac[
d][
e] * inv_det *
5914 this->gradients_quad[(f * dim +
e) * nqp + q_point];
5916 grad_out[comp][
d] = tmp;
5927 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
5929 "update_hessians"));
5931 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
5933 this->jacobian[q_point];
5937 const VectorizedArrayType inv_det =
5938 (is_face && dim == 2 && this->get_face_no() < 2) ?
5942 VectorizedArrayType tmp;
5944 for (
unsigned int comp = 0; comp < n_components; ++comp)
5945 for (
unsigned int d = 0;
d < dim; ++
d)
5948 for (
unsigned int f = 0; f < dim; ++f)
5949 for (
unsigned int e = 0;
e < dim; ++
e)
5950 tmp += t_jac[f][comp] * inv_t_jac[
d][
e] *
5951 this->gradients_quad[(f * dim +
e) * nqp + q_point];
5953 grad_out[comp][
d] = tmp * inv_det;
5962 for (
unsigned int i = 0; i < dim; ++i)
5963 for (
unsigned int j = 0; j < dim; ++j)
5965 tmp = jac_grad[0][i] * inv_t_jac[j][0] *
5966 this->values_quad[q_point];
5967 for (
unsigned int f = 1; f < dim; ++f)
5968 tmp += jac_grad[f][i] * inv_t_jac[j][f] *
5969 this->values_quad[f * nqp + q_point];
5971 grad_out[i][j] += tmp * inv_det;
5974 for (
unsigned int i = 0; i < dim; ++i)
5975 for (
unsigned int j = 0; j < dim; ++j)
5978 for (
unsigned int f = 0; f < dim; ++f)
5979 for (
unsigned int n = 0; n < dim; ++n)
5980 for (
unsigned int m = 0; m < dim; ++m)
5981 tmp += inv_t_jac[m][f] * jac_grad[f][m] *
5982 inv_t_jac[j][f] * t_jac[n][i] *
5983 this->values_quad[n * nqp + q_point];
5984 grad_out[i][j] -= tmp * inv_det;
5993 for (
unsigned int i = 0; i < dim; ++i)
5994 for (
unsigned int j = 0; j < dim; ++j)
5997 for (
unsigned int r = 0, f = dim; r < dim; ++r)
5998 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6000 tmp += jac_grad[f][i] *
6002 this->values_quad[r * nqp + q_point] +
6004 this->values_quad[k * nqp + q_point]);
6005 for (
unsigned int n = 0; n < dim; ++n)
6006 for (
unsigned int m = 0; m < dim; ++m)
6007 tmp -= jac_grad[f][m] * t_jac[n][i] *
6008 this->values_quad[n * nqp + q_point] *
6009 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6010 inv_t_jac[m][r] * inv_t_jac[j][k]);
6012 grad_out[i][j] += tmp * inv_det;
6020 return BaseClass::get_gradient(q_point);
6026 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6032 Assert(this->gradients_quad_initialized ==
true,
6036 Assert(this->jacobian !=
nullptr,
6038 "update_gradients"));
6040 VectorizedArrayType divergence;
6041 const std::size_t nqp = this->n_quadrature_points;
6043 if (this->data->element_type ==
6050 const VectorizedArrayType inv_det =
6051 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
6052 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
6053 this->jacobian[0][2][2];
6056 divergence = this->gradients_quad[q_point] * inv_det;
6057 for (
unsigned int d = 1;
d < dim; ++
d)
6059 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
6065 const VectorizedArrayType inv_det =
6067 this->jacobian[this->cell_type >
6071 Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6074 divergence = this->gradients_quad[q_point] * inv_det;
6075 for (
unsigned int d = 1;
d < dim; ++
d)
6077 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
6086 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6087 for (
unsigned int d = 1;
d < dim; ++
d)
6088 divergence += this->gradients_quad[(dim *
d +
d) * nqp + q_point] *
6089 this->jacobian[0][
d][
d];
6096 this->jacobian[q_point] :
6098 divergence = jac[0][0] * this->gradients_quad[q_point];
6099 for (
unsigned int e = 1;
e < dim; ++
e)
6100 divergence += jac[0][
e] * this->gradients_quad[
e * nqp + q_point];
6101 for (
unsigned int d = 1;
d < dim; ++
d)
6102 for (
unsigned int e = 0;
e < dim; ++
e)
6104 jac[
d][
e] * this->gradients_quad[(
d * dim +
e) * nqp + q_point];
6112 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6118 const auto grad = get_gradient(q_point);
6119 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6120 VectorizedArrayType half = Number(0.5);
6121 for (
unsigned int d = 0;
d < dim; ++
d)
6122 symmetrized[
d] = grad[
d][
d];
6128 symmetrized[2] = grad[0][1] + grad[1][0];
6129 symmetrized[2] *= half;
6132 symmetrized[3] = grad[0][1] + grad[1][0];
6133 symmetrized[3] *= half;
6134 symmetrized[4] = grad[0][2] + grad[2][0];
6135 symmetrized[4] *= half;
6136 symmetrized[5] = grad[1][2] + grad[2][1];
6137 symmetrized[5] *= half;
6147 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6149 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6151 const unsigned int q_point)
const
6155 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6161 "Computing the curl in 1d is not a useful operation"));
6164 curl[0] = grad[1][0] - grad[0][1];
6167 curl[0] = grad[2][1] - grad[1][2];
6168 curl[1] = grad[0][2] - grad[2][0];
6169 curl[2] = grad[1][0] - grad[0][1];
6179 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6184 return BaseClass::get_hessian_diagonal(q_point);
6189 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6192 const unsigned int q_point)
const
6195 Assert(this->hessians_quad_initialized ==
true,
6199 return BaseClass::get_hessian(q_point);
6203 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6207 const unsigned int q_point)
6209 if (this->data->element_type ==
6214 Assert(this->J_value !=
nullptr,
6219 this->values_quad_submitted =
true;
6222 const std::size_t nqp = this->n_quadrature_points;
6227 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6229 for (
unsigned int comp = 0; comp < n_components; ++comp)
6230 this->values_quad[comp * nqp + q_point] =
6231 val_in[comp] * weight * jac[comp][comp];
6238 this->jacobian[q_point] :
6247 const VectorizedArrayType fac =
6249 this->quadrature_weights[q_point] :
6251 this->J_value[q_point] :
6252 this->J_value[0] * this->quadrature_weights[q_point]) *
6253 ((dim == 2 && this->get_face_no() < 2) ?
6258 for (
unsigned int comp = 0; comp < n_components; ++comp)
6260 this->values_quad[comp * nqp + q_point] =
6261 val_in[0] * jac[0][comp] * fac;
6262 for (
unsigned int e = 1;
e < dim; ++
e)
6263 this->values_quad[comp * nqp + q_point] +=
6264 val_in[
e] * jac[
e][comp] * fac;
6271 BaseClass::submit_value(val_in, q_point);
6275 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6279 const unsigned int q_point)
6281 if (this->data->element_type ==
6290 Assert(this->J_value !=
nullptr,
6292 "update_gradients"));
6293 Assert(this->jacobian !=
nullptr,
6295 "update_gradients"));
6297 this->gradients_quad_submitted =
true;
6300 const std::size_t nqp = this->n_quadrature_points;
6308 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6309 for (
unsigned int d = 0;
d < dim; ++
d)
6310 for (
unsigned int comp = 0; comp < n_components; ++comp)
6311 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6312 grad_in[comp][
d] * inv_t_jac[
d][
d] * jac[comp][comp] * weight;
6323 const VectorizedArrayType fac =
6324 (!is_face) ? this->quadrature_weights[q_point] :
6325 this->J_value[0] * this->quadrature_weights[q_point] *
6326 ((dim == 2 && this->get_face_no() < 2) ?
6331 for (
unsigned int comp = 0; comp < n_components; ++comp)
6332 for (
unsigned int d = 0;
d < dim; ++
d)
6334 VectorizedArrayType tmp = 0;
6335 for (
unsigned int f = 0; f < dim; ++f)
6336 for (
unsigned int e = 0;
e < dim; ++
e)
6337 tmp += jac[f][comp] * inv_t_jac[
e][
d] * grad_in[f][
e];
6339 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6347 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
6349 this->jacobian[q_point];
6354 const VectorizedArrayType fac =
6356 this->quadrature_weights[q_point] :
6357 this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6361 VectorizedArrayType tmp;
6363 for (
unsigned int comp = 0; comp < n_components; ++comp)
6364 for (
unsigned int d = 0;
d < dim; ++
d)
6367 for (
unsigned int f = 0; f < dim; ++f)
6368 for (
unsigned int e = 0;
e < dim; ++
e)
6369 tmp += t_jac[comp][f] * inv_t_jac[
e][
d] * grad_in[f][
e];
6371 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6381 for (
unsigned int f = 0; f < dim; ++f)
6384 for (
unsigned int i = 0; i < dim; ++i)
6385 for (
unsigned int j = 0; j < dim; ++j)
6387 tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
6388 for (
unsigned int m = 0; m < dim; ++m)
6389 for (
unsigned int k = 0; k < dim; ++k)
6390 tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
6391 inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
6393 this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
6401 for (
unsigned int r = 0, f = dim; r < dim; ++r)
6402 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6404 tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
6405 for (
unsigned int j = 1; j < dim; ++j)
6406 tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
6407 for (
unsigned int i = 1; i < dim; ++i)
6408 for (
unsigned int j = 0; j < dim; ++j)
6409 tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
6410 this->values_from_gradients_quad[r * nqp + q_point] +=
6413 tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
6414 for (
unsigned int j = 1; j < dim; ++j)
6415 tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
6416 for (
unsigned int i = 1; i < dim; ++i)
6417 for (
unsigned int j = 0; j < dim; ++j)
6418 tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
6419 this->values_from_gradients_quad[k * nqp + q_point] +=
6424 for (
unsigned int n = 0; n < dim; ++n)
6427 for (
unsigned int r = 0, f = dim; r < dim; ++r)
6428 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6429 for (
unsigned int i = 0; i < dim; ++i)
6430 for (
unsigned int j = 0; j < dim; ++j)
6431 for (
unsigned int m = 0; m < dim; ++m)
6432 tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
6433 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6434 inv_t_jac[m][r] * inv_t_jac[j][k]);
6436 this->values_from_gradients_quad[n * nqp + q_point] -=
6444 BaseClass::submit_gradient(grad_in, q_point);
6450 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6455 const unsigned int q_point)
6457 if (this->data->element_type ==
6467 BaseClass::submit_gradient(grad_in, q_point);
6473 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6477 const unsigned int q_point)
6483 Assert(this->J_value !=
nullptr,
6485 "update_gradients"));
6486 Assert(this->jacobian !=
nullptr,
6488 "update_gradients"));
6490 this->gradients_quad_submitted =
true;
6493 const std::size_t nqp = this->n_quadrature_points;
6494 if (this->data->element_type ==
6501 const VectorizedArrayType fac =
6503 this->quadrature_weights[q_point] * div_in :
6505 this->J_value[q_point] :
6506 this->J_value[0] * this->quadrature_weights[q_point]) *
6509 this->jacobian[this->cell_type >
6513 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6515 for (
unsigned int d = 0;
d < dim; ++
d)
6517 this->gradients_quad[(dim *
d +
d) * nqp + q_point] = fac;
6518 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6520 this->gradients_quad[(dim *
d +
e) * nqp + q_point] =
6521 VectorizedArrayType();
6522 this->gradients_quad[(dim *
e +
d) * nqp + q_point] =
6523 VectorizedArrayType();
6526 this->divergence_is_requested =
true;
6533 const VectorizedArrayType fac =
6534 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6535 for (
unsigned int d = 0;
d < dim; ++
d)
6537 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6538 (fac * this->jacobian[0][
d][
d]);
6539 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6541 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6542 VectorizedArrayType();
6543 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6544 VectorizedArrayType();
6552 this->jacobian[q_point] :
6554 const VectorizedArrayType fac =
6556 this->J_value[q_point] :
6557 this->J_value[0] * this->quadrature_weights[q_point]) *
6559 for (
unsigned int d = 0;
d < dim; ++
d)
6561 for (
unsigned int e = 0;
e < dim; ++
e)
6562 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6571 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6576 const unsigned int q_point)
6579 this->data->element_type !=
6590 Assert(this->J_value !=
nullptr,
6592 "update_gradients"));
6593 Assert(this->jacobian !=
nullptr,
6595 "update_gradients"));
6597 this->gradients_quad_submitted =
true;
6600 const std::size_t nqp = this->n_quadrature_points;
6603 const VectorizedArrayType JxW =
6604 this->J_value[0] * this->quadrature_weights[q_point];
6605 for (
unsigned int d = 0;
d < dim; ++
d)
6606 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6608 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6609 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6611 const VectorizedArrayType
value =
6613 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6614 value * this->jacobian[0][
d][
d];
6615 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6616 value * this->jacobian[0][
e][
e];
6622 const VectorizedArrayType JxW =
6624 this->J_value[q_point] :
6625 this->J_value[0] * this->quadrature_weights[q_point];
6628 this->jacobian[q_point] :
6630 VectorizedArrayType weighted[dim][dim];
6631 for (
unsigned int i = 0; i < dim; ++i)
6633 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6634 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6636 const VectorizedArrayType
value =
6638 weighted[i][j] =
value;
6639 weighted[j][i] =
value;
6641 for (
unsigned int comp = 0; comp < dim; ++comp)
6642 for (
unsigned int d = 0;
d < dim; ++
d)
6644 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6645 for (
unsigned int e = 1;
e < dim; ++
e)
6646 new_val += jac[
e][
d] * weighted[comp][
e];
6647 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
6654 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6658 const unsigned int q_point)
6666 "Testing by the curl in 1d is not a useful operation"));
6669 grad[1][0] = curl[0];
6670 grad[0][1] = -curl[0];
6673 grad[2][1] = curl[0];
6674 grad[1][2] = -curl[0];
6675 grad[0][2] = curl[1];
6676 grad[2][0] = -curl[1];
6677 grad[1][0] = curl[2];
6678 grad[0][1] = -curl[2];
6683 submit_gradient(grad, q_point);
6690 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6694 const unsigned int dof_no,
6695 const unsigned int first_selected_component,
6696 const unsigned int quad_no,
6697 const unsigned int fe_degree,
6698 const unsigned int n_q_points,
6699 const bool is_interior_face,
6700 const unsigned int active_fe_index,
6701 const unsigned int active_quad_index,
6702 const unsigned int face_type)
6706 first_selected_component,
6718 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6725 const unsigned int first_selected_component,
6732 first_selected_component,
6738 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6747 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6759 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6762 const unsigned int dof)
const
6765 return this->values_dofs[dof];
6770 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6773 const unsigned int q_point)
const
6776 Assert(this->values_quad_initialized ==
true,
6780 return this->values_quad[q_point];
6785 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6788 const unsigned int q_point)
const
6794 Assert(this->gradients_quad_initialized ==
true,
6801 this->jacobian[q_point] :
6805 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6812 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6815 const unsigned int q_point)
const
6817 return get_gradient(q_point)[0];
6822 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6827 return BaseClass::get_normal_derivative(q_point)[0];
6832 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6835 const unsigned int q_point)
const
6837 return BaseClass::get_hessian(q_point)[0];
6842 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6847 return BaseClass::get_hessian_diagonal(q_point)[0];
6852 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6855 const unsigned int q_point)
const
6857 return BaseClass::get_laplacian(q_point)[0];
6862 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6868 this->dof_values_initialized =
true;
6871 this->values_dofs[dof] = val_in;
6876 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6879 const VectorizedArrayType val_in,
6880 const unsigned int q_point)
6887 this->values_quad_submitted =
true;
6892 const VectorizedArrayType JxW = this->J_value[q_point];
6893 this->values_quad[q_point] = val_in * JxW;
6897 const VectorizedArrayType JxW =
6898 this->J_value[0] * this->quadrature_weights[q_point];
6899 this->values_quad[q_point] = val_in * JxW;
6905 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6909 const unsigned int q_point)
6911 submit_value(val_in[0], q_point);
6916 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6920 const unsigned int q_point)
6922 submit_gradient(grad_in[0], q_point);
6927 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6930 const VectorizedArrayType grad_in,
6931 const unsigned int q_point)
6938 this->gradients_quad_submitted =
true;
6943 this->jacobian[q_point] :
6945 const VectorizedArrayType JxW =
6947 this->J_value[q_point] :
6948 this->J_value[0] * this->quadrature_weights[q_point];
6950 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6955 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6959 const unsigned int q_point)
6961 submit_gradient(grad_in[0][0], q_point);
6966 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6970 const unsigned int q_point)
6974 BaseClass::submit_normal_derivative(grad, q_point);
6979 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6983 const unsigned int q_point)
6985 BaseClass::submit_normal_derivative(grad_in, q_point);
6989 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6993 const unsigned int q_point)
6997 BaseClass::submit_hessian(hessian, q_point);
7001 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7002 inline VectorizedArrayType
7006 return BaseClass::integrate_value()[0];
7019 typename VectorizedArrayType>
7025 VectorizedArrayType>::
7027 const unsigned int fe_no,
7028 const unsigned int quad_no,
7029 const unsigned int first_selected_component,
7030 const unsigned int active_fe_index,
7031 const unsigned int active_quad_index)
7032 : BaseClass(matrix_free,
7034 first_selected_component,
7041 , dofs_per_component(this->data->dofs_per_component_on_cell)
7042 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7043 , n_q_points(this->data->n_q_points)
7045 check_template_arguments(fe_no, 0);
7055 typename VectorizedArrayType>
7061 VectorizedArrayType>::
7063 const std::pair<unsigned int, unsigned int> & range,
7064 const unsigned int dof_no,
7065 const unsigned int quad_no,
7066 const unsigned int first_selected_component)
7070 first_selected_component,
7071 matrix_free.get_cell_active_fe_index(range))
7081 typename VectorizedArrayType>
7087 VectorizedArrayType>::
7092 const unsigned int first_selected_component)
7093 : BaseClass(mapping,
7097 first_selected_component,
7099 , dofs_per_component(this->data->dofs_per_component_on_cell)
7100 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7101 , n_q_points(this->data->n_q_points)
7113 typename VectorizedArrayType>
7119 VectorizedArrayType>::
7123 const unsigned int first_selected_component)
7128 first_selected_component,
7130 , dofs_per_component(this->data->dofs_per_component_on_cell)
7131 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7132 , n_q_points(this->data->n_q_points)
7144 typename VectorizedArrayType>
7150 VectorizedArrayType>::
7153 const unsigned int first_selected_component)
7154 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7156 other.mapped_geometry->get_quadrature(),
7157 other.mapped_geometry->get_fe_values().get_update_flags(),
7158 first_selected_component,
7160 , dofs_per_component(this->data->dofs_per_component_on_cell)
7161 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7162 , n_q_points(this->data->n_q_points)
7174 typename VectorizedArrayType>
7183 , dofs_per_component(this->data->dofs_per_component_on_cell)
7184 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7185 , n_q_points(this->data->n_q_points)
7197 typename VectorizedArrayType>
7203 VectorizedArrayType> &
7209 VectorizedArrayType>::operator=(
const FEEvaluation &other)
7211 BaseClass::operator=(other);
7223 typename VectorizedArrayType>
7230 VectorizedArrayType>::
7231 check_template_arguments(
const unsigned int dof_no,
7232 const unsigned int first_selected_component)
7235 (void)first_selected_component;
7238 this->data->dofs_per_component_on_cell > 0,
7240 "There is nothing useful you can do with an FEEvaluation object with "
7241 "FE_Nothing, i.e., without DoFs! If you have passed to "
7242 "MatrixFree::reinit() a collection of finite elements also containing "
7243 "FE_Nothing, please check - before creating FEEvaluation - the category "
7244 "of the current range by calling either "
7245 "MatrixFree::get_cell_range_category(range) or "
7246 "MatrixFree::get_face_range_category(range). The returned category "
7247 "is the index of the active FE, which you can use to exclude "
7254 static_cast<unsigned int>(fe_degree) !=
7255 this->data->data.front().fe_degree) ||
7256 n_q_points != this->n_quadrature_points)
7258 std::string message =
7259 "-------------------------------------------------------\n";
7260 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7261 message +=
" Called --> FEEvaluation<dim,";
7265 message +=
",Number>(data";
7281 if (
static_cast<unsigned int>(fe_degree) ==
7282 this->data->data.front().fe_degree)
7284 proposed_dof_comp = dof_no;
7285 proposed_fe_comp = first_selected_component;
7288 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
7290 for (
unsigned int nf = 0;
7293 if (this->matrix_free
7294 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7296 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7298 proposed_dof_comp = no;
7299 proposed_fe_comp = nf;
7303 this->mapping_data->descriptor[this->active_quad_index]
7305 proposed_quad_comp = this->quad_no;
7307 for (
unsigned int no = 0;
7312 .descriptor[this->active_quad_index]
7313 .n_q_points == n_q_points)
7315 proposed_quad_comp = no;
7322 if (proposed_dof_comp != first_selected_component)
7323 message +=
"Wrong vector component selection:\n";
7325 message +=
"Wrong quadrature formula selection:\n";
7326 message +=
" Did you mean FEEvaluation<dim,";
7330 message +=
",Number>(data";
7339 std::string correct_pos;
7340 if (proposed_dof_comp != dof_no)
7341 correct_pos =
" ^ ";
7344 if (proposed_quad_comp != this->quad_no)
7345 correct_pos +=
" ^ ";
7348 if (proposed_fe_comp != first_selected_component)
7349 correct_pos +=
" ^\n";
7351 correct_pos +=
" \n";
7357 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7358 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7359 message +=
"Wrong template arguments:\n";
7360 message +=
" Did you mean FEEvaluation<dim,";
7365 message +=
",Number>(data";
7373 std::string correct_pos;
7374 if (this->data->data.front().fe_degree !=
7375 static_cast<unsigned int>(fe_degree))
7379 if (proposed_n_q_points_1d != n_q_points_1d)
7380 correct_pos +=
" ^\n";
7382 correct_pos +=
" \n";
7383 message +=
" " + correct_pos;
7385 Assert(
static_cast<unsigned int>(fe_degree) ==
7386 this->data->data.front().fe_degree &&
7387 n_q_points == this->n_quadrature_points,
7393 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7404 typename VectorizedArrayType>
7413 Assert(this->matrix_free !=
nullptr,
7414 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7415 " Integer indexing is not possible."));
7423 const unsigned int offsets =
7424 this->mapping_data->data_index_offsets[
cell_index];
7425 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7426 this->J_value = &this->mapping_data->JxW_values[offsets];
7427 if (!this->mapping_data->jacobian_gradients[0].empty())
7429 this->jacobian_gradients =
7430 this->mapping_data->jacobian_gradients[0].data() + offsets;
7431 this->jacobian_gradients_non_inverse =
7432 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
7436 VectorizedArrayType::size())
7439 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
7440 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7447 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7448 for (; i < VectorizedArrayType::size(); ++i)
7452 if (this->mapping_data->quadrature_points.empty() ==
false)
7454 &this->mapping_data->quadrature_points
7455 [this->mapping_data->quadrature_point_offsets[this->cell]];
7458 this->is_reinitialized =
true;
7459 this->dof_values_initialized =
false;
7460 this->values_quad_initialized =
false;
7461 this->gradients_quad_initialized =
false;
7462 this->hessians_quad_initialized =
false;
7473 typename VectorizedArrayType>
7480 VectorizedArrayType>
::
7481 reinit(
const std::array<
unsigned int, VectorizedArrayType::size()> &cell_ids)
7487 this->cell_ids = cell_ids;
7492 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7506 if (this->mapped_geometry ==
nullptr)
7507 this->mapped_geometry =
7508 std::make_shared<internal::MatrixFreeFunctions::
7509 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7511 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7513 auto &this_jacobian_data = mapping_storage.jacobians[0];
7514 auto &this_J_value_data = mapping_storage.JxW_values;
7515 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7516 auto &this_jacobian_gradients_non_inverse_data =
7517 mapping_storage.jacobian_gradients_non_inverse[0];
7518 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7522 if (this->mapping_data->jacobians[0].size() > 0)
7523 this_jacobian_data.resize_fast(2);
7525 if (this->mapping_data->JxW_values.size() > 0)
7526 this_J_value_data.resize_fast(1);
7528 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7529 this_jacobian_gradients_data.resize_fast(1);
7531 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7532 this_jacobian_gradients_non_inverse_data.resize_fast(1);
7534 if (this->mapping_data->quadrature_points.size() > 0)
7535 this_quadrature_points_data.resize_fast(1);
7539 if (this->mapping_data->jacobians[0].size() > 0)
7540 this_jacobian_data.resize_fast(this->n_quadrature_points);
7542 if (this->mapping_data->JxW_values.size() > 0)
7543 this_J_value_data.resize_fast(this->n_quadrature_points);
7545 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7546 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7548 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7549 this_jacobian_gradients_non_inverse_data.resize_fast(
7550 this->n_quadrature_points);
7552 if (this->mapping_data->quadrature_points.size() > 0)
7553 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7557 this->jacobian = this_jacobian_data.data();
7558 this->J_value = this_J_value_data.data();
7559 this->jacobian_gradients = this_jacobian_gradients_data.data();
7560 this->jacobian_gradients_non_inverse =
7561 this_jacobian_gradients_non_inverse_data.data();