Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_evaluation.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_fe_evaluation_h
18#define dealii_matrix_free_fe_evaluation_h
19
20
21#include <deal.II/base/config.h>
22
29
31
42
43
45
46
47
48namespace internal
49{
51
54 std::string,
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.");
62} // namespace internal
63
64template <int dim,
65 int fe_degree,
66 int n_q_points_1d = fe_degree + 1,
67 int n_components_ = 1,
68 typename Number = double,
69 typename VectorizedArrayType = VectorizedArray<Number>>
70class FEEvaluation;
71
72
73
100template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
102{
103 static_assert(
104 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
105 "Type of Number and of VectorizedArrayType do not match.");
106
107public:
108 static constexpr unsigned int dimension = dim;
109
114
123 unsigned int
125
134
140
146
151 VectorizedArrayType
152 JxW(const unsigned int q_point) const;
153
166 inverse_jacobian(const unsigned int q_point) const;
167
181 get_normal_vector(const unsigned int q_point) const;
182
189 VectorizedArrayType
191
198 void
200 const VectorizedArrayType & value) const;
201
206 template <typename T>
207 std::array<T, VectorizedArrayType::size()>
208 read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
209 &array) const;
210
215 template <typename T>
216 void
218 AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
219 const std::array<T, VectorizedArrayType::size()> & value) const;
220
225 std::array<unsigned int, VectorizedArrayType::size()>
227
232 std::array<unsigned int, VectorizedArrayType::size()>
234
235
241 const std::vector<unsigned int> &
243
252
256 unsigned int
258
262 unsigned int
264
269 unsigned int
271
276 unsigned int
278
284
285protected:
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);
303
309 const Mapping<dim> & mapping,
310 const FiniteElement<dim> &fe,
311 const Quadrature<1> & quadrature,
312 const UpdateFlags update_flags,
313 const unsigned int first_selected_component,
315 *other);
316
324
333
338
344 VectorizedArrayType *scratch_data;
345
349 const unsigned int quad_no;
350
355
362
370 (is_face ? dim - 1 : dim),
371 dim,
372 Number,
373 VectorizedArrayType> *mapping_data;
374
378 const unsigned int active_fe_index;
379
384 const unsigned int active_quad_index;
385
392 (is_face ? dim - 1 : dim),
393 dim,
394 Number,
395 VectorizedArrayType>::QuadratureDescriptor *descriptor;
396
400 const unsigned int n_quadrature_points;
401
409
415
422 const VectorizedArrayType *J_value;
423
428
433
437 const Number *quadrature_weights;
438
443 unsigned int cell;
444
450
456
461 unsigned int face_no;
462
467 unsigned int face_orientation;
468
476 unsigned int subface_index;
477
485
490 std::shared_ptr<internal::MatrixFreeFunctions::
491 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
493
494 // Make FEEvaluation objects friends for access to protected member
495 // mapped_geometry.
496 template <int, int, int, int, typename, typename>
497 friend class FEEvaluation;
498};
499
500
501
539template <int dim,
540 int n_components_,
541 typename Number,
542 bool is_face = false,
543 typename VectorizedArrayType = VectorizedArray<Number>>
545 : public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
546{
547public:
548 using number_type = Number;
552 static constexpr unsigned int dimension = dim;
553 static constexpr unsigned int n_components = n_components_;
554
591 template <typename VectorType>
592 void
593 read_dof_values(const VectorType &src, const unsigned int first_index = 0);
594
623 template <typename VectorType>
624 void
625 read_dof_values_plain(const VectorType & src,
626 const unsigned int first_index = 0);
627
659 template <typename VectorType>
660 void
662 VectorType & dst,
663 const unsigned int first_index = 0,
664 const std::bitset<VectorizedArrayType::size()> &mask =
665 std::bitset<VectorizedArrayType::size()>().flip()) const;
666
705 template <typename VectorType>
706 void
707 set_dof_values(VectorType & dst,
708 const unsigned int first_index = 0,
709 const std::bitset<VectorizedArrayType::size()> &mask =
710 std::bitset<VectorizedArrayType::size()>().flip()) const;
711
715 template <typename VectorType>
716 void
718 VectorType & dst,
719 const unsigned int first_index = 0,
720 const std::bitset<VectorizedArrayType::size()> &mask =
721 std::bitset<VectorizedArrayType::size()>().flip()) const;
722
724
746 get_dof_value(const unsigned int dof) const;
747
758 void
759 submit_dof_value(const value_type val_in, const unsigned int dof);
760
774 get_value(const unsigned int q_point) const;
775
788 void
789 submit_value(const value_type val_in, const unsigned int q_point);
790
802 get_gradient(const unsigned int q_point) const;
803
819 get_normal_derivative(const unsigned int q_point) const;
820
833 void
834 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
835
854 void
856 const unsigned int q_point);
857
870 get_hessian(const unsigned int q_point) const;
871
882 get_hessian_diagonal(const unsigned int q_point) const;
883
896 get_laplacian(const unsigned int q_point) const;
897
898#ifdef DOXYGEN
899 // doxygen does not anyhow mention functions coming from partial template
900 // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
901 // For now, hack in those functions manually only to fix documentation:
902
909 VectorizedArrayType
910 get_divergence(const unsigned int q_point) const;
911
921 get_symmetric_gradient(const unsigned int q_point) const;
922
929 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
930 get_curl(const unsigned int q_point) const;
931
947 void
948 submit_divergence(const VectorizedArrayType div_in,
949 const unsigned int q_point);
950
967 void
970 const unsigned int q_point);
971
984 void
986 const unsigned int q_point);
987
988#endif
989
1008
1010
1023 const VectorizedArrayType *
1025
1034 VectorizedArrayType *
1036
1047 const VectorizedArrayType *
1049
1060 VectorizedArrayType *
1062
1074 const VectorizedArrayType *
1076
1088 VectorizedArrayType *
1090
1103 const VectorizedArrayType *
1105
1118 VectorizedArrayType *
1120
1122
1126 unsigned int
1128
1129protected:
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);
1149
1186 const Mapping<dim> & mapping,
1187 const FiniteElement<dim> &fe,
1188 const Quadrature<1> & quadrature,
1189 const UpdateFlags update_flags,
1190 const unsigned int first_selected_component,
1192 *other);
1193
1201
1210
1217 template <typename VectorType, typename VectorOperation>
1218 void
1220 const VectorOperation & operation,
1221 const std::array<VectorType *, n_components_> &vectors,
1222 const std::array<
1224 n_components_> & vectors_sm,
1225 const std::bitset<VectorizedArrayType::size()> &mask,
1226 const bool apply_constraints = true) const;
1227
1235 template <typename VectorType, typename VectorOperation>
1236 void
1238 const VectorOperation & operation,
1239 const std::array<VectorType *, n_components_> &vectors,
1240 const std::array<
1242 n_components_> & vectors_sm,
1243 const std::bitset<VectorizedArrayType::size()> &mask) const;
1244
1252 template <typename VectorType, typename VectorOperation>
1253 void
1255 const VectorOperation & operation,
1256 const std::array<VectorType *, n_components_> &vectors) const;
1257
1270 VectorizedArrayType *values_dofs[n_components];
1271
1283 VectorizedArrayType *values_quad;
1284
1298 VectorizedArrayType *gradients_quad;
1299
1311 VectorizedArrayType *hessians_quad;
1312
1317 const unsigned int n_fe_components;
1318
1325
1332
1339
1346
1353
1360
1365 const unsigned int first_selected_component;
1366
1371 mutable std::vector<types::global_dof_index> local_dof_indices;
1372
1373private:
1378 void
1380};
1381
1382
1383
1391template <int dim,
1392 int n_components_,
1393 typename Number,
1394 bool is_face,
1395 typename VectorizedArrayType = VectorizedArray<Number>>
1397 n_components_,
1398 Number,
1399 is_face,
1400 VectorizedArrayType>
1401{
1402 static_assert(
1403 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1404 "Type of Number and of VectorizedArrayType do not match.");
1405
1406public:
1407 using number_type = Number;
1411 static constexpr unsigned int dimension = dim;
1412 static constexpr unsigned int n_components = n_components_;
1415
1416protected:
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,
1434 const unsigned int face_type = numbers::invalid_unsigned_int);
1435
1441 const Mapping<dim> & mapping,
1442 const FiniteElement<dim> &fe,
1443 const Quadrature<1> & quadrature,
1444 const UpdateFlags update_flags,
1445 const unsigned int first_selected_component,
1447 *other);
1448
1453
1459};
1460
1461
1462
1471template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1472class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1473 : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1474{
1475 static_assert(
1476 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1477 "Type of Number and of VectorizedArrayType do not match.");
1478
1479public:
1480 using number_type = Number;
1481 using value_type = VectorizedArrayType;
1483 static constexpr unsigned int dimension = dim;
1486
1491 get_dof_value(const unsigned int dof) const;
1492
1496 void
1497 submit_dof_value(const value_type val_in, const unsigned int dof);
1498
1503 get_value(const unsigned int q_point) const;
1504
1508 void
1509 submit_value(const value_type val_in, const unsigned int q_point);
1510
1514 void
1516 const unsigned int q_point);
1517
1522 get_gradient(const unsigned int q_point) const;
1523
1528 get_normal_derivative(const unsigned int q_point) const;
1529
1533 void
1534 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1535
1539 void
1541 const unsigned int q_point);
1542
1547 get_hessian(unsigned int q_point) const;
1548
1553 get_hessian_diagonal(const unsigned int q_point) const;
1554
1559 get_laplacian(const unsigned int q_point) const;
1560
1566
1567protected:
1577 const unsigned int dof_no,
1578 const unsigned int first_selected_component,
1579 const unsigned int quad_no,
1580 const unsigned int fe_degree,
1581 const unsigned int n_q_points,
1582 const bool is_interior_face = true,
1585 const unsigned int face_type = numbers::invalid_unsigned_int);
1586
1592 const Mapping<dim> & mapping,
1593 const FiniteElement<dim> &fe,
1594 const Quadrature<1> & quadrature,
1595 const UpdateFlags update_flags,
1596 const unsigned int first_selected_component,
1598 *other);
1599
1604
1610};
1611
1612
1613
1623template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1624class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1625 : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1626{
1627 static_assert(
1628 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1629 "Type of Number and of VectorizedArrayType do not match.");
1630
1631public:
1632 using number_type = Number;
1635 static constexpr unsigned int dimension = dim;
1636 static constexpr unsigned int n_components = dim;
1639
1644 get_gradient(const unsigned int q_point) const;
1645
1650 VectorizedArrayType
1651 get_divergence(const unsigned int q_point) const;
1652
1660 get_symmetric_gradient(const unsigned int q_point) const;
1661
1666 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1667 get_curl(const unsigned int q_point) const;
1668
1673 get_hessian(const unsigned int q_point) const;
1674
1679 get_hessian_diagonal(const unsigned int q_point) const;
1680
1684 void
1685 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1686
1695 void
1697 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1698 const unsigned int q_point);
1699
1708 void
1709 submit_divergence(const VectorizedArrayType div_in,
1710 const unsigned int q_point);
1711
1720 void
1723 const unsigned int q_point);
1724
1729 void
1731 const unsigned int q_point);
1732
1733protected:
1743 const unsigned int dof_no,
1744 const unsigned int first_selected_component,
1745 const unsigned int quad_no,
1746 const unsigned int dofs_per_cell,
1747 const unsigned int n_q_points,
1748 const bool is_interior_face = true,
1751 const unsigned int face_type = numbers::invalid_unsigned_int);
1752
1758 const Mapping<dim> & mapping,
1759 const FiniteElement<dim> &fe,
1760 const Quadrature<1> & quadrature,
1761 const UpdateFlags update_flags,
1762 const unsigned int first_selected_component,
1764 *other);
1765
1770
1776};
1777
1778
1787template <typename Number, bool is_face, typename VectorizedArrayType>
1788class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1789 : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1790{
1791 static_assert(
1792 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1793 "Type of Number and of VectorizedArrayType do not match.");
1794
1795public:
1796 using number_type = Number;
1797 using value_type = VectorizedArrayType;
1799 static constexpr unsigned int dimension = 1;
1802
1807 get_dof_value(const unsigned int dof) const;
1808
1812 void
1813 submit_dof_value(const value_type val_in, const unsigned int dof);
1814
1819 get_value(const unsigned int q_point) const;
1820
1824 void
1825 submit_value(const value_type val_in, const unsigned int q_point);
1826
1830 void
1831 submit_value(const gradient_type val_in, const unsigned int q_point);
1832
1837 get_gradient(const unsigned int q_point) const;
1838
1843 get_divergence(const unsigned int q_point) const;
1844
1849 get_normal_derivative(const unsigned int q_point) const;
1850
1854 void
1855 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1856
1860 void
1861 submit_gradient(const value_type grad_in, const unsigned int q_point);
1862
1866 void
1868 const unsigned int q_point);
1869
1873 void
1875 const unsigned int q_point);
1876
1880 void
1882 const unsigned int q_point);
1883
1888 get_hessian(unsigned int q_point) const;
1889
1894 get_hessian_diagonal(const unsigned int q_point) const;
1895
1900 get_laplacian(const unsigned int q_point) const;
1901
1907
1908protected:
1918 const unsigned int dof_no,
1919 const unsigned int first_selected_component,
1920 const unsigned int quad_no,
1921 const unsigned int fe_degree,
1922 const unsigned int n_q_points,
1923 const bool is_interior_face = true,
1926 const unsigned int face_type = numbers::invalid_unsigned_int);
1927
1933 const Mapping<1> & mapping,
1934 const FiniteElement<1> &fe,
1935 const Quadrature<1> & quadrature,
1936 const UpdateFlags update_flags,
1937 const unsigned int first_selected_component,
1939
1944
1950};
1951
1952
1953
2508template <int dim,
2509 int fe_degree,
2510 int n_q_points_1d,
2511 int n_components_,
2512 typename Number,
2513 typename VectorizedArrayType>
2515 n_components_,
2516 Number,
2517 false,
2518 VectorizedArrayType>
2519{
2520 static_assert(
2521 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2522 "Type of Number and of VectorizedArrayType do not match.");
2523
2524public:
2530
2534 using number_type = Number;
2535
2542
2549
2553 static constexpr unsigned int dimension = dim;
2554
2559 static constexpr unsigned int n_components = n_components_;
2560
2567 static constexpr unsigned int static_n_q_points =
2568 Utilities::pow(n_q_points_1d, dim);
2569
2577 static constexpr unsigned int static_dofs_per_component =
2578 Utilities::pow(fe_degree + 1, dim);
2579
2587 static constexpr unsigned int tensor_dofs_per_cell =
2589
2597 static constexpr unsigned int static_dofs_per_cell =
2599
2636 const unsigned int dof_no = 0,
2637 const unsigned int quad_no = 0,
2638 const unsigned int first_selected_component = 0,
2641
2650 const std::pair<unsigned int, unsigned int> & range,
2651 const unsigned int dof_no = 0,
2652 const unsigned int quad_no = 0,
2653 const unsigned int first_selected_component = 0);
2654
2682 const FiniteElement<dim> &fe,
2683 const Quadrature<1> & quadrature,
2684 const UpdateFlags update_flags,
2685 const unsigned int first_selected_component = 0);
2686
2693 const Quadrature<1> & quadrature,
2694 const UpdateFlags update_flags,
2695 const unsigned int first_selected_component = 0);
2696
2708 const FiniteElement<dim> & fe,
2710 const unsigned int first_selected_component = 0);
2711
2719
2726 FEEvaluation &
2728
2737 void
2738 reinit(const unsigned int cell_batch_index);
2739
2752 template <bool level_dof_access>
2753 void
2755
2766 void
2768
2772 static bool
2773 fast_evaluation_supported(const unsigned int given_degree,
2774 const unsigned int give_n_q_points_1d);
2775
2785 void
2787
2792 void
2793 evaluate(const bool evaluate_values,
2794 const bool evaluate_gradients,
2795 const bool evaluate_hessians = false);
2796
2809 void
2810 evaluate(const VectorizedArrayType * values_array,
2811 const EvaluationFlags::EvaluationFlags evaluation_flag);
2812
2817 void
2818 evaluate(const VectorizedArrayType *values_array,
2819 const bool evaluate_values,
2820 const bool evaluate_gradients,
2821 const bool evaluate_hessians = false);
2822
2836 template <typename VectorType>
2837 void
2838 gather_evaluate(const VectorType & input_vector,
2839 const EvaluationFlags::EvaluationFlags evaluation_flag);
2840
2844 template <typename VectorType>
2845 void
2846 gather_evaluate(const VectorType &input_vector,
2847 const bool evaluate_values,
2848 const bool evaluate_gradients,
2849 const bool evaluate_hessians = false);
2850
2861 void
2863
2864
2868 void
2869 integrate(const bool integrate_values, const bool integrate_gradients);
2870
2882 void
2884 VectorizedArrayType * values_array);
2885
2889 void
2890 integrate(const bool integrate_values,
2891 const bool integrate_gradients,
2892 VectorizedArrayType *values_array);
2893
2907 template <typename VectorType>
2908 void
2910 VectorType & output_vector);
2911
2915 template <typename VectorType>
2916 void
2917 integrate_scatter(const bool integrate_values,
2918 const bool integrate_gradients,
2919 VectorType &output_vector);
2920
2926 quadrature_point(const unsigned int q_point) const;
2927
2934 const unsigned int dofs_per_component;
2935
2942 const unsigned int dofs_per_cell;
2943
2951 const unsigned int n_q_points;
2952
2953private:
2958 void
2959 check_template_arguments(const unsigned int fe_no,
2960 const unsigned int first_selected_component);
2961};
2962
2963
2964
3000template <int dim,
3001 int fe_degree,
3002 int n_q_points_1d = fe_degree + 1,
3003 int n_components_ = 1,
3004 typename Number = double,
3005 typename VectorizedArrayType = VectorizedArray<Number>>
3007 n_components_,
3008 Number,
3009 true,
3010 VectorizedArrayType>
3011{
3012 static_assert(
3013 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
3014 "Type of Number and of VectorizedArrayType do not match.");
3015
3016public:
3022
3026 using number_type = Number;
3027
3034
3041
3045 static constexpr unsigned int dimension = dim;
3046
3051 static constexpr unsigned int n_components = n_components_;
3052
3060 static constexpr unsigned int static_n_q_points =
3061 Utilities::pow(n_q_points_1d, dim - 1);
3062
3069 static constexpr unsigned int static_n_q_points_cell =
3070 Utilities::pow(n_q_points_1d, dim);
3071
3078 static constexpr unsigned int static_dofs_per_component =
3079 Utilities::pow(fe_degree + 1, dim);
3080
3087 static constexpr unsigned int tensor_dofs_per_cell =
3089
3096 static constexpr unsigned int static_dofs_per_cell =
3098
3142 const bool is_interior_face = true,
3143 const unsigned int dof_no = 0,
3144 const unsigned int quad_no = 0,
3145 const unsigned int first_selected_component = 0,
3148 const unsigned int face_type = numbers::invalid_unsigned_int);
3149
3159 const std::pair<unsigned int, unsigned int> & range,
3160 const bool is_interior_face = true,
3161 const unsigned int dof_no = 0,
3162 const unsigned int quad_no = 0,
3163 const unsigned int first_selected_component = 0);
3164
3175 void
3176 reinit(const unsigned int face_batch_number);
3177
3185 void
3186 reinit(const unsigned int cell_batch_number, const unsigned int face_number);
3187
3191 static bool
3192 fast_evaluation_supported(const unsigned int given_degree,
3193 const unsigned int give_n_q_points_1d);
3194
3205 void
3207
3211 void
3212 evaluate(const bool evaluate_values, const bool evaluate_gradients);
3213
3226 void
3227 evaluate(const VectorizedArrayType * values_array,
3228 const EvaluationFlags::EvaluationFlags evaluation_flag);
3229
3233 void
3234 evaluate(const VectorizedArrayType *values_array,
3235 const bool evaluate_values,
3236 const bool evaluate_gradients);
3237
3249 template <typename VectorType>
3250 void
3251 gather_evaluate(const VectorType & input_vector,
3252 const EvaluationFlags::EvaluationFlags evaluation_flag);
3253
3257 template <typename VectorType>
3258 void
3259 gather_evaluate(const VectorType &input_vector,
3260 const bool evaluate_values,
3261 const bool evaluate_gradients);
3262
3272 void
3274
3278 void
3279 integrate(const bool integrate_values, const bool integrate_gradients);
3280
3289 void
3291 VectorizedArrayType * values_array);
3292
3296 void
3297 integrate(const bool integrate_values,
3298 const bool integrate_gradients,
3299 VectorizedArrayType *values_array);
3300
3312 template <typename VectorType>
3313 void
3315 VectorType & output_vector);
3316
3320 template <typename VectorType>
3321 void
3322 integrate_scatter(const bool integrate_values,
3323 const bool integrate_gradients,
3324 VectorType &output_vector);
3325
3331 quadrature_point(const unsigned int q_point) const;
3332
3339 const unsigned int dofs_per_component;
3340
3347 const unsigned int dofs_per_cell;
3348
3356 const unsigned int n_q_points;
3357
3358
3359private:
3363 std::array<unsigned int, VectorizedArrayType::size()>
3365
3369 std::array<unsigned int, VectorizedArrayType::size()>
3371};
3372
3373
3374
3375namespace internal
3376{
3377 namespace MatrixFreeFunctions
3378 {
3379 // a helper function to compute the number of DoFs of a DGP element at
3380 // compile time, depending on the degree
3381 template <int dim, int degree>
3383 {
3384 // this division is always without remainder
3385 static constexpr unsigned int value =
3386 (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
3387 };
3388
3389 // base specialization: 1d elements have 'degree+1' degrees of freedom
3390 template <int degree>
3391 struct DGP_dofs_per_component<1, degree>
3392 {
3393 static constexpr unsigned int value = degree + 1;
3394 };
3395 } // namespace MatrixFreeFunctions
3396} // namespace internal
3397
3398
3399/*----------------------- Inline functions ----------------------------------*/
3400
3401#ifndef DOXYGEN
3402
3403
3404/*----------------------- FEEvaluationBaseData ------------------------*/
3405
3406template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3410 const unsigned int dof_no,
3411 const unsigned int first_selected_component,
3412 const unsigned int quad_no_in,
3413 const unsigned int fe_degree,
3414 const unsigned int n_q_points,
3415 const bool is_interior_face,
3416 const unsigned int active_fe_index_in,
3417 const unsigned int active_quad_index_in,
3418 const unsigned int face_type)
3419 : scratch_data_array(data_in.acquire_scratch_data())
3420 , quad_no(quad_no_in)
3421 , matrix_info(&data_in)
3422 , dof_info(&data_in.get_dof_info(dof_no))
3423 , mapping_data(
3424 internal::MatrixFreeFunctions::
3425 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3426 data_in.get_mapping_info(),
3427 quad_no))
3428 , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
3429 data_in.get_dof_info(dof_no).fe_index_from_degree(
3430 first_selected_component,
3431 fe_degree) :
3432 (active_fe_index_in != numbers::invalid_unsigned_int ?
3433 active_fe_index_in :
3434 0))
3435 , active_quad_index(
3436 fe_degree != numbers::invalid_unsigned_int ?
3437 (mapping_data->quad_index_from_n_q_points(n_q_points)) :
3438 (active_quad_index_in != numbers::invalid_unsigned_int ?
3439 active_quad_index_in :
3440 std::min<unsigned int>(active_fe_index,
3441 mapping_data->descriptor.size() - 1)))
3442 , descriptor(
3443 &mapping_data->descriptor
3444 [is_face ?
3445 (active_quad_index * std::max<unsigned int>(1, dim - 1) +
3446 (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
3447 active_quad_index])
3448 , n_quadrature_points(descriptor->n_q_points)
3449 , data(&data_in.get_shape_info(
3450 dof_no,
3451 quad_no_in,
3452 dof_info->component_to_base_index[first_selected_component],
3453 active_fe_index,
3454 active_quad_index))
3455 , jacobian(nullptr)
3456 , J_value(nullptr)
3457 , normal_vectors(nullptr)
3458 , normal_x_jacobian(nullptr)
3459 , quadrature_weights(descriptor->quadrature_weights.begin())
3461 , is_interior_face(is_interior_face)
3462 , dof_access_index(
3463 is_face ?
3464 (is_interior_face ?
3465 internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3466 internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3467 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3468 , cell_type(internal::MatrixFreeFunctions::general)
3469{
3470 Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3471 AssertDimension(matrix_info->get_task_info().vectorization_length,
3472 VectorizedArrayType::size());
3473 AssertDimension(n_quadrature_points, descriptor->n_q_points);
3474}
3475
3476
3477
3478template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3481 const Mapping<dim> & mapping,
3482 const FiniteElement<dim> &fe,
3483 const Quadrature<1> & quadrature,
3484 const UpdateFlags update_flags,
3485 const unsigned int first_selected_component,
3487 *other)
3488 : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3489 , quad_no(numbers::invalid_unsigned_int)
3490 , active_fe_index(numbers::invalid_unsigned_int)
3491 , active_quad_index(numbers::invalid_unsigned_int)
3492 , descriptor(nullptr)
3493 , n_quadrature_points(
3494 Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3495 , matrix_info(nullptr)
3496 , dof_info(nullptr)
3497 , mapping_data(nullptr)
3498 ,
3499 // select the correct base element from the given FE component
3500 data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
3501 Quadrature<dim - is_face>(quadrature),
3502 fe,
3503 fe.component_to_base_index(first_selected_component).first))
3504 , jacobian(nullptr)
3505 , J_value(nullptr)
3506 , normal_vectors(nullptr)
3507 , normal_x_jacobian(nullptr)
3508 , quadrature_weights(nullptr)
3509 , cell(0)
3510 , cell_type(internal::MatrixFreeFunctions::general)
3511 , is_interior_face(true)
3512 , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3513{
3514 Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3516 if (other != nullptr &&
3517 other->mapped_geometry->get_quadrature() == quadrature)
3518 mapped_geometry = other->mapped_geometry;
3519 else
3520 mapped_geometry =
3521 std::make_shared<internal::MatrixFreeFunctions::
3522 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3523 mapping, quadrature, update_flags);
3524 cell = 0;
3525
3526 mapping_data = &mapped_geometry->get_data_storage();
3527 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3528 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3529}
3530
3531
3532
3533template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3537 &other)
3538 : scratch_data_array(other.matrix_info == nullptr ?
3539 new AlignedVector<VectorizedArrayType>() :
3540 other.matrix_info->acquire_scratch_data())
3541 , quad_no(other.quad_no)
3542 , active_fe_index(other.active_fe_index)
3543 , active_quad_index(other.active_quad_index)
3544 , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
3545 , n_quadrature_points(other.n_quadrature_points)
3546 , matrix_info(other.matrix_info)
3547 , dof_info(other.dof_info)
3548 , mapping_data(other.mapping_data)
3549 , data(other.matrix_info == nullptr ?
3550 new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
3551 *other.data) :
3552 other.data)
3553 , jacobian(nullptr)
3554 , J_value(nullptr)
3555 , normal_vectors(nullptr)
3556 , normal_x_jacobian(nullptr)
3557 , quadrature_weights(other.descriptor == nullptr ?
3558 nullptr :
3559 descriptor->quadrature_weights.begin())
3561 , cell_type(internal::MatrixFreeFunctions::general)
3562 , is_interior_face(other.is_interior_face)
3563 , dof_access_index(other.dof_access_index)
3564{
3565 // Create deep copy of mapped geometry for use in parallel...
3566 if (other.mapped_geometry.get() != nullptr)
3567 {
3568 mapped_geometry = std::make_shared<
3569 internal::MatrixFreeFunctions::
3570 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3571 other.mapped_geometry->get_fe_values().get_mapping(),
3572 other.mapped_geometry->get_quadrature(),
3573 other.mapped_geometry->get_fe_values().get_update_flags());
3574 mapping_data = &mapped_geometry->get_data_storage();
3575 cell = 0;
3576
3577 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3578 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3579 }
3580}
3581
3582
3583
3584template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3588{
3589 AssertDimension(quad_no, other.quad_no);
3590 AssertDimension(active_fe_index, other.active_fe_index);
3591 AssertDimension(active_quad_index, other.active_quad_index);
3592
3593 // release old memory
3594 if (matrix_info == nullptr)
3595 {
3596 delete data;
3597 delete scratch_data_array;
3598 }
3599 else
3600 {
3601 matrix_info->release_scratch_data(scratch_data_array);
3602 }
3603
3604 matrix_info = other.matrix_info;
3605 dof_info = other.dof_info;
3606 descriptor = other.descriptor;
3607 mapping_data = other.mapping_data;
3608 if (other.matrix_info == nullptr)
3609 {
3611 *other.data);
3612 scratch_data_array = new AlignedVector<VectorizedArrayType>();
3613 }
3614 else
3615 {
3616 data = other.data;
3617 scratch_data_array = matrix_info->acquire_scratch_data();
3618 }
3619
3620 quadrature_weights =
3621 (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr);
3624 is_interior_face = other.is_interior_face;
3625 dof_access_index = other.dof_access_index;
3626
3627 // Create deep copy of mapped geometry for use in parallel...
3628 if (other.mapped_geometry.get() != nullptr)
3629 {
3630 mapped_geometry = std::make_shared<
3631 internal::MatrixFreeFunctions::
3632 MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3633 other.mapped_geometry->get_fe_values().get_mapping(),
3634 other.mapped_geometry->get_quadrature(),
3635 other.mapped_geometry->get_fe_values().get_update_flags());
3636 cell = 0;
3637 mapping_data = &mapped_geometry->get_data_storage();
3638 jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3639 J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3640 }
3641
3642 return *this;
3643}
3644
3645
3646
3647template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3650{
3651 if (matrix_info != nullptr)
3652 {
3653 try
3654 {
3655 matrix_info->release_scratch_data(scratch_data_array);
3656 }
3657 catch (...)
3658 {}
3659 }
3660 else
3661 {
3662 delete scratch_data_array;
3663 delete data;
3664 data = nullptr;
3665 }
3666 scratch_data_array = nullptr;
3667}
3668
3669
3670
3671template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3672inline unsigned int
3675{
3676 if (matrix_info == nullptr)
3677 return 0;
3678 else
3679 {
3680 AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3681 return this->mapping_data->data_index_offsets[cell];
3682 }
3683}
3684
3685
3686
3687template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3690 const
3691{
3693 return cell_type;
3694}
3695
3696
3697
3698template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3701 get_shape_info() const
3702{
3703 Assert(data != nullptr, ExcInternalError());
3704 return *data;
3705}
3706
3707
3708
3709template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3712 const
3713{
3714 Assert(dof_info != nullptr,
3715 ExcMessage(
3716 "FEEvaluation was not initialized with a MatrixFree object!"));
3717 return *dof_info;
3718}
3719
3720
3721
3722template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3725 get_normal_vector(const unsigned int q_point) const
3726{
3727 AssertIndexRange(q_point, n_quadrature_points);
3728 Assert(normal_vectors != nullptr,
3730 "update_normal_vectors"));
3731 if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3732 return normal_vectors[0];
3733 else
3734 return normal_vectors[q_point];
3735}
3736
3737
3738
3739template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3740inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3742 const unsigned int q_point) const
3743{
3744 AssertIndexRange(q_point, n_quadrature_points);
3745 Assert(J_value != nullptr,
3747 "update_values|update_gradients"));
3748 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3749 {
3750 Assert(this->quadrature_weights != nullptr, ExcInternalError());
3751 return J_value[0] * this->quadrature_weights[q_point];
3752 }
3753 else
3754 return J_value[q_point];
3755}
3756
3757
3758
3759template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3762 inverse_jacobian(const unsigned int q_point) const
3763{
3764 AssertIndexRange(q_point, n_quadrature_points);
3765 Assert(this->jacobian != nullptr,
3767 "update_gradients"));
3768 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3769 return jacobian[0];
3770 else
3771 return jacobian[q_point];
3772}
3773
3774
3775
3776template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3777inline std::array<unsigned int, VectorizedArrayType::size()>
3779 const
3780{
3781 Assert(this->matrix_info != nullptr, ExcNotInitialized());
3782
3783 const unsigned int n_lanes = VectorizedArrayType::size();
3784 std::array<unsigned int, n_lanes> cells;
3785
3786 // initialize array
3787 for (unsigned int i = 0; i < n_lanes; ++i)
3789
3790 if ((is_face == false) ||
3791 (is_face &&
3792 this->dof_access_index ==
3794 this->is_interior_face))
3795 {
3796 // cell or interior face face (element-centric loop)
3797 for (unsigned int i = 0; i < n_lanes; ++i)
3798 cells[i] = cell * n_lanes + i;
3799 }
3800 else if (is_face &&
3801 this->dof_access_index ==
3803 this->is_interior_face == false)
3804 {
3805 // exterior face (element-centric loop): for this case, we need to
3806 // look into the FaceInfo field that collects information from both
3807 // sides of a face once for the global mesh, and pick the face id that
3808 // is not the local one (cell_this).
3809 for (unsigned int i = 0; i < n_lanes; i++)
3810 {
3811 // compute actual (non vectorized) cell ID
3812 const unsigned int cell_this = this->cell * n_lanes + i;
3813 // compute face ID
3814 unsigned int face_index =
3815 this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3816 this->face_no,
3817 i);
3818
3819 if (face_index == numbers::invalid_unsigned_int)
3820 continue; // invalid face ID: no neighbor on boundary
3821
3822 // get cell ID on both sides of face
3823 auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3824 .cells_interior[face_index % n_lanes];
3825 auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3826 .cells_exterior[face_index % n_lanes];
3827
3828 // compare the IDs with the given cell ID
3829 if (cell_m == cell_this)
3830 cells[i] = cell_p; // neighbor has the other ID
3831 else if (cell_p == cell_this)
3832 cells[i] = cell_m;
3833 }
3834 }
3835 else if (is_face)
3836 {
3837 // face-centric faces
3838 const unsigned int *cells_ =
3839 is_interior_face ?
3840 &this->matrix_info->get_face_info(cell).cells_interior[0] :
3841 &this->matrix_info->get_face_info(cell).cells_exterior[0];
3842 for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3843 if (cells_[i] != numbers::invalid_unsigned_int)
3844 cells[i] = cells_[i];
3845 }
3846
3847 return cells;
3848}
3849
3850
3851namespace internal
3852{
3853 template <int dim,
3854 typename Number,
3855 bool is_face,
3856 typename VectorizedArrayType,
3857 typename VectorizedArrayType2,
3858 typename GlobalVectorType,
3859 typename FU>
3860 inline void
3861 process_cell_data(
3864 GlobalVectorType & array,
3865 VectorizedArrayType2 & out,
3866 const FU & fu)
3867 {
3868 (void)matrix_info;
3869 Assert(matrix_info != nullptr, ExcNotImplemented());
3870 AssertDimension(array.size(),
3871 matrix_info->get_task_info().cell_partition_data.back());
3872
3873 // 1) collect ids of cell
3874 const auto cells = phi.get_cell_ids();
3875
3876 // 2) actually gather values
3877 for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3878 if (cells[i] != numbers::invalid_unsigned_int)
3879 fu(out[i],
3880 array[cells[i] / VectorizedArrayType::size()]
3881 [cells[i] % VectorizedArrayType::size()]);
3882 }
3883} // namespace internal
3884
3885
3886
3887template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3888std::array<unsigned int, VectorizedArrayType::size()>
3891{
3892 const unsigned int v_len = VectorizedArrayType::size();
3893 std::array<unsigned int, VectorizedArrayType::size()> cells;
3894
3895 // initialize array
3896 for (unsigned int i = 0; i < v_len; ++i)
3898
3899 if (is_face &&
3900 this->dof_access_index ==
3902 this->is_interior_face == false)
3903 {
3904 // cell-based face-loop: plus face
3905 for (unsigned int i = 0; i < v_len; i++)
3906 {
3907 // compute actual (non vectorized) cell ID
3908 const unsigned int cell_this = this->cell * v_len + i;
3909 // compute face ID
3910 unsigned int fn =
3911 this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3912 this->face_no,
3913 i);
3914
3916 continue; // invalid face ID: no neighbor on boundary
3917
3918 // get cell ID on both sides of face
3919 auto cell_m = this->matrix_info->get_face_info(fn / v_len)
3920 .cells_interior[fn % v_len];
3921 auto cell_p = this->matrix_info->get_face_info(fn / v_len)
3922 .cells_exterior[fn % v_len];
3923
3924 // compare the IDs with the given cell ID
3925 if (cell_m == cell_this)
3926 cells[i] = cell_p; // neighbor has the other ID
3927 else if (cell_p == cell_this)
3928 cells[i] = cell_m;
3929 }
3930 }
3931 else
3932 {
3933 for (unsigned int i = 0; i < v_len; ++i)
3934 cells[i] = cell * v_len + i;
3935 }
3936
3937 return cells;
3938}
3939
3940
3941
3942template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3943inline VectorizedArrayType
3945 const AlignedVector<VectorizedArrayType> &array) const
3946{
3947 VectorizedArrayType out = Number(1.);
3948 internal::process_cell_data(
3949 *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3950 local = global;
3951 });
3952 return out;
3953}
3954
3955
3956
3957template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3958inline void
3961 const VectorizedArrayType & in) const
3962{
3963 internal::process_cell_data(
3964 *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3965 global = local;
3966 });
3967}
3968
3969
3970
3971template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3972template <typename T>
3973inline std::array<T, VectorizedArrayType::size()>
3975 const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
3976{
3977 std::array<T, VectorizedArrayType::size()> out;
3978 internal::process_cell_data(
3979 *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3980 local = global;
3981 });
3982 return out;
3983}
3984
3985
3986
3987template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3988template <typename T>
3989inline void
3991 AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
3992 const std::array<T, VectorizedArrayType::size()> & in) const
3993{
3994 internal::process_cell_data(
3995 *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3996 global = local;
3997 });
3998}
3999
4000
4001
4002template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4003inline const std::vector<unsigned int> &
4006{
4007 return data->lexicographic_numbering;
4008}
4009
4010
4011
4012template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4015 get_scratch_data() const
4016{
4018 const_cast<VectorizedArrayType *>(scratch_data),
4019 scratch_data_array->end() - scratch_data);
4020}
4021
4022
4023
4024template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4025inline unsigned int
4028{
4029 return this->quad_no;
4030}
4031
4032
4033
4034template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4035inline unsigned int
4038{
4039 if (is_face && this->dof_access_index ==
4041 return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
4042 else
4043 return this->cell;
4044}
4045
4046
4047
4048template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4049inline unsigned int
4052{
4053 return active_fe_index;
4054}
4055
4056
4057
4058template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4059inline unsigned int
4062{
4063 return active_quad_index;
4064}
4065
4066
4067
4068template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4071 get_matrix_free() const
4072{
4073 Assert(matrix_info != nullptr,
4074 ExcMessage(
4075 "FEEvaluation was not initialized with a MatrixFree object!"));
4076 return *matrix_info;
4077}
4078
4079
4080/*----------------------- FEEvaluationBase ----------------------------------*/
4081
4082template <int dim,
4083 int n_components_,
4084 typename Number,
4085 bool is_face,
4086 typename VectorizedArrayType>
4087inline FEEvaluationBase<dim,
4088 n_components_,
4089 Number,
4090 is_face,
4091 VectorizedArrayType>::
4092 FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
4093 const unsigned int dof_no,
4094 const unsigned int first_selected_component,
4095 const unsigned int quad_no_in,
4096 const unsigned int fe_degree,
4097 const unsigned int n_q_points,
4098 const bool is_interior_face,
4099 const unsigned int active_fe_index,
4100 const unsigned int active_quad_index,
4101 const unsigned int face_type)
4102 : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
4103 data_in,
4104 dof_no,
4105 first_selected_component,
4106 quad_no_in,
4107 fe_degree,
4108 n_q_points,
4109 is_interior_face,
4110 active_fe_index,
4111 active_quad_index,
4112 face_type)
4113 , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
4114 , dof_values_initialized(false)
4115 , values_quad_initialized(false)
4116 , gradients_quad_initialized(false)
4117 , hessians_quad_initialized(false)
4118 , values_quad_submitted(false)
4119 , gradients_quad_submitted(false)
4120 , first_selected_component(first_selected_component)
4121{
4122 set_data_pointers();
4123 Assert(
4124 this->dof_info->start_components.back() == 1 ||
4125 static_cast<int>(n_components_) <=
4126 static_cast<int>(
4127 this->dof_info->start_components
4128 [this->dof_info->component_to_base_index[first_selected_component] +
4129 1]) -
4130 first_selected_component,
4131 ExcMessage(
4132 "You tried to construct a vector-valued evaluator with " +
4133 std::to_string(n_components) +
4134 " components. However, "
4135 "the current base element has only " +
4137 this->dof_info->start_components
4138 [this->dof_info->component_to_base_index[first_selected_component] +
4139 1] -
4140 first_selected_component) +
4141 " components left when starting from local element index " +
4143 first_selected_component -
4144 this->dof_info->start_components
4145 [this->dof_info->component_to_base_index[first_selected_component]]) +
4146 " (global index " + std::to_string(first_selected_component) + ")"));
4147
4148 // do not check for correct dimensions of data fields here, should be done
4149 // in derived classes
4150}
4151
4152
4153
4154template <int dim,
4155 int n_components_,
4156 typename Number,
4157 bool is_face,
4158 typename VectorizedArrayType>
4159inline FEEvaluationBase<dim,
4160 n_components_,
4161 Number,
4162 is_face,
4163 VectorizedArrayType>::
4164 FEEvaluationBase(
4165 const Mapping<dim> & mapping,
4166 const FiniteElement<dim> &fe,
4167 const Quadrature<1> & quadrature,
4168 const UpdateFlags update_flags,
4169 const unsigned int first_selected_component,
4171 *other)
4172 : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
4173 mapping,
4174 fe,
4175 quadrature,
4176 update_flags,
4177 first_selected_component,
4178 other)
4179 , n_fe_components(n_components_)
4180 , dof_values_initialized(false)
4181 , values_quad_initialized(false)
4182 , gradients_quad_initialized(false)
4183 , hessians_quad_initialized(false)
4184 , values_quad_submitted(false)
4185 , gradients_quad_submitted(false)
4186 // keep the number of the selected component within the current base element
4187 // for reading dof values
4188 , first_selected_component(first_selected_component)
4189{
4190 set_data_pointers();
4191
4192 const unsigned int base_element_number =
4193 fe.component_to_base_index(first_selected_component).first;
4194 Assert(fe.element_multiplicity(base_element_number) == 1 ||
4195 fe.element_multiplicity(base_element_number) -
4196 first_selected_component >=
4197 n_components_,
4198 ExcMessage("The underlying element must at least contain as many "
4199 "components as requested by this class"));
4200 (void)base_element_number;
4201}
4202
4203
4204
4205template <int dim,
4206 int n_components_,
4207 typename Number,
4208 bool is_face,
4209 typename VectorizedArrayType>
4210inline FEEvaluationBase<dim,
4211 n_components_,
4212 Number,
4213 is_face,
4214 VectorizedArrayType>::
4215 FEEvaluationBase(const FEEvaluationBase<dim,
4216 n_components_,
4217 Number,
4218 is_face,
4219 VectorizedArrayType> &other)
4220 : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(other)
4221 , n_fe_components(other.n_fe_components)
4222 , dof_values_initialized(false)
4223 , values_quad_initialized(false)
4224 , gradients_quad_initialized(false)
4225 , hessians_quad_initialized(false)
4226 , values_quad_submitted(false)
4227 , gradients_quad_submitted(false)
4228 , first_selected_component(other.first_selected_component)
4229{
4230 set_data_pointers();
4231}
4232
4233
4234
4235template <int dim,
4236 int n_components_,
4237 typename Number,
4238 bool is_face,
4239 typename VectorizedArrayType>
4240inline FEEvaluationBase<dim,
4241 n_components_,
4242 Number,
4243 is_face,
4244 VectorizedArrayType> &
4246operator=(const FEEvaluationBase<dim,
4247 n_components_,
4248 Number,
4249 is_face,
4250 VectorizedArrayType> &other)
4251{
4253 operator=(other);
4254 AssertDimension(n_fe_components, other.n_fe_components);
4255 AssertDimension(first_selected_component, other.first_selected_component);
4256
4257 return *this;
4258}
4259
4260
4261
4262template <int dim,
4263 int n_components_,
4264 typename Number,
4265 bool is_face,
4266 typename VectorizedArrayType>
4267inline void
4270{
4271 Assert(this->scratch_data_array != nullptr, ExcInternalError());
4272
4273 const unsigned int tensor_dofs_per_component =
4274 Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4275 const unsigned int dofs_per_component =
4276 this->data->dofs_per_component_on_cell;
4277 const unsigned int n_quadrature_points = this->n_quadrature_points;
4278
4279 const unsigned int shift =
4280 std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4281 n_components_ * 3 +
4282 2 * n_quadrature_points;
4283 const unsigned int allocated_size =
4284 shift + n_components_ * dofs_per_component +
4285 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4286 this->scratch_data_array->resize_fast(allocated_size);
4287
4288 // set the pointers to the correct position in the data array
4289 for (unsigned int c = 0; c < n_components_; ++c)
4290 {
4291 values_dofs[c] =
4292 this->scratch_data_array->begin() + c * dofs_per_component;
4293 }
4294 values_quad =
4295 this->scratch_data_array->begin() + n_components * dofs_per_component;
4296 gradients_quad = this->scratch_data_array->begin() +
4297 n_components * (dofs_per_component + n_quadrature_points);
4298 hessians_quad =
4299 this->scratch_data_array->begin() +
4300 n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4301 this->scratch_data =
4302 this->scratch_data_array->begin() + n_components_ * dofs_per_component +
4303 (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4304}
4305
4306
4307
4308namespace internal
4309{
4310 // allows to select between block vectors and non-block vectors, which
4311 // allows to use a unified interface for extracting blocks on block vectors
4312 // and doing nothing on usual vectors
4313 template <typename VectorType, bool>
4314 struct BlockVectorSelector
4315 {};
4316
4317 template <typename VectorType>
4318 struct BlockVectorSelector<VectorType, true>
4319 {
4320 using BaseVectorType = typename VectorType::BlockType;
4321
4322 static BaseVectorType *
4323 get_vector_component(VectorType &vec, const unsigned int component)
4324 {
4325 AssertIndexRange(component, vec.n_blocks());
4326 return &vec.block(component);
4327 }
4328 };
4329
4330 template <typename VectorType>
4331 struct BlockVectorSelector<VectorType, false>
4332 {
4333 using BaseVectorType = VectorType;
4334
4335 static BaseVectorType *
4336 get_vector_component(VectorType &vec, const unsigned int component)
4337 {
4338 // FEEvaluation allows to combine several vectors from a scalar
4339 // FiniteElement into a "vector-valued" FEEvaluation object with
4340 // multiple components. These components can be extracted with the other
4341 // get_vector_component functions. If we do not get a vector of vectors
4342 // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
4343 // must make sure that we do not duplicate the components in input
4344 // and/or duplicate the resulting integrals. In such a case, we should
4345 // only get the zeroth component in the vector contained set nullptr for
4346 // the others which allows us to catch unintended use in
4347 // read_write_operation.
4348 if (component == 0)
4349 return &vec;
4350 else
4351 return nullptr;
4352 }
4353 };
4354
4355 template <typename VectorType>
4356 struct BlockVectorSelector<std::vector<VectorType>, false>
4357 {
4358 using BaseVectorType = VectorType;
4359
4360 static BaseVectorType *
4361 get_vector_component(std::vector<VectorType> &vec,
4362 const unsigned int component)
4363 {
4364 AssertIndexRange(component, vec.size());
4365 return &vec[component];
4366 }
4367 };
4368
4369 template <typename VectorType>
4370 struct BlockVectorSelector<std::vector<VectorType *>, false>
4371 {
4372 using BaseVectorType = VectorType;
4373
4374 static BaseVectorType *
4375 get_vector_component(std::vector<VectorType *> &vec,
4376 const unsigned int component)
4377 {
4378 AssertIndexRange(component, vec.size());
4379 return vec[component];
4380 }
4381 };
4382} // namespace internal
4383
4384
4385
4386template <int dim,
4387 int n_components_,
4388 typename Number,
4389 bool is_face,
4390 typename VectorizedArrayType>
4391template <typename VectorType, typename VectorOperation>
4392inline void
4395 const VectorOperation & operation,
4396 const std::array<VectorType *, n_components_> &src,
4397 const std::array<
4399 n_components_> & src_sm,
4400 const std::bitset<VectorizedArrayType::size()> &mask,
4401 const bool apply_constraints) const
4402{
4403 // Case 1: No MatrixFree object given, simple case because we do not need to
4404 // process constraints and need not care about vectorization -> go to
4405 // separate function
4406 if (this->matrix_info == nullptr)
4407 {
4408 read_write_operation_global(operation, src);
4409 return;
4410 }
4411
4412 Assert(this->dof_info != nullptr, ExcNotInitialized());
4413 Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
4414 if (n_fe_components == 1)
4415 for (unsigned int comp = 0; comp < n_components; ++comp)
4416 {
4417 Assert(src[comp] != nullptr,
4418 ExcMessage("The finite element underlying this FEEvaluation "
4419 "object is scalar, but you requested " +
4420 std::to_string(n_components) +
4421 " components via the template argument in "
4422 "FEEvaluation. In that case, you must pass an "
4423 "std::vector<VectorType> or a BlockVector to " +
4424 "read_dof_values and distribute_local_to_global."));
4425 internal::check_vector_compatibility(*src[comp], *this->dof_info);
4426 }
4427 else
4428 {
4429 internal::check_vector_compatibility(*src[0], *this->dof_info);
4430 }
4431
4432 // Case 2: contiguous indices which use reduced storage of indices and can
4433 // use vectorized load/store operations -> go to separate function
4435 this->cell,
4436 this->dof_info->index_storage_variants[this->dof_access_index].size());
4437 if (this->dof_info->index_storage_variants
4438 [is_face ? this->dof_access_index :
4440 [this->cell] >=
4442 {
4443 read_write_operation_contiguous(operation, src, src_sm, mask);
4444 return;
4445 }
4446
4447 // Case 3: standard operation with one index per degree of freedom -> go on
4448 // here
4449 constexpr unsigned int n_lanes = VectorizedArrayType::size();
4450 Assert(mask.count() == n_lanes,
4451 ExcNotImplemented("Masking currently not implemented for "
4452 "non-contiguous DoF storage"));
4453
4454 std::integral_constant<bool,
4455 internal::is_vectorizable<VectorType, Number>::value>
4456 vector_selector;
4457
4458 const unsigned int dofs_per_component =
4459 this->data->dofs_per_component_on_cell;
4460 if (this->dof_info->index_storage_variants
4461 [is_face ? this->dof_access_index :
4463 [this->cell] ==
4465 {
4466 const unsigned int *dof_indices =
4467 this->dof_info->dof_indices_interleaved.data() +
4468 this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4469 .first +
4470 this->dof_info
4471 ->component_dof_indices_offset[this->active_fe_index]
4472 [this->first_selected_component] *
4473 n_lanes;
4474 if (n_components == 1 || n_fe_components == 1)
4475 for (unsigned int i = 0; i < dofs_per_component;
4476 ++i, dof_indices += n_lanes)
4477 for (unsigned int comp = 0; comp < n_components; ++comp)
4478 operation.process_dof_gather(dof_indices,
4479 *src[comp],
4480 0,
4481 values_dofs[comp][i],
4482 vector_selector);
4483 else
4484 for (unsigned int comp = 0; comp < n_components; ++comp)
4485 for (unsigned int i = 0; i < dofs_per_component;
4486 ++i, dof_indices += n_lanes)
4487 operation.process_dof_gather(
4488 dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4489 return;
4490 }
4491
4492 const unsigned int * dof_indices[n_lanes];
4493 VectorizedArrayType **values_dofs =
4494 const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
4495
4496 // Assign the appropriate cell ids for face/cell case and get the pointers
4497 // to the dof indices of the cells on all lanes
4498 unsigned int cells_copied[n_lanes];
4499 const unsigned int *cells;
4500 unsigned int n_vectorization_actual =
4501 this->dof_info
4502 ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4503 bool has_constraints = false;
4504 const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4505 if (is_face)
4506 {
4507 if (this->dof_access_index ==
4509 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4510 cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4511 cells =
4512 this->dof_access_index ==
4514 &cells_copied[0] :
4515 (this->is_interior_face ?
4516 &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
4517 &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
4518 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4519 {
4520 Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4522 const std::pair<unsigned int, unsigned int> *my_index_start =
4523 &this->dof_info->row_starts[cells[v] * n_fe_components +
4524 first_selected_component];
4525
4526 // check whether any of the SIMD lanes has constraints, i.e., the
4527 // constraint indicator which is the second entry of row_starts
4528 // increments on this cell
4529 if (my_index_start[n_components_read].second !=
4530 my_index_start[0].second)
4531 has_constraints = true;
4532
4533 dof_indices[v] =
4534 this->dof_info->dof_indices.data() + my_index_start[0].first;
4535 }
4536 for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4537 dof_indices[v] = nullptr;
4538 }
4539 else
4540 {
4541 AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
4542 this->dof_info->row_starts.size());
4543 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4544 {
4545 const std::pair<unsigned int, unsigned int> *my_index_start =
4546 &this->dof_info
4547 ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4548 first_selected_component];
4549 if (my_index_start[n_components_read].second !=
4550 my_index_start[0].second)
4551 has_constraints = true;
4552 Assert(my_index_start[n_components_read].first ==
4553 my_index_start[0].first ||
4554 my_index_start[0].first < this->dof_info->dof_indices.size(),
4555 ExcIndexRange(0,
4556 my_index_start[0].first,
4557 this->dof_info->dof_indices.size()));
4558 dof_indices[v] =
4559 this->dof_info->dof_indices.data() + my_index_start[0].first;
4560 }
4561 for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4562 dof_indices[v] = nullptr;
4563 }
4564
4565 // Case where we have no constraints throughout the whole cell: Can go
4566 // through the list of DoFs directly
4567 if (!has_constraints)
4568 {
4569 if (n_vectorization_actual < n_lanes)
4570 for (unsigned int comp = 0; comp < n_components; ++comp)
4571 for (unsigned int i = 0; i < dofs_per_component; ++i)
4572 operation.process_empty(values_dofs[comp][i]);
4573 if (n_components == 1 || n_fe_components == 1)
4574 {
4575 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4576 for (unsigned int i = 0; i < dofs_per_component; ++i)
4577 for (unsigned int comp = 0; comp < n_components; ++comp)
4578 operation.process_dof(dof_indices[v][i],
4579 *src[comp],
4580 values_dofs[comp][i][v]);
4581 }
4582 else
4583 {
4584 for (unsigned int comp = 0; comp < n_components; ++comp)
4585 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4586 for (unsigned int i = 0; i < dofs_per_component; ++i)
4587 operation.process_dof(
4588 dof_indices[v][comp * dofs_per_component + i],
4589 *src[0],
4590 values_dofs[comp][i][v]);
4591 }
4592 return;
4593 }
4594
4595 // In the case where there are some constraints to be resolved, loop over
4596 // all vector components that are filled and then over local dofs. ind_local
4597 // holds local number on cell, index iterates over the elements of
4598 // index_local_to_global and dof_indices points to the global indices stored
4599 // in index_local_to_global
4600 if (n_vectorization_actual < n_lanes)
4601 for (unsigned int comp = 0; comp < n_components; ++comp)
4602 for (unsigned int i = 0; i < dofs_per_component; ++i)
4603 operation.process_empty(values_dofs[comp][i]);
4604 for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4605 {
4606 const unsigned int cell_index =
4607 is_face ? cells[v] : this->cell * n_lanes + v;
4608 const unsigned int cell_dof_index =
4609 cell_index * n_fe_components + first_selected_component;
4610 const unsigned int n_components_read =
4611 n_fe_components > 1 ? n_components : 1;
4612 unsigned int index_indicators =
4613 this->dof_info->row_starts[cell_dof_index].second;
4614 unsigned int next_index_indicators =
4615 this->dof_info->row_starts[cell_dof_index + 1].second;
4616
4617 // For read_dof_values_plain, redirect the dof_indices field to the
4618 // unconstrained indices
4619 if (apply_constraints == false &&
4620 this->dof_info->row_starts[cell_dof_index].second !=
4621 this->dof_info->row_starts[cell_dof_index + n_components_read]
4622 .second)
4623 {
4624 Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
4627 dof_indices[v] =
4628 this->dof_info->plain_dof_indices.data() +
4629 this->dof_info
4630 ->component_dof_indices_offset[this->active_fe_index]
4631 [this->first_selected_component] +
4632 this->dof_info->row_starts_plain_indices[cell_index];
4633 next_index_indicators = index_indicators;
4634 }
4635
4636 if (n_components == 1 || n_fe_components == 1)
4637 {
4638 unsigned int ind_local = 0;
4639 for (; index_indicators != next_index_indicators; ++index_indicators)
4640 {
4641 const std::pair<unsigned short, unsigned short> indicator =
4642 this->dof_info->constraint_indicator[index_indicators];
4643 // run through values up to next constraint
4644 for (unsigned int j = 0; j < indicator.first; ++j)
4645 for (unsigned int comp = 0; comp < n_components; ++comp)
4646 operation.process_dof(dof_indices[v][j],
4647 *src[comp],
4648 values_dofs[comp][ind_local + j][v]);
4649
4650 ind_local += indicator.first;
4651 dof_indices[v] += indicator.first;
4652
4653 // constrained case: build the local value as a linear
4654 // combination of the global value according to constraints
4655 Number value[n_components];
4656 for (unsigned int comp = 0; comp < n_components; ++comp)
4657 operation.pre_constraints(values_dofs[comp][ind_local][v],
4658 value[comp]);
4659
4660 const Number *data_val =
4661 this->matrix_info->constraint_pool_begin(indicator.second);
4662 const Number *end_pool =
4663 this->matrix_info->constraint_pool_end(indicator.second);
4664 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4665 for (unsigned int comp = 0; comp < n_components; ++comp)
4666 operation.process_constraint(*dof_indices[v],
4667 *data_val,
4668 *src[comp],
4669 value[comp]);
4670
4671 for (unsigned int comp = 0; comp < n_components; ++comp)
4672 operation.post_constraints(value[comp],
4673 values_dofs[comp][ind_local][v]);
4674 ind_local++;
4675 }
4676
4677 AssertIndexRange(ind_local, dofs_per_component + 1);
4678
4679 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4680 for (unsigned int comp = 0; comp < n_components; ++comp)
4681 operation.process_dof(*dof_indices[v],
4682 *src[comp],
4683 values_dofs[comp][ind_local][v]);
4684 }
4685 else
4686 {
4687 // case with vector-valued finite elements where all components are
4688 // included in one single vector. Assumption: first come all entries
4689 // to the first component, then all entries to the second one, and
4690 // so on. This is ensured by the way MatrixFree reads out the
4691 // indices.
4692 for (unsigned int comp = 0; comp < n_components; ++comp)
4693 {
4694 unsigned int ind_local = 0;
4695
4696 // check whether there is any constraint on the current cell
4697 for (; index_indicators != next_index_indicators;
4698 ++index_indicators)
4699 {
4700 const std::pair<unsigned short, unsigned short> indicator =
4701 this->dof_info->constraint_indicator[index_indicators];
4702
4703 // run through values up to next constraint
4704 for (unsigned int j = 0; j < indicator.first; ++j)
4705 operation.process_dof(dof_indices[v][j],
4706 *src[0],
4707 values_dofs[comp][ind_local + j][v]);
4708 ind_local += indicator.first;
4709 dof_indices[v] += indicator.first;
4710
4711 // constrained case: build the local value as a linear
4712 // combination of the global value according to constraints
4713 Number value;
4714 operation.pre_constraints(values_dofs[comp][ind_local][v],
4715 value);
4716
4717 const Number *data_val =
4718 this->matrix_info->constraint_pool_begin(indicator.second);
4719 const Number *end_pool =
4720 this->matrix_info->constraint_pool_end(indicator.second);
4721
4722 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4723 operation.process_constraint(*dof_indices[v],
4724 *data_val,
4725 *src[0],
4726 value);
4727
4728 operation.post_constraints(value,
4729 values_dofs[comp][ind_local][v]);
4730 ind_local++;
4731 }
4732
4733 AssertIndexRange(ind_local, dofs_per_component + 1);
4734
4735 // get the dof values past the last constraint
4736 for (; ind_local < dofs_per_component;
4737 ++dof_indices[v], ++ind_local)
4738 {
4739 AssertIndexRange(*dof_indices[v], src[0]->size());
4740 operation.process_dof(*dof_indices[v],
4741 *src[0],
4742 values_dofs[comp][ind_local][v]);
4743 }
4744
4745 if (apply_constraints == true && comp + 1 < n_components)
4746 next_index_indicators =
4747 this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4748 }
4749 }
4750 }
4751}
4752
4753
4754
4755template <int dim,
4756 int n_components_,
4757 typename Number,
4758 bool is_face,
4759 typename VectorizedArrayType>
4760template <typename VectorType, typename VectorOperation>
4761inline void
4764 const VectorOperation & operation,
4765 const std::array<VectorType *, n_components_> &src) const
4766{
4767 Assert(!local_dof_indices.empty(), ExcNotInitialized());
4768
4769 unsigned int index =
4770 first_selected_component * this->data->dofs_per_component_on_cell;
4771 for (unsigned int comp = 0; comp < n_components; ++comp)
4772 {
4773 for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4774 ++i, ++index)
4775 {
4776 operation.process_empty(values_dofs[comp][i]);
4777 operation.process_dof_global(
4778 local_dof_indices[this->data->lexicographic_numbering[index]],
4779 *src[0],
4780 values_dofs[comp][i][0]);
4781 }
4782 }
4783}
4784
4785
4786
4787template <int dim,
4788 int n_components_,
4789 typename Number,
4790 bool is_face,
4791 typename VectorizedArrayType>
4792template <typename VectorType, typename VectorOperation>
4793inline void
4796 const VectorOperation & operation,
4797 const std::array<VectorType *, n_components_> &src,
4798 const std::array<
4800 n_components_> & vectors_sm,
4801 const std::bitset<VectorizedArrayType::size()> &mask) const
4802{
4803 // This functions processes the functions read_dof_values,
4804 // distribute_local_to_global, and set_dof_values with the same code for
4805 // contiguous cell indices (DG case). The distinction between these three
4806 // cases is made by the input VectorOperation that either reads values from
4807 // a vector and puts the data into the local data field or write local data
4808 // into the vector. Certain operations are no-ops for the given use case.
4809
4810 std::integral_constant<bool,
4811 internal::is_vectorizable<VectorType, Number>::value>
4812 vector_selector;
4814 is_face ? this->dof_access_index :
4816 const unsigned int n_lanes = mask.count();
4817
4818 const std::vector<unsigned int> &dof_indices_cont =
4819 this->dof_info->dof_indices_contiguous[ind];
4820
4821 // Simple case: We have contiguous storage, so we can simply copy out the
4822 // data
4823 if ((this->dof_info->index_storage_variants[ind][this->cell] ==
4825 interleaved_contiguous &&
4826 n_lanes == VectorizedArrayType::size()) &&
4827 !(is_face &&
4828 this->dof_access_index ==
4830 this->is_interior_face == false))
4831 {
4832 const unsigned int dof_index =
4833 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4834 this->dof_info->component_dof_indices_offset[this->active_fe_index]
4835 [first_selected_component] *
4836 VectorizedArrayType::size();
4837 if (n_components == 1 || n_fe_components == 1)
4838 for (unsigned int comp = 0; comp < n_components; ++comp)
4839 operation.process_dofs_vectorized(
4840 this->data->dofs_per_component_on_cell,
4841 dof_index,
4842 *src[comp],
4843 values_dofs[comp],
4844 vector_selector);
4845 else
4846 operation.process_dofs_vectorized(
4847 this->data->dofs_per_component_on_cell * n_components,
4848 dof_index,
4849 *src[0],
4850 values_dofs[0],
4851 vector_selector);
4852 return;
4853 }
4854
4855 std::array<unsigned int, VectorizedArrayType::size()> cells =
4856 this->get_cell_or_face_ids();
4857
4858 // More general case: Must go through the components one by one and apply
4859 // some transformations
4860 const unsigned int n_filled_lanes =
4861 this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4862
4863 const bool is_ecl =
4864 this->dof_access_index ==
4866 this->is_interior_face == false;
4867
4868 if (vectors_sm[0] != nullptr)
4869 {
4870 const auto compute_vector_ptrs = [&](const unsigned int comp) {
4871 std::array<typename VectorType::value_type *,
4872 VectorizedArrayType::size()>
4873 vector_ptrs = {};
4874
4875 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4876 {
4879 Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
4881 ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
4882 Assert(cells[v] <
4883 this->dof_info->dof_indices_contiguous_sm[ind].size(),
4885 cells[v],
4886 0,
4887 this->dof_info->dof_indices_contiguous_sm[ind].size()));
4888
4889 const auto &temp =
4890 this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
4891
4892 if (temp.first != numbers::invalid_unsigned_int)
4893 vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
4894 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
4895 this->dof_info->component_dof_indices_offset
4896 [this->active_fe_index][this->first_selected_component]);
4897 else
4898 vector_ptrs[v] = nullptr;
4899 }
4900 for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
4901 ++v)
4902 vector_ptrs[v] = nullptr;
4903
4904 return vector_ptrs;
4905 };
4906
4907 if (n_filled_lanes == VectorizedArrayType::size() &&
4908 n_lanes == VectorizedArrayType::size() && !is_ecl)
4909 {
4910 if (n_components == 1 || n_fe_components == 1)
4911 {
4912 for (unsigned int comp = 0; comp < n_components; ++comp)
4913 {
4914 auto vector_ptrs = compute_vector_ptrs(comp);
4915 operation.process_dofs_vectorized_transpose(
4916 this->data->dofs_per_component_on_cell,
4917 vector_ptrs,
4918 values_dofs[comp],
4919 vector_selector);
4920 }
4921 }
4922 else
4923 {
4924 auto vector_ptrs = compute_vector_ptrs(0);
4925 operation.process_dofs_vectorized_transpose(
4926 this->data->dofs_per_component_on_cell * n_components,
4927 vector_ptrs,
4928 &values_dofs[0][0],
4929 vector_selector);
4930 }
4931 }
4932 else
4933 for (unsigned int comp = 0; comp < n_components; ++comp)
4934 {
4935 auto vector_ptrs = compute_vector_ptrs(
4936 (n_components == 1 || n_fe_components == 1) ? comp : 0);
4937
4938 for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4939 ++i)
4940 operation.process_empty(values_dofs[comp][i]);
4941
4942 if (n_components == 1 || n_fe_components == 1)
4943 {
4944 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4945 if (mask[v] == true)
4946 for (unsigned int i = 0;
4947 i < this->data->dofs_per_component_on_cell;
4948 ++i)
4949 operation.process_dof(vector_ptrs[v][i],
4950 values_dofs[comp][i][v]);
4951 }
4952 else
4953 {
4954 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4955 if (mask[v] == true)
4956 for (unsigned int i = 0;
4957 i < this->data->dofs_per_component_on_cell;
4958 ++i)
4959 operation.process_dof(
4960 vector_ptrs[v]
4961 [i + comp * this->data
4962 ->dofs_per_component_on_cell],
4963 values_dofs[comp][i][v]);
4964 }
4965 }
4966 return;
4967 }
4968
4969 unsigned int dof_indices[VectorizedArrayType::size()];
4970
4971 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4972 {
4974 dof_indices[v] =
4975 dof_indices_cont[cells[v]] +
4976 this->dof_info
4977 ->component_dof_indices_offset[this->active_fe_index]
4978 [this->first_selected_component] *
4979 this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
4980 }
4981
4982 for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4983 dof_indices[v] = numbers::invalid_unsigned_int;
4984
4985 // In the case with contiguous cell indices, we know that there are no
4986 // constraints and that the indices within each element are contiguous
4987 if (n_filled_lanes == VectorizedArrayType::size() &&
4988 n_lanes == VectorizedArrayType::size() && !is_ecl)
4989 {
4990 if (this->dof_info->index_storage_variants[ind][this->cell] ==
4992 contiguous)
4993 {
4994 if (n_components == 1 || n_fe_components == 1)
4995 for (unsigned int comp = 0; comp < n_components; ++comp)
4996 operation.process_dofs_vectorized_transpose(
4997 this->data->dofs_per_component_on_cell,
4998 dof_indices,
4999 *src[comp],
5000 values_dofs[comp],
5001 vector_selector);
5002 else
5003 operation.process_dofs_vectorized_transpose(
5004 this->data->dofs_per_component_on_cell * n_components,
5005 dof_indices,
5006 *src[0],
5007 &values_dofs[0][0],
5008 vector_selector);
5009 }
5010 else if (this->dof_info->index_storage_variants[ind][this->cell] ==
5012 interleaved_contiguous_strided)
5013 {
5014 if (n_components == 1 || n_fe_components == 1)
5015 for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5016 ++i)
5017 {
5018 for (unsigned int comp = 0; comp < n_components; ++comp)
5019 operation.process_dof_gather(dof_indices,
5020 *src[comp],
5021 i * VectorizedArrayType::size(),
5022 values_dofs[comp][i],
5023 vector_selector);
5024 }
5025 else
5026 for (unsigned int comp = 0; comp < n_components; ++comp)
5027 for (unsigned int i = 0;
5028 i < this->data->dofs_per_component_on_cell;
5029 ++i)
5030 {
5031 operation.process_dof_gather(
5032 dof_indices,
5033 *src[0],
5034 (comp * this->data->dofs_per_component_on_cell + i) *
5035 VectorizedArrayType::size(),
5036 values_dofs[comp][i],
5037 vector_selector);
5038 }
5039 }
5040 else
5041 {
5042 Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
5044 IndexStorageVariants::interleaved_contiguous_mixed_strides,
5046 const unsigned int *offsets =
5047 &this->dof_info->dof_indices_interleave_strides
5048 [ind][VectorizedArrayType::size() * this->cell];
5049 if (n_components == 1 || n_fe_components == 1)
5050 for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5051 ++i)
5052 {
5053 for (unsigned int comp = 0; comp < n_components; ++comp)
5054 operation.process_dof_gather(dof_indices,
5055 *src[comp],
5056 0,
5057 values_dofs[comp][i],
5058 vector_selector);
5060 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5061 dof_indices[v] += offsets[v];
5062 }
5063 else
5064 for (unsigned int comp = 0; comp < n_components; ++comp)
5065 for (unsigned int i = 0;
5066 i < this->data->dofs_per_component_on_cell;
5067 ++i)
5068 {
5069 operation.process_dof_gather(dof_indices,
5070 *src[0],
5071 0,
5072 values_dofs[comp][i],
5073 vector_selector);
5075 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5076 dof_indices[v] += offsets[v];
5077 }
5078 }
5079 }
5080 else
5081 for (unsigned int comp = 0; comp < n_components; ++comp)
5082 {
5083 for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5084 ++i)
5085 operation.process_empty(values_dofs[comp][i]);
5086 if (this->dof_info->index_storage_variants[ind][this->cell] ==
5088 contiguous)
5089 {
5090 if (n_components == 1 || n_fe_components == 1)
5091 {
5092 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5093 if (mask[v] == true)
5094 for (unsigned int i = 0;
5095 i < this->data->dofs_per_component_on_cell;
5096 ++i)
5097 operation.process_dof(dof_indices[v] + i,
5098 *src[comp],
5099 values_dofs[comp][i][v]);
5100 }
5101 else
5102 {
5103 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5104 if (mask[v] == true)
5105 for (unsigned int i = 0;
5106 i < this->data->dofs_per_component_on_cell;
5107 ++i)
5108 operation.process_dof(
5109 dof_indices[v] + i +
5110 comp * this->data->dofs_per_component_on_cell,
5111 *src[0],
5112 values_dofs[comp][i][v]);
5113 }
5114 }
5115 else
5116 {
5117 const unsigned int *offsets =
5118 &this->dof_info->dof_indices_interleave_strides
5119 [ind][VectorizedArrayType::size() * this->cell];
5120 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5121 AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
5122 if (n_components == 1 || n_fe_components == 1)
5123 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5124 {
5125 if (mask[v] == true)
5126 for (unsigned int i = 0;
5127 i < this->data->dofs_per_component_on_cell;
5128 ++i)
5129 operation.process_dof(dof_indices[v] + i * offsets[v],
5130 *src[comp],
5131 values_dofs[comp][i][v]);
5132 }
5133 else
5134 {
5135 for (unsigned int v = 0; v < n_filled_lanes; ++v)
5136 if (mask[v] == true)
5137 for (unsigned int i = 0;
5138 i < this->data->dofs_per_component_on_cell;
5139 ++i)
5140 operation.process_dof(
5141 dof_indices[v] +
5142 (i + comp * this->data->dofs_per_component_on_cell) *
5143 offsets[v],
5144 *src[0],
5145 values_dofs[comp][i][v]);
5146 }
5147 }
5148 }
5149}
5150
5151namespace internal
5152{
5153 template <typename Number,
5154 typename VectorType,
5155 typename std::enable_if<!IsBlockVector<VectorType>::value,
5156 VectorType>::type * = nullptr>
5157 decltype(std::declval<VectorType>().begin())
5158 get_beginning(VectorType &vec)
5159 {
5160 return vec.begin();
5161 }
5162
5163 template <typename Number,
5164 typename VectorType,
5165 typename std::enable_if<IsBlockVector<VectorType>::value,
5166 VectorType>::type * = nullptr>
5167 typename VectorType::value_type *
5168 get_beginning(VectorType &)
5169 {
5170 return nullptr;
5171 }
5172
5173 template <typename VectorType,
5174 typename std::enable_if<has_shared_vector_data<VectorType>::value,
5175 VectorType>::type * = nullptr>
5176 const std::vector<ArrayView<const typename VectorType::value_type>> *
5177 get_shared_vector_data(VectorType & vec,
5178 const bool is_valid_mode_for_sm,
5179 const unsigned int active_fe_index,
5181 {
5182 // note: no hp is supported
5183 if (is_valid_mode_for_sm &&
5184 dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
5185 .size() > 0 &&
5186 active_fe_index == 0)
5187 return &vec.shared_vector_data();
5188 else
5189 return nullptr;
5190 }
5191
5192 template <typename VectorType,
5193 typename std::enable_if<!has_shared_vector_data<VectorType>::value,
5194 VectorType>::type * = nullptr>
5195 const std::vector<ArrayView<const typename VectorType::value_type>> *
5196 get_shared_vector_data(VectorType &,
5197 const bool,
5198 const unsigned int,
5200 {
5201 return nullptr;
5202 }
5203
5204 template <int n_components, typename VectorType>
5205 std::pair<
5206 std::array<typename internal::BlockVectorSelector<
5207 typename std::remove_const<VectorType>::type,
5209 value>::BaseVectorType *,
5210 n_components>,
5211 std::array<
5212 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
5213 typename std::remove_const<VectorType>::type,
5215 BaseVectorType::value_type>> *,
5216 n_components>>
5217 get_vector_data(VectorType & src,
5218 const unsigned int first_index,
5219 const bool is_valid_mode_for_sm,
5220 const unsigned int active_fe_index,
5222 {
5223 // select between block vectors and non-block vectors. Note that the number
5224 // of components is checked in the internal data
5225 std::pair<
5226 std::array<typename internal::BlockVectorSelector<
5227 typename std::remove_const<VectorType>::type,
5229 value>::BaseVectorType *,
5230 n_components>,
5231 std::array<
5232 const std::vector<
5233 ArrayView<const typename internal::BlockVectorSelector<
5234 typename std::remove_const<VectorType>::type,
5236 value>::BaseVectorType::value_type>> *,
5237 n_components>>
5238 src_data;
5239
5240 for (unsigned int d = 0; d < n_components; ++d)
5241 src_data.first[d] = internal::BlockVectorSelector<
5242 typename std::remove_const<VectorType>::type,
5243 IsBlockVector<typename std::remove_const<VectorType>::type>::value>::
5244 get_vector_component(
5245 const_cast<typename std::remove_const<VectorType>::type &>(src),
5246 d + first_index);
5247
5248 for (unsigned int d = 0; d < n_components; ++d)
5249 src_data.second[d] = get_shared_vector_data(*src_data.first[d],
5250 is_valid_mode_for_sm,
5251 active_fe_index,
5252 dof_info);
5253
5254 return src_data;
5255 }
5256} // namespace internal
5257
5258
5259
5260template <int dim,
5261 int n_components_,
5262 typename Number,
5263 bool is_face,
5264 typename VectorizedArrayType>
5265template <typename VectorType>
5266inline void
5268 read_dof_values(const VectorType &src, const unsigned int first_index)
5269{
5270 const auto src_data = internal::get_vector_data<n_components_>(
5271 src,
5272 first_index,
5273 this->dof_access_index ==
5275 this->active_fe_index,
5276 this->dof_info);
5277
5279 read_write_operation(reader,
5280 src_data.first,
5281 src_data.second,
5282 std::bitset<VectorizedArrayType::size()>().flip(),
5283 true);
5284
5285# ifdef DEBUG
5286 dof_values_initialized = true;
5287# endif
5288}
5289
5290
5291
5292template <int dim,
5293 int n_components_,
5294 typename Number,
5295 bool is_face,
5296 typename VectorizedArrayType>
5297template <typename VectorType>
5298inline void
5300 read_dof_values_plain(const VectorType &src, const unsigned int first_index)
5301{
5302 const auto src_data = internal::get_vector_data<n_components_>(
5303 src,
5304 first_index,
5305 this->dof_access_index ==
5307 this->active_fe_index,
5308 this->dof_info);
5309
5311 read_write_operation(reader,
5312 src_data.first,
5313 src_data.second,
5314 std::bitset<VectorizedArrayType::size()>().flip(),
5315 false);
5316
5317# ifdef DEBUG
5318 dof_values_initialized = true;
5319# endif
5320}
5321
5322
5323
5324template <int dim,
5325 int n_components_,
5326 typename Number,
5327 bool is_face,
5328 typename VectorizedArrayType>
5329template <typename VectorType>
5330inline void
5333 VectorType & dst,
5334 const unsigned int first_index,
5335 const std::bitset<VectorizedArrayType::size()> &mask) const
5336{
5337# ifdef DEBUG
5338 Assert(dof_values_initialized == true,
5340# endif
5341
5342 const auto dst_data = internal::get_vector_data<n_components_>(
5343 dst,
5344 first_index,
5345 this->dof_access_index ==
5347 this->active_fe_index,
5348 this->dof_info);
5349
5351 distributor;
5352 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
5353}
5354
5355
5356
5357template <int dim,
5358 int n_components_,
5359 typename Number,
5360 bool is_face,
5361 typename VectorizedArrayType>
5362template <typename VectorType>
5363inline void
5365 set_dof_values(VectorType & dst,
5366 const unsigned int first_index,
5367 const std::bitset<VectorizedArrayType::size()> &mask) const
5368{
5369# ifdef DEBUG
5370 Assert(dof_values_initialized == true,
5372# endif
5373
5374 const auto dst_data = internal::get_vector_data<n_components_>(
5375 dst,
5376 first_index,
5377 this->dof_access_index ==
5379 this->active_fe_index,
5380 this->dof_info);
5381
5383 read_write_operation(setter, dst_data.first, dst_data.second, mask);
5384}
5385
5386
5387
5388template <int dim,
5389 int n_components_,
5390 typename Number,
5391 bool is_face,
5392 typename VectorizedArrayType>
5393template <typename VectorType>
5394inline void
5397 VectorType & dst,
5398 const unsigned int first_index,
5399 const std::bitset<VectorizedArrayType::size()> &mask) const
5400{
5401# ifdef DEBUG
5402 Assert(dof_values_initialized == true,
5404# endif
5405
5406 const auto dst_data = internal::get_vector_data<n_components_>(
5407 dst,
5408 first_index,
5409 this->dof_access_index ==
5411 this->active_fe_index,
5412 this->dof_info);
5413
5415 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
5416}
5417
5418
5419
5420/*------------------------------ access to data fields ----------------------*/
5421
5422
5423
5424template <int dim,
5425 int n_components,
5426 typename Number,
5427 bool is_face,
5428 typename VectorizedArrayType>
5429inline const VectorizedArrayType *
5431 begin_dof_values() const
5432{
5433 return &values_dofs[0][0];
5434}
5435
5436
5437
5438template <int dim,
5439 int n_components,
5440 typename Number,
5441 bool is_face,
5442 typename VectorizedArrayType>
5443inline VectorizedArrayType *
5446{
5447# ifdef DEBUG
5448 dof_values_initialized = true;
5449# endif
5450 return &values_dofs[0][0];
5451}
5452
5453
5454
5455template <int dim,
5456 int n_components,
5457 typename Number,
5458 bool is_face,
5459 typename VectorizedArrayType>
5460inline const VectorizedArrayType *
5462 begin_values() const
5463{
5464# ifdef DEBUG
5465 Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
5466# endif
5467 return values_quad;
5468}
5469
5470
5471
5472template <int dim,
5473 int n_components,
5474 typename Number,
5475 bool is_face,
5476 typename VectorizedArrayType>
5477inline VectorizedArrayType *
5480{
5481# ifdef DEBUG
5482 values_quad_initialized = true;
5483 values_quad_submitted = true;
5484# endif
5485 return values_quad;
5486}
5487
5488
5489
5490template <int dim,
5491 int n_components,
5492 typename Number,
5493 bool is_face,
5494 typename VectorizedArrayType>
5495inline const VectorizedArrayType *
5497 begin_gradients() const
5498{
5499# ifdef DEBUG
5500 Assert(gradients_quad_initialized || gradients_quad_submitted,
5502# endif
5503 return gradients_quad;
5504}
5505
5506
5507
5508template <int dim,
5509 int n_components,
5510 typename Number,
5511 bool is_face,
5512 typename VectorizedArrayType>
5513inline VectorizedArrayType *
5516{
5517# ifdef DEBUG
5518 gradients_quad_submitted = true;
5519 gradients_quad_initialized = true;
5520# endif
5521 return gradients_quad;
5522}
5523
5524
5525
5526template <int dim,
5527 int n_components,
5528 typename Number,
5529 bool is_face,
5530 typename VectorizedArrayType>
5531inline const VectorizedArrayType *
5533 begin_hessians() const
5534{
5535# ifdef DEBUG
5536 Assert(hessians_quad_initialized, ExcNotInitialized());
5537# endif
5538 return hessians_quad;
5539}
5540
5541
5542
5543template <int dim,
5544 int n_components,
5545 typename Number,
5546 bool is_face,
5547 typename VectorizedArrayType>
5548inline VectorizedArrayType *
5551{
5552# ifdef DEBUG
5553 hessians_quad_initialized = true;
5554# endif
5555 return hessians_quad;
5556}
5557
5558
5559
5560template <int dim,
5561 int n_components,
5562 typename Number,
5563 bool is_face,
5564 typename VectorizedArrayType>
5565inline unsigned int
5568{
5569 return first_selected_component;
5570}
5571
5572
5573
5574template <int dim,
5575 int n_components_,
5576 typename Number,
5577 bool is_face,
5578 typename VectorizedArrayType>
5581 get_dof_value(const unsigned int dof) const
5582{
5583 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5585 for (unsigned int comp = 0; comp < n_components; comp++)
5586 return_value[comp] = this->values_dofs[comp][dof];
5587 return return_value;
5588}
5589
5590
5591
5592template <int dim,
5593 int n_components_,
5594 typename Number,
5595 bool is_face,
5596 typename VectorizedArrayType>
5599 get_value(const unsigned int q_point) const
5600{
5601# ifdef DEBUG
5602 Assert(this->values_quad_initialized == true,
5604# endif
5605
5606 AssertIndexRange(q_point, this->n_quadrature_points);
5607 const std::size_t nqp = this->n_quadrature_points;
5609 for (unsigned int comp = 0; comp < n_components; comp++)
5610 return_value[comp] = values_quad[comp * nqp + q_point];
5611 return return_value;
5612}
5613
5614
5615
5616template <int dim,
5617 int n_components_,
5618 typename Number,
5619 bool is_face,
5620 typename VectorizedArrayType>
5624 get_gradient(const unsigned int q_point) const
5625{
5626# ifdef DEBUG
5627 Assert(this->gradients_quad_initialized == true,
5629# endif
5630
5631 AssertIndexRange(q_point, this->n_quadrature_points);
5632 Assert(this->jacobian != nullptr,
5634 "update_gradients"));
5635 const std::size_t nqp = this->n_quadrature_points;
5637
5638 // Cartesian cell
5639 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5640 {
5641 for (unsigned int d = 0; d < dim; ++d)
5642 for (unsigned int comp = 0; comp < n_components; comp++)
5643 grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] *
5644 this->jacobian[0][d][d];
5645 }
5646 // cell with general/affine Jacobian
5647 else
5648 {
5650 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5651 q_point :
5652 0];
5653 for (unsigned int comp = 0; comp < n_components; comp++)
5654 for (unsigned int d = 0; d < dim; ++d)
5655 {
5656 grad_out[comp][d] =
5657 jac[d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5658 for (unsigned int e = 1; e < dim; ++e)
5659 grad_out[comp][d] +=
5660 jac[d][e] * gradients_quad[(comp * dim + e) * nqp + q_point];
5661 }
5662 }
5663 return grad_out;
5664}
5665
5666
5667
5668template <int dim,
5669 int n_components_,
5670 typename Number,
5671 bool is_face,
5672 typename VectorizedArrayType>
5675 get_normal_derivative(const unsigned int q_point) const
5676{
5677 AssertIndexRange(q_point, this->n_quadrature_points);
5678# ifdef DEBUG
5679 Assert(this->gradients_quad_initialized == true,
5681# endif
5682
5683 Assert(this->normal_x_jacobian != nullptr,
5685 "update_gradients"));
5686
5687 const std::size_t nqp = this->n_quadrature_points;
5689
5690 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5691 for (unsigned int comp = 0; comp < n_components; comp++)
5692 grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5693 (this->normal_x_jacobian[0][dim - 1]);
5694 else
5695 {
5696 const std::size_t index =
5697 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5698 for (unsigned int comp = 0; comp < n_components; comp++)
5699 {
5700 grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5701 this->normal_x_jacobian[index][0];
5702 for (unsigned int d = 1; d < dim; ++d)
5703 grad_out[comp] += gradients_quad[(comp * dim + d) * nqp + q_point] *
5704 this->normal_x_jacobian[index][d];
5705 }
5706 }
5707 return grad_out;
5708}
5709
5710
5711
5712namespace internal
5713{
5714 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5715 // store the lower diagonal because of symmetry
5716 template <typename VectorizedArrayType>
5717 inline void
5718 hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
5719 const VectorizedArrayType *const hessians,
5720 const unsigned int,
5721 VectorizedArrayType (&tmp)[1][1])
5722 {
5723 tmp[0][0] = jac[0][0] * hessians[0];
5724 }
5725
5726 template <typename VectorizedArrayType>
5727 inline void
5728 hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
5729 const VectorizedArrayType *const hessians,
5730 const unsigned int nqp,
5731 VectorizedArrayType (&tmp)[2][2])
5732 {
5733 for (unsigned int d = 0; d < 2; ++d)
5734 {
5735 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
5736 tmp[1][d] =
5737 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
5738 }
5739 }
5740
5741 template <typename VectorizedArrayType>
5742 inline void
5743 hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
5744 const VectorizedArrayType *const hessians,
5745 const unsigned int nqp,
5746 VectorizedArrayType (&tmp)[3][3])
5747 {
5748 for (unsigned int d = 0; d < 3; ++d)
5749 {
5750 tmp[0][d] =
5751 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
5752 jac[d][2] * hessians[4 * nqp]);
5753 tmp[1][d] =
5754 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
5755 jac[d][2] * hessians[5 * nqp]);
5756 tmp[2][d] =
5757 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
5758 jac[d][2] * hessians[2 * nqp]);
5759 }
5760 }
5761} // namespace internal
5762
5763
5764
5765template <int dim,
5766 int n_components_,
5767 typename Number,
5768 bool is_face,
5769 typename VectorizedArrayType>
5772 get_hessian(const unsigned int q_point) const
5773{
5774 Assert(!is_face, ExcNotImplemented());
5775# ifdef DEBUG
5776 Assert(this->hessians_quad_initialized == true,
5778# endif
5779 AssertIndexRange(q_point, this->n_quadrature_points);
5780
5781 Assert(this->jacobian != nullptr,
5783 "update_hessian"));
5785 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5786 0 :
5787 q_point];
5788
5790
5791 const std::size_t nqp = this->n_quadrature_points;
5792 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5793
5794 // Cartesian cell
5795 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5796 {
5797 for (unsigned int comp = 0; comp < n_components; comp++)
5798 {
5799 for (unsigned int d = 0; d < dim; ++d)
5800 hessian_out[comp][d][d] =
5801 hessians_quad[(comp * hdim + d) * nqp + q_point] *
5802 (jac[d][d] * jac[d][d]);
5803 switch (dim)
5804 {
5805 case 1:
5806 break;
5807 case 2:
5808 hessian_out[comp][0][1] =
5809 hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5810 (jac[0][0] * jac[1][1]);
5811 break;
5812 case 3:
5813 hessian_out[comp][0][1] =
5814 hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5815 (jac[0][0] * jac[1][1]);
5816 hessian_out[comp][0][2] =
5817 hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5818 (jac[0][0] * jac[2][2]);
5819 hessian_out[comp][1][2] =
5820 hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5821 (jac[1][1] * jac[2][2]);
5822 break;
5823 default:
5824 Assert(false, ExcNotImplemented());
5825 }
5826 for (unsigned int d = 0; d < dim; ++d)
5827 for (unsigned int e = d + 1; e < dim; ++e)
5828 hessian_out[comp][e][d] = hessian_out[comp][d][e];
5829 }
5830 }
5831 // cell with general Jacobian, but constant within the cell
5832 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5833 {
5834 for (unsigned int comp = 0; comp < n_components; comp++)
5835 {
5836 VectorizedArrayType tmp[dim][dim];
5837 internal::hessian_unit_times_jac(
5838 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5839
5840 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5841 for (unsigned int d = 0; d < dim; ++d)
5842 for (unsigned int e = d; e < dim; ++e)
5843 {
5844 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5845 for (unsigned int f = 1; f < dim; ++f)
5846 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5847 }
5848
5849 // no J' * grad(u) part here because the Jacobian is constant
5850 // throughout the cell and hence, its derivative is zero
5851
5852 // take symmetric part
5853 for (unsigned int d = 0; d < dim; ++d)
5854 for (unsigned int e = d + 1; e < dim; ++e)
5855 hessian_out[comp][e][d] = hessian_out[comp][d][e];
5856 }
5857 }
5858 // cell with general Jacobian
5859 else
5860 {
5861 const auto &jac_grad =
5862 this->mapping_data->jacobian_gradients
5863 [1 - this->is_interior_face]
5864 [this->mapping_data->data_index_offsets[this->cell] + q_point];
5865 for (unsigned int comp = 0; comp < n_components; comp++)
5866 {
5867 // compute laplacian before the gradient because it needs to access
5868 // unscaled gradient data
5869 VectorizedArrayType tmp[dim][dim];
5870 internal::hessian_unit_times_jac(
5871 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5872
5873 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5874 for (unsigned int d = 0; d < dim; ++d)
5875 for (unsigned int e = d; e < dim; ++e)
5876 {
5877 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5878 for (unsigned int f = 1; f < dim; ++f)
5879 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5880 }
5881
5882 // add diagonal part of J' * grad(u)
5883 for (unsigned int d = 0; d < dim; ++d)
5884 for (unsigned int e = 0; e < dim; ++e)
5885 hessian_out[comp][d][d] +=
5886 jac_grad[d][e] *
5887 gradients_quad[(comp * dim + e) * nqp + q_point];
5888
5889 // add off-diagonal part of J' * grad(u)
5890 for (unsigned int d = 0, count = dim; d < dim; ++d)
5891 for (unsigned int e = d + 1; e < dim; ++e, ++count)
5892 for (unsigned int f = 0; f < dim; ++f)
5893 hessian_out[comp][d][e] +=
5894 jac_grad[count][f] *
5895 gradients_quad[(comp * dim + f) * nqp + q_point];
5896
5897 // take symmetric part
5898 for (unsigned int d = 0; d < dim; ++d)
5899 for (unsigned int e = d + 1; e < dim; ++e)
5900 hessian_out[comp][e][d] = hessian_out[comp][d][e];
5901 }
5902 }
5903 return hessian_out;
5904}
5905
5906
5907
5908template <int dim,
5909 int n_components_,
5910 typename Number,
5911 bool is_face,
5912 typename VectorizedArrayType>
5915 get_hessian_diagonal(const unsigned int q_point) const
5916{
5917 Assert(!is_face, ExcNotImplemented());
5918# ifdef DEBUG
5919 Assert(this->hessians_quad_initialized == true,
5921# endif
5922 AssertIndexRange(q_point, this->n_quadrature_points);
5923
5924 Assert(this->jacobian != nullptr, ExcNotImplemented());
5926 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5927 0 :
5928 q_point];
5929
5930 const std::size_t nqp = this->n_quadrature_points;
5931 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5933
5934 // Cartesian cell
5935 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5936 {
5937 for (unsigned int comp = 0; comp < n_components; comp++)
5938 for (unsigned int d = 0; d < dim; ++d)
5939 hessian_out[comp][d] =
5940 hessians_quad[(comp * hdim + d) * nqp + q_point] *
5941 (jac[d][d] * jac[d][d]);
5942 }
5943 // cell with general Jacobian, but constant within the cell
5944 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5945 {
5946 for (unsigned int comp = 0; comp < n_components; comp++)
5947 {
5948 // compute laplacian before the gradient because it needs to access
5949 // unscaled gradient data
5950 VectorizedArrayType tmp[dim][dim];
5951 internal::hessian_unit_times_jac(
5952 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5953
5954 // compute only the trace part of hessian, J * tmp = J *
5955 // hess_unit(u) * J^T
5956 for (unsigned int d = 0; d < dim; ++d)
5957 {
5958 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5959 for (unsigned int f = 1; f < dim; ++f)
5960 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5961 }
5962 }
5963 }
5964 // cell with general Jacobian
5965 else
5966 {
5967 const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5968 &jac_grad =
5969 this->mapping_data->jacobian_gradients
5970 [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5971 for (unsigned int comp = 0; comp < n_components; comp++)
5972 {
5973 // compute laplacian before the gradient because it needs to access
5974 // unscaled gradient data
5975 VectorizedArrayType tmp[dim][dim];
5976 internal::hessian_unit_times_jac(
5977 jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5978
5979 // compute only the trace part of hessian, J * tmp = J *
5980 // hess_unit(u) * J^T
5981 for (unsigned int d = 0; d < dim; ++d)
5982 {
5983 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5984 for (unsigned int f = 1; f < dim; ++f)
5985 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5986 }
5987
5988 for (unsigned int d = 0; d < dim; ++d)
5989 for (unsigned int e = 0; e < dim; ++e)
5990 hessian_out[comp][d] +=
5991 jac_grad[d][e] *
5992 gradients_quad[(comp * dim + e) * nqp + q_point];
5993 }
5994 }
5995 return hessian_out;
5996}
5997
5998
5999
6000template <int dim,
6001 int n_components_,
6002 typename Number,
6003 bool is_face,
6004 typename VectorizedArrayType>
6007 get_laplacian(const unsigned int q_point) const
6008{
6009 Assert(is_face == false, ExcNotImplemented());
6010# ifdef DEBUG
6011 Assert(this->hessians_quad_initialized == true,
6013# endif
6014 AssertIndexRange(q_point, this->n_quadrature_points);
6015
6017 const auto hess_diag = get_hessian_diagonal(q_point);
6018 for (unsigned int comp = 0; comp < n_components; ++comp)
6019 {
6020 laplacian_out[comp] = hess_diag[comp][0];
6021 for (unsigned int d = 1; d < dim; ++d)
6022 laplacian_out[comp] += hess_diag[comp][d];
6023 }
6024 return laplacian_out;
6025}
6026
6027
6028
6029template <int dim,
6030 int n_components_,
6031 typename Number,
6032 bool is_face,
6033 typename VectorizedArrayType>
6034inline DEAL_II_ALWAYS_INLINE void
6037 const unsigned int dof)
6038{
6039# ifdef DEBUG
6040 this->dof_values_initialized = true;
6041# endif
6042 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6043 for (unsigned int comp = 0; comp < n_components; comp++)
6044 this->values_dofs[comp][dof] = val_in[comp];
6045}
6046
6047
6048
6049template <int dim,
6050 int n_components_,
6051 typename Number,
6052 bool is_face,
6053 typename VectorizedArrayType>
6054inline DEAL_II_ALWAYS_INLINE void
6057 const unsigned int q_point)
6058{
6060 AssertIndexRange(q_point, this->n_quadrature_points);
6061 Assert(this->J_value != nullptr,
6063 "update_values"));
6064# ifdef DEBUG
6065 this->values_quad_submitted = true;
6066# endif
6067
6068 const std::size_t nqp = this->n_quadrature_points;
6069 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6070 {
6071 const VectorizedArrayType JxW =
6072 this->J_value[0] * this->quadrature_weights[q_point];
6073 for (unsigned int comp = 0; comp < n_components; ++comp)
6074 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6075 }
6076 else
6077 {
6078 const VectorizedArrayType JxW = this->J_value[q_point];
6079 for (unsigned int comp = 0; comp < n_components; ++comp)
6080 values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6081 }
6082}
6083
6084
6085
6086template <int dim,
6087 int n_components_,
6088 typename Number,
6089 bool is_face,
6090 typename VectorizedArrayType>
6091inline DEAL_II_ALWAYS_INLINE void
6094 const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
6095 const unsigned int q_point)
6096{
6098 AssertIndexRange(q_point, this->n_quadrature_points);
6099 Assert(this->J_value != nullptr,
6101 "update_gradients"));
6102 Assert(this->jacobian != nullptr,
6104 "update_gradients"));
6105# ifdef DEBUG
6106 this->gradients_quad_submitted = true;
6107# endif
6108
6109 const std::size_t nqp = this->n_quadrature_points;
6110 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6111 {
6112 const VectorizedArrayType JxW =
6113 this->J_value[0] * this->quadrature_weights[q_point];
6114 for (unsigned int d = 0; d < dim; ++d)
6115 {
6116 const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
6117 for (unsigned int comp = 0; comp < n_components; comp++)
6118 gradients_quad[(comp * dim + d) * nqp + q_point] =
6119 grad_in[comp][d] * factor;
6120 }
6121 }
6122 else
6123 {
6125 this->cell_type > internal::MatrixFreeFunctions::affine ?
6126 this->jacobian[q_point] :
6127 this->jacobian[0];
6128 const VectorizedArrayType JxW =
6129 this->cell_type > internal::MatrixFreeFunctions::affine ?
6130 this->J_value[q_point] :
6131 this->J_value[0] * this->quadrature_weights[q_point];
6132 for (unsigned int comp = 0; comp < n_components; ++comp)
6133 for (unsigned int d = 0; d < dim; ++d)
6134 {
6135 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
6136 for (unsigned int e = 1; e < dim; ++e)
6137 new_val += (jac[e][d] * grad_in[comp][e]);
6138 gradients_quad[(comp * dim + d) * nqp + q_point] = new_val * JxW;
6139 }
6140 }
6141}
6142
6143
6144
6145template <int dim,
6146 int n_components_,
6147 typename Number,
6148 bool is_face,
6149 typename VectorizedArrayType>
6150inline DEAL_II_ALWAYS_INLINE void
6154 const unsigned int q_point)
6155{
6156 AssertIndexRange(q_point, this->n_quadrature_points);
6157 Assert(this->normal_x_jacobian != nullptr,
6159 "update_gradients"));
6160# ifdef DEBUG
6161 this->gradients_quad_submitted = true;
6162# endif
6163
6164 const std::size_t nqp = this->n_quadrature_points;
6165 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
6166 for (unsigned int comp = 0; comp < n_components; comp++)
6167 {
6168 for (unsigned int d = 0; d < dim - 1; ++d)
6169 gradients_quad[(comp * dim + d) * nqp + q_point] =
6170 VectorizedArrayType();
6171 gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
6172 grad_in[comp] *
6173 (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
6174 this->quadrature_weights[q_point]);
6175 }
6176 else
6177 {
6178 const unsigned int index =
6179 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
6181 this->normal_x_jacobian[index];
6182 for (unsigned int comp = 0; comp < n_components; comp++)
6183 {
6184 VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
6185 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6186 factor = factor * this->quadrature_weights[q_point];
6187 for (unsigned int d = 0; d < dim; ++d)
6188 gradients_quad[(comp * dim + d) * nqp + q_point] = factor * jac[d];
6189 }
6190 }
6191}
6192
6193
6194
6195template <int dim,
6196 int n_components_,
6197 typename Number,
6198 bool is_face,
6199 typename VectorizedArrayType>
6202 integrate_value() const
6203{
6205# ifdef DEBUG
6206 Assert(this->values_quad_submitted == true,
6208# endif
6209
6211 const std::size_t nqp = this->n_quadrature_points;
6212 for (unsigned int q = 0; q < nqp; ++q)
6213 for (unsigned int comp = 0; comp < n_components; ++comp)
6214 return_value[comp] += this->values_quad[comp * nqp + q];
6215 return (return_value);
6216}
6217
6218
6219
6220/*----------------------- FEEvaluationAccess --------------------------------*/
6221
6222
6223template <int dim,
6224 int n_components_,
6225 typename Number,
6226 bool is_face,
6227 typename VectorizedArrayType>
6228inline FEEvaluationAccess<dim,
6229 n_components_,
6230 Number,
6231 is_face,
6232 VectorizedArrayType>::
6233 FEEvaluationAccess(
6235 const unsigned int dof_no,
6236 const unsigned int first_selected_component,
6237 const unsigned int quad_no_in,
6238 const unsigned int fe_degree,
6239 const unsigned int n_q_points,
6240 const bool is_interior_face,
6241 const unsigned int active_fe_index,
6242 const unsigned int active_quad_index,
6243 const unsigned int face_type)
6244 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6245 data_in,
6246 dof_no,
6247 first_selected_component,
6248 quad_no_in,
6249 fe_degree,
6250 n_q_points,
6251 is_interior_face,
6252 active_fe_index,
6253 active_quad_index,
6254 face_type)
6255{}
6256
6257
6258
6259template <int dim,
6260 int n_components_,
6261 typename Number,
6262 bool is_face,
6263 typename VectorizedArrayType>
6264inline FEEvaluationAccess<dim,
6265 n_components_,
6266 Number,
6267 is_face,
6268 VectorizedArrayType>::
6269 FEEvaluationAccess(
6270 const Mapping<dim> & mapping,
6271 const FiniteElement<dim> &fe,
6272 const Quadrature<1> & quadrature,
6273 const UpdateFlags update_flags,
6274 const unsigned int first_selected_component,
6276 *other)
6277 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6278 mapping,
6279 fe,
6280 quadrature,
6281 update_flags,
6282 first_selected_component,
6283 other)
6284{}
6285
6286
6287
6288template <int dim,
6289 int n_components_,
6290 typename Number,
6291 bool is_face,
6292 typename VectorizedArrayType>
6293inline FEEvaluationAccess<dim,
6294 n_components_,
6295 Number,
6296 is_face,
6297 VectorizedArrayType>::
6298 FEEvaluationAccess(const FEEvaluationAccess<dim,
6299 n_components_,
6300 Number,
6301 is_face,
6302 VectorizedArrayType> &other)
6303 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
6304 other)
6305{}
6306
6307
6308
6309template <int dim,
6310 int n_components_,
6311 typename Number,
6312 bool is_face,
6313 typename VectorizedArrayType>
6314inline FEEvaluationAccess<dim,
6315 n_components_,
6316 Number,
6317 is_face,
6318 VectorizedArrayType> &
6321 n_components_,
6322 Number,
6323 is_face,
6324 VectorizedArrayType> &other)
6325{
6326 this->FEEvaluationBase<dim,
6327 n_components_,
6328 Number,
6329 is_face,
6330 VectorizedArrayType>::operator=(other);
6331 return *this;
6332}
6333
6334
6335
6336/*-------------------- FEEvaluationAccess scalar ----------------------------*/
6337
6338
6339template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6343 const unsigned int dof_no,
6344 const unsigned int first_selected_component,
6345 const unsigned int quad_no_in,
6346 const unsigned int fe_degree,
6347 const unsigned int n_q_points,
6348 const bool is_interior_face,
6349 const unsigned int active_fe_index,
6350 const unsigned int active_quad_index,
6351 const unsigned int face_type)
6352 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
6353 data_in,
6354 dof_no,
6355 first_selected_component,
6356 quad_no_in,
6357 fe_degree,
6358 n_q_points,
6359 is_interior_face,
6360 active_fe_index,
6361 active_quad_index,
6362 face_type)
6363{}
6364
6365
6366
6367template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6370 const Mapping<dim> & mapping,
6371 const FiniteElement<dim> &fe,
6372 const Quadrature<1> & quadrature,
6373 const UpdateFlags update_flags,
6374 const unsigned int first_selected_component,
6376 *other)
6377 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
6378 mapping,
6379 fe,
6380 quadrature,
6381 update_flags,
6382 first_selected_component,
6383 other)
6384{}
6385
6386
6387
6388template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6392 &other)
6393 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
6394{}
6395
6396
6397
6398template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6402{
6404 operator=(other);
6405 return *this;
6406}
6407
6408
6409
6410template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6411inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6413 const unsigned int dof) const
6414{
6415 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6416 return this->values_dofs[0][dof];
6417}
6418
6419
6420
6421template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6422inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6424 const unsigned int q_point) const
6425{
6426# ifdef DEBUG
6427 Assert(this->values_quad_initialized == true,
6429# endif
6430 AssertIndexRange(q_point, this->n_quadrature_points);
6431 return this->values_quad[q_point];
6432}
6433
6434
6435
6436template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6437inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6439 get_normal_derivative(const unsigned int q_point) const
6440{
6441 return BaseClass::get_normal_derivative(q_point)[0];
6442}
6443
6444
6445
6446template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6449 const unsigned int q_point) const
6450{
6451 // could use the base class gradient, but that involves too many expensive
6452 // initialization operations on tensors
6453
6454# ifdef DEBUG
6455 Assert(this->gradients_quad_initialized == true,
6457# endif
6458 AssertIndexRange(q_point, this->n_quadrature_points);
6459
6460 Assert(this->jacobian != nullptr,
6462 "update_gradients"));
6463
6465
6466 const std::size_t nqp = this->n_quadrature_points;
6467 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6468 {
6469 for (unsigned int d = 0; d < dim; ++d)
6470 grad_out[d] =
6471 this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
6472 }
6473 // cell with general/affine Jacobian
6474 else
6475 {
6477 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
6478 q_point :
6479 0];
6480 for (unsigned int d = 0; d < dim; ++d)
6481 {
6482 grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
6483 for (unsigned int e = 1; e < dim; ++e)
6484 grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
6485 }
6486 }
6487 return grad_out;
6488}
6489
6490
6491
6492template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6495 const unsigned int q_point) const
6496{
6497 return BaseClass::get_hessian(q_point)[0];
6498}
6499
6500
6501
6502template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6505 get_hessian_diagonal(const unsigned int q_point) const
6506{
6507 return BaseClass::get_hessian_diagonal(q_point)[0];
6508}
6509
6510
6511
6512template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6513inline VectorizedArrayType
6515 const unsigned int q_point) const
6516{
6517 return BaseClass::get_laplacian(q_point)[0];
6518}
6519
6520
6521
6522template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6523inline void DEAL_II_ALWAYS_INLINE
6525 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6526{
6527# ifdef DEBUG
6528 this->dof_values_initialized = true;
6529 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6530# endif
6531 this->values_dofs[0][dof] = val_in;
6532}
6533
6534
6535
6536template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6537inline void DEAL_II_ALWAYS_INLINE
6539 const VectorizedArrayType val_in,
6540 const unsigned int q_point)
6541{
6543 AssertIndexRange(q_point, this->n_quadrature_points);
6544 Assert(this->J_value != nullptr,
6546 "update_value"));
6547# ifdef DEBUG
6548 this->values_quad_submitted = true;
6549# endif
6550
6551 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6552 {
6553 const VectorizedArrayType JxW =
6554 this->J_value[0] * this->quadrature_weights[q_point];
6555 this->values_quad[q_point] = val_in * JxW;
6556 }
6557 else // if (this->cell_type < internal::MatrixFreeFunctions::general)
6558 {
6559 this->values_quad[q_point] = val_in * this->J_value[q_point];
6560 }
6561}
6562
6563
6564
6565template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6566inline DEAL_II_ALWAYS_INLINE void
6569 const unsigned int q_point)
6570{
6571 submit_value(val_in[0], q_point);
6572}
6573
6574
6575
6576template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6577inline DEAL_II_ALWAYS_INLINE void
6579 submit_normal_derivative(const VectorizedArrayType grad_in,
6580 const unsigned int q_point)
6581{
6583 grad[0] = grad_in;
6584 BaseClass::submit_normal_derivative(grad, q_point);
6585}
6586
6587
6588
6589template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6590inline DEAL_II_ALWAYS_INLINE void
6593 const unsigned int q_point)
6594{
6596 AssertIndexRange(q_point, this->n_quadrature_points);
6597 Assert(this->J_value != nullptr,
6599 "update_gradients"));
6600 Assert(this->jacobian != nullptr,
6602 "update_gradients"));
6603# ifdef DEBUG
6604 this->gradients_quad_submitted = true;
6605# endif
6606
6607 const std::size_t nqp = this->n_quadrature_points;
6608 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6609 {
6610 const VectorizedArrayType JxW =
6611 this->J_value[0] * this->quadrature_weights[q_point];
6612 for (unsigned int d = 0; d < dim; ++d)
6613 this->gradients_quad[d * nqp + q_point] =
6614 (grad_in[d] * this->jacobian[0][d][d] * JxW);
6615 }
6616 // general/affine cell type
6617 else
6618 {
6620 this->cell_type > internal::MatrixFreeFunctions::affine ?
6621 this->jacobian[q_point] :
6622 this->jacobian[0];
6623 const VectorizedArrayType JxW =
6624 this->cell_type > internal::MatrixFreeFunctions::affine ?
6625 this->J_value[q_point] :
6626 this->J_value[0] * this->quadrature_weights[q_point];
6627 for (unsigned int d = 0; d < dim; ++d)
6628 {
6629 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
6630 for (unsigned int e = 1; e < dim; ++e)
6631 new_val += jac[e][d] * grad_in[e];
6632 this->gradients_quad[d * nqp + q_point] = new_val * JxW;
6633 }
6634 }
6635}
6636
6637
6638
6639template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6640inline VectorizedArrayType
6642 integrate_value() const
6643{
6644 return BaseClass::integrate_value()[0];
6645}
6646
6647
6648
6649/*----------------- FEEvaluationAccess vector-valued ------------------------*/
6650
6651
6652template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6656 const unsigned int dof_no,
6657 const unsigned int first_selected_component,
6658 const unsigned int quad_no_in,
6659 const unsigned int fe_degree,
6660 const unsigned int n_q_points,
6661 const bool is_interior_face,
6662 const unsigned int active_fe_index,
6663 const unsigned int active_quad_index,
6664 const unsigned int face_type)
6665 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
6666 data_in,
6667 dof_no,
6668 first_selected_component,
6669 quad_no_in,
6670 fe_degree,
6671 n_q_points,
6672 is_interior_face,
6673 active_fe_index,
6674 active_quad_index,
6675 face_type)
6676{}
6677
6678
6679
6680template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6683 const Mapping<dim> & mapping,
6684 const FiniteElement<dim> &fe,
6685 const Quadrature<1> & quadrature,
6686 const UpdateFlags update_flags,
6687 const unsigned int first_selected_component,
6689 *other)
6690 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
6691 mapping,
6692 fe,
6693 quadrature,
6694 update_flags,
6695 first_selected_component,
6696 other)
6697{}
6698
6699
6700
6701template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6705 &other)
6706 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
6707{}
6708
6709
6710
6711template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6715 &other)
6716{
6718 operator=(other);
6719 return *this;
6720}
6721
6722
6723
6724template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6727 get_gradient(const unsigned int q_point) const
6728{
6729 return BaseClass::get_gradient(q_point);
6730}
6731
6732
6733
6734template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6735inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6737 get_divergence(const unsigned int q_point) const
6738{
6739# ifdef DEBUG
6740 Assert(this->gradients_quad_initialized == true,
6742# endif
6743 AssertIndexRange(q_point, this->n_quadrature_points);
6744 Assert(this->jacobian != nullptr,
6746 "update_gradients"));
6747
6748 VectorizedArrayType divergence;
6749 const std::size_t nqp = this->n_quadrature_points;
6750
6751 // Cartesian cell
6752 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6753 {
6754 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6755 for (unsigned int d = 1; d < dim; ++d)
6756 divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
6757 this->jacobian[0][d][d];
6758 }
6759 // cell with general/constant Jacobian
6760 else
6761 {
6763 this->cell_type == internal::MatrixFreeFunctions::general ?
6764 this->jacobian[q_point] :
6765 this->jacobian[0];
6766 divergence = jac[0][0] * this->gradients_quad[q_point];
6767 for (unsigned int e = 1; e < dim; ++e)
6768 divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
6769 for (unsigned int d = 1; d < dim; ++d)
6770 for (unsigned int e = 0; e < dim; ++e)
6771 divergence +=
6772 jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
6773 }
6774 return divergence;
6775}
6776
6777
6778
6779template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6782 get_symmetric_gradient(const unsigned int q_point) const
6783{
6784 // copy from generic function into dim-specialization function
6785 const auto grad = get_gradient(q_point);
6786 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6787 VectorizedArrayType half = Number(0.5);
6788 for (unsigned int d = 0; d < dim; ++d)
6789 symmetrized[d] = grad[d][d];
6790 switch (dim)
6791 {
6792 case 1:
6793 break;
6794 case 2:
6795 symmetrized[2] = grad[0][1] + grad[1][0];
6796 symmetrized[2] *= half;
6797 break;
6798 case 3:
6799 symmetrized[3] = grad[0][1] + grad[1][0];
6800 symmetrized[3] *= half;
6801 symmetrized[4] = grad[0][2] + grad[2][0];
6802 symmetrized[4] *= half;
6803 symmetrized[5] = grad[1][2] + grad[2][1];
6804 symmetrized[5] *= half;
6805 break;
6806 default:
6807 Assert(false, ExcNotImplemented());
6808 }
6810}
6811
6812
6813
6814template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6816 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6818 const unsigned int q_point) const
6819{
6820 // copy from generic function into dim-specialization function
6821 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6822 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6823 switch (dim)
6824 {
6825 case 1:
6826 Assert(false,
6827 ExcMessage(
6828 "Computing the curl in 1d is not a useful operation"));
6829 break;
6830 case 2:
6831 curl[0] = grad[1][0] - grad[0][1];
6832 break;
6833 case 3:
6834 curl[0] = grad[2][1] - grad[1][2];
6835 curl[1] = grad[0][2] - grad[2][0];
6836 curl[2] = grad[1][0] - grad[0][1];
6837 break;
6838 default:
6839 Assert(false, ExcNotImplemented());
6840 }
6841 return curl;
6842}
6843
6844
6845
6846template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6849 get_hessian_diagonal(const unsigned int q_point) const
6850{
6851 return BaseClass::get_hessian_diagonal(q_point);
6852}
6853
6854
6855
6856template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6859 const unsigned int q_point) const
6860{
6861# ifdef DEBUG
6862 Assert(this->hessians_quad_initialized == true,
6864# endif
6865 AssertIndexRange(q_point, this->n_quadrature_points);
6866 return BaseClass::get_hessian(q_point);
6867}
6868
6869
6870
6871template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6872inline DEAL_II_ALWAYS_INLINE void
6875 const unsigned int q_point)
6876{
6877 BaseClass::submit_gradient(grad_in, q_point);
6878}
6879
6880
6881
6882template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6883inline DEAL_II_ALWAYS_INLINE void
6886 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6887 const unsigned int q_point)
6888{
6889 BaseClass::submit_gradient(grad_in, q_point);
6890}
6891
6892
6893
6894template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6895inline DEAL_II_ALWAYS_INLINE void
6897 submit_divergence(const VectorizedArrayType div_in,
6898 const unsigned int q_point)
6899{
6901 AssertIndexRange(q_point, this->n_quadrature_points);
6902 Assert(this->J_value != nullptr,
6904 "update_gradients"));
6905 Assert(this->jacobian != nullptr,
6907 "update_gradients"));
6908# ifdef DEBUG
6909 this->gradients_quad_submitted = true;
6910# endif
6911
6912 const std::size_t nqp = this->n_quadrature_points;
6913 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6914 {
6915 const VectorizedArrayType fac =
6916 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6917 for (unsigned int d = 0; d < dim; ++d)
6918 {
6919 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6920 (fac * this->jacobian[0][d][d]);
6921 for (unsigned int e = d + 1; e < dim; ++e)
6922 {
6923 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6924 VectorizedArrayType();
6925 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6926 VectorizedArrayType();
6927 }
6928 }
6929 }
6930 else
6931 {
6933 this->cell_type == internal::MatrixFreeFunctions::general ?
6934 this->jacobian[q_point] :
6935 this->jacobian[0];
6936 const VectorizedArrayType fac =
6937 (this->cell_type == internal::MatrixFreeFunctions::general ?
6938 this->J_value[q_point] :
6939 this->J_value[0] * this->quadrature_weights[q_point]) *
6940 div_in;
6941 for (unsigned int d = 0; d < dim; ++d)
6942 {
6943 for (unsigned int e = 0; e < dim; ++e)
6944 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6945 jac[d][e] * fac;
6946 }
6947 }
6948}
6949
6950
6951
6952template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6953inline DEAL_II_ALWAYS_INLINE void
6957 const unsigned int q_point)
6958{
6959 // could have used base class operator, but that involves some overhead
6960 // which is inefficient. it is nice to have the symmetric tensor because
6961 // that saves some operations
6963 AssertIndexRange(q_point, this->n_quadrature_points);
6964 Assert(this->J_value != nullptr,
6966 "update_gradients"));
6967 Assert(this->jacobian != nullptr,
6969 "update_gradients"));
6970# ifdef DEBUG
6971 this->gradients_quad_submitted = true;
6972# endif
6973
6974 const std::size_t nqp = this->n_quadrature_points;
6975 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6976 {
6977 const VectorizedArrayType JxW =
6978 this->J_value[0] * this->quadrature_weights[q_point];
6979 for (unsigned int d = 0; d < dim; ++d)
6980 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6981 (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6982 for (unsigned int e = 0, counter = dim; e < dim; ++e)
6983 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6984 {
6985 const VectorizedArrayType value =
6986 sym_grad.access_raw_entry(counter) * JxW;
6987 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6988 value * this->jacobian[0][d][d];
6989 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6990 value * this->jacobian[0][e][e];
6991 }
6992 }
6993 // general/affine cell type
6994 else
6995 {
6996 const VectorizedArrayType JxW =
6997 this->cell_type == internal::MatrixFreeFunctions::general ?
6998 this->J_value[q_point] :
6999 this->J_value[0] * this->quadrature_weights[q_point];
7001 this->cell_type == internal::MatrixFreeFunctions::general ?
7002 this->jacobian[q_point] :
7003 this->jacobian[0];
7004 VectorizedArrayType weighted[dim][dim];
7005 for (unsigned int i = 0; i < dim; ++i)
7006 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
7007 for (unsigned int i = 0, counter = dim; i < dim; ++i)
7008 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
7009 {
7010 const VectorizedArrayType value =
7011 sym_grad.access_raw_entry(counter) * JxW;
7012 weighted[i][j] = value;
7013 weighted[j][i] = value;
7014 }
7015 for (unsigned int comp = 0; comp < dim; ++comp)
7016 for (unsigned int d = 0; d < dim; ++d)
7017 {
7018 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
7019 for (unsigned int e = 1; e < dim; ++e)
7020 new_val += jac[e][d] * weighted[comp][e];
7021 this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
7022 }
7023 }
7024}
7025
7026
7027
7028template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7029inline DEAL_II_ALWAYS_INLINE void
7032 const unsigned int q_point)
7033{
7035 switch (dim)
7036 {
7037 case 1:
7038 Assert(false,
7039 ExcMessage(
7040 "Testing by the curl in 1d is not a useful operation"));
7041 break;
7042 case 2:
7043 grad[1][0] = curl[0];
7044 grad[0][1] = -curl[0];
7045 break;
7046 case 3:
7047 grad[2][1] = curl[0];
7048 grad[1][2] = -curl[0];
7049 grad[0][2] = curl[1];
7050 grad[2][0] = -curl[1];
7051 grad[1][0] = curl[2];
7052 grad[0][1] = -curl[2];
7053 break;
7054 default:
7055 Assert(false, ExcNotImplemented());
7056 }
7057 submit_gradient(grad, q_point);
7058}
7059
7060
7061/*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
7062
7063
7064template <typename Number, bool is_face, typename VectorizedArrayType>
7067 const unsigned int dof_no,
7068 const unsigned int first_selected_component,
7069 const unsigned int quad_no_in,
7070 const unsigned int fe_degree,
7071 const unsigned int n_q_points,
7072 const bool is_interior_face,
7073 const unsigned int active_fe_index,
7074 const unsigned int active_quad_index,
7075 const unsigned int face_type)
7076 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
7077 data_in,
7078 dof_no,
7079 first_selected_component,
7080 quad_no_in,
7081 fe_degree,
7082 n_q_points,
7083 is_interior_face,
7084 active_fe_index,
7085 active_quad_index,
7086 face_type)
7087{}
7088
7089
7090
7091template <typename Number, bool is_face, typename VectorizedArrayType>
7094 const Mapping<1> & mapping,
7095 const FiniteElement<1> &fe,
7096 const Quadrature<1> & quadrature,
7097 const UpdateFlags update_flags,
7098 const unsigned int first_selected_component,
7100 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
7101 mapping,
7102 fe,
7103 quadrature,
7104 update_flags,
7105 first_selected_component,
7106 other)
7107{}
7108
7109
7110
7111template <typename Number, bool is_face, typename VectorizedArrayType>
7115 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
7116{}
7117
7118
7119
7120template <typename Number, bool is_face, typename VectorizedArrayType>
7124{
7126 other);
7127 return *this;
7128}
7129
7130
7131
7132template <typename Number, bool is_face, typename VectorizedArrayType>
7133inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7135 const unsigned int dof) const
7136{
7137 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7138 return this->values_dofs[0][dof];
7139}
7140
7141
7142
7143template <typename Number, bool is_face, typename VectorizedArrayType>
7144inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7146 const unsigned int q_point) const
7147{
7148# ifdef DEBUG
7149 Assert(this->values_quad_initialized == true,
7151# endif
7152 AssertIndexRange(q_point, this->n_quadrature_points);
7153 return this->values_quad[q_point];
7154}
7155
7156
7157
7158template <typename Number, bool is_face, typename VectorizedArrayType>
7161 const unsigned int q_point) const
7162{
7163 // could use the base class gradient, but that involves too many inefficient
7164 // initialization operations on tensors
7165
7166# ifdef DEBUG
7167 Assert(this->gradients_quad_initialized == true,
7169# endif
7170 AssertIndexRange(q_point, this->n_quadrature_points);
7171
7173 this->cell_type == internal::MatrixFreeFunctions::general ?
7174 this->jacobian[q_point] :
7175 this->jacobian[0];
7176
7178 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
7179
7180 return grad_out;
7181}
7182
7183
7184
7185template <typename Number, bool is_face, typename VectorizedArrayType>
7186inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7188 const unsigned int q_point) const
7189{
7190 return get_gradient(q_point)[0];
7191}
7192
7193
7194
7195template <typename Number, bool is_face, typename VectorizedArrayType>
7196inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7198 get_normal_derivative(const unsigned int q_point) const
7199{
7200 return BaseClass::get_normal_derivative(q_point)[0];
7201}
7202
7203
7204
7205template <typename Number, bool is_face, typename VectorizedArrayType>
7208 const unsigned int q_point) const
7209{
7210 return BaseClass::get_hessian(q_point)[0];
7211}
7212
7213
7214
7215template <typename Number, bool is_face, typename VectorizedArrayType>
7218 get_hessian_diagonal(const unsigned int q_point) const
7219{
7220 return BaseClass::get_hessian_diagonal(q_point)[0];
7221}
7222
7223
7224
7225template <typename Number, bool is_face, typename VectorizedArrayType>
7226inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7228 const unsigned int q_point) const
7229{
7230 return BaseClass::get_laplacian(q_point)[0];
7231}
7232
7233
7234
7235template <typename Number, bool is_face, typename VectorizedArrayType>
7238 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
7239{
7240# ifdef DEBUG
7241 this->dof_values_initialized = true;
7242 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7243# endif
7244 this->values_dofs[0][dof] = val_in;
7245}
7246
7247
7248
7249template <typename Number, bool is_face, typename VectorizedArrayType>
7250inline DEAL_II_ALWAYS_INLINE void
7252 const VectorizedArrayType val_in,
7253 const unsigned int q_point)
7254{
7256 AssertIndexRange(q_point, this->n_quadrature_points);
7257# ifdef DEBUG
7258 this->values_quad_submitted = true;
7259# endif
7260
7261 if (this->cell_type == internal::MatrixFreeFunctions::general)
7262 {
7263 const VectorizedArrayType JxW = this->J_value[q_point];
7264 this->values_quad[q_point] = val_in * JxW;
7265 }
7266 else // if (this->cell_type == internal::MatrixFreeFunctions::general)
7267 {
7268 const VectorizedArrayType JxW =
7269 this->J_value[0] * this->quadrature_weights[q_point];
7270 this->values_quad[q_point] = val_in * JxW;
7271 }
7272}
7273
7274
7275
7276template <typename Number, bool is_face, typename VectorizedArrayType>
7277inline DEAL_II_ALWAYS_INLINE void
7280 const unsigned int q_point)
7281{
7282 submit_value(val_in[0], q_point);
7283}
7284
7285
7286
7287template <typename Number, bool is_face, typename VectorizedArrayType>
7288inline DEAL_II_ALWAYS_INLINE void
7291 const unsigned int q_point)
7292{
7293 submit_gradient(grad_in[0], q_point);
7294}
7295
7296
7297
7298template <typename Number, bool is_face, typename VectorizedArrayType>
7299inline DEAL_II_ALWAYS_INLINE void
7301 const VectorizedArrayType grad_in,
7302 const unsigned int q_point)
7303{
7305 AssertIndexRange(q_point, this->n_quadrature_points);
7306# ifdef DEBUG
7307 this->gradients_quad_submitted = true;
7308# endif
7309
7311 this->cell_type == internal::MatrixFreeFunctions::general ?
7312 this->jacobian[q_point] :
7313 this->jacobian[0];
7314 const VectorizedArrayType JxW =
7315 this->cell_type == internal::MatrixFreeFunctions::general ?
7316 this->J_value[q_point] :
7317 this->J_value[0] * this->quadrature_weights[q_point];
7318
7319 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7320}
7321
7322
7323
7324template <typename Number, bool is_face, typename VectorizedArrayType>
7325inline DEAL_II_ALWAYS_INLINE void
7328 const unsigned int q_point)
7329{
7330 submit_gradient(grad_in[0][0], q_point);
7331}
7332
7333
7334
7335template <typename Number, bool is_face, typename VectorizedArrayType>
7336inline DEAL_II_ALWAYS_INLINE void
7338 submit_normal_derivative(const VectorizedArrayType grad_in,
7339 const unsigned int q_point)
7340{
7342 grad[0] = grad_in;
7343 BaseClass::submit_normal_derivative(grad, q_point);
7344}
7345
7346
7347
7348template <typename Number, bool is_face, typename VectorizedArrayType>
7349inline DEAL_II_ALWAYS_INLINE void
7352 const unsigned int q_point)
7353{
7354 BaseClass::submit_normal_derivative(grad_in, q_point);
7355}
7356
7357
7358
7359template <typename Number, bool is_face, typename VectorizedArrayType>
7360inline VectorizedArrayType
7362 integrate_value() const
7363{
7364 return BaseClass::integrate_value()[0];
7365}
7366
7367
7368
7369/*-------------------------- FEEvaluation -----------------------------------*/
7370
7371
7372template <int dim,
7373 int fe_degree,
7374 int n_q_points_1d,
7375 int n_components_,
7376 typename Number,
7377 typename VectorizedArrayType>
7378inline FEEvaluation<dim,
7379 fe_degree,
7380 n_q_points_1d,
7381 n_components_,
7382 Number,
7383 VectorizedArrayType>::
7384 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
7385 const unsigned int fe_no,
7386 const unsigned int quad_no,
7387 const unsigned int first_selected_component,
7388 const unsigned int active_fe_index,
7389 const unsigned int active_quad_index)
7390 : BaseClass(data_in,
7391 fe_no,
7392 first_selected_component,
7393 quad_no,
7394 fe_degree,
7395 static_n_q_points,
7396 true /*note: this is not a face*/,
7397 active_fe_index,
7398 active_quad_index)
7399 , dofs_per_component(this->data->dofs_per_component_on_cell)
7400 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7401 , n_q_points(this->data->n_q_points)
7402{
7403 check_template_arguments(fe_no, 0);
7404}
7405
7406
7407
7408template <int dim,
7409 int fe_degree,
7410 int n_q_points_1d,
7411 int n_components_,
7412 typename Number,
7413 typename VectorizedArrayType>
7414inline FEEvaluation<dim,
7415 fe_degree,
7416 n_q_points_1d,
7417 n_components_,
7418 Number,
7419 VectorizedArrayType>::
7420 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7421 const std::pair<unsigned int, unsigned int> & range,
7422 const unsigned int dof_no,
7423 const unsigned int quad_no,
7424 const unsigned int first_selected_component)
7425 : FEEvaluation(matrix_free,
7426 dof_no,
7427 quad_no,
7428 first_selected_component,
7429 matrix_free.get_cell_active_fe_index(range))
7430{}
7431
7432
7433
7434template <int dim,
7435 int fe_degree,
7436 int n_q_points_1d,
7437 int n_components_,
7438 typename Number,
7439 typename VectorizedArrayType>
7440inline FEEvaluation<dim,
7441 fe_degree,
7442 n_q_points_1d,
7443 n_components_,
7444 Number,
7445 VectorizedArrayType>::
7446 FEEvaluation(const Mapping<dim> & mapping,
7447 const FiniteElement<dim> &fe,
7448 const Quadrature<1> & quadrature,
7449 const UpdateFlags update_flags,
7450 const unsigned int first_selected_component)
7451 : BaseClass(mapping,
7452 fe,
7453 quadrature,
7454 update_flags,
7455 first_selected_component,
7456 nullptr)
7457 , dofs_per_component(this->data->dofs_per_component_on_cell)
7458 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7459 , n_q_points(this->data->n_q_points)
7460{
7461 check_template_arguments(numbers::invalid_unsigned_int, 0);
7462}
7463
7464
7465
7466template <int dim,
7467 int fe_degree,
7468 int n_q_points_1d,
7469 int n_components_,
7470 typename Number,
7471 typename VectorizedArrayType>
7472inline FEEvaluation<dim,
7473 fe_degree,
7474 n_q_points_1d,
7475 n_components_,
7476 Number,
7477 VectorizedArrayType>::
7478 FEEvaluation(const FiniteElement<dim> &fe,
7479 const Quadrature<1> & quadrature,
7480 const UpdateFlags update_flags,
7481 const unsigned int first_selected_component)
7482 : BaseClass(StaticMappingQ1<dim>::mapping,
7483 fe,
7484 quadrature,
7485 update_flags,
7486 first_selected_component,
7487 nullptr)
7488 , dofs_per_component(this->data->dofs_per_component_on_cell)
7489 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7490 , n_q_points(this->data->n_q_points)
7491{
7492 check_template_arguments(numbers::invalid_unsigned_int, 0);
7493}
7494
7495
7496
7497template <int dim,
7498 int fe_degree,
7499 int n_q_points_1d,
7500 int n_components_,
7501 typename Number,
7502 typename VectorizedArrayType>
7503inline FEEvaluation<dim,
7504 fe_degree,
7505 n_q_points_1d,
7506 n_components_,
7507 Number,
7508 VectorizedArrayType>::
7509 FEEvaluation(
7510 const FiniteElement<dim> & fe,
7512 const unsigned int first_selected_component)
7513 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7514 fe,
7515 other.mapped_geometry->get_quadrature(),
7516 other.mapped_geometry->get_fe_values().get_update_flags(),
7517 first_selected_component,
7518 &other)
7519 , dofs_per_component(this->data->dofs_per_component_on_cell)
7520 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7521 , n_q_points(this->data->n_q_points)
7522{
7523 check_template_arguments(numbers::invalid_unsigned_int, 0);
7524}
7525
7526
7527
7528template <int dim,
7529 int fe_degree,
7530 int n_q_points_1d,
7531 int n_components_,
7532 typename Number,
7533 typename VectorizedArrayType>
7534inline FEEvaluation<dim,
7535 fe_degree,
7536 n_q_points_1d,
7537 n_components_,
7538 Number,
7539 VectorizedArrayType>::FEEvaluation(const FEEvaluation
7540 &other)
7541 : BaseClass(other)
7542 , dofs_per_component(this->data->dofs_per_component_on_cell)
7543 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7544 , n_q_points(this->data->n_q_points)
7545{
7546 check_template_arguments(numbers::invalid_unsigned_int, 0);
7547}
7548
7549
7550
7551template <int dim,
7552 int fe_degree,
7553 int n_q_points_1d,
7554 int n_components_,
7555 typename Number,
7556 typename VectorizedArrayType>
7557inline FEEvaluation<dim,
7558 fe_degree,
7559 n_q_points_1d,
7560 n_components_,
7561 Number,
7562 VectorizedArrayType> &
7563FEEvaluation<dim,
7564 fe_degree,
7565 n_q_points_1d,
7566 n_components_,
7567 Number,
7568 VectorizedArrayType>::operator=(const FEEvaluation &other)
7569{
7570 BaseClass::operator=(other);
7571 check_template_arguments(numbers::invalid_unsigned_int, 0);
7572 return *this;
7573}
7574
7575
7576
7577template <int dim,
7578 int fe_degree,
7579 int n_q_points_1d,
7580 int n_components_,
7581 typename Number,
7582 typename VectorizedArrayType>
7583inline void
7584FEEvaluation<dim,
7585 fe_degree,
7586 n_q_points_1d,
7587 n_components_,
7588 Number,
7589 VectorizedArrayType>::
7590 check_template_arguments(const unsigned int dof_no,
7591 const unsigned int first_selected_component)
7592{
7593 (void)dof_no;
7594 (void)first_selected_component;
7595
7596# ifdef DEBUG
7597 // print error message when the dimensions do not match. Propose a possible
7598 // fix
7599 if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7600 static_cast<unsigned int>(fe_degree) !=
7601 this->data->data.front().fe_degree) ||
7602 n_q_points != this->n_quadrature_points)
7603 {
7604 std::string message =
7605 "-------------------------------------------------------\n";
7606 message += "Illegal arguments in constructor/wrong template arguments!\n";
7607 message += " Called --> FEEvaluation<dim,";
7608 message += Utilities::int_to_string(fe_degree) + ",";
7609 message += Utilities::int_to_string(n_q_points_1d);
7610 message += "," + Utilities::int_to_string(n_components);
7611 message += ",Number>(data";
7612 if (first_selected_component != numbers::invalid_unsigned_int)
7613 {
7614 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7615 message += Utilities::int_to_string(this->quad_no) + ", ";
7616 message += Utilities::int_to_string(first_selected_component);
7617 }
7618 message += ")\n";
7619
7620 // check whether some other vector component has the correct number of
7621 // points
7622 unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7623 proposed_fe_comp = numbers::invalid_unsigned_int,
7624 proposed_quad_comp = numbers::invalid_unsigned_int;
7625 if (dof_no != numbers::invalid_unsigned_int)
7626 {
7627 if (static_cast<unsigned int>(fe_degree) ==
7628 this->data->data.front().fe_degree)
7629 {
7630 proposed_dof_comp = dof_no;
7631 proposed_fe_comp = first_selected_component;
7632 }
7633 else
7634 for (unsigned int no = 0; no < this->matrix_info->n_components();
7635 ++no)
7636 for (unsigned int nf = 0;
7637 nf < this->matrix_info->n_base_elements(no);
7638 ++nf)
7639 if (this->matrix_info
7640 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7641 .data.front()
7642 .fe_degree == static_cast<unsigned int>(fe_degree))
7643 {
7644 proposed_dof_comp = no;
7645 proposed_fe_comp = nf;
7646 break;
7647 }
7648 if (n_q_points ==
7649 this->mapping_data->descriptor[this->active_quad_index]
7650 .n_q_points)
7651 proposed_quad_comp = this->quad_no;
7652 else
7653 for (unsigned int no = 0;
7654 no < this->matrix_info->get_mapping_info().cell_data.size();
7655 ++no)
7656 if (this->matrix_info->get_mapping_info()
7657 .cell_data[no]
7658 .descriptor[this->active_quad_index]
7659 .n_q_points == n_q_points)
7660 {
7661 proposed_quad_comp = no;
7662 break;
7663 }
7664 }
7665 if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7666 proposed_quad_comp != numbers::invalid_unsigned_int)
7667 {
7668 if (proposed_dof_comp != first_selected_component)
7669 message += "Wrong vector component selection:\n";
7670 else
7671 message += "Wrong quadrature formula selection:\n";
7672 message += " Did you mean FEEvaluation<dim,";
7673 message += Utilities::int_to_string(fe_degree) + ",";
7674 message += Utilities::int_to_string(n_q_points_1d);
7675 message += "," + Utilities::int_to_string(n_components);
7676 message += ",Number>(data";
7677 if (dof_no != numbers::invalid_unsigned_int)
7678 {
7679 message +=
7680 ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7681 message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7682 message += Utilities::int_to_string(proposed_fe_comp);
7683 }
7684 message += ")?\n";
7685 std::string correct_pos;
7686 if (proposed_dof_comp != dof_no)
7687 correct_pos = " ^ ";
7688 else
7689 correct_pos = " ";
7690 if (proposed_quad_comp != this->quad_no)
7691 correct_pos += " ^ ";
7692 else
7693 correct_pos += " ";
7694 if (proposed_fe_comp != first_selected_component)
7695 correct_pos += " ^\n";
7696 else
7697 correct_pos += " \n";
7698 message += " " +
7699 correct_pos;
7700 }
7701 // ok, did not find the numbers specified by the template arguments in
7702 // the given list. Suggest correct template arguments
7703 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7704 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7705 message += "Wrong template arguments:\n";
7706 message += " Did you mean FEEvaluation<dim,";
7707 message +=
7708 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7709 message += Utilities::int_to_string(proposed_n_q_points_1d);
7710 message += "," + Utilities::int_to_string(n_components);
7711 message += ",Number>(data";
7712 if (dof_no != numbers::invalid_unsigned_int)
7713 {
7714 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7715 message += Utilities::int_to_string(this->quad_no);
7716 message += ", " + Utilities::int_to_string(first_selected_component);
7717 }
7718 message += ")?\n";
7719 std::string correct_pos;
7720 if (this->data->data.front().fe_degree !=
7721 static_cast<unsigned int>(fe_degree))
7722 correct_pos = " ^";
7723 else
7724 correct_pos = " ";
7725 if (proposed_n_q_points_1d != n_q_points_1d)
7726 correct_pos += " ^\n";
7727 else
7728 correct_pos += " \n";
7729 message += " " + correct_pos;
7730
7731 Assert(static_cast<unsigned int>(fe_degree) ==
7732 this->data->data.front().fe_degree &&