17 #ifndef dealii_matrix_free_fe_evaluation_h 18 #define dealii_matrix_free_fe_evaluation_h 55 <<
"You are requesting information from an FEEvaluation/FEFaceEvaluation " 56 <<
"object for which this kind of information has not been computed. What " 57 <<
"information these objects compute is determined by the update_* flags you " 58 <<
"pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, " 59 <<
"the operation you are attempting requires the <" << arg1
60 <<
"> flag to be set, but it was apparently not specified " 61 <<
"upon initialization.");
66 int n_q_points_1d = fe_degree + 1,
67 int n_components_ = 1,
100 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
104 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
105 "Type of Number and of VectorizedArrayType do not match.");
108 static constexpr
unsigned int dimension = dim;
124 get_mapping_data_index_offset()
const;
133 get_cell_type()
const;
139 get_shape_info()
const;
145 get_dof_info()
const;
152 JxW(
const unsigned int q_point)
const;
166 inverse_jacobian(
const unsigned int q_point)
const;
181 get_normal_vector(
const unsigned int q_point)
const;
200 const VectorizedArrayType &
value)
const;
206 template <
typename T>
207 std::array<T, VectorizedArrayType::size()>
208 read_cell_data(
const AlignedVector<std::array<
T, VectorizedArrayType::size()>>
215 template <
typename T>
219 const std::array<
T, VectorizedArrayType::size()> &
value)
const;
225 std::array<unsigned int, VectorizedArrayType::size()>
226 get_cell_ids()
const;
232 std::array<unsigned int, VectorizedArrayType::size()>
233 get_cell_or_face_ids()
const;
241 const std::vector<unsigned int> &
242 get_internal_dof_numbering()
const;
251 get_scratch_data()
const;
257 get_quadrature_index()
const;
263 get_current_cell_index()
const;
270 get_active_fe_index()
const;
277 get_active_quadrature_index()
const;
283 get_matrix_free()
const;
294 const unsigned int dof_no,
295 const unsigned int first_selected_component,
296 const unsigned int quad_no,
297 const unsigned int fe_degree,
298 const unsigned int n_q_points,
299 const bool is_interior_face,
300 const unsigned int active_fe_index,
301 const unsigned int active_quad_index,
302 const unsigned int face_type);
313 const unsigned int first_selected_component,
370 (is_face ? dim - 1 : dim),
373 VectorizedArrayType> *mapping_data;
392 (is_face ? dim - 1 : dim),
395 VectorizedArrayType>::QuadratureDescriptor *descriptor;
490 std::shared_ptr<internal::MatrixFreeFunctions::
491 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
496 template <
int,
int,
int,
int,
typename,
typename>
542 bool is_face =
false,
552 static constexpr
unsigned int dimension = dim;
553 static constexpr
unsigned int n_components = n_components_;
591 template <
typename VectorType>
593 read_dof_values(
const VectorType &src,
const unsigned int first_index = 0);
623 template <
typename VectorType>
626 const unsigned int first_index = 0);
659 template <
typename VectorType>
661 distribute_local_to_global(
663 const unsigned int first_index = 0,
664 const std::bitset<VectorizedArrayType::size()> &mask =
665 std::bitset<VectorizedArrayType::size()>().flip())
const;
705 template <
typename VectorType>
708 const unsigned int first_index = 0,
709 const std::bitset<VectorizedArrayType::size()> &mask =
710 std::bitset<VectorizedArrayType::size()>().flip())
const;
715 template <
typename VectorType>
717 set_dof_values_plain(
719 const unsigned int first_index = 0,
720 const std::bitset<VectorizedArrayType::size()> &mask =
721 std::bitset<VectorizedArrayType::size()>().flip())
const;
746 get_dof_value(
const unsigned int dof)
const;
759 submit_dof_value(
const value_type val_in,
const unsigned int dof);
774 get_value(
const unsigned int q_point)
const;
789 submit_value(
const value_type val_in,
const unsigned int q_point);
802 get_gradient(
const unsigned int q_point)
const;
819 get_normal_derivative(
const unsigned int q_point)
const;
834 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point);
855 submit_normal_derivative(
const value_type grad_in,
856 const unsigned int q_point);
870 get_hessian(
const unsigned int q_point)
const;
882 get_hessian_diagonal(
const unsigned int q_point)
const;
896 get_laplacian(
const unsigned int q_point)
const;
910 get_divergence(
const unsigned int q_point)
const;
921 get_symmetric_gradient(
const unsigned int q_point)
const;
930 get_curl(
const unsigned int q_point)
const;
948 submit_divergence(
const VectorizedArrayType div_in,
949 const unsigned int q_point);
968 submit_symmetric_gradient(
970 const unsigned int q_point);
986 const unsigned int q_point);
1007 integrate_value()
const;
1023 const VectorizedArrayType *
1024 begin_dof_values()
const;
1034 VectorizedArrayType *
1047 const VectorizedArrayType *
1048 begin_values()
const;
1060 VectorizedArrayType *
1074 const VectorizedArrayType *
1075 begin_gradients()
const;
1088 VectorizedArrayType *
1103 const VectorizedArrayType *
1104 begin_hessians()
const;
1118 VectorizedArrayType *
1127 get_first_selected_component()
const;
1140 const unsigned int dof_no,
1141 const unsigned int first_selected_component,
1142 const unsigned int quad_no,
1143 const unsigned int fe_degree,
1144 const unsigned int n_q_points,
1145 const bool is_interior_face,
1146 const unsigned int active_fe_index,
1147 const unsigned int active_quad_index,
1148 const unsigned int face_type);
1190 const unsigned int first_selected_component,
1217 template <
typename VectorType,
typename VectorOperation>
1219 read_write_operation(
1221 const std::array<VectorType *, n_components_> &vectors,
1224 n_components_> & vectors_sm,
1225 const std::bitset<VectorizedArrayType::size()> &mask,
1226 const bool apply_constraints =
true)
const;
1235 template <
typename VectorType,
typename VectorOperation>
1237 read_write_operation_contiguous(
1239 const std::array<VectorType *, n_components_> &vectors,
1242 n_components_> & vectors_sm,
1243 const std::bitset<VectorizedArrayType::size()> &mask)
const;
1252 template <
typename VectorType,
typename VectorOperation>
1254 read_write_operation_global(
1256 const std::array<VectorType *, n_components_> &vectors)
const;
1270 VectorizedArrayType *values_dofs[n_components];
1379 set_data_pointers();
1400 VectorizedArrayType>
1403 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1404 "Type of Number and of VectorizedArrayType do not match.");
1411 static constexpr
unsigned int dimension = dim;
1412 static constexpr
unsigned int n_components = n_components_;
1426 const unsigned int dof_no,
1427 const unsigned int first_selected_component,
1428 const unsigned int quad_no,
1429 const unsigned int fe_degree,
1430 const unsigned int n_q_points,
1431 const bool is_interior_face =
true,
1445 const unsigned int first_selected_component,
1471 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
1476 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1477 "Type of Number and of VectorizedArrayType do not match.");
1483 static constexpr
unsigned int dimension = dim;
1490 get_dof_value(
const unsigned int dof)
const;
1495 submit_dof_value(
const value_type val_in,
const unsigned int dof);
1500 get_value(
const unsigned int q_point)
const;
1505 submit_value(
const value_type val_in,
const unsigned int q_point);
1511 const unsigned int q_point);
1516 get_gradient(
const unsigned int q_point)
const;
1521 get_normal_derivative(
const unsigned int q_point)
const;
1526 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point);
1531 submit_normal_derivative(
const value_type grad_in,
1532 const unsigned int q_point);
1537 get_hessian(
unsigned int q_point)
const;
1542 get_hessian_diagonal(
const unsigned int q_point)
const;
1547 get_laplacian(
const unsigned int q_point)
const;
1552 integrate_value()
const;
1564 const unsigned int dof_no,
1565 const unsigned int first_selected_component,
1566 const unsigned int quad_no,
1567 const unsigned int fe_degree,
1568 const unsigned int n_q_points,
1569 const bool is_interior_face =
true,
1583 const unsigned int first_selected_component,
1610 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
1615 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1616 "Type of Number and of VectorizedArrayType do not match.");
1622 static constexpr
unsigned int dimension = dim;
1623 static constexpr
unsigned int n_components = dim;
1630 get_gradient(
const unsigned int q_point)
const;
1637 get_divergence(
const unsigned int q_point)
const;
1646 get_symmetric_gradient(
const unsigned int q_point)
const;
1653 get_curl(
const unsigned int q_point)
const;
1658 get_hessian(
const unsigned int q_point)
const;
1663 get_hessian_diagonal(
const unsigned int q_point)
const;
1668 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point);
1681 const unsigned int q_point);
1692 submit_divergence(
const VectorizedArrayType div_in,
1693 const unsigned int q_point);
1704 submit_symmetric_gradient(
1706 const unsigned int q_point);
1714 const unsigned int q_point);
1726 const unsigned int dof_no,
1727 const unsigned int first_selected_component,
1728 const unsigned int quad_no,
1729 const unsigned int dofs_per_cell,
1730 const unsigned int n_q_points,
1731 const bool is_interior_face =
true,
1745 const unsigned int first_selected_component,
1770 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
1775 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1776 "Type of Number and of VectorizedArrayType do not match.");
1782 static constexpr
unsigned int dimension = 1;
1789 get_dof_value(
const unsigned int dof)
const;
1794 submit_dof_value(
const value_type val_in,
const unsigned int dof);
1799 get_value(
const unsigned int q_point)
const;
1804 submit_value(
const value_type val_in,
const unsigned int q_point);
1809 submit_value(
const gradient_type val_in,
const unsigned int q_point);
1814 get_gradient(
const unsigned int q_point)
const;
1819 get_normal_derivative(
const unsigned int q_point)
const;
1824 submit_gradient(
const gradient_type grad_in,
const unsigned int q_point);
1829 submit_gradient(
const value_type grad_in,
const unsigned int q_point);
1834 submit_normal_derivative(
const value_type grad_in,
1835 const unsigned int q_point);
1841 const unsigned int q_point);
1846 get_hessian(
unsigned int q_point)
const;
1851 get_hessian_diagonal(
const unsigned int q_point)
const;
1856 get_laplacian(
const unsigned int q_point)
const;
1861 integrate_value()
const;
1873 const unsigned int dof_no,
1874 const unsigned int first_selected_component,
1875 const unsigned int quad_no,
1876 const unsigned int fe_degree,
1877 const unsigned int n_q_points,
1878 const bool is_interior_face =
true,
1892 const unsigned int first_selected_component,
2465 typename VectorizedArrayType>
2470 VectorizedArrayType>
2473 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2474 "Type of Number and of VectorizedArrayType do not match.");
2505 static constexpr
unsigned int dimension = dim;
2511 static constexpr
unsigned int n_components = n_components_;
2519 static constexpr
unsigned int static_n_q_points =
2529 static constexpr
unsigned int static_dofs_per_component =
2539 static constexpr
unsigned int tensor_dofs_per_cell =
2540 static_dofs_per_component * n_components;
2549 static constexpr
unsigned int static_dofs_per_cell =
2550 static_dofs_per_component * n_components;
2588 const unsigned int dof_no = 0,
2589 const unsigned int quad_no = 0,
2590 const unsigned int first_selected_component = 0,
2602 const std::pair<unsigned int, unsigned int> & range,
2603 const unsigned int dof_no = 0,
2604 const unsigned int quad_no = 0,
2605 const unsigned int first_selected_component = 0);
2637 const unsigned int first_selected_component = 0);
2647 const unsigned int first_selected_component = 0);
2662 const unsigned int first_selected_component = 0);
2690 reinit(
const unsigned int cell_batch_index);
2704 template <
bool level_dof_access>
2738 evaluate(
const bool evaluate_values,
2739 const bool evaluate_gradients,
2740 const bool evaluate_hessians =
false);
2755 evaluate(
const VectorizedArrayType * values_array,
2763 evaluate(
const VectorizedArrayType *values_array,
2764 const bool evaluate_values,
2765 const bool evaluate_gradients,
2766 const bool evaluate_hessians =
false);
2781 template <
typename VectorType>
2783 gather_evaluate(
const VectorType & input_vector,
2789 template <
typename VectorType>
2791 gather_evaluate(
const VectorType &input_vector,
2792 const bool evaluate_values,
2793 const bool evaluate_gradients,
2794 const bool evaluate_hessians =
false);
2814 integrate(
const bool integrate_values,
const bool integrate_gradients);
2829 VectorizedArrayType * values_array);
2835 integrate(
const bool integrate_values,
2836 const bool integrate_gradients,
2837 VectorizedArrayType *values_array);
2852 template <
typename VectorType>
2860 template <
typename VectorType>
2862 integrate_scatter(
const bool integrate_values,
2863 const bool integrate_gradients,
2871 quadrature_point(
const unsigned int q_point)
const;
2904 check_template_arguments(
const unsigned int fe_no,
2905 const unsigned int first_selected_component);
2947 int n_q_points_1d = fe_degree + 1,
2948 int n_components_ = 1,
2949 typename Number =
double,
2955 VectorizedArrayType>
2958 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2959 "Type of Number and of VectorizedArrayType do not match.");
2990 static constexpr
unsigned int dimension = dim;
2996 static constexpr
unsigned int n_components = n_components_;
3005 static constexpr
unsigned int static_n_q_points =
3014 static constexpr
unsigned int static_n_q_points_cell =
3023 static constexpr
unsigned int static_dofs_per_component =
3032 static constexpr
unsigned int tensor_dofs_per_cell =
3033 static_dofs_per_component * n_components;
3041 static constexpr
unsigned int static_dofs_per_cell =
3042 static_dofs_per_component * n_components;
3087 const bool is_interior_face =
true,
3088 const unsigned int dof_no = 0,
3089 const unsigned int quad_no = 0,
3090 const unsigned int first_selected_component = 0,
3104 const std::pair<unsigned int, unsigned int> & range,
3105 const bool is_interior_face =
true,
3106 const unsigned int dof_no = 0,
3107 const unsigned int quad_no = 0,
3108 const unsigned int first_selected_component = 0);
3121 reinit(
const unsigned int face_batch_number);
3131 reinit(
const unsigned int cell_batch_number,
const unsigned int face_number);
3150 evaluate(
const bool evaluate_values,
const bool evaluate_gradients);
3165 evaluate(
const VectorizedArrayType * values_array,
3172 evaluate(
const VectorizedArrayType *values_array,
3173 const bool evaluate_values,
3174 const bool evaluate_gradients);
3187 template <
typename VectorType>
3189 gather_evaluate(
const VectorType & input_vector,
3195 template <
typename VectorType>
3197 gather_evaluate(
const VectorType &input_vector,
3198 const bool evaluate_values,
3199 const bool evaluate_gradients);
3217 integrate(
const bool integrate_values,
const bool integrate_gradients);
3229 VectorizedArrayType * values_array);
3235 integrate(
const bool integrate_values,
3236 const bool integrate_gradients,
3237 VectorizedArrayType *values_array);
3250 template <
typename VectorType>
3258 template <
typename VectorType>
3260 integrate_scatter(
const bool integrate_values,
3261 const bool integrate_gradients,
3269 quadrature_point(
const unsigned int q_point)
const;
3301 std::array<unsigned int, VectorizedArrayType::size()>
3302 compute_face_no_data();
3307 std::array<unsigned int, VectorizedArrayType::size()>
3308 compute_face_orientations();
3315 namespace MatrixFreeFunctions
3319 template <
int dim,
int degree>
3328 template <
int degree>
3331 static constexpr
unsigned int value = degree + 1;
3344 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3348 const unsigned int dof_no,
3349 const unsigned int first_selected_component,
3350 const unsigned int quad_no_in,
3351 const unsigned int fe_degree,
3352 const unsigned int n_q_points,
3353 const bool is_interior_face,
3354 const unsigned int active_fe_index_in,
3355 const unsigned int active_quad_index_in,
3356 const unsigned int face_type)
3357 : scratch_data_array(data_in.acquire_scratch_data())
3358 , quad_no(quad_no_in)
3359 , matrix_info(&data_in)
3360 , dof_info(&data_in.get_dof_info(dof_no))
3364 data_in.get_mapping_info(),
3367 data_in.get_dof_info(dof_no).fe_index_from_degree(
3368 first_selected_component,
3371 active_fe_index_in :
3373 , active_quad_index(
3375 (mapping_data->quad_index_from_n_q_points(n_q_points)) :
3377 active_quad_index_in :
3379 mapping_data->descriptor.size() - 1)))
3381 &mapping_data->descriptor
3383 (active_quad_index *
std::
max<unsigned
int>(1, dim - 1) +
3386 , n_quadrature_points(descriptor->n_q_points)
3387 , data(&data_in.get_shape_info(
3390 dof_info->component_to_base_index[first_selected_component],
3395 , normal_vectors(nullptr)
3396 , normal_x_jacobian(nullptr)
3397 , quadrature_weights(descriptor->quadrature_weights.
begin())
3399 , is_interior_face(is_interior_face)
3410 VectorizedArrayType::size());
3416 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3423 const unsigned int first_selected_component,
3430 , descriptor(
nullptr)
3431 , n_quadrature_points(
3432 Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.
size()))
3433 , matrix_info(
nullptr)
3435 , mapping_data(
nullptr)
3444 , normal_vectors(
nullptr)
3445 , normal_x_jacobian(
nullptr)
3446 , quadrature_weights(
nullptr)
3449 , is_interior_face(
true)
3454 if (other !=
nullptr &&
3461 mapping, quadrature, update_flags);
3465 jacobian = mapped_geometry->get_data_storage().
jacobians[0].begin();
3466 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3471 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3476 : scratch_data_array(other.
matrix_info ==
nullptr ?
3493 , normal_vectors(
nullptr)
3494 , normal_x_jacobian(
nullptr)
3495 , quadrature_weights(other.
descriptor ==
nullptr ?
3497 descriptor->quadrature_weights.begin())
3506 mapped_geometry = std::make_shared<
3507 internal::MatrixFreeFunctions::
3508 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3512 mapping_data = &mapped_geometry->get_data_storage();
3515 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3516 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3522 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3532 if (matrix_info ==
nullptr)
3535 delete scratch_data_array;
3539 matrix_info->release_scratch_data(scratch_data_array);
3555 scratch_data_array = matrix_info->acquire_scratch_data();
3558 quadrature_weights =
3559 (descriptor !=
nullptr ? descriptor->quadrature_weights.begin() :
nullptr);
3568 mapped_geometry = std::make_shared<
3569 internal::MatrixFreeFunctions::
3570 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3575 mapping_data = &mapped_geometry->get_data_storage();
3576 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3577 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3585 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3589 if (matrix_info !=
nullptr)
3593 matrix_info->release_scratch_data(scratch_data_array);
3600 delete scratch_data_array;
3604 scratch_data_array =
nullptr;
3609 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3614 if (matrix_info ==
nullptr)
3619 return this->mapping_data->data_index_offsets[cell];
3625 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3636 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3647 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3652 Assert(dof_info !=
nullptr,
3654 "FEEvaluation was not initialized with a MatrixFree object!"));
3660 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3666 Assert(normal_vectors !=
nullptr,
3668 "update_normal_vectors"));
3670 return normal_vectors[0];
3672 return normal_vectors[q_point];
3677 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3680 const unsigned int q_point)
const 3683 Assert(J_value !=
nullptr,
3685 "update_values|update_gradients"));
3689 return J_value[0] * this->quadrature_weights[q_point];
3692 return J_value[q_point];
3697 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3703 Assert(this->jacobian !=
nullptr,
3705 "update_gradients"));
3709 return jacobian[q_point];
3714 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3715 inline std::array<unsigned int, VectorizedArrayType::size()>
3721 const unsigned int n_lanes = VectorizedArrayType::size();
3722 std::array<unsigned int, n_lanes> cells;
3725 for (
unsigned int i = 0; i < n_lanes; ++i)
3728 if ((is_face ==
false) ||
3730 this->dof_access_index ==
3732 this->is_interior_face))
3735 for (
unsigned int i = 0; i < n_lanes; ++i)
3736 cells[i] = cell * n_lanes + i;
3739 this->dof_access_index ==
3741 this->is_interior_face ==
false)
3747 for (
unsigned int i = 0; i < n_lanes; i++)
3750 const unsigned int cell_this = this->cell * n_lanes + i;
3752 unsigned int face_index =
3753 this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3761 auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3762 .cells_interior[face_index % n_lanes];
3763 auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3764 .cells_exterior[face_index % n_lanes];
3767 if (cell_m == cell_this)
3769 else if (cell_p == cell_this)
3776 const unsigned int *cells_ =
3778 &this->matrix_info->get_face_info(cell).cells_interior[0] :
3779 &this->matrix_info->get_face_info(cell).cells_exterior[0];
3780 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3782 cells[i] = cells_[i];
3794 typename VectorizedArrayType,
3795 typename VectorizedArrayType2,
3796 typename GlobalVectorType,
3802 GlobalVectorType & array,
3803 VectorizedArrayType2 & out,
3815 for (
unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3818 array[cells[i] / VectorizedArrayType::size()]
3819 [cells[i] % VectorizedArrayType::size()]);
3825 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3826 std::array<unsigned int, VectorizedArrayType::size()>
3830 const unsigned int v_len = VectorizedArrayType::size();
3831 std::array<unsigned int, VectorizedArrayType::size()> cells;
3834 for (
unsigned int i = 0; i < v_len; ++i)
3838 this->dof_access_index ==
3840 this->is_interior_face ==
false)
3843 for (
unsigned int i = 0; i < v_len; i++)
3846 const unsigned int cell_this = this->cell * v_len + i;
3849 this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3857 auto cell_m = this->matrix_info->get_face_info(fn / v_len)
3858 .cells_interior[fn % v_len];
3859 auto cell_p = this->matrix_info->get_face_info(fn / v_len)
3860 .cells_exterior[fn % v_len];
3863 if (cell_m == cell_this)
3865 else if (cell_p == cell_this)
3871 for (
unsigned int i = 0; i < v_len; ++i)
3872 cells[i] = cell * v_len + i;
3880 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3881 inline VectorizedArrayType
3885 VectorizedArrayType out = Number(1.);
3886 internal::process_cell_data(
3887 *
this, this->matrix_info, array, out, [](
auto &local,
const auto &global) {
3895 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3899 const VectorizedArrayType & in)
const 3901 internal::process_cell_data(
3902 *
this, this->matrix_info, array, in, [](
const auto &local,
auto &global) {
3909 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3910 template <
typename T>
3911 inline std::array<T, VectorizedArrayType::size()>
3913 const AlignedVector<std::array<
T, VectorizedArrayType::size()>> &array)
const 3915 std::array<T, VectorizedArrayType::size()> out;
3916 internal::process_cell_data(
3917 *
this, this->matrix_info, array, out, [](
auto &local,
const auto &global) {
3925 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3926 template <
typename T>
3929 AlignedVector<std::array<
T, VectorizedArrayType::size()>> &array,
3930 const std::array<
T, VectorizedArrayType::size()> & in)
const 3932 internal::process_cell_data(
3933 *
this, this->matrix_info, array, in, [](
const auto &local,
auto &global) {
3940 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3941 inline const std::vector<unsigned int> &
3945 return data->lexicographic_numbering;
3950 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3956 const_cast<VectorizedArrayType *
>(scratch_data),
3957 scratch_data_array->end() - scratch_data);
3962 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3967 return this->quad_no;
3972 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3977 if (is_face && this->dof_access_index ==
3986 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
3991 return active_fe_index;
3996 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4001 return active_quad_index;
4006 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
4011 Assert(matrix_info !=
nullptr,
4013 "FEEvaluation was not initialized with a MatrixFree object!"));
4014 return *matrix_info;
4024 typename VectorizedArrayType>
4029 VectorizedArrayType>::
4031 const unsigned int dof_no,
4032 const unsigned int first_selected_component,
4033 const unsigned int quad_no_in,
4034 const unsigned int fe_degree,
4035 const unsigned int n_q_points,
4036 const bool is_interior_face,
4037 const unsigned int active_fe_index,
4038 const unsigned int active_quad_index,
4039 const unsigned int face_type)
4043 first_selected_component,
4052 , dof_values_initialized(
false)
4053 , values_quad_initialized(
false)
4054 , gradients_quad_initialized(
false)
4055 , hessians_quad_initialized(
false)
4056 , values_quad_submitted(
false)
4057 , gradients_quad_submitted(
false)
4058 , first_selected_component(first_selected_component)
4060 set_data_pointers();
4062 this->dof_info->start_components.back() == 1 ||
4063 static_cast<int>(n_components_) <=
4065 this->dof_info->start_components
4066 [this->dof_info->component_to_base_index[first_selected_component] +
4068 first_selected_component,
4070 "You tried to construct a vector-valued evaluator with " +
4072 " components. However, " 4073 "the current base element has only " +
4075 this->dof_info->start_components
4076 [this->dof_info->component_to_base_index[first_selected_component] +
4078 first_selected_component) +
4079 " components left when starting from local element index " +
4081 first_selected_component -
4083 [this->dof_info->component_to_base_index[first_selected_component]]) +
4084 " (global index " +
std::to_string(first_selected_component) +
")"));
4096 typename VectorizedArrayType>
4101 VectorizedArrayType>::
4107 const unsigned int first_selected_component,
4115 first_selected_component,
4117 , n_fe_components(n_components_)
4118 , dof_values_initialized(
false)
4119 , values_quad_initialized(
false)
4120 , gradients_quad_initialized(
false)
4121 , hessians_quad_initialized(
false)
4122 , values_quad_submitted(
false)
4123 , gradients_quad_submitted(
false)
4126 , first_selected_component(first_selected_component)
4128 set_data_pointers();
4130 const unsigned int base_element_number =
4131 fe.component_to_base_index(first_selected_component).first;
4132 Assert(fe.element_multiplicity(base_element_number) == 1 ||
4133 fe.element_multiplicity(base_element_number) -
4134 first_selected_component >=
4136 ExcMessage(
"The underlying element must at least contain as many " 4137 "components as requested by this class"));
4138 (void)base_element_number;
4147 typename VectorizedArrayType>
4152 VectorizedArrayType>::
4157 VectorizedArrayType> &other)
4159 , n_fe_components(other.n_fe_components)
4160 , dof_values_initialized(
false)
4161 , values_quad_initialized(
false)
4162 , gradients_quad_initialized(
false)
4163 , hessians_quad_initialized(
false)
4164 , values_quad_submitted(
false)
4165 , gradients_quad_submitted(
false)
4166 , first_selected_component(other.first_selected_component)
4168 set_data_pointers();
4177 typename VectorizedArrayType>
4182 VectorizedArrayType> &
4188 VectorizedArrayType> &other)
4193 AssertDimension(first_selected_component, other.first_selected_component);
4204 typename VectorizedArrayType>
4211 const unsigned int tensor_dofs_per_component =
4212 Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4213 const unsigned int dofs_per_component =
4214 this->data->dofs_per_component_on_cell;
4215 const unsigned int n_quadrature_points = this->n_quadrature_points;
4217 const unsigned int shift =
4218 std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4220 2 * n_quadrature_points;
4221 const unsigned int allocated_size =
4222 shift + n_components_ * dofs_per_component +
4223 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4224 this->scratch_data_array->resize_fast(allocated_size);
4227 for (
unsigned int c = 0; c < n_components_; ++c)
4230 this->scratch_data_array->begin() + c * dofs_per_component;
4233 this->scratch_data_array->begin() + n_components * dofs_per_component;
4234 gradients_quad = this->scratch_data_array->begin() +
4235 n_components * (dofs_per_component + n_quadrature_points);
4237 this->scratch_data_array->begin() +
4238 n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4239 this->scratch_data =
4240 this->scratch_data_array->begin() + n_components_ * dofs_per_component +
4241 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4251 template <
typename VectorType,
bool>
4252 struct BlockVectorSelector
4255 template <
typename VectorType>
4256 struct BlockVectorSelector<VectorType, true>
4258 using BaseVectorType =
typename VectorType::BlockType;
4260 static BaseVectorType *
4261 get_vector_component(VectorType &vec,
const unsigned int component)
4264 return &vec.block(component);
4268 template <
typename VectorType>
4269 struct BlockVectorSelector<VectorType, false>
4271 using BaseVectorType = VectorType;
4273 static BaseVectorType *
4274 get_vector_component(VectorType &vec,
const unsigned int component)
4293 template <
typename VectorType>
4294 struct BlockVectorSelector<std::vector<VectorType>, false>
4296 using BaseVectorType = VectorType;
4298 static BaseVectorType *
4299 get_vector_component(std::vector<VectorType> &vec,
4300 const unsigned int component)
4303 return &vec[component];
4307 template <
typename VectorType>
4308 struct BlockVectorSelector<std::vector<VectorType *>, false>
4310 using BaseVectorType = VectorType;
4312 static BaseVectorType *
4313 get_vector_component(std::vector<VectorType *> &vec,
4314 const unsigned int component)
4317 return vec[component];
4328 typename VectorizedArrayType>
4329 template <
typename VectorType,
typename VectorOperation>
4334 const std::array<VectorType *, n_components_> &src,
4337 n_components_> & src_sm,
4338 const std::bitset<VectorizedArrayType::size()> &mask,
4339 const bool apply_constraints)
const 4344 if (this->matrix_info ==
nullptr)
4346 read_write_operation_global(operation, src);
4352 if (n_fe_components == 1)
4353 for (
unsigned int comp = 0; comp < n_components; ++comp)
4355 Assert(src[comp] !=
nullptr,
4356 ExcMessage(
"The finite element underlying this FEEvaluation " 4357 "object is scalar, but you requested " +
4359 " components via the template argument in " 4360 "FEEvaluation. In that case, you must pass an " 4361 "std::vector<VectorType> or a BlockVector to " +
4362 "read_dof_values and distribute_local_to_global."));
4374 this->dof_info->index_storage_variants[this->dof_access_index].size());
4375 if (this->dof_info->index_storage_variants
4376 [is_face ? this->dof_access_index :
4381 read_write_operation_contiguous(operation, src, src_sm, mask);
4387 constexpr
unsigned int n_lanes = VectorizedArrayType::size();
4388 Assert(mask.count() == n_lanes,
4390 "non-contiguous DoF storage"));
4392 std::integral_constant<
bool,
4396 const unsigned int dofs_per_component =
4397 this->data->dofs_per_component_on_cell;
4398 if (this->dof_info->index_storage_variants
4399 [is_face ? this->dof_access_index :
4404 const unsigned int *dof_indices =
4405 this->dof_info->dof_indices_interleaved.data() +
4406 this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4409 ->component_dof_indices_offset[this->active_fe_index]
4410 [this->first_selected_component] *
4412 if (n_components == 1 || n_fe_components == 1)
4413 for (
unsigned int i = 0; i < dofs_per_component;
4414 ++i, dof_indices += n_lanes)
4415 for (
unsigned int comp = 0; comp < n_components; ++comp)
4416 operation.process_dof_gather(dof_indices,
4419 values_dofs[comp][i],
4422 for (
unsigned int comp = 0; comp < n_components; ++comp)
4423 for (
unsigned int i = 0; i < dofs_per_component;
4424 ++i, dof_indices += n_lanes)
4425 operation.process_dof_gather(
4426 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4430 const unsigned int * dof_indices[n_lanes];
4431 VectorizedArrayType **values_dofs =
4432 const_cast<VectorizedArrayType **
>(&this->values_dofs[0]);
4436 unsigned int cells_copied[n_lanes];
4437 const unsigned int *cells;
4438 unsigned int n_vectorization_actual =
4440 ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4441 bool has_constraints =
false;
4442 const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4445 if (this->dof_access_index ==
4447 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4448 cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4450 this->dof_access_index ==
4453 (this->is_interior_face ?
4454 &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
4455 &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
4456 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4458 Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4460 const std::pair<unsigned int, unsigned int> *my_index_start =
4461 &this->dof_info->row_starts[cells[v] * n_fe_components +
4462 first_selected_component];
4467 if (my_index_start[n_components_read].
second !=
4468 my_index_start[0].
second)
4469 has_constraints =
true;
4472 this->dof_info->dof_indices.data() + my_index_start[0].first;
4474 for (
unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4475 dof_indices[v] =
nullptr;
4480 this->dof_info->row_starts.size());
4481 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4483 const std::pair<unsigned int, unsigned int> *my_index_start =
4485 ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4486 first_selected_component];
4487 if (my_index_start[n_components_read].
second !=
4488 my_index_start[0].
second)
4489 has_constraints =
true;
4491 my_index_start[0].
first ||
4492 my_index_start[0].first < this->dof_info->dof_indices.size(),
4494 my_index_start[0].
first,
4495 this->dof_info->dof_indices.size()));
4497 this->dof_info->dof_indices.data() + my_index_start[0].first;
4499 for (
unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4500 dof_indices[v] =
nullptr;
4505 if (!has_constraints)
4507 if (n_vectorization_actual < n_lanes)
4508 for (
unsigned int comp = 0; comp < n_components; ++comp)
4509 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4510 operation.process_empty(values_dofs[comp][i]);
4511 if (n_components == 1 || n_fe_components == 1)
4513 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4514 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4515 for (
unsigned int comp = 0; comp < n_components; ++comp)
4516 operation.process_dof(dof_indices[v][i],
4518 values_dofs[comp][i][v]);
4522 for (
unsigned int comp = 0; comp < n_components; ++comp)
4523 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4524 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4525 operation.process_dof(
4526 dof_indices[v][comp * dofs_per_component + i],
4528 values_dofs[comp][i][v]);
4538 if (n_vectorization_actual < n_lanes)
4539 for (
unsigned int comp = 0; comp < n_components; ++comp)
4540 for (
unsigned int i = 0; i < dofs_per_component; ++i)
4541 operation.process_empty(values_dofs[comp][i]);
4542 for (
unsigned int v = 0; v < n_vectorization_actual; ++v)
4544 const unsigned int cell_index =
4545 is_face ? cells[v] : this->cell * n_lanes + v;
4546 const unsigned int cell_dof_index =
4547 cell_index * n_fe_components + first_selected_component;
4548 const unsigned int n_components_read =
4549 n_fe_components > 1 ? n_components : 1;
4550 unsigned int index_indicators =
4551 this->dof_info->row_starts[cell_dof_index].second;
4552 unsigned int next_index_indicators =
4553 this->dof_info->row_starts[cell_dof_index + 1].second;
4557 if (apply_constraints ==
false &&
4558 this->dof_info->row_starts[cell_dof_index].second !=
4559 this->dof_info->row_starts[cell_dof_index + n_components_read]
4562 Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
4566 this->dof_info->plain_dof_indices.data() +
4568 ->component_dof_indices_offset[this->active_fe_index]
4569 [this->first_selected_component] +
4570 this->dof_info->row_starts_plain_indices[cell_index];
4571 next_index_indicators = index_indicators;
4574 if (n_components == 1 || n_fe_components == 1)
4576 unsigned int ind_local = 0;
4577 for (; index_indicators != next_index_indicators; ++index_indicators)
4579 const std::pair<unsigned short, unsigned short> indicator =
4580 this->dof_info->constraint_indicator[index_indicators];
4582 for (
unsigned int j = 0; j < indicator.first; ++j)
4583 for (
unsigned int comp = 0; comp < n_components; ++comp)
4584 operation.process_dof(dof_indices[v][j],
4586 values_dofs[comp][ind_local + j][v]);
4588 ind_local += indicator.first;
4589 dof_indices[v] += indicator.first;
4593 Number
value[n_components];
4594 for (
unsigned int comp = 0; comp < n_components; ++comp)
4595 operation.pre_constraints(values_dofs[comp][ind_local][v],
4598 const Number *data_val =
4599 this->matrix_info->constraint_pool_begin(indicator.second);
4600 const Number *end_pool =
4601 this->matrix_info->constraint_pool_end(indicator.second);
4602 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4603 for (
unsigned int comp = 0; comp < n_components; ++comp)
4604 operation.process_constraint(*dof_indices[v],
4609 for (
unsigned int comp = 0; comp < n_components; ++comp)
4610 operation.post_constraints(value[comp],
4611 values_dofs[comp][ind_local][v]);
4617 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4618 for (
unsigned int comp = 0; comp < n_components; ++comp)
4619 operation.process_dof(*dof_indices[v],
4621 values_dofs[comp][ind_local][v]);
4630 for (
unsigned int comp = 0; comp < n_components; ++comp)
4632 unsigned int ind_local = 0;
4635 for (; index_indicators != next_index_indicators;
4638 const std::pair<unsigned short, unsigned short> indicator =
4639 this->dof_info->constraint_indicator[index_indicators];
4642 for (
unsigned int j = 0; j < indicator.first; ++j)
4643 operation.process_dof(dof_indices[v][j],
4645 values_dofs[comp][ind_local + j][v]);
4646 ind_local += indicator.first;
4647 dof_indices[v] += indicator.first;
4652 operation.pre_constraints(values_dofs[comp][ind_local][v],
4655 const Number *data_val =
4656 this->matrix_info->constraint_pool_begin(indicator.second);
4657 const Number *end_pool =
4658 this->matrix_info->constraint_pool_end(indicator.second);
4660 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4661 operation.process_constraint(*dof_indices[v],
4666 operation.post_constraints(value,
4667 values_dofs[comp][ind_local][v]);
4674 for (; ind_local < dofs_per_component;
4675 ++dof_indices[v], ++ind_local)
4678 operation.process_dof(*dof_indices[v],
4680 values_dofs[comp][ind_local][v]);
4683 if (apply_constraints ==
true && comp + 1 < n_components)
4684 next_index_indicators =
4685 this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4697 typename VectorizedArrayType>
4698 template <
typename VectorType,
typename VectorOperation>
4703 const std::array<VectorType *, n_components_> &src)
const 4707 unsigned int index =
4708 first_selected_component * this->data->dofs_per_component_on_cell;
4709 for (
unsigned int comp = 0; comp < n_components; ++comp)
4711 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4714 operation.process_empty(values_dofs[comp][i]);
4715 operation.process_dof_global(
4716 local_dof_indices[this->data->lexicographic_numbering[index]],
4718 values_dofs[comp][i][0]);
4729 typename VectorizedArrayType>
4730 template <
typename VectorType,
typename VectorOperation>
4735 const std::array<VectorType *, n_components_> &src,
4738 n_components_> & vectors_sm,
4739 const std::bitset<VectorizedArrayType::size()> &mask)
const 4748 std::integral_constant<
bool,
4752 is_face ? this->dof_access_index :
4754 const unsigned int n_lanes = mask.count();
4756 const std::vector<unsigned int> &dof_indices_cont =
4757 this->dof_info->dof_indices_contiguous[ind];
4761 if ((this->dof_info->index_storage_variants[ind][this->cell] ==
4763 interleaved_contiguous &&
4764 n_lanes == VectorizedArrayType::size()) &&
4766 this->dof_access_index ==
4768 this->is_interior_face ==
false))
4770 const unsigned int dof_index =
4771 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4772 this->dof_info->component_dof_indices_offset[this->active_fe_index]
4773 [first_selected_component] *
4774 VectorizedArrayType::size();
4775 if (n_components == 1 || n_fe_components == 1)
4776 for (
unsigned int comp = 0; comp < n_components; ++comp)
4777 operation.process_dofs_vectorized(
4778 this->data->dofs_per_component_on_cell,
4784 operation.process_dofs_vectorized(
4785 this->data->dofs_per_component_on_cell * n_components,
4793 std::array<unsigned int, VectorizedArrayType::size()> cells =
4794 this->get_cell_or_face_ids();
4798 const unsigned int n_filled_lanes =
4799 this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4802 this->dof_access_index ==
4804 this->is_interior_face ==
false;
4806 if (vectors_sm[0] !=
nullptr)
4808 const auto compute_vector_ptrs = [&](
const unsigned int comp) {
4809 std::array<
typename VectorType::value_type *,
4810 VectorizedArrayType::size()>
4813 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4817 Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
4819 ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
4821 this->dof_info->dof_indices_contiguous_sm[ind].size(),
4825 this->dof_info->dof_indices_contiguous_sm[ind].size()));
4828 this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
4831 vector_ptrs[v] =
const_cast<typename VectorType::value_type *
>(
4832 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
4833 this->dof_info->component_dof_indices_offset
4834 [this->active_fe_index][this->first_selected_component]);
4836 vector_ptrs[v] =
nullptr;
4838 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
4840 vector_ptrs[v] =
nullptr;
4845 if (n_filled_lanes == VectorizedArrayType::size() &&
4846 n_lanes == VectorizedArrayType::size() && !is_ecl)
4848 if (n_components == 1 || n_fe_components == 1)
4850 for (
unsigned int comp = 0; comp < n_components; ++comp)
4852 auto vector_ptrs = compute_vector_ptrs(comp);
4853 operation.process_dofs_vectorized_transpose(
4854 this->data->dofs_per_component_on_cell,
4862 auto vector_ptrs = compute_vector_ptrs(0);
4863 operation.process_dofs_vectorized_transpose(
4864 this->data->dofs_per_component_on_cell * n_components,
4871 for (
unsigned int comp = 0; comp < n_components; ++comp)
4873 auto vector_ptrs = compute_vector_ptrs(
4874 (n_components == 1 || n_fe_components == 1) ? comp : 0);
4876 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4878 operation.process_empty(values_dofs[comp][i]);
4880 if (n_components == 1 || n_fe_components == 1)
4882 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4883 if (mask[v] ==
true)
4884 for (
unsigned int i = 0;
4885 i < this->data->dofs_per_component_on_cell;
4887 operation.process_dof(vector_ptrs[v][i],
4888 values_dofs[comp][i][v]);
4892 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4893 if (mask[v] ==
true)
4894 for (
unsigned int i = 0;
4895 i < this->data->dofs_per_component_on_cell;
4897 operation.process_dof(
4899 [i + comp * this->data
4900 ->dofs_per_component_on_cell],
4901 values_dofs[comp][i][v]);
4907 unsigned int dof_indices[VectorizedArrayType::size()];
4909 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
4913 dof_indices_cont[cells[v]] +
4915 ->component_dof_indices_offset[this->active_fe_index]
4916 [this->first_selected_component] *
4917 this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
4920 for (
unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4925 if (n_filled_lanes == VectorizedArrayType::size() &&
4926 n_lanes == VectorizedArrayType::size() && !is_ecl)
4928 if (this->dof_info->index_storage_variants[ind][this->cell] ==
4932 if (n_components == 1 || n_fe_components == 1)
4933 for (
unsigned int comp = 0; comp < n_components; ++comp)
4934 operation.process_dofs_vectorized_transpose(
4935 this->data->dofs_per_component_on_cell,
4941 operation.process_dofs_vectorized_transpose(
4942 this->data->dofs_per_component_on_cell * n_components,
4948 else if (this->dof_info->index_storage_variants[ind][this->cell] ==
4950 interleaved_contiguous_strided)
4952 if (n_components == 1 || n_fe_components == 1)
4953 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4956 for (
unsigned int comp = 0; comp < n_components; ++comp)
4957 operation.process_dof_gather(dof_indices,
4959 i * VectorizedArrayType::size(),
4960 values_dofs[comp][i],
4964 for (
unsigned int comp = 0; comp < n_components; ++comp)
4965 for (
unsigned int i = 0;
4966 i < this->data->dofs_per_component_on_cell;
4969 operation.process_dof_gather(
4972 (comp * this->data->dofs_per_component_on_cell + i) *
4973 VectorizedArrayType::size(),
4974 values_dofs[comp][i],
4980 Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
4982 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4984 const unsigned int *offsets =
4985 &this->dof_info->dof_indices_interleave_strides
4986 [ind][VectorizedArrayType::size() * this->cell];
4987 if (n_components == 1 || n_fe_components == 1)
4988 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4991 for (
unsigned int comp = 0; comp < n_components; ++comp)
4992 operation.process_dof_gather(dof_indices,
4995 values_dofs[comp][i],
4998 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4999 dof_indices[v] += offsets[v];
5002 for (
unsigned int comp = 0; comp < n_components; ++comp)
5003 for (
unsigned int i = 0;
5004 i < this->data->dofs_per_component_on_cell;
5007 operation.process_dof_gather(dof_indices,
5010 values_dofs[comp][i],
5013 for (
unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5014 dof_indices[v] += offsets[v];
5019 for (
unsigned int comp = 0; comp < n_components; ++comp)
5021 for (
unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5023 operation.process_empty(values_dofs[comp][i]);
5024 if (this->dof_info->index_storage_variants[ind][this->cell] ==
5028 if (n_components == 1 || n_fe_components == 1)
5030 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5031 if (mask[v] ==
true)
5032 for (
unsigned int i = 0;
5033 i < this->data->dofs_per_component_on_cell;
5035 operation.process_dof(dof_indices[v] + i,
5037 values_dofs[comp][i][v]);
5041 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5042 if (mask[v] ==
true)
5043 for (
unsigned int i = 0;
5044 i < this->data->dofs_per_component_on_cell;
5046 operation.process_dof(
5047 dof_indices[v] + i +
5048 comp * this->data->dofs_per_component_on_cell,
5050 values_dofs[comp][i][v]);
5055 const unsigned int *offsets =
5056 &this->dof_info->dof_indices_interleave_strides
5057 [ind][VectorizedArrayType::size() * this->cell];
5058 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5060 if (n_components == 1 || n_fe_components == 1)
5061 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5063 if (mask[v] ==
true)
5064 for (
unsigned int i = 0;
5065 i < this->data->dofs_per_component_on_cell;
5067 operation.process_dof(dof_indices[v] + i * offsets[v],
5069 values_dofs[comp][i][v]);
5073 for (
unsigned int v = 0; v < n_filled_lanes; ++v)
5074 if (mask[v] ==
true)
5075 for (
unsigned int i = 0;
5076 i < this->data->dofs_per_component_on_cell;
5078 operation.process_dof(
5080 (i + comp * this->data->dofs_per_component_on_cell) *
5083 values_dofs[comp][i][v]);
5091 template <
typename Number,
5093 typename std::enable_if<!IsBlockVector<VectorType>::value,
5094 VectorType>::type * =
nullptr>
5095 decltype(std::declval<VectorType>().
begin())
5096 get_beginning(VectorType &vec)
5101 template <
typename Number,
5102 typename VectorType,
5103 typename std::enable_if<IsBlockVector<VectorType>::value,
5104 VectorType>::type * =
nullptr>
5105 typename VectorType::value_type *
5106 get_beginning(VectorType &)
5111 template <
typename VectorType,
5112 typename std::enable_if<has_shared_vector_data<VectorType>::value,
5113 VectorType>::type * =
nullptr>
5114 const std::vector<ArrayView<const typename VectorType::value_type>> *
5115 get_shared_vector_data(VectorType & vec,
5116 const bool is_valid_mode_for_sm,
5117 const unsigned int active_fe_index,
5121 if (is_valid_mode_for_sm &&
5124 active_fe_index == 0)
5125 return &vec.shared_vector_data();
5130 template <
typename VectorType,
5131 typename std::enable_if<!has_shared_vector_data<VectorType>::value,
5132 VectorType>::type * =
nullptr>
5133 const std::vector<ArrayView<const typename VectorType::value_type>> *
5134 get_shared_vector_data(VectorType &,
5142 template <
int n_components,
typename VectorType>
5144 std::array<
typename internal::BlockVectorSelector<
5145 typename std::remove_const<VectorType>::type,
5147 value>::BaseVectorType *,
5150 const std::vector<
ArrayView<
const typename internal::BlockVectorSelector<
5151 typename std::remove_const<VectorType>::type,
5153 BaseVectorType::value_type>> *,
5155 get_vector_data(VectorType & src,
5156 const unsigned int first_index,
5157 const bool is_valid_mode_for_sm,
5158 const unsigned int active_fe_index,
5164 std::array<
typename internal::BlockVectorSelector<
5165 typename std::remove_const<VectorType>::type,
5167 value>::BaseVectorType *,
5171 ArrayView<
const typename internal::BlockVectorSelector<
5172 typename std::remove_const<VectorType>::type,
5174 value>::BaseVectorType::value_type>> *,
5178 for (
unsigned int d = 0;
d < n_components; ++
d)
5179 src_data.first[
d] = internal::BlockVectorSelector<
5180 typename std::remove_const<VectorType>::type,
5182 get_vector_component(
5183 const_cast<typename std::remove_const<VectorType>::type &
>(src),
5186 for (
unsigned int d = 0;
d < n_components; ++
d)
5187 src_data.second[
d] = get_shared_vector_data(*src_data.first[
d],
5188 is_valid_mode_for_sm,
5202 typename VectorizedArrayType>
5203 template <
typename VectorType>
5208 const auto src_data = internal::get_vector_data<n_components_>(
5211 this->dof_access_index ==
5213 this->active_fe_index,
5217 read_write_operation(reader,
5220 std::bitset<VectorizedArrayType::size()>().flip(),
5224 dof_values_initialized =
true;
5234 typename VectorizedArrayType>
5235 template <
typename VectorType>
5240 const auto src_data = internal::get_vector_data<n_components_>(
5243 this->dof_access_index ==
5245 this->active_fe_index,
5249 read_write_operation(reader,
5252 std::bitset<VectorizedArrayType::size()>().flip(),
5256 dof_values_initialized =
true;
5266 typename VectorizedArrayType>
5267 template <
typename VectorType>
5272 const unsigned int first_index,
5273 const std::bitset<VectorizedArrayType::size()> &mask)
const 5276 Assert(dof_values_initialized ==
true,
5280 const auto dst_data = internal::get_vector_data<n_components_>(
5283 this->dof_access_index ==
5285 this->active_fe_index,
5290 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
5299 typename VectorizedArrayType>
5300 template <
typename VectorType>
5304 const unsigned int first_index,
5305 const std::bitset<VectorizedArrayType::size()> &mask)
const 5308 Assert(dof_values_initialized ==
true,
5312 const auto dst_data = internal::get_vector_data<n_components_>(
5315 this->dof_access_index ==
5317 this->active_fe_index,
5321 read_write_operation(setter, dst_data.first, dst_data.second, mask);
5330 typename VectorizedArrayType>
5331 template <
typename VectorType>
5336 const unsigned int first_index,
5337 const std::bitset<VectorizedArrayType::size()> &mask)
const 5340 Assert(dof_values_initialized ==
true,
5344 const auto dst_data = internal::get_vector_data<n_components_>(
5347 this->dof_access_index ==
5349 this->active_fe_index,
5353 read_write_operation(setter, dst_data.first, dst_data.second, mask,
false);
5366 typename VectorizedArrayType>
5367 inline const VectorizedArrayType *
5371 return &values_dofs[0][0];
5380 typename VectorizedArrayType>
5381 inline VectorizedArrayType *
5386 dof_values_initialized =
true;
5388 return &values_dofs[0][0];
5397 typename VectorizedArrayType>
5398 inline const VectorizedArrayType *
5414 typename VectorizedArrayType>
5415 inline VectorizedArrayType *
5420 values_quad_initialized =
true;
5421 values_quad_submitted =
true;
5432 typename VectorizedArrayType>
5433 inline const VectorizedArrayType *
5438 Assert(gradients_quad_initialized || gradients_quad_submitted,
5441 return gradients_quad;
5450 typename VectorizedArrayType>
5451 inline VectorizedArrayType *
5456 gradients_quad_submitted =
true;
5457 gradients_quad_initialized =
true;
5459 return gradients_quad;
5468 typename VectorizedArrayType>
5469 inline const VectorizedArrayType *
5476 return hessians_quad;
5485 typename VectorizedArrayType>
5486 inline VectorizedArrayType *
5491 hessians_quad_initialized =
true;
5493 return hessians_quad;
5502 typename VectorizedArrayType>
5507 return first_selected_component;
5516 typename VectorizedArrayType>
5523 for (
unsigned int comp = 0; comp < n_components; comp++)
5524 return_value[comp] = this->values_dofs[comp][dof];
5525 return return_value;
5534 typename VectorizedArrayType>
5537 get_value(
const unsigned int q_point)
const 5540 Assert(this->values_quad_initialized ==
true,
5545 const std::size_t nqp = this->n_quadrature_points;
5547 for (
unsigned int comp = 0; comp < n_components; comp++)
5548 return_value[comp] = values_quad[comp * nqp + q_point];
5549 return return_value;
5558 typename VectorizedArrayType>
5565 Assert(this->gradients_quad_initialized ==
true,
5570 Assert(this->jacobian !=
nullptr,
5572 "update_gradients"));
5573 const std::size_t nqp = this->n_quadrature_points;
5579 for (
unsigned int d = 0;
d < dim; ++
d)
5580 for (
unsigned int comp = 0; comp < n_components; comp++)
5581 grad_out[comp][
d] = gradients_quad[(comp * dim +
d) * nqp + q_point] *
5582 this->jacobian[0][
d][
d];
5591 for (
unsigned int comp = 0; comp < n_components; comp++)
5592 for (
unsigned int d = 0;
d < dim; ++
d)
5595 jac[
d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5596 for (
unsigned int e = 1;
e < dim; ++
e)
5597 grad_out[comp][
d] +=
5598 jac[
d][
e] * gradients_quad[(comp * dim +
e) * nqp + q_point];
5610 typename VectorizedArrayType>
5617 Assert(this->gradients_quad_initialized ==
true,
5621 Assert(this->normal_x_jacobian !=
nullptr,
5623 "update_gradients"));
5625 const std::size_t nqp = this->n_quadrature_points;
5629 for (
unsigned int comp = 0; comp < n_components; comp++)
5630 grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5631 (this->normal_x_jacobian[0][dim - 1]);
5634 const std::size_t index =
5636 for (
unsigned int comp = 0; comp < n_components; comp++)
5638 grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5639 this->normal_x_jacobian[index][0];
5640 for (
unsigned int d = 1;
d < dim; ++
d)
5641 grad_out[comp] += gradients_quad[(comp * dim +
d) * nqp + q_point] *
5642 this->normal_x_jacobian[index][
d];
5654 template <
typename VectorizedArrayType>
5657 const VectorizedArrayType *
const hessians,
5659 VectorizedArrayType (&tmp)[1][1])
5661 tmp[0][0] = jac[0][0] * hessians[0];
5664 template <
typename VectorizedArrayType>
5667 const VectorizedArrayType *
const hessians,
5668 const unsigned int nqp,
5669 VectorizedArrayType (&tmp)[2][2])
5671 for (
unsigned int d = 0;
d < 2; ++
d)
5673 tmp[0][
d] = (jac[
d][0] * hessians[0] + jac[
d][1] * hessians[2 * nqp]);
5675 (jac[
d][0] * hessians[2 * nqp] + jac[
d][1] * hessians[1 * nqp]);
5679 template <
typename VectorizedArrayType>
5682 const VectorizedArrayType *
const hessians,
5683 const unsigned int nqp,
5684 VectorizedArrayType (&tmp)[3][3])
5686 for (
unsigned int d = 0;
d < 3; ++
d)
5689 (jac[
d][0] * hessians[0 * nqp] + jac[
d][1] * hessians[3 * nqp] +
5690 jac[
d][2] * hessians[4 * nqp]);
5692 (jac[
d][0] * hessians[3 * nqp] + jac[
d][1] * hessians[1 * nqp] +
5693 jac[
d][2] * hessians[5 * nqp]);
5695 (jac[
d][0] * hessians[4 * nqp] + jac[
d][1] * hessians[5 * nqp] +
5696 jac[
d][2] * hessians[2 * nqp]);
5707 typename VectorizedArrayType>
5714 Assert(this->hessians_quad_initialized ==
true,
5719 Assert(this->jacobian !=
nullptr,
5729 const std::size_t nqp = this->n_quadrature_points;
5730 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
5735 for (
unsigned int comp = 0; comp < n_components; comp++)
5737 for (
unsigned int d = 0;
d < dim; ++
d)
5738 hessian_out[comp][
d][
d] =
5739 hessians_quad[(comp * hdim +
d) * nqp + q_point] *
5740 (jac[
d][
d] * jac[
d][
d]);
5746 hessian_out[comp][0][1] =
5747 hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5748 (jac[0][0] * jac[1][1]);
5751 hessian_out[comp][0][1] =
5752 hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5753 (jac[0][0] * jac[1][1]);
5754 hessian_out[comp][0][2] =
5755 hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5756 (jac[0][0] * jac[2][2]);
5757 hessian_out[comp][1][2] =
5758 hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5759 (jac[1][1] * jac[2][2]);
5764 for (
unsigned int d = 0;
d < dim; ++
d)
5765 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5766 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5772 for (
unsigned int comp = 0; comp < n_components; comp++)
5774 VectorizedArrayType tmp[dim][dim];
5775 internal::hessian_unit_times_jac(
5776 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5779 for (
unsigned int d = 0;
d < dim; ++
d)
5780 for (
unsigned int e =
d;
e < dim; ++
e)
5782 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
5783 for (
unsigned int f = 1; f < dim; ++f)
5784 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
5791 for (
unsigned int d = 0;
d < dim; ++
d)
5792 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5793 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5799 const auto &jac_grad =
5800 this->mapping_data->jacobian_gradients
5801 [1 - this->is_interior_face]
5802 [this->mapping_data->data_index_offsets[this->cell] + q_point];
5803 for (
unsigned int comp = 0; comp < n_components; comp++)
5807 VectorizedArrayType tmp[dim][dim];
5808 internal::hessian_unit_times_jac(
5809 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5812 for (
unsigned int d = 0;
d < dim; ++
d)
5813 for (
unsigned int e =
d;
e < dim; ++
e)
5815 hessian_out[comp][
d][
e] = jac[
d][0] * tmp[0][
e];
5816 for (
unsigned int f = 1; f < dim; ++f)
5817 hessian_out[comp][
d][
e] += jac[
d][f] * tmp[f][
e];
5821 for (
unsigned int d = 0;
d < dim; ++
d)
5822 for (
unsigned int e = 0;
e < dim; ++
e)
5823 hessian_out[comp][
d][
d] +=
5825 gradients_quad[(comp * dim +
e) * nqp + q_point];
5828 for (
unsigned int d = 0, count = dim;
d < dim; ++
d)
5829 for (
unsigned int e =
d + 1;
e < dim; ++
e, ++count)
5830 for (
unsigned int f = 0; f < dim; ++f)
5831 hessian_out[comp][
d][
e] +=
5832 jac_grad[count][f] *
5833 gradients_quad[(comp * dim + f) * nqp + q_point];
5836 for (
unsigned int d = 0;
d < dim; ++
d)
5837 for (
unsigned int e =
d + 1;
e < dim; ++
e)
5838 hessian_out[comp][
e][
d] = hessian_out[comp][
d][
e];
5850 typename VectorizedArrayType>
5857 Assert(this->hessians_quad_initialized ==
true,
5868 const std::size_t nqp = this->n_quadrature_points;
5869 constexpr
unsigned int hdim = (dim * (dim + 1)) / 2;
5875 for (
unsigned int comp = 0; comp < n_components; comp++)
5876 for (
unsigned int d = 0;
d < dim; ++
d)
5877 hessian_out[comp][
d] =
5878 hessians_quad[(comp * hdim +
d) * nqp + q_point] *
5879 (jac[
d][
d] * jac[
d][
d]);
5884 for (
unsigned int comp = 0; comp < n_components; comp++)
5888 VectorizedArrayType tmp[dim][dim];
5889 internal::hessian_unit_times_jac(
5890 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5894 for (
unsigned int d = 0;
d < dim; ++
d)
5896 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
5897 for (
unsigned int f = 1; f < dim; ++f)
5898 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
5907 this->mapping_data->jacobian_gradients
5908 [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5909 for (
unsigned int comp = 0; comp < n_components; comp++)
5913 VectorizedArrayType tmp[dim][dim];
5914 internal::hessian_unit_times_jac(
5915 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5919 for (
unsigned int d = 0;
d < dim; ++
d)
5921 hessian_out[comp][
d] = jac[
d][0] * tmp[0][
d];
5922 for (
unsigned int f = 1; f < dim; ++f)
5923 hessian_out[comp][
d] += jac[
d][f] * tmp[f][
d];
5926 for (
unsigned int d = 0;
d < dim; ++
d)
5927 for (
unsigned int e = 0;
e < dim; ++
e)
5928 hessian_out[comp][
d] +=
5930 gradients_quad[(comp * dim +
e) * nqp + q_point];
5942 typename VectorizedArrayType>
5949 Assert(this->hessians_quad_initialized ==
true,
5955 const auto hess_diag = get_hessian_diagonal(q_point);
5956 for (
unsigned int comp = 0; comp < n_components; ++comp)
5958 laplacian_out[comp] = hess_diag[comp][0];
5959 for (
unsigned int d = 1;
d < dim; ++
d)
5960 laplacian_out[comp] += hess_diag[comp][
d];
5962 return laplacian_out;
5971 typename VectorizedArrayType>
5975 const unsigned int dof)
5978 this->dof_values_initialized =
true;
5981 for (
unsigned int comp = 0; comp < n_components; comp++)
5982 this->values_dofs[comp][dof] = val_in[comp];
5991 typename VectorizedArrayType>
5995 const unsigned int q_point)
5999 Assert(this->J_value !=
nullptr,
6003 this->values_quad_submitted =
true;
6006 const std::size_t nqp = this->n_quadrature_points;
6009 const VectorizedArrayType JxW =
6010 this->J_value[0] * this->quadrature_weights[q_point];
6011 for (
unsigned int comp = 0; comp < n_components; ++comp)
6012 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6016 const VectorizedArrayType JxW = this->J_value[q_point];
6017 for (
unsigned int comp = 0; comp < n_components; ++comp)
6018 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6028 typename VectorizedArrayType>
6033 const unsigned int q_point)
6037 Assert(this->J_value !=
nullptr,
6039 "update_gradients"));
6040 Assert(this->jacobian !=
nullptr,
6042 "update_gradients"));
6044 this->gradients_quad_submitted =
true;
6047 const std::size_t nqp = this->n_quadrature_points;
6050 const VectorizedArrayType JxW =
6051 this->J_value[0] * this->quadrature_weights[q_point];
6052 for (
unsigned int d = 0;
d < dim; ++
d)
6054 const VectorizedArrayType factor = this->jacobian[0][
d][
d] * JxW;
6055 for (
unsigned int comp = 0; comp < n_components; comp++)
6056 gradients_quad[(comp * dim +
d) * nqp + q_point] =
6057 grad_in[comp][
d] * factor;
6064 this->jacobian[q_point] :
6066 const VectorizedArrayType JxW =
6068 this->J_value[q_point] :
6069 this->J_value[0] * this->quadrature_weights[q_point];
6070 for (
unsigned int comp = 0; comp < n_components; ++comp)
6071 for (
unsigned int d = 0;
d < dim; ++
d)
6073 VectorizedArrayType new_val = jac[0][
d] * grad_in[comp][0];
6074 for (
unsigned int e = 1;
e < dim; ++
e)
6075 new_val += (jac[
e][
d] * grad_in[comp][
e]);
6076 gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val * JxW;
6087 typename VectorizedArrayType>
6092 const unsigned int q_point)
6095 Assert(this->normal_x_jacobian !=
nullptr,
6097 "update_gradients"));
6099 this->gradients_quad_submitted =
true;
6102 const std::size_t nqp = this->n_quadrature_points;
6104 for (
unsigned int comp = 0; comp < n_components; comp++)
6106 for (
unsigned int d = 0;
d < dim - 1; ++
d)
6107 gradients_quad[(comp * dim +
d) * nqp + q_point] =
6108 VectorizedArrayType();
6109 gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
6111 (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
6112 this->quadrature_weights[q_point]);
6116 const unsigned int index =
6119 this->normal_x_jacobian[index];
6120 for (
unsigned int comp = 0; comp < n_components; comp++)
6122 VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
6124 factor = factor * this->quadrature_weights[q_point];
6125 for (
unsigned int d = 0;
d < dim; ++
d)
6126 gradients_quad[(comp * dim +
d) * nqp + q_point] = factor * jac[
d];
6137 typename VectorizedArrayType>
6144 Assert(this->values_quad_submitted ==
true,
6149 const std::size_t nqp = this->n_quadrature_points;
6150 for (
unsigned int q = 0; q < nqp; ++q)
6151 for (
unsigned int comp = 0; comp < n_components; ++comp)
6152 return_value[comp] += this->values_quad[comp * nqp + q];
6153 return (return_value);
6165 typename VectorizedArrayType>
6170 VectorizedArrayType>::
6173 const unsigned int dof_no,
6174 const unsigned int first_selected_component,
6175 const unsigned int quad_no_in,
6176 const unsigned int fe_degree,
6177 const unsigned int n_q_points,
6178 const bool is_interior_face,
6179 const unsigned int active_fe_index,
6180 const unsigned int active_quad_index,
6181 const unsigned int face_type)
6185 first_selected_component,
6201 typename VectorizedArrayType>
6206 VectorizedArrayType>::
6212 const unsigned int first_selected_component,
6220 first_selected_component,
6230 typename VectorizedArrayType>
6235 VectorizedArrayType>::
6240 VectorizedArrayType> &other)
6251 typename VectorizedArrayType>
6256 VectorizedArrayType> &
6262 VectorizedArrayType> &other)
6268 VectorizedArrayType>::operator=(other);
6277 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6281 const unsigned int dof_no,
6282 const unsigned int first_selected_component,
6283 const unsigned int quad_no_in,
6284 const unsigned int fe_degree,
6285 const unsigned int n_q_points,
6286 const bool is_interior_face,
6287 const unsigned int active_fe_index,
6288 const unsigned int active_quad_index,
6289 const unsigned int face_type)
6293 first_selected_component,
6305 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6312 const unsigned int first_selected_component,
6320 first_selected_component,
6326 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6336 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6348 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6351 const unsigned int dof)
const 6354 return this->values_dofs[0][dof];
6359 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6362 const unsigned int q_point)
const 6365 Assert(this->values_quad_initialized ==
true,
6369 return this->values_quad[q_point];
6374 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6379 return BaseClass::get_normal_derivative(q_point)[0];
6384 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6387 const unsigned int q_point)
const 6393 Assert(this->gradients_quad_initialized ==
true,
6398 Assert(this->jacobian !=
nullptr,
6400 "update_gradients"));
6404 const std::size_t nqp = this->n_quadrature_points;
6407 for (
unsigned int d = 0;
d < dim; ++
d)
6409 this->gradients_quad[
d * nqp + q_point] * this->jacobian[0][
d][
d];
6418 for (
unsigned int d = 0;
d < dim; ++
d)
6420 grad_out[
d] = jac[
d][0] * this->gradients_quad[q_point];
6421 for (
unsigned int e = 1;
e < dim; ++
e)
6422 grad_out[
d] += jac[
d][
e] * this->gradients_quad[
e * nqp + q_point];
6430 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6433 const unsigned int q_point)
const 6435 return BaseClass::get_hessian(q_point)[0];
6440 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6445 return BaseClass::get_hessian_diagonal(q_point)[0];
6450 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6451 inline VectorizedArrayType
6453 const unsigned int q_point)
const 6455 return BaseClass::get_laplacian(q_point)[0];
6460 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6466 this->dof_values_initialized =
true;
6469 this->values_dofs[0][dof] = val_in;
6474 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6477 const VectorizedArrayType val_in,
6478 const unsigned int q_point)
6482 Assert(this->J_value !=
nullptr,
6486 this->values_quad_submitted =
true;
6491 const VectorizedArrayType JxW =
6492 this->J_value[0] * this->quadrature_weights[q_point];
6493 this->values_quad[q_point] = val_in * JxW;
6497 this->values_quad[q_point] = val_in * this->J_value[q_point];
6503 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6507 const unsigned int q_point)
6509 submit_value(val_in[0], q_point);
6514 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6518 const unsigned int q_point)
6522 BaseClass::submit_normal_derivative(grad, q_point);
6527 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6531 const unsigned int q_point)
6535 Assert(this->J_value !=
nullptr,
6537 "update_gradients"));
6538 Assert(this->jacobian !=
nullptr,
6540 "update_gradients"));
6542 this->gradients_quad_submitted =
true;
6545 const std::size_t nqp = this->n_quadrature_points;
6548 const VectorizedArrayType JxW =
6549 this->J_value[0] * this->quadrature_weights[q_point];
6550 for (
unsigned int d = 0;
d < dim; ++
d)
6551 this->gradients_quad[
d * nqp + q_point] =
6552 (grad_in[
d] * this->jacobian[0][
d][
d] * JxW);
6559 this->jacobian[q_point] :
6561 const VectorizedArrayType JxW =
6563 this->J_value[q_point] :
6564 this->J_value[0] * this->quadrature_weights[q_point];
6565 for (
unsigned int d = 0;
d < dim; ++
d)
6567 VectorizedArrayType new_val = jac[0][
d] * grad_in[0];
6568 for (
unsigned int e = 1;
e < dim; ++
e)
6569 new_val += jac[
e][
d] * grad_in[
e];
6570 this->gradients_quad[
d * nqp + q_point] = new_val * JxW;
6577 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6578 inline VectorizedArrayType
6582 return BaseClass::integrate_value()[0];
6590 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6594 const unsigned int dof_no,
6595 const unsigned int first_selected_component,
6596 const unsigned int quad_no_in,
6597 const unsigned int fe_degree,
6598 const unsigned int n_q_points,
6599 const bool is_interior_face,
6600 const unsigned int active_fe_index,
6601 const unsigned int active_quad_index,
6602 const unsigned int face_type)
6606 first_selected_component,
6618 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6625 const unsigned int first_selected_component,
6633 first_selected_component,
6639 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6649 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6662 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6667 return BaseClass::get_gradient(q_point);
6672 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6678 Assert(this->gradients_quad_initialized ==
true,
6682 Assert(this->jacobian !=
nullptr,
6684 "update_gradients"));
6686 VectorizedArrayType divergence;
6687 const std::size_t nqp = this->n_quadrature_points;
6692 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6693 for (
unsigned int d = 1;
d < dim; ++
d)
6694 divergence += this->gradients_quad[(dim *
d +
d) * nqp + q_point] *
6695 this->jacobian[0][
d][
d];
6702 this->jacobian[q_point] :
6704 divergence = jac[0][0] * this->gradients_quad[q_point];
6705 for (
unsigned int e = 1;
e < dim; ++
e)
6706 divergence += jac[0][
e] * this->gradients_quad[
e * nqp + q_point];
6707 for (
unsigned int d = 1;
d < dim; ++
d)
6708 for (
unsigned int e = 0;
e < dim; ++
e)
6710 jac[
d][
e] * this->gradients_quad[(
d * dim +
e) * nqp + q_point];
6717 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6723 const auto grad = get_gradient(q_point);
6724 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6725 VectorizedArrayType half = Number(0.5);
6726 for (
unsigned int d = 0;
d < dim; ++
d)
6727 symmetrized[
d] = grad[
d][
d];
6733 symmetrized[2] = grad[0][1] + grad[1][0];
6734 symmetrized[2] *= half;
6737 symmetrized[3] = grad[0][1] + grad[1][0];
6738 symmetrized[3] *= half;
6739 symmetrized[4] = grad[0][2] + grad[2][0];
6740 symmetrized[4] *= half;
6741 symmetrized[5] = grad[1][2] + grad[2][1];
6742 symmetrized[5] *= half;
6752 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6756 const unsigned int q_point)
const 6766 "Computing the curl in 1d is not a useful operation"));
6769 curl[0] = grad[1][0] - grad[0][1];
6772 curl[0] = grad[2][1] - grad[1][2];
6773 curl[1] = grad[0][2] - grad[2][0];
6774 curl[2] = grad[1][0] - grad[0][1];
6784 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6789 return BaseClass::get_hessian_diagonal(q_point);
6794 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6797 const unsigned int q_point)
const 6800 Assert(this->hessians_quad_initialized ==
true,
6804 return BaseClass::get_hessian(q_point);
6809 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6813 const unsigned int q_point)
6815 BaseClass::submit_gradient(grad_in, q_point);
6820 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6825 const unsigned int q_point)
6827 BaseClass::submit_gradient(grad_in, q_point);
6832 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6836 const unsigned int q_point)
6840 Assert(this->J_value !=
nullptr,
6842 "update_gradients"));
6843 Assert(this->jacobian !=
nullptr,
6845 "update_gradients"));
6847 this->gradients_quad_submitted =
true;
6850 const std::size_t nqp = this->n_quadrature_points;
6853 const VectorizedArrayType fac =
6854 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6855 for (
unsigned int d = 0;
d < dim; ++
d)
6857 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6858 (fac * this->jacobian[0][
d][
d]);
6859 for (
unsigned int e = d + 1;
e < dim; ++
e)
6861 this->gradients_quad[(d * dim +
e) * nqp + q_point] =
6862 VectorizedArrayType();
6863 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6864 VectorizedArrayType();
6872 this->jacobian[q_point] :
6874 const VectorizedArrayType fac =
6876 this->J_value[q_point] :
6877 this->J_value[0] * this->quadrature_weights[q_point]) *
6879 for (
unsigned int d = 0;
d < dim; ++
d)
6881 for (
unsigned int e = 0;
e < dim; ++
e)
6882 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6890 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6895 const unsigned int q_point)
6902 Assert(this->J_value !=
nullptr,
6904 "update_gradients"));
6905 Assert(this->jacobian !=
nullptr,
6907 "update_gradients"));
6909 this->gradients_quad_submitted =
true;
6912 const std::size_t nqp = this->n_quadrature_points;
6915 const VectorizedArrayType JxW =
6916 this->J_value[0] * this->quadrature_weights[q_point];
6917 for (
unsigned int d = 0;
d < dim; ++
d)
6918 this->gradients_quad[(
d * dim +
d) * nqp + q_point] =
6920 for (
unsigned int e = 0, counter = dim;
e < dim; ++
e)
6921 for (
unsigned int d =
e + 1;
d < dim; ++
d, ++counter)
6923 const VectorizedArrayType
value =
6925 this->gradients_quad[(
e * dim +
d) * nqp + q_point] =
6926 value * this->jacobian[0][
d][
d];
6927 this->gradients_quad[(
d * dim +
e) * nqp + q_point] =
6928 value * this->jacobian[0][
e][
e];
6934 const VectorizedArrayType JxW =
6936 this->J_value[q_point] :
6937 this->J_value[0] * this->quadrature_weights[q_point];
6940 this->jacobian[q_point] :
6942 VectorizedArrayType weighted[dim][dim];
6943 for (
unsigned int i = 0; i < dim; ++i)
6945 for (
unsigned int i = 0, counter = dim; i < dim; ++i)
6946 for (
unsigned int j = i + 1; j < dim; ++j, ++counter)
6948 const VectorizedArrayType
value =
6950 weighted[i][j] =
value;
6951 weighted[j][i] =
value;
6953 for (
unsigned int comp = 0; comp < dim; ++comp)
6954 for (
unsigned int d = 0;
d < dim; ++
d)
6956 VectorizedArrayType new_val = jac[0][
d] * weighted[comp][0];
6957 for (
unsigned int e = 1;
e < dim; ++
e)
6958 new_val += jac[
e][
d] * weighted[comp][
e];
6959 this->gradients_quad[(comp * dim +
d) * nqp + q_point] = new_val;
6966 template <
int dim,
typename Number,
bool is_face,
typename VectorizedArrayType>
6970 const unsigned int q_point)
6978 "Testing by the curl in 1d is not a useful operation"));
6981 grad[1][0] = curl[0];
6982 grad[0][1] = -curl[0];
6985 grad[2][1] = curl[0];
6986 grad[1][2] = -curl[0];
6987 grad[0][2] = curl[1];
6988 grad[2][0] = -curl[1];
6989 grad[1][0] = curl[2];
6990 grad[0][1] = -curl[2];
6995 submit_gradient(grad, q_point);
7002 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7005 const unsigned int dof_no,
7006 const unsigned int first_selected_component,
7007 const unsigned int quad_no_in,
7008 const unsigned int fe_degree,
7009 const unsigned int n_q_points,
7010 const bool is_interior_face,
7011 const unsigned int active_fe_index,
7012 const unsigned int active_quad_index,
7013 const unsigned int face_type)
7017 first_selected_component,
7029 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7036 const unsigned int first_selected_component,
7043 first_selected_component,
7049 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7058 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7070 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7073 const unsigned int dof)
const 7076 return this->values_dofs[0][dof];
7081 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7084 const unsigned int q_point)
const 7087 Assert(this->values_quad_initialized ==
true,
7091 return this->values_quad[q_point];
7096 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7099 const unsigned int q_point)
const 7105 Assert(this->gradients_quad_initialized ==
true,
7112 this->jacobian[q_point] :
7116 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
7123 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7128 return BaseClass::get_normal_derivative(q_point)[0];
7133 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7136 const unsigned int q_point)
const 7138 return BaseClass::get_hessian(q_point)[0];
7143 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7148 return BaseClass::get_hessian_diagonal(q_point)[0];
7153 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7156 const unsigned int q_point)
const 7158 return BaseClass::get_laplacian(q_point)[0];
7163 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7169 this->dof_values_initialized =
true;
7172 this->values_dofs[0][dof] = val_in;
7177 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7180 const VectorizedArrayType val_in,
7181 const unsigned int q_point)
7186 this->values_quad_submitted =
true;
7191 const VectorizedArrayType JxW = this->J_value[q_point];
7192 this->values_quad[q_point] = val_in * JxW;
7196 const VectorizedArrayType JxW =
7197 this->J_value[0] * this->quadrature_weights[q_point];
7198 this->values_quad[q_point] = val_in * JxW;
7204 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7208 const unsigned int q_point)
7210 submit_value(val_in[0], q_point);
7215 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7219 const unsigned int q_point)
7221 submit_gradient(grad_in[0], q_point);
7226 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7229 const VectorizedArrayType grad_in,
7230 const unsigned int q_point)
7235 this->gradients_quad_submitted =
true;
7240 this->jacobian[q_point] :
7242 const VectorizedArrayType JxW =
7244 this->J_value[q_point] :
7245 this->J_value[0] * this->quadrature_weights[q_point];
7247 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7252 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7256 const unsigned int q_point)
7260 BaseClass::submit_normal_derivative(grad, q_point);
7265 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7269 const unsigned int q_point)
7271 BaseClass::submit_normal_derivative(grad_in, q_point);
7276 template <
typename Number,
bool is_face,
typename VectorizedArrayType>
7277 inline VectorizedArrayType
7281 return BaseClass::integrate_value()[0];
7294 typename VectorizedArrayType>
7300 VectorizedArrayType>::
7302 const unsigned int fe_no,
7303 const unsigned int quad_no,
7304 const unsigned int first_selected_component,
7305 const unsigned int active_fe_index,
7306 const unsigned int active_quad_index)
7307 : BaseClass(data_in,
7309 first_selected_component,
7316 , dofs_per_component(this->data->dofs_per_component_on_cell)
7317 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7318 , n_q_points(this->data->n_q_points)
7320 check_template_arguments(fe_no, 0);
7330 typename VectorizedArrayType>
7336 VectorizedArrayType>::
7338 const std::pair<unsigned int, unsigned int> & range,
7339 const unsigned int dof_no,
7340 const unsigned int quad_no,
7341 const unsigned int first_selected_component)
7345 first_selected_component,
7356 typename VectorizedArrayType>
7362 VectorizedArrayType>::
7367 const unsigned int first_selected_component)
7368 : BaseClass(mapping,
7372 first_selected_component,
7374 , dofs_per_component(this->data->dofs_per_component_on_cell)
7375 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7376 , n_q_points(this->data->n_q_points)
7388 typename VectorizedArrayType>
7394 VectorizedArrayType>::
7398 const unsigned int first_selected_component)
7403 first_selected_component,
7405 , dofs_per_component(this->data->dofs_per_component_on_cell)
7406 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7407 , n_q_points(this->data->n_q_points)
7419 typename VectorizedArrayType>
7425 VectorizedArrayType>::
7429 const unsigned int first_selected_component)
7434 first_selected_component,
7436 , dofs_per_component(this->data->dofs_per_component_on_cell)
7437 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7438 , n_q_points(this->data->n_q_points)
7450 typename VectorizedArrayType>
7459 , dofs_per_component(this->data->dofs_per_component_on_cell)
7460 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7461 , n_q_points(this->data->n_q_points)
7473 typename VectorizedArrayType>
7479 VectorizedArrayType> &
7485 VectorizedArrayType>::operator=(
const FEEvaluation &other)
7487 BaseClass::operator=(other);
7499 typename VectorizedArrayType>
7506 VectorizedArrayType>::
7507 check_template_arguments(
const unsigned int dof_no,
7508 const unsigned int first_selected_component)
7511 (void)first_selected_component;
7517 static_cast<unsigned int>(fe_degree) !=
7518 this->data->data.front().fe_degree) ||
7519 n_q_points != this->n_quadrature_points)
7521 std::string message =
7522 "-------------------------------------------------------\n";
7523 message +=
"Illegal arguments in constructor/wrong template arguments!\n";
7524 message +=
" Called --> FEEvaluation<dim,";
7528 message +=
",Number>(data";
7544 if (static_cast<unsigned int>(fe_degree) ==
7545 this->data->data.front().fe_degree)
7547 proposed_dof_comp = dof_no;
7548 proposed_fe_comp = first_selected_component;
7551 for (
unsigned int no = 0; no < this->matrix_info->n_components();
7553 for (
unsigned int nf = 0;
7554 nf < this->matrix_info->n_base_elements(no);
7556 if (this->matrix_info
7557 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7559 .fe_degree ==
static_cast<unsigned int>(fe_degree))
7561 proposed_dof_comp = no;
7562 proposed_fe_comp = nf;
7566 this->mapping_data->
descriptor[this->active_quad_index]
7568 proposed_quad_comp = this->quad_no;
7570 for (
unsigned int no = 0;
7571 no < this->matrix_info->get_mapping_info().cell_data.size();
7573 if (this->matrix_info->get_mapping_info()
7575 .descriptor[this->active_quad_index]
7576 .n_q_points == n_q_points)
7578 proposed_quad_comp = no;
7585 if (proposed_dof_comp != first_selected_component)
7586 message +=
"Wrong vector component selection:\n";
7588 message +=
"Wrong quadrature formula selection:\n";
7589 message +=
" Did you mean FEEvaluation<dim,";
7593 message +=
",Number>(data";
7602 std::string correct_pos;
7603 if (proposed_dof_comp != dof_no)
7604 correct_pos =
" ^ ";
7607 if (proposed_quad_comp != this->quad_no)
7608 correct_pos +=
" ^ ";
7611 if (proposed_fe_comp != first_selected_component)
7612 correct_pos +=
" ^\n";
7614 correct_pos +=
" \n";
7620 const unsigned int proposed_n_q_points_1d =
static_cast<unsigned int>(
7621 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7622 message +=
"Wrong template arguments:\n";
7623 message +=
" Did you mean FEEvaluation<dim,";
7628 message +=
",Number>(data";
7636 std::string correct_pos;
7637 if (this->data->data.front().fe_degree !=
7638 static_cast<unsigned int>(fe_degree))
7642 if (proposed_n_q_points_1d != n_q_points_1d)
7643 correct_pos +=
" ^\n";
7645 correct_pos +=
" \n";
7646 message +=
" " + correct_pos;
7648 Assert(static_cast<unsigned int>(fe_degree) ==
7649 this->data->data.front().fe_degree &&
7650 n_q_points == this->n_quadrature_points,
7656 this->mapping_data->
descriptor[this->active_quad_index].n_q_points);
7667 typename VectorizedArrayType>
7674 VectorizedArrayType>
::reinit(
const unsigned int cell_index)
7676 Assert(this->mapped_geometry ==
nullptr,
7677 ExcMessage(
"FEEvaluation was initialized without a matrix-free object." 7678 " Integer indexing is not possible"));
7679 if (this->mapped_geometry !=
nullptr)
7684 this->cell = cell_index;
7686 this->matrix_info->get_mapping_info().get_cell_type(cell_index);
7688 const unsigned int offsets =
7689 this->mapping_data->data_index_offsets[cell_index];
7690 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7691 this->J_value = &this->mapping_data->JxW_values[offsets];
7694 this->dof_values_initialized =
false;
7695 this->values_quad_initialized =
false;
7696 this->gradients_quad_initialized =
false;
7697 this->hessians_quad_initialized =
false;
7708 typename VectorizedArrayType>
7709 template <
bool level_dof_access>
7716 VectorizedArrayType>
:: 7719 Assert(this->matrix_info ==
nullptr,
7720 ExcMessage(
"Cannot use initialization from cell iterator if " 7721 "initialized from MatrixFree object. Use variant for " 7722 "on the fly computation with arguments as for FEValues " 7725 this->mapped_geometry->reinit(
7727 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7728 if (level_dof_access)
7729 cell->get_mg_dof_indices(this->local_dof_indices);
7731 cell->get_dof_indices(this->local_dof_indices);
7741 typename VectorizedArrayType>
7748 VectorizedArrayType>
:: 7751 Assert(this->matrix_info == 0,
7752 ExcMessage(
"Cannot use initialization from cell iterator if " 7753 "initialized from MatrixFree object. Use variant for " 7754 "on the fly computation with arguments as for FEValues " 7757 this->mapped_geometry->reinit(cell);
7767 typename VectorizedArrayType>
7774 VectorizedArrayType>::quadrature_point(
const unsigned int q)
const 7776 if (this->matrix_info ==
nullptr)
7778 Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7781 "update_quadrature_points"));
7785 Assert(this->mapping_data->quadrature_point_offsets.empty() ==
false,
7787 "update_quadrature_points"));
7793 &this->mapping_data->quadrature_points
7794 [this->mapping_data->quadrature_point_offsets[this->cell]];
7806 for (
unsigned int d = 0;
d < dim; ++
d)
7807 point[
d] += jac[
d][
d] * static_cast<Number>(
7808 this->descriptor->quadrature.point(q)[
d]);
7810 for (
unsigned int d = 0;
d < dim; ++
d)
7811 for (
unsigned int e = 0;
e < dim; ++
e)
7812 point[
d] += jac[
d][
e] * static_cast<Number>(
7813 this->descriptor->quadrature.point(q)[
e]);
7817 return quadrature_points[q];
7827 typename VectorizedArrayType>
7834 VectorizedArrayType>::evaluate(
const bool evaluate_values,
7835 const bool evaluate_gradients,
7836 const bool evaluate_hessians)
7839 Assert(this->dof_values_initialized ==
true,
7842 evaluate(this->values_dofs[0],
7854 typename VectorizedArrayType>
7861 VectorizedArrayType>::
7865 Assert(this->dof_values_initialized ==
true,
7868 evaluate(this->values_dofs[0], evaluation_flags);
7878 typename VectorizedArrayType>
7885 VectorizedArrayType>::evaluate(
const VectorizedArrayType
7887 const bool evaluate_values,
7888 const bool evaluate_gradients,
7889 const bool evaluate_hessians)
7894 EvaluationFlags::nothing) |
7896 EvaluationFlags::nothing);
7898 evaluate(values_array, flag);
7908 typename VectorizedArrayType>
7915 VectorizedArrayType>::
7916 evaluate(
const VectorizedArrayType * values_array,
7924 const_cast<VectorizedArrayType *>(values_array),
7926 this->gradients_quad,
7927 this->hessians_quad,
7928 this->scratch_data);
7934 const_cast<VectorizedArrayType *>(values_array),
7936 this->gradients_quad,
7937 this->hessians_quad,
7938 this->scratch_data);
7942 this->values_quad_initialized =
true;
7944 this->gradients_quad_initialized =
true;
7946 this->hessians_quad_initialized =
true;
7957 typename VectorizedArrayType>
7958 template <
typename VectorType>
7966 VectorizedArrayType>::gather_evaluate(
const VectorType &input_vector,
7967 const bool evaluate_values,
7968 const bool evaluate_gradients,
7969 const bool evaluate_hessians)
7974 EvaluationFlags::nothing) |
7976 EvaluationFlags::nothing);
7978 gather_evaluate(input_vector, flag);
7987 template <
typename Number,
7988 typename VectorizedArrayType,
7991 typename std::enable_if<
7993 std::is_same<decltype(std::declval<VectorType>().
begin()),
7995 VectorType>::type * =
nullptr>
7997 try_gather_evaluate_inplace(
7999 const VectorType & input_vector,
8000 const unsigned int cell,
8001 const unsigned int active_fe_index,
8002 const unsigned int first_selected_component,
8010 if (std::is_same<typename VectorType::value_type, Number>::value &&
8014 interleaved_contiguous &&
8015 reinterpret_cast<std::size_t>(
8016 input_vector.begin() +
8019 [cell * VectorizedArrayType::size()]) %
8020 sizeof(VectorizedArrayType) ==
8023 const VectorizedArrayType *vec_values =
8024 reinterpret_cast<const VectorizedArrayType *
>(
8025 input_vector.begin() +
8028 [cell * VectorizedArrayType::size()] +
8030 [first_selected_component] *
8031 VectorizedArrayType::size());
8033 phi->evaluate(vec_values, evaluation_flag);
8044 template <
typename Number,
8045 typename VectorizedArrayType,
8046 typename VectorType,
8048 typename std::enable_if<
8050 !std::is_same<decltype(std::declval<VectorType>().
begin()),
8052 VectorType>::type * =
nullptr>
8054 try_gather_evaluate_inplace(T,
8072 typename VectorizedArrayType,
8073 typename VectorType,
8074 typename std::enable_if<
8076 std::is_same<decltype(std::declval<VectorType>().
begin()),
8078 VectorType>::type * =
nullptr>
8080 try_integrate_scatter_inplace(
8081 VectorType & destination,
8082 const unsigned int cell,
8083 const unsigned int n_components,
8084 const unsigned int active_fe_index,
8085 const unsigned int first_selected_component,
8087 VectorizedArrayType * values_quad,
8088 VectorizedArrayType * gradients_quad,
8089 VectorizedArrayType * scratch_data,
8098 if (std::is_same<typename VectorType::value_type, Number>::value &&
8102 interleaved_contiguous &&
8103 reinterpret_cast<std::size_t>(
8104 destination.begin() +
8107 [cell * VectorizedArrayType::size()]) %
8108 sizeof(VectorizedArrayType) ==
8111 VectorizedArrayType *vec_values =
8112 reinterpret_cast<VectorizedArrayType *
>(
8113 destination.begin() +
8116 [cell * VectorizedArrayType::size()] +
8118 [first_selected_component] *
8119 VectorizedArrayType::size());
8154 typename VectorizedArrayType,
8155 typename VectorType,
8156 typename std::enable_if<
8158 !std::is_same<decltype(std::declval<VectorType>().
begin()),
8160 VectorType>::type * =
nullptr>
8162 try_integrate_scatter_inplace(
8169 const VectorizedArrayType *,
8170 const VectorizedArrayType *,
8171 const VectorizedArrayType *,
8186 typename VectorizedArrayType>
8187 template <
typename VectorType>
8194 VectorizedArrayType>::
8195 gather_evaluate(
const VectorType & input_vector,
8198 if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
8202 this->active_fe_index,
8203 this->first_selected_component,
8205 evaluation_flag) ==
false)
8207 this->read_dof_values(input_vector);
8208 evaluate(this->begin_dof_values(), evaluation_flag);
8219 typename VectorizedArrayType>
8226 VectorizedArrayType>::integrate(
const bool integrate_values,
8227 const bool integrate_gradients)
8229 integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8232 this->dof_values_initialized =
true;
8243 typename VectorizedArrayType>
8250 VectorizedArrayType>::
8253 integrate(integration_flag, this->values_dofs[0]);
8256 this->dof_values_initialized =
true;
8267 typename VectorizedArrayType>
8274 VectorizedArrayType>::integrate(
const bool integrate_values,
8275 const bool integrate_gradients,
8276 VectorizedArrayType *values_array)
8282 integrate(flag, values_array);
8292 typename VectorizedArrayType>
8299 VectorizedArrayType>::
8301 VectorizedArrayType * values_array)
8305 Assert(this->values_quad_submitted ==
true,
8308 Assert(this->gradients_quad_submitted ==
true,
8311 Assert(this->matrix_info !=
nullptr ||
8312 this->mapped_geometry->is_initialized(),
8317 ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
8319 "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8328 this->gradients_quad,
8338 this->gradients_quad,
8343 this->dof_values_initialized =
true;
8354 typename VectorizedArrayType>
8355 template <
typename VectorType>