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.");
2033 const unsigned int dof_no = 0,
2034 const unsigned int quad_no = 0,
2047 const std::pair<unsigned int, unsigned int> & range,
2048 const unsigned int dof_no = 0,
2049 const unsigned int quad_no = 0,
2159 template <
bool level_dof_access>
2181 const unsigned int given_n_q_points_1d);
2201 const bool evaluate_gradients,
2202 const bool evaluate_hessians =
false);
2226 const bool evaluate_values,
2227 const bool evaluate_gradients,
2228 const bool evaluate_hessians =
false);
2243 template <
typename VectorType>
2251 template <
typename VectorType>
2254 const bool evaluate_values,
2255 const bool evaluate_gradients,
2256 const bool evaluate_hessians =
false);
2275 integrate(
const bool integrate_values,
const bool integrate_gradients);
2290 VectorizedArrayType * values_array,
2291 const bool sum_into_values =
false);
2298 const bool integrate_gradients,
2299 VectorizedArrayType *values_array);
2314 template <
typename VectorType>
2317 VectorType & output_vector);
2322 template <
typename VectorType>
2325 const bool integrate_gradients,
2326 VectorType &output_vector);
2410 int n_q_points_1d = fe_degree + 1,
2411 int n_components_ = 1,
2412 typename Number =
double,
2418 VectorizedArrayType>
2421 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2422 "Type of Number and of VectorizedArrayType do not match.");
2561 const unsigned int dof_no = 0,
2562 const unsigned int quad_no = 0,
2577 const std::pair<unsigned int, unsigned int> & range,
2579 const unsigned int dof_no = 0,
2580 const unsigned int quad_no = 0,
2594 reinit(
const unsigned int face_batch_number);
2604 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
2611 const unsigned int given_n_q_points_1d);
2630 evaluate(
const bool evaluate_values,
const bool evaluate_gradients);
2653 const bool evaluate_values,
2654 const bool evaluate_gradients);
2667 template <
typename VectorType>
2675 template <
typename VectorType>
2678 const bool evaluate_values,
2679 const bool evaluate_gradients);
2697 integrate(
const bool integrate_values,
const bool integrate_gradients);
2709 VectorizedArrayType * values_array);
2716 const bool integrate_gradients,
2717 VectorizedArrayType *values_array);
2730 template <
typename VectorType>
2733 VectorType & output_vector);
2738 template <
typename VectorType>
2741 const bool integrate_gradients,
2742 VectorType &output_vector);
2804 namespace MatrixFreeFunctions
2808 template <
int dim,
int degree>
2817 template <
int degree>
2820 static constexpr
unsigned int value = degree + 1;
2836 template <
bool is_face,
2839 typename VectorizedArrayType>
2842 extract_initialization_data(
2844 const unsigned int dof_no,
2845 const unsigned int first_selected_component,
2846 const unsigned int quad_no,
2847 const unsigned int fe_degree,
2848 const unsigned int n_q_points,
2849 const unsigned int active_fe_index_given,
2850 const unsigned int active_quad_index_given,
2851 const unsigned int face_type)
2854 InitializationData init_data;
2858 &internal::MatrixFreeFunctions::
2859 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2867 active_fe_index_given :
2872 active_quad_index_given :
2903 typename VectorizedArrayType>
2908 VectorizedArrayType>::
2911 const unsigned int dof_no,
2912 const unsigned int first_selected_component,
2913 const unsigned int quad_no,
2914 const unsigned int fe_degree,
2915 const unsigned int n_q_points,
2916 const bool is_interior_face,
2917 const unsigned int active_fe_index,
2918 const unsigned int active_quad_index,
2919 const unsigned int face_type)
2921 internal::extract_initialization_data<is_face>(matrix_free,
2923 first_selected_component,
2932 first_selected_component)
2933 , scratch_data_array(matrix_free.acquire_scratch_data())
2934 , matrix_free(&matrix_free)
2936 this->set_data_pointers(scratch_data_array, n_components_);
2938 this->dof_info->start_components.back() == 1 ||
2939 static_cast<int>(n_components_) <=
2941 this->dof_info->start_components
2942 [this->dof_info->component_to_base_index[first_selected_component] +
2944 first_selected_component,
2946 "You tried to construct a vector-valued evaluator with " +
2948 " components. However, "
2949 "the current base element has only " +
2951 this->dof_info->start_components
2952 [this->dof_info->component_to_base_index[first_selected_component] +
2954 first_selected_component) +
2955 " components left when starting from local element index " +
2957 first_selected_component -
2958 this->dof_info->start_components
2959 [this->dof_info->component_to_base_index[first_selected_component]]) +
2960 " (global index " +
std::to_string(first_selected_component) +
")"));
2972 typename VectorizedArrayType>
2977 VectorizedArrayType>::
2983 const unsigned int first_selected_component,
2987 other->mapped_geometry->get_quadrature() == quadrature ?
2988 other->mapped_geometry :
2989 std::make_shared<
internal::MatrixFreeFunctions::
2990 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2995 first_selected_component)
2996 , scratch_data_array(new
AlignedVector<VectorizedArrayType>())
2997 , matrix_free(nullptr)
2999 const unsigned int base_element_number =
3003 first_selected_component >=
3005 ExcMessage(
"The underlying element must at least contain as many "
3006 "components as requested by this class"));
3007 (void)base_element_number;
3012 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
3016 this->set_data_pointers(scratch_data_array, n_components_);
3025 typename VectorizedArrayType>
3030 VectorizedArrayType>::
3035 VectorizedArrayType> &other)
3037 , scratch_data_array(other.matrix_free == nullptr ?
3039 other.matrix_free->acquire_scratch_data())
3040 , matrix_free(other.matrix_free)
3042 if (other.matrix_free ==
nullptr)
3050 this->mapped_geometry =
3051 std::make_shared<internal::MatrixFreeFunctions::
3052 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3056 this->mapping_data = &this->mapped_geometry->get_data_storage();
3060 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3062 this->mapped_geometry->get_data_storage().JxW_values.begin();
3063 this->jacobian_gradients =
3064 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3065 this->jacobian_gradients_non_inverse =
3066 this->mapped_geometry->get_data_storage()
3067 .jacobian_gradients_non_inverse[0]
3070 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3073 this->set_data_pointers(scratch_data_array, n_components_);
3082 typename VectorizedArrayType>
3087 VectorizedArrayType> &
3093 VectorizedArrayType> &other)
3096 if (matrix_free ==
nullptr)
3099 delete scratch_data_array;
3108 matrix_free = other.matrix_free;
3110 if (other.matrix_free ==
nullptr)
3119 this->mapped_geometry =
3120 std::make_shared<internal::MatrixFreeFunctions::
3121 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3126 this->mapping_data = &this->mapped_geometry->get_data_storage();
3128 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3130 this->mapped_geometry->get_data_storage().JxW_values.begin();
3131 this->jacobian_gradients =
3132 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3133 this->jacobian_gradients_non_inverse =
3134 this->mapped_geometry->get_data_storage()
3135 .jacobian_gradients_non_inverse[0]
3138 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3145 this->set_data_pointers(scratch_data_array, n_components_);
3156 typename VectorizedArrayType>
3161 VectorizedArrayType>::~FEEvaluationBase()
3163 if (matrix_free !=
nullptr)
3174 delete scratch_data_array;
3185 typename VectorizedArrayType>
3190 Assert(matrix_free !=
nullptr,
3192 "FEEvaluation was not initialized with a MatrixFree object!"));
3193 return *matrix_free;
3202 template <
typename VectorType,
bool>
3203 struct ConstBlockVectorSelector;
3205 template <
typename VectorType>
3206 struct ConstBlockVectorSelector<VectorType, true>
3208 using BaseVectorType =
const typename VectorType::BlockType;
3211 template <
typename VectorType>
3212 struct ConstBlockVectorSelector<VectorType, false>
3214 using BaseVectorType =
typename VectorType::BlockType;
3220 template <
typename VectorType,
bool>
3221 struct BlockVectorSelector;
3223 template <
typename VectorType>
3224 struct BlockVectorSelector<VectorType, true>
3226 using BaseVectorType =
typename ConstBlockVectorSelector<
3228 std::is_const<VectorType>::value>::BaseVectorType;
3230 static BaseVectorType *
3231 get_vector_component(VectorType &vec,
const unsigned int component)
3234 return &vec.block(component);
3238 template <
typename VectorType>
3239 struct BlockVectorSelector<VectorType, false>
3241 using BaseVectorType = VectorType;
3243 static BaseVectorType *
3244 get_vector_component(VectorType &vec,
const unsigned int component)
3263 template <
typename VectorType>
3264 struct BlockVectorSelector<std::vector<VectorType>, false>
3266 using BaseVectorType = VectorType;
3268 static BaseVectorType *
3269 get_vector_component(std::vector<VectorType> &vec,
3270 const unsigned int component)
3273 return &vec[component];
3277 template <
typename VectorType>
3278 struct BlockVectorSelector<const std::vector<VectorType>, false>
3280 using BaseVectorType =
const VectorType;
3282 static const BaseVectorType *
3283 get_vector_component(
const std::vector<VectorType> &vec,
3284 const unsigned int component)
3287 return &vec[component];
3291 template <
typename VectorType>
3292 struct BlockVectorSelector<std::vector<VectorType *>, false>
3294 using BaseVectorType = VectorType;
3296 static BaseVectorType *
3297 get_vector_component(std::vector<VectorType *> &vec,
3298 const unsigned int component)
3301 return vec[component];
3305 template <
typename VectorType>
3306 struct BlockVectorSelector<const std::vector<VectorType *>, false>
3308 using BaseVectorType =
const VectorType;
3310 static const BaseVectorType *
3311 get_vector_component(
const std::vector<VectorType *> &vec,
3312 const unsigned int component)
3315 return vec[component];
3326 typename VectorizedArrayType>
3327 template <
typename VectorType,
typename VectorOperation>
3332 const std::array<VectorType *, n_components_> &src,
3335 n_components_> & src_sm,
3336 const std::bitset<VectorizedArrayType::size()> &mask,
3337 const bool apply_constraints)
const
3342 if (this->matrix_free ==
nullptr)
3344 read_write_operation_global(operation, src);
3351 if (this->n_fe_components == 1)
3352 for (
unsigned int comp = 0; comp < n_components; ++comp)
3354 Assert(src[comp] !=
nullptr,
3355 ExcMessage(
"The finite element underlying this FEEvaluation "
3356 "object is scalar, but you requested " +
3358 " components via the template argument in "
3359 "FEEvaluation. In that case, you must pass an "
3360 "std::vector<VectorType> or a BlockVector to " +
3361 "read_dof_values and distribute_local_to_global."));
3381 [is_face ? this->dof_access_index :
3384 IndexStorageVariants::contiguous)
3386 read_write_operation_contiguous(operation, src, src_sm, mask);
3393 constexpr
unsigned int n_lanes = VectorizedArrayType::size();
3395 std::array<
unsigned int, VectorizedArrayType::size()> cells =
3396 this->get_cell_ids();
3398 const bool masking_is_active =
mask.count() < n_lanes;
3399 if (masking_is_active)
3400 for (
unsigned int v = 0; v < n_lanes; ++v)
3401 if (mask[v] ==
false)
3404 bool has_hn_constraints =
false;
3406 if (is_face ==
false)
3412 [this->first_selected_component])
3413 for (
unsigned int v = 0; v < n_lanes; ++v)
3418 has_hn_constraints =
true;
3421 std::integral_constant<
bool,
3422 internal::is_vectorizable<VectorType, Number>::value>
3425 const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
3427 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3428 std::array<VectorizedArrayType *, n_components> values_dofs;
3429 for (
unsigned int c = 0; c < n_components; ++c)
3430 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3431 c * dofs_per_component;
3435 [is_face ? this->dof_access_index :
3438 IndexStorageVariants::interleaved &&
3439 use_vectorized_path)
3441 const unsigned int *dof_indices =
3443 dof_info.
row_starts[this->cell * this->n_fe_components * n_lanes]
3447 [this->first_selected_component] *
3449 if (n_components == 1 || this->n_fe_components == 1)
3450 for (
unsigned int i = 0; i < dofs_per_component;
3451 ++i, dof_indices += n_lanes)
3452 for (
unsigned int comp = 0; comp < n_components; ++comp)
3453 operation.process_dof_gather(dof_indices,
3456 values_dofs[comp][i],
3459 for (
unsigned int comp = 0; comp < n_components; ++comp)
3460 for (
unsigned int i = 0; i < dofs_per_component;
3461 ++i, dof_indices += n_lanes)
3462 operation.process_dof_gather(
3463 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3469 std::array<const unsigned int *, n_lanes> dof_indices;
3470 dof_indices.fill(
nullptr);
3475 bool has_constraints =
false;
3476 const unsigned int n_components_read =
3477 this->n_fe_components > 1 ? n_components : 1;
3481 for (
unsigned int v = 0; v < n_lanes; ++v)
3487 const std::pair<unsigned int, unsigned int> *my_index_start =
3488 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3489 this->first_selected_component];
3494 if (my_index_start[n_components_read].
second !=
3495 my_index_start[0].
second)
3496 has_constraints =
true;
3499 dof_info.
dof_indices.data() + my_index_start[0].first;
3504 for (
unsigned int v = 0; v < n_lanes; ++v)
3509 const std::pair<unsigned int, unsigned int> *my_index_start =
3510 &dof_info.
row_starts[cells[v] * this->n_fe_components +
3511 this->first_selected_component];
3512 if (my_index_start[n_components_read].
second !=
3513 my_index_start[0].
second)
3514 has_constraints =
true;
3522 [this->active_fe_index][this->first_selected_component])
3523 has_hn_constraints =
true;
3526 my_index_start[0].
first ||
3529 my_index_start[0].
first,
3532 dof_info.
dof_indices.data() + my_index_start[0].first;
3536 if (std::count_if(cells.begin(), cells.end(), [](
const auto i) {
3537 return i != numbers::invalid_unsigned_int;
3539 for (
unsigned int comp = 0; comp < n_components; ++comp)
3540 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3541 operation.process_empty(values_dofs[comp][i]);
3545 if (!has_constraints && apply_constraints)
3547 if (n_components == 1 || this->n_fe_components == 1)
3549 for (
unsigned int v = 0; v < n_lanes; ++v)
3554 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3555 for (
unsigned int comp = 0; comp < n_components; ++comp)
3556 operation.process_dof(dof_indices[v][i],
3558 values_dofs[comp][i][v]);
3563 for (
unsigned int comp = 0; comp < n_components; ++comp)
3564 for (
unsigned int v = 0; v < n_lanes; ++v)
3569 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3570 operation.process_dof(
3571 dof_indices[v][comp * dofs_per_component + i],
3573 values_dofs[comp][i][v]);
3585 for (
unsigned int v = 0; v < n_lanes; ++v)
3591 const unsigned int cell_dof_index =
3592 cell_index * this->n_fe_components + this->first_selected_component;
3593 const unsigned int n_components_read =
3594 this->n_fe_components > 1 ? n_components : 1;
3595 unsigned int index_indicators =
3597 unsigned int next_index_indicators =
3598 dof_info.
row_starts[cell_dof_index + 1].second;
3602 if (apply_constraints ==
false &&
3603 (dof_info.
row_starts[cell_dof_index].second !=
3604 dof_info.
row_starts[cell_dof_index + n_components_read].second ||
3611 [this->active_fe_index][this->first_selected_component])))
3620 [this->first_selected_component] +
3622 next_index_indicators = index_indicators;
3625 if (n_components == 1 || this->n_fe_components == 1)
3627 unsigned int ind_local = 0;
3628 for (; index_indicators != next_index_indicators; ++index_indicators)
3630 const std::pair<unsigned short, unsigned short> indicator =
3633 for (
unsigned int j = 0; j < indicator.first; ++j)
3634 for (
unsigned int comp = 0; comp < n_components; ++comp)
3635 operation.process_dof(dof_indices[v][j],
3637 values_dofs[comp][ind_local + j][v]);
3639 ind_local += indicator.first;
3640 dof_indices[v] += indicator.first;
3644 Number
value[n_components];
3645 for (
unsigned int comp = 0; comp < n_components; ++comp)
3646 operation.pre_constraints(values_dofs[comp][ind_local][v],
3649 const Number *data_val =
3651 const Number *end_pool =
3653 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3654 for (
unsigned int comp = 0; comp < n_components; ++comp)
3655 operation.process_constraint(*dof_indices[v],
3660 for (
unsigned int comp = 0; comp < n_components; ++comp)
3661 operation.post_constraints(value[comp],
3662 values_dofs[comp][ind_local][v]);
3668 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3669 for (
unsigned int comp = 0; comp < n_components; ++comp)
3670 operation.process_dof(*dof_indices[v],
3672 values_dofs[comp][ind_local][v]);
3681 for (
unsigned int comp = 0; comp < n_components; ++comp)
3683 unsigned int ind_local = 0;
3686 for (; index_indicators != next_index_indicators;
3689 const std::pair<unsigned short, unsigned short> indicator =
3693 for (
unsigned int j = 0; j < indicator.first; ++j)
3694 operation.process_dof(dof_indices[v][j],
3696 values_dofs[comp][ind_local + j][v]);
3697 ind_local += indicator.first;
3698 dof_indices[v] += indicator.first;
3703 operation.pre_constraints(values_dofs[comp][ind_local][v],
3706 const Number *data_val =
3708 const Number *end_pool =
3711 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3712 operation.process_constraint(*dof_indices[v],
3717 operation.post_constraints(value,
3718 values_dofs[comp][ind_local][v]);
3725 for (; ind_local < dofs_per_component;
3726 ++dof_indices[v], ++ind_local)
3729 operation.process_dof(*dof_indices[v],
3731 values_dofs[comp][ind_local][v]);
3734 if (apply_constraints ==
true && comp + 1 < n_components)
3735 next_index_indicators =
3736 dof_info.
row_starts[cell_dof_index + comp + 2].second;
3748 typename VectorizedArrayType>
3749 template <
typename VectorType,
typename VectorOperation>
3754 const std::array<VectorType *, n_components_> &src)
const
3758 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3759 unsigned int index = this->first_selected_component * dofs_per_component;
3760 for (
unsigned int comp = 0; comp < n_components; ++comp)
3762 for (
unsigned int i = 0; i < dofs_per_component; ++i, ++
index)
3764 operation.process_empty(
3765 this->values_dofs[comp * dofs_per_component + i]);
3766 operation.process_dof_global(
3767 local_dof_indices[this->data->lexicographic_numbering[index]],
3769 this->values_dofs[comp * dofs_per_component + i][0]);
3780 typename VectorizedArrayType>
3781 template <
typename VectorType,
typename VectorOperation>
3786 const std::array<VectorType *, n_components_> &src,
3789 n_components_> & vectors_sm,
3790 const std::bitset<VectorizedArrayType::size()> &mask)
const
3799 std::integral_constant<
bool,
3800 internal::is_vectorizable<VectorType, Number>::value>
3803 is_face ? this->dof_access_index :
3805 const unsigned int n_lanes =
mask.count();
3808 const std::vector<unsigned int> & dof_indices_cont =
3811 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3812 std::array<VectorizedArrayType *, n_components> values_dofs;
3813 for (
unsigned int c = 0; c < n_components; ++c)
3814 values_dofs[c] =
const_cast<VectorizedArrayType *
>(this->values_dofs) +
3815 c * dofs_per_component;
3823 interleaved_contiguous &&
3824 n_lanes == VectorizedArrayType::size()) &&
3826 this->dof_access_index ==
3828 this->is_interior_face() ==
false) &&
3829 !(!is_face && !this->is_interior_face()))
3831 const unsigned int dof_index =
3832 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3835 [this->first_selected_component] *
3836 VectorizedArrayType::size();
3837 if (n_components == 1 || this->n_fe_components == 1)
3838 for (
unsigned int comp = 0; comp < n_components; ++comp)
3839 operation.process_dofs_vectorized(dofs_per_component,
3845 operation.process_dofs_vectorized(dofs_per_component * n_components,
3853 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
3854 this->get_cell_or_face_ids();
3858 const unsigned int n_filled_lanes =
3862 (this->dof_access_index ==
3864 this->is_interior_face() ==
false) ||
3865 (!is_face && !this->is_interior_face());
3867 if (vectors_sm[0] !=
nullptr)
3869 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
3870 std::array<
typename VectorType::value_type *,
3871 VectorizedArrayType::size()>
3874 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3876 if (mask[v] ==
false)
3878 vector_ptrs[v] =
nullptr;
3898 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
3899 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3901 [this->active_fe_index][this->first_selected_component]);
3903 vector_ptrs[v] =
nullptr;
3905 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3907 vector_ptrs[v] =
nullptr;
3912 if (n_filled_lanes == VectorizedArrayType::size() &&
3913 n_lanes == VectorizedArrayType::size() && !is_ecl)
3915 if (n_components == 1 || this->n_fe_components == 1)
3917 for (
unsigned int comp = 0; comp < n_components; ++comp)
3919 auto vector_ptrs = compute_vector_ptrs(comp);
3920 operation.process_dofs_vectorized_transpose(
3929 auto vector_ptrs = compute_vector_ptrs(0);
3930 operation.process_dofs_vectorized_transpose(dofs_per_component *
3938 for (
unsigned int comp = 0; comp < n_components; ++comp)
3940 auto vector_ptrs = compute_vector_ptrs(
3941 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3943 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3944 operation.process_empty(values_dofs[comp][i]);
3946 if (n_components == 1 || this->n_fe_components == 1)
3948 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3949 if (mask[v] ==
true)
3950 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3951 operation.process_dof(vector_ptrs[v][i],
3952 values_dofs[comp][i][v]);
3956 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3957 if (mask[v] ==
true)
3958 for (
unsigned int i = 0; i < dofs_per_component; ++i)
3959 operation.process_dof(
3960 vector_ptrs[v][i + comp * dofs_per_component],
3961 values_dofs[comp][i][v]);
3967 unsigned int dof_indices[VectorizedArrayType::size()];
3969 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
3973 if (mask[v] ==
true)
3975 dof_indices_cont[cells[v]] +
3978 [this->first_selected_component] *
3984 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3989 if (n_filled_lanes == VectorizedArrayType::size() &&
3990 n_lanes == VectorizedArrayType::size() && !is_ecl)
3996 if (n_components == 1 || this->n_fe_components == 1)
3997 for (
unsigned int comp = 0; comp < n_components; ++comp)
3998 operation.process_dofs_vectorized_transpose(dofs_per_component,
4004 operation.process_dofs_vectorized_transpose(dofs_per_component *
4013 interleaved_contiguous_strided)
4015 if (n_components == 1 || this->n_fe_components == 1)
4016 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4018 for (
unsigned int comp = 0; comp < n_components; ++comp)
4019 operation.process_dof_gather(dof_indices,
4021 i * VectorizedArrayType::size(),
4022 values_dofs[comp][i],
4026 for (
unsigned int comp = 0; comp < n_components; ++comp)
4027 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4029 operation.process_dof_gather(dof_indices,
4031 (comp * dofs_per_component + i) *
4032 VectorizedArrayType::size(),
4033 values_dofs[comp][i],
4041 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4043 const unsigned int *offsets =
4045 [ind][VectorizedArrayType::size() * this->cell];
4046 if (n_components == 1 || this->n_fe_components == 1)
4047 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4049 for (
unsigned int comp = 0; comp < n_components; ++comp)
4050 operation.process_dof_gather(dof_indices,
4053 values_dofs[comp][i],
4056 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4057 dof_indices[v] += offsets[v];
4060 for (
unsigned int comp = 0; comp < n_components; ++comp)
4061 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4063 operation.process_dof_gather(dof_indices,
4066 values_dofs[comp][i],
4069 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4070 dof_indices[v] += offsets[v];
4075 for (
unsigned int comp = 0; comp < n_components; ++comp)
4077 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4078 operation.process_empty(values_dofs[comp][i]);
4083 if (n_components == 1 || this->n_fe_components == 1)
4085 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4086 if (mask[v] ==
true)
4087 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4088 operation.process_dof(dof_indices[v] + i,
4090 values_dofs[comp][i][v]);
4094 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4095 if (mask[v] ==
true)
4096 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4097 operation.process_dof(dof_indices[v] + i +
4098 comp * dofs_per_component,
4100 values_dofs[comp][i][v]);
4105 const unsigned int *offsets =
4107 [ind][VectorizedArrayType::size() * this->cell];
4108 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4110 if (n_components == 1 || this->n_fe_components == 1)
4111 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4113 if (mask[v] ==
true)
4114 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4115 operation.process_dof(dof_indices[v] + i * offsets[v],
4117 values_dofs[comp][i][v]);
4121 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4122 if (mask[v] ==
true)
4123 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4124 operation.process_dof(dof_indices[v] +
4125 (i + comp * dofs_per_component) *
4128 values_dofs[comp][i][v]);
4138 typename VectorType,
4139 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4140 decltype(std::declval<VectorType>().
begin())
4141 get_beginning(VectorType &vec)
4148 typename VectorType,
4149 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * =
nullptr>
4150 typename VectorType::value_type *
4151 get_beginning(VectorType &)
4156 template <
typename VectorType,
4157 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4159 const std::vector<ArrayView<const typename VectorType::value_type>> *
4160 get_shared_vector_data(VectorType * vec,
4161 const bool is_valid_mode_for_sm,
4162 const unsigned int active_fe_index,
4166 if (is_valid_mode_for_sm &&
4169 active_fe_index == 0)
4170 return &vec->shared_vector_data();
4175 template <
typename VectorType,
4176 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4178 const std::vector<ArrayView<const typename VectorType::value_type>> *
4179 get_shared_vector_data(VectorType *,
4187 template <
int n_components,
typename VectorType>
4189 std::array<
typename internal::BlockVectorSelector<
4194 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
4198 get_vector_data(VectorType & src,
4199 const unsigned int first_index,
4200 const bool is_valid_mode_for_sm,
4201 const unsigned int active_fe_index,
4207 std::array<
typename internal::BlockVectorSelector<
4213 ArrayView<
const typename internal::BlockVectorSelector<
4219 for (
unsigned int d = 0;
d < n_components; ++
d)
4220 src_data.first[
d] = internal::BlockVectorSelector<
4226 for (
unsigned int d = 0;
d < n_components; ++
d)
4227 src_data.second[
d] = get_shared_vector_data(
4228 const_cast<typename internal::BlockVectorSelector<
4229 typename std::remove_const<VectorType>::type,
4231 BaseVectorType *>(src_data.first[
d]),
4232 is_valid_mode_for_sm,
4246 typename VectorizedArrayType>
4251 if (this->dof_info ==
nullptr ||
4253 this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4254 this->dof_info->hanging_node_constraint_masks_comp
4255 [this->active_fe_index][this->first_selected_component] ==
false)
4258 constexpr
unsigned int n_lanes = VectorizedArrayType::size();
4259 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4262 bool hn_available =
false;
4264 const std::array<
unsigned int, VectorizedArrayType::size()> &cells =
4265 this->get_cell_ids();
4267 for (
unsigned int v = 0; v < n_lanes; ++v)
4279 constraint_mask[v] =
mask;
4285 if (hn_available ==
false)
4290 this->data->data.front().fe_degree,
4291 this->get_shape_info(),
4303 typename VectorizedArrayType>
4304 template <
typename VectorType>
4308 const unsigned int first_index,
4309 const std::bitset<VectorizedArrayType::size()> &mask)
4311 const auto src_data = internal::get_vector_data<n_components_>(
4314 this->dof_access_index ==
4316 this->active_fe_index,
4320 read_write_operation(reader, src_data.first, src_data.second, mask,
true);
4322 apply_hanging_node_constraints(
false);
4325 this->dof_values_initialized =
true;
4335 typename VectorizedArrayType>
4336 template <
typename VectorType>
4340 const unsigned int first_index,
4341 const std::bitset<VectorizedArrayType::size()> &mask)
4343 const auto src_data = internal::get_vector_data<n_components_>(
4346 this->dof_access_index ==
4348 this->active_fe_index,
4352 read_write_operation(reader, src_data.first, src_data.second, mask,
false);
4355 this->dof_values_initialized =
true;
4365 typename VectorizedArrayType>
4366 template <
typename VectorType>
4371 const unsigned int first_index,
4372 const std::bitset<VectorizedArrayType::size()> &mask)
const
4375 Assert(this->dof_values_initialized ==
true,
4379 apply_hanging_node_constraints(
true);
4381 const auto dst_data = internal::get_vector_data<n_components_>(
4384 this->dof_access_index ==
4386 this->active_fe_index,
4391 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4400 typename VectorizedArrayType>
4401 template <
typename VectorType>
4405 const unsigned int first_index,
4406 const std::bitset<VectorizedArrayType::size()> &mask)
const
4409 Assert(this->dof_values_initialized ==
true,
4413 const auto dst_data = internal::get_vector_data<n_components_>(
4416 this->dof_access_index ==
4418 this->active_fe_index,
4422 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4431 typename VectorizedArrayType>
4432 template <
typename VectorType>
4437 const unsigned int first_index,
4438 const std::bitset<VectorizedArrayType::size()> &mask)
const
4441 Assert(this->dof_values_initialized ==
true,
4445 const auto dst_data = internal::get_vector_data<n_components_>(
4448 this->dof_access_index ==
4450 this->active_fe_index,
4454 read_write_operation(setter, dst_data.first, dst_data.second, mask,
false);
4467 typename VectorizedArrayType>
4473 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4475 for (
unsigned int comp = 0; comp < n_components; ++comp)
4476 return_value[comp] = this->values_dofs[comp * dofs + dof];
4477 return return_value;
4486 typename VectorizedArrayType>
4489 get_value(
const unsigned int q_point)
const
4492 Assert(this->values_quad_initialized ==
true,
4497 const std::size_t nqp = this->n_quadrature_points;
4499 for (
unsigned int comp = 0; comp < n_components; ++comp)
4500 return_value[comp] = this->values_quad[comp * nqp + q_point];
4501 return return_value;
4510 typename VectorizedArrayType>
4517 Assert(this->gradients_quad_initialized ==
true,
4522 Assert(this->jacobian !=
nullptr,
4524 "update_gradients"));
4525 const std::size_t nqp = this->n_quadrature_points;
4531 for (
unsigned int d = 0;
d < dim; ++
d)
4532 for (
unsigned int comp = 0; comp < n_components; ++comp)
4534 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
4535 this->jacobian[0][
d][
d];
4544 for (
unsigned int comp = 0; comp < n_components; ++comp)
4545 for (
unsigned int d = 0;
d < dim; ++
d)
4548 jac[
d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4549 for (
unsigned int e = 1;
e < dim; ++
e)
4550 grad_out[comp][
d] +=
4552 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4564 typename VectorizedArrayType>
4571 Assert(this->gradients_quad_initialized ==
true,
4575 Assert(this->normal_x_jacobian !=
nullptr,
4577 "update_gradients"));
4579 const std::size_t nqp = this->n_quadrature_points;
4583 for (
unsigned int comp = 0; comp < n_components; ++comp)
4585 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4586 (this->normal_x_jacobian[0][dim - 1]);
4589 const std::size_t
index =
4591 for (
unsigned int comp = 0; comp < n_components; ++comp)
4593 grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4594 this->normal_x_jacobian[
index][0];
4595 for (
unsigned int d = 1;
d < dim; ++
d)
4597 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
4598 this->normal_x_jacobian[
index][
d];
4610 template <
typename VectorizedArrayType>
4613 const VectorizedArrayType *
const hessians,
4615 VectorizedArrayType (&tmp)[1][1])
4617 tmp[0][0] = jac[0][0] *
hessians[0];
4620 template <
typename VectorizedArrayType>
4623 const VectorizedArrayType *
const hessians,
4624 const unsigned int nqp,
4625 VectorizedArrayType (&tmp)[2][2])
4627 for (
unsigned int d = 0;
d < 2; ++
d)
4635 template <
typename VectorizedArrayType>
4638 const VectorizedArrayType *
const hessians,
4639 const unsigned int nqp,
4640 VectorizedArrayType (&tmp)[3][3])
4642 for (
unsigned int d = 0;
d < 3; ++
d)
4663 typename VectorizedArrayType>
4669 Assert(this->hessians_quad_initialized ==
true,
4674 Assert(this->jacobian !=
nullptr,
4684 const std::size_t nqp = this->n_quadrature_points;
4685 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4690 for (
unsigned int comp = 0; comp < n_components; ++comp)
4692 for (
unsigned int d = 0;
d < dim; ++
d)
4693 hessian_out[comp][
d][
d] =
4694 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4695 (jac[
d][
d] * jac[
d][
d]);
4701 hessian_out[comp][0][1] =
4702 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4703 (jac[0][0] * jac[1][1]);
4706 hessian_out[comp][0][1] =
4707 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4708 (jac[0][0] * jac[1][1]);
4709 hessian_out[comp][0][2] =
4710 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4711 (jac[0][0] * jac[2][2]);
4712 hessian_out[comp][1][2] =
4713 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4714 (jac[1][1] * jac[2][2]);
4719 for (
unsigned int d = 0;
d < dim; ++
d)
4720 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4721 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4727 for (
unsigned int comp = 0; comp < n_components; ++comp)
4729 VectorizedArrayType tmp[dim][dim];
4730 internal::hessian_unit_times_jac(
4731 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4734 for (
unsigned int d = 0;
d < dim; ++
d)
4735 for (
unsigned int e =
d;
e < dim; ++
e)
4737 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4738 for (
unsigned int f = 1; f < dim; ++f)
4739 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4746 for (
unsigned int d = 0;
d < dim; ++
d)
4747 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4748 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4754 const auto &jac_grad = this->jacobian_gradients[q_point];
4755 for (
unsigned int comp = 0; comp < n_components; ++comp)
4757 VectorizedArrayType tmp[dim][dim];
4758 internal::hessian_unit_times_jac(
4759 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4762 for (
unsigned int d = 0;
d < dim; ++
d)
4763 for (
unsigned int e =
d;
e < dim; ++
e)
4765 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
4766 for (
unsigned int f = 1; f < dim; ++f)
4767 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
4771 for (
unsigned int d = 0;
d < dim; ++
d)
4772 for (
unsigned int e = 0;
e < dim; ++
e)
4773 hessian_out[comp][
d][
d] +=
4775 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4778 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
4779 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++count)
4780 for (
unsigned int f = 0; f < dim; ++f)
4781 hessian_out[comp][
d][
e] +=
4782 jac_grad[count][f] *
4783 this->gradients_quad[(comp * dim + f) * nqp + q_point];
4786 for (
unsigned int d = 0;
d < dim; ++
d)
4787 for (
unsigned int e =
d + 1;
e < dim; ++
e)
4788 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
4800 typename VectorizedArrayType>
4807 Assert(this->hessians_quad_initialized ==
true,
4818 const std::size_t nqp = this->n_quadrature_points;
4819 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
4825 for (
unsigned int comp = 0; comp < n_components; ++comp)
4826 for (
unsigned int d = 0;
d < dim; ++
d)
4827 hessian_out[comp][
d] =
4828 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] *
4829 (jac[
d][
d] * jac[
d][
d]);
4834 for (
unsigned int comp = 0; comp < n_components; ++comp)
4838 VectorizedArrayType tmp[dim][dim];
4839 internal::hessian_unit_times_jac(
4840 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4844 for (
unsigned int d = 0;
d < dim; ++
d)
4846 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4847 for (
unsigned int f = 1; f < dim; ++f)
4848 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4855 const auto &jac_grad = this->jacobian_gradients[q_point];
4856 for (
unsigned int comp = 0; comp < n_components; ++comp)
4860 VectorizedArrayType tmp[dim][dim];
4861 internal::hessian_unit_times_jac(
4862 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4866 for (
unsigned int d = 0;
d < dim; ++
d)
4868 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
4869 for (
unsigned int f = 1; f < dim; ++f)
4870 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
4873 for (
unsigned int d = 0;
d < dim; ++
d)
4874 for (
unsigned int e = 0;
e < dim; ++
e)
4875 hessian_out[comp][
d] +=
4877 this->gradients_quad[(comp * dim +
e) * nqp + q_point];
4889 typename VectorizedArrayType>
4896 Assert(this->hessians_quad_initialized ==
true,
4902 const auto hess_diag = get_hessian_diagonal(q_point);
4903 for (
unsigned int comp = 0; comp < n_components; ++comp)
4905 laplacian_out[comp] = hess_diag[comp][0];
4906 for (
unsigned int d = 1;
d < dim; ++
d)
4907 laplacian_out[comp] += hess_diag[comp][
d];
4909 return laplacian_out;
4918 typename VectorizedArrayType>
4922 const unsigned int dof)
4925 this->dof_values_initialized =
true;
4927 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4929 for (
unsigned int comp = 0; comp < n_components; ++comp)
4930 this->values_dofs[comp * dofs + dof] = val_in[comp];
4939 typename VectorizedArrayType>
4943 const unsigned int q_point)
4949 Assert(this->J_value !=
nullptr,
4953 this->values_quad_submitted =
true;
4956 const std::size_t nqp = this->n_quadrature_points;
4959 const VectorizedArrayType JxW =
4960 this->J_value[0] * this->quadrature_weights[q_point];
4961 for (
unsigned int comp = 0; comp < n_components; ++comp)
4962 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4966 const VectorizedArrayType JxW = this->J_value[q_point];
4967 for (
unsigned int comp = 0; comp < n_components; ++comp)
4968 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4978 typename VectorizedArrayType>
4983 const unsigned int q_point)
4989 Assert(this->J_value !=
nullptr,
4991 "update_gradients"));
4992 Assert(this->jacobian !=
nullptr,
4994 "update_gradients"));
4996 this->gradients_quad_submitted =
true;
4999 const std::size_t nqp = this->n_quadrature_points;
5002 const VectorizedArrayType JxW =
5003 this->J_value[0] * this->quadrature_weights[q_point];
5004 std::array<VectorizedArrayType, dim> jac;
5005 for (
unsigned int d = 0;
d < dim; ++
d)
5006 jac[
d] = this->jacobian[0][
d][
d];
5007 for (
unsigned int d = 0;
d < dim; ++
d)
5009 const VectorizedArrayType factor = jac[
d] * JxW;
5010 for (
unsigned int comp = 0; comp < n_components; ++comp)
5011 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5012 grad_in[comp][
d] * factor;
5019 this->jacobian[q_point] :
5021 const VectorizedArrayType JxW =
5023 this->J_value[q_point] :
5024 this->J_value[0] * this->quadrature_weights[q_point];
5025 for (
unsigned int comp = 0; comp < n_components; ++comp)
5026 for (
unsigned int d = 0;
d < dim; ++
d)
5028 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
5029 for (
unsigned int e = 1;
e < dim; ++
e)
5030 new_val += (jac[
e][
d] * grad_in[comp][
e]);
5031 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5043 typename VectorizedArrayType>
5048 const unsigned int q_point)
5051 Assert(this->normal_x_jacobian !=
nullptr,
5053 "update_gradients"));
5055 this->gradients_quad_submitted =
true;
5058 const std::size_t nqp = this->n_quadrature_points;
5061 const VectorizedArrayType JxW_jac = this->J_value[0] *
5062 this->quadrature_weights[q_point] *
5063 this->normal_x_jacobian[0][dim - 1];
5064 for (
unsigned int comp = 0; comp < n_components; ++comp)
5066 for (
unsigned int d = 0;
d < dim - 1; ++
d)
5067 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5068 VectorizedArrayType();
5069 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5070 grad_in[comp] * JxW_jac;
5075 const unsigned int index =
5078 this->normal_x_jacobian[
index];
5079 const VectorizedArrayType JxW =
5081 this->J_value[index] * this->quadrature_weights[q_point] :
5082 this->J_value[index];
5083 for (
unsigned int comp = 0; comp < n_components; ++comp)
5085 for (
unsigned int d = 0;
d < dim; ++
d)
5086 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
5087 (grad_in[comp] * JxW) * jac[
d];
5098 typename VectorizedArrayType>
5104 const unsigned int q_point)
5110 Assert(this->J_value !=
nullptr,
5112 "update_hessians"));
5113 Assert(this->jacobian !=
nullptr,
5115 "update_hessians"));
5117 this->hessians_quad_submitted =
true;
5121 const std::size_t nqp = this->n_quadrature_points;
5122 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
5125 const VectorizedArrayType JxW =
5126 this->J_value[0] * this->quadrature_weights[q_point];
5129 for (
unsigned int d = 0;
d < dim; ++
d)
5131 const auto jac_d = this->jacobian[0][
d][
d];
5132 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5133 for (
unsigned int comp = 0; comp < n_components; ++comp)
5134 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5135 hessian_in[comp][
d][
d] * factor;
5139 for (
unsigned int d = 1, off_dia = dim;
d < dim; ++
d)
5140 for (
unsigned int e = 0;
e <
d; ++
e, ++off_dia)
5142 const auto jac_d = this->jacobian[0][
d][
d];
5143 const auto jac_e = this->jacobian[0][
e][
e];
5144 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5145 for (
unsigned int comp = 0; comp < n_components; ++comp)
5146 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5147 (hessian_in[comp][
d][
e] + hessian_in[comp][
e][
d]) * factor;
5154 const VectorizedArrayType JxW =
5155 this->J_value[0] * this->quadrature_weights[q_point];
5156 for (
unsigned int comp = 0; comp < n_components; ++comp)
5159 VectorizedArrayType tmp[dim][dim];
5160 for (
unsigned int i = 0; i < dim; ++i)
5161 for (
unsigned int j = 0; j < dim; ++j)
5163 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5164 for (
unsigned int k = 1; k < dim; ++k)
5165 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5169 VectorizedArrayType tmp2[dim][dim];
5170 for (
unsigned int i = 0; i < dim; ++i)
5171 for (
unsigned int j = 0; j < dim; ++j)
5173 tmp2[i][j] = jac[0][i] * tmp[0][j];
5174 for (
unsigned int k = 1; k < dim; ++k)
5175 tmp2[i][j] += jac[k][i] * tmp[k][j];
5179 for (
unsigned int d = 0;
d < dim; ++
d)
5180 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5184 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5185 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5186 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5187 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5193 const VectorizedArrayType JxW = this->J_value[q_point];
5194 const auto &jac_grad = this->jacobian_gradients[q_point];
5195 for (
unsigned int comp = 0; comp < n_components; ++comp)
5198 VectorizedArrayType tmp[dim][dim];
5199 for (
unsigned int i = 0; i < dim; ++i)
5200 for (
unsigned int j = 0; j < dim; ++j)
5202 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5203 for (
unsigned int k = 1; k < dim; ++k)
5204 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5208 VectorizedArrayType tmp2[dim][dim];
5209 for (
unsigned int i = 0; i < dim; ++i)
5210 for (
unsigned int j = 0; j < dim; ++j)
5212 tmp2[i][j] = jac[0][i] * tmp[0][j];
5213 for (
unsigned int k = 1; k < dim; ++k)
5214 tmp2[i][j] += jac[k][i] * tmp[k][j];
5218 for (
unsigned int d = 0;
d < dim; ++
d)
5219 this->hessians_quad[(comp * hdim +
d) * nqp + q_point] =
5223 for (
unsigned int d = 0, off_diag = dim;
d < dim; ++
d)
5224 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++off_diag)
5225 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5226 (tmp2[
d][
e] + tmp2[
e][
d]) * JxW;
5229 for (
unsigned int d = 0;
d < dim; ++
d)
5231 VectorizedArrayType
sum = 0;
5232 for (
unsigned int e = 0;
e < dim; ++
e)
5233 sum += hessian_in[comp][
e][
e] * jac_grad[
e][
d];
5234 for (
unsigned int e = 0, count = dim;
e < dim; ++
e)
5235 for (
unsigned int f =
e + 1; f < dim; ++f, ++count)
5236 sum += (hessian_in[comp][
e][f] + hessian_in[comp][f][
e]) *
5238 this->gradients_from_hessians_quad[(comp * dim +
d) * nqp +
5239 q_point] =
sum * JxW;
5251 typename VectorizedArrayType>
5258 Assert(this->values_quad_submitted ==
true,
5263 const std::size_t nqp = this->n_quadrature_points;
5264 for (
unsigned int q = 0; q < nqp; ++q)
5265 for (
unsigned int comp = 0; comp < n_components; ++comp)
5266 return_value[comp] += this->values_quad[comp * nqp + q];
5267 return (return_value);
5279 typename VectorizedArrayType>
5284 VectorizedArrayType>::
5287 const unsigned int dof_no,
5288 const unsigned int first_selected_component,
5289 const unsigned int quad_no,
5290 const unsigned int fe_degree,
5291 const unsigned int n_q_points,
5292 const bool is_interior_face,
5293 const unsigned int active_fe_index,
5294 const unsigned int active_quad_index,
5295 const unsigned int face_type)
5296 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5299 first_selected_component,
5315 typename VectorizedArrayType>
5320 VectorizedArrayType>::
5326 const unsigned int first_selected_component,
5328 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5333 first_selected_component,
5343 typename VectorizedArrayType>
5348 VectorizedArrayType>::
5353 VectorizedArrayType> &other)
5354 :
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5364 typename VectorizedArrayType>
5369 VectorizedArrayType> &
5375 VectorizedArrayType> &other)
5381 VectorizedArrayType>::operator=(other);
5390 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5394 const unsigned int dof_no,
5395 const unsigned int first_selected_component,
5396 const unsigned int quad_no,
5397 const unsigned int fe_degree,
5398 const unsigned int n_q_points,
5399 const bool is_interior_face,
5400 const unsigned int active_fe_index,
5401 const unsigned int active_quad_index,
5402 const unsigned int face_type)
5406 first_selected_component,
5418 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5425 const unsigned int first_selected_component,
5432 first_selected_component,
5438 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5448 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5454 ->
FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5461 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5464 const unsigned int dof)
const
5467 return this->values_dofs[dof];
5472 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5475 const unsigned int q_point)
const
5478 Assert(this->values_quad_initialized ==
true,
5482 return this->values_quad[q_point];
5487 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5492 return BaseClass::get_normal_derivative(q_point)[0];
5497 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5500 const unsigned int q_point)
const
5506 Assert(this->gradients_quad_initialized ==
true,
5511 Assert(this->jacobian !=
nullptr,
5513 "update_gradients"));
5517 const std::size_t nqp = this->n_quadrature_points;
5520 for (
unsigned int d = 0;
d < dim; ++
d)
5522 this->gradients_quad[
d * nqp + q_point] * this->jacobian[0][
d][
d];
5531 for (
unsigned int d = 0;
d < dim; ++
d)
5533 grad_out[
d] = jac[
d][0] * this->gradients_quad[q_point];
5534 for (
unsigned int e = 1;
e < dim; ++
e)
5535 grad_out[
d] += jac[
d][
e] * this->gradients_quad[
e * nqp + q_point];
5543 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5546 const unsigned int q_point)
const
5548 return BaseClass::get_hessian(q_point)[0];
5553 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5558 return BaseClass::get_hessian_diagonal(q_point)[0];
5563 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5564 inline VectorizedArrayType
5566 const unsigned int q_point)
const
5568 return BaseClass::get_laplacian(q_point)[0];
5573 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5579 this->dof_values_initialized =
true;
5582 this->values_dofs[dof] = val_in;
5587 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5590 const VectorizedArrayType val_in,
5591 const unsigned int q_point)
5597 Assert(this->J_value !=
nullptr,
5601 this->values_quad_submitted =
true;
5606 const VectorizedArrayType JxW =
5607 this->J_value[0] * this->quadrature_weights[q_point];
5608 this->values_quad[q_point] = val_in * JxW;
5612 this->values_quad[q_point] = val_in * this->J_value[q_point];
5618 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5622 const unsigned int q_point)
5624 submit_value(val_in[0], q_point);
5629 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5633 const unsigned int q_point)
5637 BaseClass::submit_normal_derivative(grad, q_point);
5642 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5646 const unsigned int q_point)
5652 Assert(this->J_value !=
nullptr,
5654 "update_gradients"));
5655 Assert(this->jacobian !=
nullptr,
5657 "update_gradients"));
5659 this->gradients_quad_submitted =
true;
5662 const std::size_t nqp = this->n_quadrature_points;
5665 const VectorizedArrayType JxW =
5666 this->J_value[0] * this->quadrature_weights[q_point];
5670 std::array<VectorizedArrayType, dim> jac;
5671 for (
unsigned int d = 0;
d < dim; ++
d)
5672 jac[
d] = this->jacobian[0][
d][
d];
5674 for (
unsigned int d = 0;
d < dim; ++
d)
5675 this->gradients_quad[
d * nqp + q_point] = grad_in[
d] * jac[
d] * JxW;
5682 this->jacobian[q_point] :
5684 const VectorizedArrayType JxW =
5686 this->J_value[q_point] :
5687 this->J_value[0] * this->quadrature_weights[q_point];
5688 for (
unsigned int d = 0;
d < dim; ++
d)
5690 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
5691 for (
unsigned int e = 1;
e < dim; ++
e)
5692 new_val += jac[
e][
d] * grad_in[
e];
5693 this->gradients_quad[
d * nqp + q_point] = new_val * JxW;
5700 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5704 const unsigned int q_point)
5708 BaseClass::submit_hessian(hessian, q_point);
5713 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5714 inline VectorizedArrayType
5718 return BaseClass::integrate_value()[0];
5726 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5730 const unsigned int dof_no,
5731 const unsigned int first_selected_component,
5732 const unsigned int quad_no,
5733 const unsigned int fe_degree,
5734 const unsigned int n_q_points,
5735 const bool is_interior_face,
5736 const unsigned int active_fe_index,
5737 const unsigned int active_quad_index,
5738 const unsigned int face_type)
5742 first_selected_component,
5754 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5761 const unsigned int first_selected_component,
5768 first_selected_component,
5774 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5784 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5796 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5799 const unsigned int q_point)
const
5801 if (this->data->element_type ==
5806 Assert(this->values_quad_initialized ==
true,
5811 Assert(this->J_value !=
nullptr,
5814 const std::size_t nqp = this->n_quadrature_points;
5822 const VectorizedArrayType inv_det =
5823 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5824 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5825 this->jacobian[0][2][2];
5828 for (
unsigned int comp = 0; comp < n_components; ++comp)
5829 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5830 jac[comp][comp] * inv_det;
5837 this->jacobian[q_point] :
5845 const VectorizedArrayType inv_det =
5846 (is_face && dim == 2 && this->get_face_no() < 2) ?
5850 for (
unsigned int comp = 0; comp < n_components; ++comp)
5853 this->values_quad[q_point] * jac[comp][0] * inv_det;
5854 for (
unsigned int e = 1;
e < dim; ++
e)
5856 this->values_quad[
e * nqp + q_point] * jac[comp][
e] * inv_det;
5864 return BaseClass::get_value(q_point);
5868 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
5873 if (this->data->element_type ==
5878 Assert(this->gradients_quad_initialized ==
true,
5883 Assert(this->jacobian !=
nullptr,
5885 "update_gradients"));
5886 const std::size_t nqp = this->n_quadrature_points;
5896 const VectorizedArrayType inv_det =
5897 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5898 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5899 this->jacobian[0][2][2];
5902 for (
unsigned int d = 0;
d < dim; ++
d)
5903 for (
unsigned int comp = 0; comp < n_components; ++comp)
5905 this->gradients_quad[(comp * dim +
d) * nqp + q_point] *
5906 inv_t_jac[
d][
d] * jac[comp][comp] * inv_det;
5916 const VectorizedArrayType inv_det =
5917 (is_face && dim == 2 && this->get_face_no() < 2) ?
5921 VectorizedArrayType tmp;
5923 for (
unsigned int comp = 0; comp < n_components; ++comp)
5924 for (
unsigned int d = 0;
d < dim; ++
d)
5927 for (
unsigned int f = 0; f < dim; ++f)
5928 for (
unsigned int e = 0;
e < dim; ++
e)
5929 tmp += jac[comp][f] * inv_t_jac[
d][
e] * inv_det *
5930 this->gradients_quad[(f * dim +
e) * nqp + q_point];
5932 grad_out[comp][
d] = tmp;
5943 Assert(this->jacobian_gradients_non_inverse !=
nullptr,
5945 "update_hessians"));
5947 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
5949 this->jacobian[q_point];
5953 const VectorizedArrayType inv_det =
5954 (is_face && dim == 2 && this->get_face_no() < 2) ?
5958 VectorizedArrayType tmp;
5960 for (
unsigned int comp = 0; comp < n_components; ++comp)
5961 for (
unsigned int d = 0;
d < dim; ++
d)
5964 for (
unsigned int f = 0; f < dim; ++f)
5965 for (
unsigned int e = 0;
e < dim; ++
e)
5966 tmp += t_jac[f][comp] * inv_t_jac[
d][
e] *
5967 this->gradients_quad[(f * dim +
e) * nqp + q_point];
5969 grad_out[comp][
d] = tmp * inv_det;
5978 for (
unsigned int i = 0; i < dim; ++i)
5979 for (
unsigned int j = 0; j < dim; ++j)
5981 tmp = jac_grad[0][i] * inv_t_jac[j][0] *
5982 this->values_quad[q_point];
5983 for (
unsigned int f = 1; f < dim; ++f)
5984 tmp += jac_grad[f][i] * inv_t_jac[j][f] *
5985 this->values_quad[f * nqp + q_point];
5987 grad_out[i][j] += tmp * inv_det;
5990 for (
unsigned int i = 0; i < dim; ++i)
5991 for (
unsigned int j = 0; j < dim; ++j)
5994 for (
unsigned int f = 0; f < dim; ++f)
5995 for (
unsigned int n = 0; n < dim; ++n)
5996 for (
unsigned int m = 0; m < dim; ++m)
5997 tmp += inv_t_jac[m][f] * jac_grad[f][m] *
5998 inv_t_jac[j][f] * t_jac[n][i] *
5999 this->values_quad[n * nqp + q_point];
6000 grad_out[i][j] -= tmp * inv_det;
6009 for (
unsigned int i = 0; i < dim; ++i)
6010 for (
unsigned int j = 0; j < dim; ++j)
6013 for (
unsigned int r = 0, f = dim; r < dim; ++r)
6014 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6016 tmp += jac_grad[f][i] *
6018 this->values_quad[r * nqp + q_point] +
6020 this->values_quad[k * nqp + q_point]);
6021 for (
unsigned int n = 0; n < dim; ++n)
6022 for (
unsigned int m = 0; m < dim; ++m)
6023 tmp -= jac_grad[f][m] * t_jac[n][i] *
6024 this->values_quad[n * nqp + q_point] *
6025 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6026 inv_t_jac[m][r] * inv_t_jac[j][k]);
6028 grad_out[i][j] += tmp * inv_det;
6036 return BaseClass::get_gradient(q_point);
6042 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6048 Assert(this->gradients_quad_initialized ==
true,
6052 Assert(this->jacobian !=
nullptr,
6054 "update_gradients"));
6056 VectorizedArrayType divergence;
6057 const std::size_t nqp = this->n_quadrature_points;
6059 if (this->data->element_type ==
6066 const VectorizedArrayType inv_det =
6067 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
6068 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
6069 this->jacobian[0][2][2];
6072 divergence = this->gradients_quad[q_point] * inv_det;
6073 for (
unsigned int d = 1;
d < dim; ++
d)
6075 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
6081 const VectorizedArrayType inv_det =
6083 this->jacobian[this->cell_type >
6087 Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6090 divergence = this->gradients_quad[q_point] * inv_det;
6091 for (
unsigned int d = 1;
d < dim; ++
d)
6093 this->gradients_quad[(dim *
d +
d) * nqp + q_point] * inv_det;
6102 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6103 for (
unsigned int d = 1;
d < dim; ++
d)
6104 divergence += this->gradients_quad[(dim *
d +
d) * nqp + q_point] *
6105 this->jacobian[0][
d][
d];
6112 this->jacobian[q_point] :
6114 divergence = jac[0][0] * this->gradients_quad[q_point];
6115 for (
unsigned int e = 1;
e < dim; ++
e)
6116 divergence += jac[0][
e] * this->gradients_quad[
e * nqp + q_point];
6117 for (
unsigned int d = 1;
d < dim; ++
d)
6118 for (
unsigned int e = 0;
e < dim; ++
e)
6120 jac[
d][
e] * this->gradients_quad[(
d * dim +
e) * nqp + q_point];
6128 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6134 const auto grad = get_gradient(q_point);
6135 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6136 VectorizedArrayType half = Number(0.5);
6137 for (
unsigned int d = 0;
d < dim; ++
d)
6138 symmetrized[
d] = grad[
d][
d];
6144 symmetrized[2] = grad[0][1] + grad[1][0];
6145 symmetrized[2] *= half;
6148 symmetrized[3] = grad[0][1] + grad[1][0];
6149 symmetrized[3] *= half;
6150 symmetrized[4] = grad[0][2] + grad[2][0];
6151 symmetrized[4] *= half;
6152 symmetrized[5] = grad[1][2] + grad[2][1];
6153 symmetrized[5] *= half;
6163 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6165 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6167 const unsigned int q_point)
const
6171 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6177 "Computing the curl in 1d is not a useful operation"));
6180 curl[0] = grad[1][0] - grad[0][1];
6183 curl[0] = grad[2][1] - grad[1][2];
6184 curl[1] = grad[0][2] - grad[2][0];
6185 curl[2] = grad[1][0] - grad[0][1];
6195 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6200 return BaseClass::get_hessian_diagonal(q_point);
6205 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6208 const unsigned int q_point)
const
6211 Assert(this->hessians_quad_initialized ==
true,
6215 return BaseClass::get_hessian(q_point);
6219 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6223 const unsigned int q_point)
6225 if (this->data->element_type ==
6230 Assert(this->J_value !=
nullptr,
6235 this->values_quad_submitted =
true;
6238 const std::size_t nqp = this->n_quadrature_points;
6243 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6245 for (
unsigned int comp = 0; comp < n_components; ++comp)
6246 this->values_quad[comp * nqp + q_point] =
6247 val_in[comp] * weight * jac[comp][comp];
6254 this->jacobian[q_point] :
6263 const VectorizedArrayType fac =
6265 this->quadrature_weights[q_point] :
6267 this->J_value[q_point] :
6268 this->J_value[0] * this->quadrature_weights[q_point]) *
6269 ((dim == 2 && this->get_face_no() < 2) ?
6274 for (
unsigned int comp = 0; comp < n_components; ++comp)
6276 this->values_quad[comp * nqp + q_point] =
6277 val_in[0] * jac[0][comp] * fac;
6278 for (
unsigned int e = 1;
e < dim; ++
e)
6279 this->values_quad[comp * nqp + q_point] +=
6280 val_in[
e] * jac[
e][comp] * fac;
6287 BaseClass::submit_value(val_in, q_point);
6291 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6295 const unsigned int q_point)
6297 if (this->data->element_type ==
6306 Assert(this->J_value !=
nullptr,
6308 "update_gradients"));
6309 Assert(this->jacobian !=
nullptr,
6311 "update_gradients"));
6313 this->gradients_quad_submitted =
true;
6316 const std::size_t nqp = this->n_quadrature_points;
6324 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6325 for (
unsigned int d = 0;
d < dim; ++
d)
6326 for (
unsigned int comp = 0; comp < n_components; ++comp)
6327 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6328 grad_in[comp][
d] * inv_t_jac[
d][
d] * jac[comp][comp] * weight;
6339 const VectorizedArrayType fac =
6340 (!is_face) ? this->quadrature_weights[q_point] :
6341 this->J_value[0] * this->quadrature_weights[q_point] *
6342 ((dim == 2 && this->get_face_no() < 2) ?
6347 for (
unsigned int comp = 0; comp < n_components; ++comp)
6348 for (
unsigned int d = 0;
d < dim; ++
d)
6350 VectorizedArrayType tmp = 0;
6351 for (
unsigned int f = 0; f < dim; ++f)
6352 for (
unsigned int e = 0;
e < dim; ++
e)
6353 tmp += jac[f][comp] * inv_t_jac[
e][
d] * grad_in[f][
e];
6355 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6363 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
6365 this->jacobian[q_point];
6370 const VectorizedArrayType fac =
6372 this->quadrature_weights[q_point] :
6373 this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6377 VectorizedArrayType tmp;
6379 for (
unsigned int comp = 0; comp < n_components; ++comp)
6380 for (
unsigned int d = 0;
d < dim; ++
d)
6383 for (
unsigned int f = 0; f < dim; ++f)
6384 for (
unsigned int e = 0;
e < dim; ++
e)
6385 tmp += t_jac[comp][f] * inv_t_jac[
e][
d] * grad_in[f][
e];
6387 this->gradients_quad[(comp * dim +
d) * nqp + q_point] =
6397 for (
unsigned int f = 0; f < dim; ++f)
6400 for (
unsigned int i = 0; i < dim; ++i)
6401 for (
unsigned int j = 0; j < dim; ++j)
6403 tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
6404 for (
unsigned int m = 0; m < dim; ++m)
6405 for (
unsigned int k = 0; k < dim; ++k)
6406 tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
6407 inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
6409 this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
6417 for (
unsigned int r = 0, f = dim; r < dim; ++r)
6418 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6420 tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
6421 for (
unsigned int j = 1; j < dim; ++j)
6422 tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
6423 for (
unsigned int i = 1; i < dim; ++i)
6424 for (
unsigned int j = 0; j < dim; ++j)
6425 tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
6426 this->values_from_gradients_quad[r * nqp + q_point] +=
6429 tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
6430 for (
unsigned int j = 1; j < dim; ++j)
6431 tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
6432 for (
unsigned int i = 1; i < dim; ++i)
6433 for (
unsigned int j = 0; j < dim; ++j)
6434 tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
6435 this->values_from_gradients_quad[k * nqp + q_point] +=
6440 for (
unsigned int n = 0; n < dim; ++n)
6443 for (
unsigned int r = 0, f = dim; r < dim; ++r)
6444 for (
unsigned int k = r + 1; k < dim; ++k, ++f)
6445 for (
unsigned int i = 0; i < dim; ++i)
6446 for (
unsigned int j = 0; j < dim; ++j)
6447 for (
unsigned int m = 0; m < dim; ++m)
6448 tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
6449 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6450 inv_t_jac[m][r] * inv_t_jac[j][k]);
6452 this->values_from_gradients_quad[n * nqp + q_point] -=
6460 BaseClass::submit_gradient(grad_in, q_point);
6466 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6471 const unsigned int q_point)
6473 if (this->data->element_type ==
6483 BaseClass::submit_gradient(grad_in, q_point);
6489 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6493 const unsigned int q_point)
6499 Assert(this->J_value !=
nullptr,
6501 "update_gradients"));
6502 Assert(this->jacobian !=
nullptr,
6504 "update_gradients"));
6506 this->gradients_quad_submitted =
true;
6509 const std::size_t nqp = this->n_quadrature_points;
6510 if (this->data->element_type ==
6517 const VectorizedArrayType fac =
6519 this->quadrature_weights[q_point] * div_in :
6521 this->J_value[q_point] :
6522 this->J_value[0] * this->quadrature_weights[q_point]) *
6525 this->jacobian[this->cell_type >
6529 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6531 for (
unsigned int d = 0;
d < dim; ++
d)
6533 this->gradients_quad[(dim *
d +
d) * nqp + q_point] = fac;
6534 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6536 this->gradients_quad[(dim *
d +
e) * nqp + q_point] =
6537 VectorizedArrayType();
6538 this->gradients_quad[(dim *
e +
d) * nqp + q_point] =
6539 VectorizedArrayType();
6542 this->divergence_is_requested =
true;
6549 const VectorizedArrayType fac =
6550 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6551 for (
unsigned int d = 0;
d < dim; ++
d)
6553 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6554 (fac * this->jacobian[0][
d][
d]);
6555 for (
unsigned int e =
d + 1;
e < dim; ++
e)
6557 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6558 VectorizedArrayType();
6559 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6560 VectorizedArrayType();
6568 this->jacobian[q_point] :
6570 const VectorizedArrayType fac =
6572 this->J_value[q_point] :
6573 this->J_value[0] * this->quadrature_weights[q_point]) *
6575 for (
unsigned int d = 0;
d < dim; ++
d)
6577 for (
unsigned int e = 0;
e < dim; ++
e)
6578 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6587 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6592 const unsigned int q_point)
6595 this->data->element_type !=
6606 Assert(this->J_value !=
nullptr,
6608 "update_gradients"));
6609 Assert(this->jacobian !=
nullptr,
6611 "update_gradients"));
6613 this->gradients_quad_submitted =
true;
6616 const std::size_t nqp = this->n_quadrature_points;
6619 const VectorizedArrayType JxW =
6620 this->J_value[0] * this->quadrature_weights[q_point];
6621 for (
unsigned int d = 0;
d < dim; ++
d)
6622 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6624 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6625 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6627 const VectorizedArrayType
value =
6629 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6630 value * this->jacobian[0][
d][
d];
6631 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6632 value * this->jacobian[0][
e][
e];
6638 const VectorizedArrayType JxW =
6640 this->J_value[q_point] :
6641 this->J_value[0] * this->quadrature_weights[q_point];
6644 this->jacobian[q_point] :
6646 VectorizedArrayType weighted[dim][dim];
6647 for (
unsigned int i = 0; i < dim; ++i)
6649 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6650 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6652 const VectorizedArrayType
value =
6654 weighted[i][j] =
value;
6655 weighted[j][i] =
value;
6657 for (
unsigned int comp = 0; comp < dim; ++comp)
6658 for (
unsigned int d = 0;
d < dim; ++
d)
6660 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6661 for (
unsigned int e = 1;
e < dim; ++
e)
6662 new_val += jac[
e][
d] * weighted[comp][
e];
6663 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
6670 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6674 const unsigned int q_point)
6682 "Testing by the curl in 1d is not a useful operation"));
6685 grad[1][0] = curl[0];
6686 grad[0][1] = -curl[0];
6689 grad[2][1] = curl[0];
6690 grad[1][2] = -curl[0];
6691 grad[0][2] = curl[1];
6692 grad[2][0] = -curl[1];
6693 grad[1][0] = curl[2];
6694 grad[0][1] = -curl[2];
6699 submit_gradient(grad, q_point);
6706 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6710 const unsigned int dof_no,
6711 const unsigned int first_selected_component,
6712 const unsigned int quad_no,
6713 const unsigned int fe_degree,
6714 const unsigned int n_q_points,
6715 const bool is_interior_face,
6716 const unsigned int active_fe_index,
6717 const unsigned int active_quad_index,
6718 const unsigned int face_type)
6722 first_selected_component,
6734 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6741 const unsigned int first_selected_component,
6748 first_selected_component,
6754 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6763 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6775 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6778 const unsigned int dof)
const
6781 return this->values_dofs[dof];
6786 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6789 const unsigned int q_point)
const
6792 Assert(this->values_quad_initialized ==
true,
6796 return this->values_quad[q_point];
6801 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6804 const unsigned int q_point)
const
6810 Assert(this->gradients_quad_initialized ==
true,
6817 this->jacobian[q_point] :
6821 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6828 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6831 const unsigned int q_point)
const
6833 return get_gradient(q_point)[0];
6838 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6843 return BaseClass::get_normal_derivative(q_point)[0];
6848 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6851 const unsigned int q_point)
const
6853 return BaseClass::get_hessian(q_point)[0];
6858 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6863 return BaseClass::get_hessian_diagonal(q_point)[0];
6868 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6871 const unsigned int q_point)
const
6873 return BaseClass::get_laplacian(q_point)[0];
6878 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6884 this->dof_values_initialized =
true;
6887 this->values_dofs[dof] = val_in;
6892 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6895 const VectorizedArrayType val_in,
6896 const unsigned int q_point)
6903 this->values_quad_submitted =
true;
6908 const VectorizedArrayType JxW = this->J_value[q_point];
6909 this->values_quad[q_point] = val_in * JxW;
6913 const VectorizedArrayType JxW =
6914 this->J_value[0] * this->quadrature_weights[q_point];
6915 this->values_quad[q_point] = val_in * JxW;
6921 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6925 const unsigned int q_point)
6927 submit_value(val_in[0], q_point);
6932 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6936 const unsigned int q_point)
6938 submit_gradient(grad_in[0], q_point);
6943 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6946 const VectorizedArrayType grad_in,
6947 const unsigned int q_point)
6954 this->gradients_quad_submitted =
true;
6959 this->jacobian[q_point] :
6961 const VectorizedArrayType JxW =
6963 this->J_value[q_point] :
6964 this->J_value[0] * this->quadrature_weights[q_point];
6966 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6971 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6975 const unsigned int q_point)
6977 submit_gradient(grad_in[0][0], q_point);
6982 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6986 const unsigned int q_point)
6990 BaseClass::submit_normal_derivative(grad, q_point);
6995 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
6999 const unsigned int q_point)
7001 BaseClass::submit_normal_derivative(grad_in, q_point);
7005 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7009 const unsigned int q_point)
7013 BaseClass::submit_hessian(hessian, q_point);
7017 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7018 inline VectorizedArrayType
7022 return BaseClass::integrate_value()[0];
7035 typename VectorizedArrayType>
7041 VectorizedArrayType>::
7043 const unsigned int fe_no,
7044 const unsigned int quad_no,
7045 const unsigned int first_selected_component,
7046 const unsigned int active_fe_index,
7047 const unsigned int active_quad_index)
7048 : BaseClass(matrix_free,
7050 first_selected_component,
7057 , dofs_per_component(this->data->dofs_per_component_on_cell)
7058 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7059 , n_q_points(this->data->n_q_points)
7061 check_template_arguments(fe_no, 0);
7071 typename VectorizedArrayType>
7077 VectorizedArrayType>::
7079 const std::pair<unsigned int, unsigned int> & range,
7080 const unsigned int dof_no,
7081 const unsigned int quad_no,
7082 const unsigned int first_selected_component)
7086 first_selected_component,
7087 matrix_free.get_cell_active_fe_index(range))
7097 typename VectorizedArrayType>
7103 VectorizedArrayType>::
7108 const unsigned int first_selected_component)
7109 : BaseClass(mapping,
7113 first_selected_component,
7115 , dofs_per_component(this->data->dofs_per_component_on_cell)
7116 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7117 , n_q_points(this->data->n_q_points)
7129 typename VectorizedArrayType>
7135 VectorizedArrayType>::
7139 const unsigned int first_selected_component)
7144 first_selected_component,
7146 , dofs_per_component(this->data->dofs_per_component_on_cell)
7147 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7148 , n_q_points(this->data->n_q_points)
7160 typename VectorizedArrayType>
7166 VectorizedArrayType>::
7169 const unsigned int first_selected_component)
7170 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7172 other.mapped_geometry->get_quadrature(),
7173 other.mapped_geometry->get_fe_values().get_update_flags(),
7174 first_selected_component,
7176 , dofs_per_component(this->data->dofs_per_component_on_cell)
7177 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7178 , n_q_points(this->data->n_q_points)
7190 typename VectorizedArrayType>
7199 , dofs_per_component(this->data->dofs_per_component_on_cell)
7200 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7201 , n_q_points(this->data->n_q_points)
7213 typename VectorizedArrayType>
7219 VectorizedArrayType> &
7225 VectorizedArrayType>::operator=(
const FEEvaluation &other)
7227 BaseClass::operator=(other);
7239 typename VectorizedArrayType>
7246 VectorizedArrayType>::
7247 check_template_arguments(
const unsigned int dof_no,
7248 const unsigned int first_selected_component)
7251 (void)first_selected_component;
7254 this->data->dofs_per_component_on_cell > 0,
7256 "There is nothing useful you can do with an FEEvaluation object with "
7257 "FE_Nothing, i.e., without DoFs! If you have passed to "
7258 "MatrixFree::reinit() a collection of finite elements also containing "
7259 "FE_Nothing, please check - before creating FEEvaluation - the category "
7260 "of the current range by calling either "
7261 "MatrixFree::get_cell_range_category(range) or "
7262 "MatrixFree::get_face_range_category(range). The returned category "
7263 "is the index of the active FE, which you can use to exclude "
7270 static_cast<unsigned int>(fe_degree) !=
7271 this->data->data.front().fe_degree) ||
7272 n_q_points != this->n_quadrature_points)
7274 std::string message =
7275 "-------------------------------------------------------\n";
7276 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7277 message +=
" Called --> FEEvaluation<dim,";
7281 message +=
",Number>(data";
7297 if (
static_cast<unsigned int>(fe_degree) ==
7298 this->data->data.front().fe_degree)
7300 proposed_dof_comp = dof_no;
7301 proposed_fe_comp = first_selected_component;
7304 for (
unsigned int no = 0; no < this->matrix_free->
n_components();
7306 for (
unsigned int nf = 0;
7309 if (this->matrix_free
7310 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7312 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7314 proposed_dof_comp = no;
7315 proposed_fe_comp = nf;
7319 this->mapping_data->descriptor[this->active_quad_index]
7321 proposed_quad_comp = this->quad_no;
7323 for (
unsigned int no = 0;
7328 .descriptor[this->active_quad_index]
7329 .n_q_points == n_q_points)
7331 proposed_quad_comp = no;
7338 if (proposed_dof_comp != first_selected_component)
7339 message +=
"Wrong vector component selection:\n";
7341 message +=
"Wrong quadrature formula selection:\n";
7342 message +=
" Did you mean FEEvaluation<dim,";
7346 message +=
",Number>(data";
7355 std::string correct_pos;
7356 if (proposed_dof_comp != dof_no)
7357 correct_pos =
" ^ ";
7360 if (proposed_quad_comp != this->quad_no)
7361 correct_pos +=
" ^ ";
7364 if (proposed_fe_comp != first_selected_component)
7365 correct_pos +=
" ^\n";
7367 correct_pos +=
" \n";
7373 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7374 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7375 message +=
"Wrong template arguments:\n";
7376 message +=
" Did you mean FEEvaluation<dim,";
7381 message +=
",Number>(data";
7389 std::string correct_pos;
7390 if (this->data->data.front().fe_degree !=
7391 static_cast<unsigned int>(fe_degree))
7395 if (proposed_n_q_points_1d != n_q_points_1d)
7396 correct_pos +=
" ^\n";
7398 correct_pos +=
" \n";
7399 message +=
" " + correct_pos;
7401 Assert(
static_cast<unsigned int>(fe_degree) ==
7402 this->data->data.front().fe_degree &&
7403 n_q_points == this->n_quadrature_points,
7409 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7420 typename VectorizedArrayType>
7429 Assert(this->matrix_free !=
nullptr,
7430 ExcMessage(
"FEEvaluation was initialized without a matrix-free object."
7431 " Integer indexing is not possible."));
7439 const unsigned int offsets =
7440 this->mapping_data->data_index_offsets[
cell_index];
7441 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7442 this->J_value = &this->mapping_data->JxW_values[offsets];
7443 if (!this->mapping_data->jacobian_gradients[0].empty())
7445 this->jacobian_gradients =
7446 this->mapping_data->jacobian_gradients[0].data() + offsets;
7447 this->jacobian_gradients_non_inverse =
7448 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
7452 VectorizedArrayType::size())
7455 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
7456 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7463 this->cell_ids[i] =
cell_index * VectorizedArrayType::size() + i;
7464 for (; i < VectorizedArrayType::size(); ++i)
7468 if (this->mapping_data->quadrature_points.empty() ==
false)
7470 &this->mapping_data->quadrature_points
7471 [this->mapping_data->quadrature_point_offsets[this->cell]];
7474 this->is_reinitialized =
true;
7475 this->dof_values_initialized =
false;
7476 this->values_quad_initialized =
false;
7477 this->gradients_quad_initialized =
false;
7478 this->hessians_quad_initialized =
false;
7489 typename VectorizedArrayType>
7496 VectorizedArrayType>
::
7497 reinit(
const std::array<
unsigned int, VectorizedArrayType::size()> &cell_ids)
7503 this->cell_ids = cell_ids;
7508 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7522 if (this->mapped_geometry ==
nullptr)
7523 this->mapped_geometry =
7524 std::make_shared<internal::MatrixFreeFunctions::
7525 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7527 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7529 auto &this_jacobian_data = mapping_storage.jacobians[0];
7530 auto &this_J_value_data = mapping_storage.JxW_values;
7531 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7532 auto &this_jacobian_gradients_non_inverse_data =
7533 mapping_storage.jacobian_gradients_non_inverse[0];
7534 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7538 if (this->mapping_data->jacobians[0].size() > 0)
7539 this_jacobian_data.resize_fast(2);
7541 if (this->mapping_data->JxW_values.size() > 0)
7542 this_J_value_data.resize_fast(1);
7544 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7545 this_jacobian_gradients_data.resize_fast(1);
7547 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7548 this_jacobian_gradients_non_inverse_data.resize_fast(1);
7550 if (this->mapping_data->quadrature_points.size() > 0)
7551 this_quadrature_points_data.resize_fast(1);
7555 if (this->mapping_data->jacobians[0].size() > 0)
7556 this_jacobian_data.resize_fast(this->n_quadrature_points);
7558 if (this->mapping_data->JxW_values.size() > 0)
7559 this_J_value_data.resize_fast(this->n_quadrature_points);
7561 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7562 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7564 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7565 this_jacobian_gradients_non_inverse_data.resize_fast(
7566 this->n_quadrature_points);
7568 if (this->mapping_data->quadrature_points.size() > 0)
7569 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7573 this->jacobian = this_jacobian_data.data();
7574 this->J_value = this_J_value_data.data();
7575 this->jacobian_gradients = this_jacobian_gradients_data.data();
7576 this->jacobian_gradients_non_inverse =
7577 this_jacobian_gradients_non_inverse_data.data();
7581 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)