Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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 
48 namespace 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 
64 template <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>>
71 
72 
73 
100 template <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 
107 public:
108  static constexpr unsigned int dimension = dim;
109 
114 
123  unsigned int
124  get_mapping_data_index_offset() const;
125 
133  get_cell_type() const;
134 
139  get_shape_info() const;
140 
145  get_dof_info() const;
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
190  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
191 
198  void
199  set_cell_data(AlignedVector<VectorizedArrayType> &array,
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
217  set_cell_data(
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()>
226  get_cell_ids() const;
227 
232  std::array<unsigned int, VectorizedArrayType::size()>
233  get_cell_or_face_ids() const;
234 
235 
241  const std::vector<unsigned int> &
242  get_internal_dof_numbering() const;
243 
251  get_scratch_data() const;
252 
256  unsigned int
257  get_quadrature_index() const;
258 
262  unsigned int
263  get_current_cell_index() const;
264 
269  unsigned int
270  get_active_fe_index() const;
271 
276  unsigned int
277  get_active_quadrature_index() const;
278 
283  get_matrix_free() const;
284 
285 protected:
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 
332  operator=(const FEEvaluationBaseData &other);
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 
539 template <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 {
547 public:
548  using number_type = Number;
550  using gradient_type =
552  using hessian_type =
554  static constexpr unsigned int dimension = dim;
555  static constexpr unsigned int n_components = n_components_;
556 
593  template <typename VectorType>
594  void
595  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
596 
625  template <typename VectorType>
626  void
627  read_dof_values_plain(const VectorType & src,
628  const unsigned int first_index = 0);
629 
661  template <typename VectorType>
662  void
663  distribute_local_to_global(
664  VectorType & dst,
665  const unsigned int first_index = 0,
666  const std::bitset<VectorizedArrayType::size()> &mask =
667  std::bitset<VectorizedArrayType::size()>().flip()) const;
668 
707  template <typename VectorType>
708  void
709  set_dof_values(VectorType & dst,
710  const unsigned int first_index = 0,
711  const std::bitset<VectorizedArrayType::size()> &mask =
712  std::bitset<VectorizedArrayType::size()>().flip()) const;
713 
717  template <typename VectorType>
718  void
719  set_dof_values_plain(
720  VectorType & dst,
721  const unsigned int first_index = 0,
722  const std::bitset<VectorizedArrayType::size()> &mask =
723  std::bitset<VectorizedArrayType::size()>().flip()) const;
724 
726 
747  value_type
748  get_dof_value(const unsigned int dof) const;
749 
760  void
761  submit_dof_value(const value_type val_in, const unsigned int dof);
762 
775  value_type
776  get_value(const unsigned int q_point) const;
777 
790  void
791  submit_value(const value_type val_in, const unsigned int q_point);
792 
804  get_gradient(const unsigned int q_point) const;
805 
820  value_type
821  get_normal_derivative(const unsigned int q_point) const;
822 
835  void
836  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
837 
856  void
857  submit_normal_derivative(const value_type grad_in,
858  const unsigned int q_point);
859 
872  void
873  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
874 
887  get_hessian(const unsigned int q_point) const;
888 
899  get_hessian_diagonal(const unsigned int q_point) const;
900 
912  value_type
913  get_laplacian(const unsigned int q_point) const;
914 
915 #ifdef DOXYGEN
916  // doxygen does not anyhow mention functions coming from partial template
917  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
918  // For now, hack in those functions manually only to fix documentation:
919 
926  VectorizedArrayType
927  get_divergence(const unsigned int q_point) const;
928 
938  get_symmetric_gradient(const unsigned int q_point) const;
939 
947  get_curl(const unsigned int q_point) const;
948 
964  void
965  submit_divergence(const VectorizedArrayType div_in,
966  const unsigned int q_point);
967 
984  void
985  submit_symmetric_gradient(
987  const unsigned int q_point);
988 
1001  void
1002  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1003  const unsigned int q_point);
1004 
1005 #endif
1006 
1023  value_type
1024  integrate_value() const;
1025 
1027 
1040  const VectorizedArrayType *
1041  begin_dof_values() const;
1042 
1051  VectorizedArrayType *
1052  begin_dof_values();
1053 
1064  const VectorizedArrayType *
1065  begin_values() const;
1066 
1077  VectorizedArrayType *
1078  begin_values();
1079 
1091  const VectorizedArrayType *
1092  begin_gradients() const;
1093 
1105  VectorizedArrayType *
1106  begin_gradients();
1107 
1120  const VectorizedArrayType *
1121  begin_hessians() const;
1122 
1135  VectorizedArrayType *
1136  begin_hessians();
1137 
1139 
1143  unsigned int
1144  get_first_selected_component() const;
1145 
1146 protected:
1157  const unsigned int dof_no,
1158  const unsigned int first_selected_component,
1159  const unsigned int quad_no,
1160  const unsigned int fe_degree,
1161  const unsigned int n_q_points,
1162  const bool is_interior_face,
1163  const unsigned int active_fe_index,
1164  const unsigned int active_quad_index,
1165  const unsigned int face_type);
1166 
1204  const Mapping<dim> & mapping,
1205  const FiniteElement<dim> &fe,
1206  const Quadrature<1> & quadrature,
1207  const UpdateFlags update_flags,
1208  const unsigned int first_selected_component,
1210  *other);
1211 
1218  FEEvaluationBase(const FEEvaluationBase &other);
1219 
1227  operator=(const FEEvaluationBase &other);
1228 
1235  template <typename VectorType, typename VectorOperation>
1236  void
1237  read_write_operation(
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,
1244  const bool apply_constraints = true) const;
1245 
1253  template <typename VectorType, typename VectorOperation>
1254  void
1255  read_write_operation_contiguous(
1256  const VectorOperation & operation,
1257  const std::array<VectorType *, n_components_> &vectors,
1258  const std::array<
1260  n_components_> & vectors_sm,
1261  const std::bitset<VectorizedArrayType::size()> &mask) const;
1262 
1270  template <typename VectorType, typename VectorOperation>
1271  void
1272  read_write_operation_global(
1273  const VectorOperation & operation,
1274  const std::array<VectorType *, n_components_> &vectors) const;
1275 
1279  template <bool transpose>
1280  void
1281  apply_hanging_node_constraints() const;
1282 
1295  VectorizedArrayType *values_dofs[n_components];
1296 
1308  VectorizedArrayType *values_quad;
1309 
1323  VectorizedArrayType *gradients_quad;
1324 
1338  VectorizedArrayType *gradients_from_hessians_quad;
1339 
1351  VectorizedArrayType *hessians_quad;
1352 
1357  const unsigned int n_fe_components;
1358 
1365 
1372 
1379 
1386 
1393 
1400 
1407 
1412  const unsigned int first_selected_component;
1413 
1418  mutable std::vector<types::global_dof_index> local_dof_indices;
1419 
1420 private:
1425  void
1426  set_data_pointers();
1427 };
1428 
1429 
1430 
1438 template <int dim,
1439  int n_components_,
1440  typename Number,
1441  bool is_face,
1442  typename VectorizedArrayType = VectorizedArray<Number>>
1444  n_components_,
1445  Number,
1446  is_face,
1447  VectorizedArrayType>
1448 {
1449  static_assert(
1450  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1451  "Type of Number and of VectorizedArrayType do not match.");
1452 
1453 public:
1454  using number_type = Number;
1456  using gradient_type =
1458  static constexpr unsigned int dimension = dim;
1459  static constexpr unsigned int n_components = n_components_;
1460  using BaseClass =
1462 
1463 protected:
1473  const unsigned int dof_no,
1474  const unsigned int first_selected_component,
1475  const unsigned int quad_no,
1476  const unsigned int fe_degree,
1477  const unsigned int n_q_points,
1478  const bool is_interior_face = true,
1479  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1480  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1481  const unsigned int face_type = numbers::invalid_unsigned_int);
1482 
1488  const Mapping<dim> & mapping,
1489  const FiniteElement<dim> &fe,
1490  const Quadrature<1> & quadrature,
1491  const UpdateFlags update_flags,
1492  const unsigned int first_selected_component,
1494  *other);
1495 
1499  FEEvaluationAccess(const FEEvaluationAccess &other);
1500 
1505  operator=(const FEEvaluationAccess &other);
1506 };
1507 
1508 
1509 
1518 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1519 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1520  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1521 {
1522  static_assert(
1523  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1524  "Type of Number and of VectorizedArrayType do not match.");
1525 
1526 public:
1527  using number_type = Number;
1528  using value_type = VectorizedArrayType;
1531  static constexpr unsigned int dimension = dim;
1532  using BaseClass =
1534 
1538  value_type
1539  get_dof_value(const unsigned int dof) const;
1540 
1544  void
1545  submit_dof_value(const value_type val_in, const unsigned int dof);
1546 
1550  value_type
1551  get_value(const unsigned int q_point) const;
1552 
1556  void
1557  submit_value(const value_type val_in, const unsigned int q_point);
1558 
1562  void
1563  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
1564  const unsigned int q_point);
1565 
1570  get_gradient(const unsigned int q_point) const;
1571 
1575  value_type
1576  get_normal_derivative(const unsigned int q_point) const;
1577 
1581  void
1582  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1583 
1587  void
1588  submit_normal_derivative(const value_type grad_in,
1589  const unsigned int q_point);
1590 
1594  hessian_type
1595  get_hessian(unsigned int q_point) const;
1596 
1601  get_hessian_diagonal(const unsigned int q_point) const;
1602 
1606  void
1607  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1608 
1612  value_type
1613  get_laplacian(const unsigned int q_point) const;
1614 
1618  value_type
1619  integrate_value() const;
1620 
1621 protected:
1631  const unsigned int dof_no,
1632  const unsigned int first_selected_component,
1633  const unsigned int quad_no,
1634  const unsigned int fe_degree,
1635  const unsigned int n_q_points,
1636  const bool is_interior_face = true,
1637  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1638  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1639  const unsigned int face_type = numbers::invalid_unsigned_int);
1640 
1646  const Mapping<dim> & mapping,
1647  const FiniteElement<dim> &fe,
1648  const Quadrature<1> & quadrature,
1649  const UpdateFlags update_flags,
1650  const unsigned int first_selected_component,
1652  *other);
1653 
1657  FEEvaluationAccess(const FEEvaluationAccess &other);
1658 
1663  operator=(const FEEvaluationAccess &other);
1664 };
1665 
1666 
1667 
1677 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1678 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1679  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1680 {
1681  static_assert(
1682  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1683  "Type of Number and of VectorizedArrayType do not match.");
1684 
1685 public:
1686  using number_type = Number;
1689  static constexpr unsigned int dimension = dim;
1690  static constexpr unsigned int n_components = dim;
1691  using BaseClass =
1693 
1698  get_gradient(const unsigned int q_point) const;
1699 
1704  VectorizedArrayType
1705  get_divergence(const unsigned int q_point) const;
1706 
1714  get_symmetric_gradient(const unsigned int q_point) const;
1715 
1721  get_curl(const unsigned int q_point) const;
1722 
1727  get_hessian(const unsigned int q_point) const;
1728 
1733  get_hessian_diagonal(const unsigned int q_point) const;
1734 
1738  void
1739  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1740 
1749  void
1750  submit_gradient(
1751  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1752  const unsigned int q_point);
1753 
1762  void
1763  submit_divergence(const VectorizedArrayType div_in,
1764  const unsigned int q_point);
1765 
1774  void
1775  submit_symmetric_gradient(
1777  const unsigned int q_point);
1778 
1783  void
1784  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1785  const unsigned int q_point);
1786 
1787 protected:
1797  const unsigned int dof_no,
1798  const unsigned int first_selected_component,
1799  const unsigned int quad_no,
1800  const unsigned int dofs_per_cell,
1801  const unsigned int n_q_points,
1802  const bool is_interior_face = true,
1803  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1804  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1805  const unsigned int face_type = numbers::invalid_unsigned_int);
1806 
1812  const Mapping<dim> & mapping,
1813  const FiniteElement<dim> &fe,
1814  const Quadrature<1> & quadrature,
1815  const UpdateFlags update_flags,
1816  const unsigned int first_selected_component,
1818  *other);
1819 
1823  FEEvaluationAccess(const FEEvaluationAccess &other);
1824 
1829  operator=(const FEEvaluationAccess &other);
1830 };
1831 
1832 
1841 template <typename Number, bool is_face, typename VectorizedArrayType>
1842 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1843  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1844 {
1845  static_assert(
1846  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1847  "Type of Number and of VectorizedArrayType do not match.");
1848 
1849 public:
1850  using number_type = Number;
1851  using value_type = VectorizedArrayType;
1854  static constexpr unsigned int dimension = 1;
1855  using BaseClass =
1857 
1861  value_type
1862  get_dof_value(const unsigned int dof) const;
1863 
1867  void
1868  submit_dof_value(const value_type val_in, const unsigned int dof);
1869 
1873  value_type
1874  get_value(const unsigned int q_point) const;
1875 
1879  void
1880  submit_value(const value_type val_in, const unsigned int q_point);
1881 
1885  void
1886  submit_value(const gradient_type val_in, const unsigned int q_point);
1887 
1892  get_gradient(const unsigned int q_point) const;
1893 
1897  value_type
1898  get_divergence(const unsigned int q_point) const;
1899 
1903  value_type
1904  get_normal_derivative(const unsigned int q_point) const;
1905 
1909  void
1910  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1911 
1915  void
1916  submit_gradient(const value_type grad_in, const unsigned int q_point);
1917 
1921  void
1922  submit_gradient(const Tensor<2, 1, VectorizedArrayType> grad_in,
1923  const unsigned int q_point);
1924 
1928  void
1929  submit_normal_derivative(const value_type grad_in,
1930  const unsigned int q_point);
1931 
1935  void
1936  submit_normal_derivative(const gradient_type grad_in,
1937  const unsigned int q_point);
1938 
1942  hessian_type
1943  get_hessian(unsigned int q_point) const;
1944 
1949  get_hessian_diagonal(const unsigned int q_point) const;
1950 
1954  void
1955  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1956 
1960  value_type
1961  get_laplacian(const unsigned int q_point) const;
1962 
1966  value_type
1967  integrate_value() const;
1968 
1969 protected:
1978  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
1979  const unsigned int dof_no,
1980  const unsigned int first_selected_component,
1981  const unsigned int quad_no,
1982  const unsigned int fe_degree,
1983  const unsigned int n_q_points,
1984  const bool is_interior_face = true,
1985  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1986  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1987  const unsigned int face_type = numbers::invalid_unsigned_int);
1988 
1994  const Mapping<1> & mapping,
1995  const FiniteElement<1> &fe,
1996  const Quadrature<1> & quadrature,
1997  const UpdateFlags update_flags,
1998  const unsigned int first_selected_component,
2000 
2004  FEEvaluationAccess(const FEEvaluationAccess &other);
2005 
2010  operator=(const FEEvaluationAccess &other);
2011 };
2012 
2013 
2014 
2569 template <int dim,
2570  int fe_degree,
2571  int n_q_points_1d,
2572  int n_components_,
2573  typename Number,
2574  typename VectorizedArrayType>
2575 class FEEvaluation : public FEEvaluationAccess<dim,
2576  n_components_,
2577  Number,
2578  false,
2579  VectorizedArrayType>
2580 {
2581  static_assert(
2582  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2583  "Type of Number and of VectorizedArrayType do not match.");
2584 
2585 public:
2589  using BaseClass =
2591 
2595  using number_type = Number;
2596 
2603 
2610 
2614  static constexpr unsigned int dimension = dim;
2615 
2620  static constexpr unsigned int n_components = n_components_;
2621 
2628  static constexpr unsigned int static_n_q_points =
2629  Utilities::pow(n_q_points_1d, dim);
2630 
2638  static constexpr unsigned int static_dofs_per_component =
2639  Utilities::pow(fe_degree + 1, dim);
2640 
2648  static constexpr unsigned int tensor_dofs_per_cell =
2649  static_dofs_per_component * n_components;
2650 
2658  static constexpr unsigned int static_dofs_per_cell =
2659  static_dofs_per_component * n_components;
2660 
2695  FEEvaluation(
2697  const unsigned int dof_no = 0,
2698  const unsigned int quad_no = 0,
2699  const unsigned int first_selected_component = 0,
2700  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2701  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2702 
2711  const std::pair<unsigned int, unsigned int> & range,
2712  const unsigned int dof_no = 0,
2713  const unsigned int quad_no = 0,
2714  const unsigned int first_selected_component = 0);
2715 
2744  FEEvaluation(const Mapping<dim> & mapping,
2745  const FiniteElement<dim> &fe,
2746  const Quadrature<1> & quadrature,
2747  const UpdateFlags update_flags,
2748  const unsigned int first_selected_component = 0);
2749 
2755  FEEvaluation(const FiniteElement<dim> &fe,
2756  const Quadrature<1> & quadrature,
2757  const UpdateFlags update_flags,
2758  const unsigned int first_selected_component = 0);
2759 
2770  FEEvaluation(
2771  const FiniteElement<dim> & fe,
2773  const unsigned int first_selected_component = 0);
2774 
2781  FEEvaluation(const FEEvaluation &other);
2782 
2789  FEEvaluation &
2790  operator=(const FEEvaluation &other);
2791 
2800  void
2801  reinit(const unsigned int cell_batch_index);
2802 
2815  template <bool level_dof_access>
2816  void
2818 
2829  void
2830  reinit(const typename Triangulation<dim>::cell_iterator &cell);
2831 
2835  static bool
2836  fast_evaluation_supported(const unsigned int given_degree,
2837  const unsigned int give_n_q_points_1d);
2838 
2848  void
2849  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2850 
2855  DEAL_II_DEPRECATED_EARLY void
2856  evaluate(const bool evaluate_values,
2857  const bool evaluate_gradients,
2858  const bool evaluate_hessians = false);
2859 
2872  void
2873  evaluate(const VectorizedArrayType * values_array,
2874  const EvaluationFlags::EvaluationFlags evaluation_flag);
2875 
2880  DEAL_II_DEPRECATED_EARLY void
2881  evaluate(const VectorizedArrayType *values_array,
2882  const bool evaluate_values,
2883  const bool evaluate_gradients,
2884  const bool evaluate_hessians = false);
2885 
2899  template <typename VectorType>
2900  void
2901  gather_evaluate(const VectorType & input_vector,
2902  const EvaluationFlags::EvaluationFlags evaluation_flag);
2903 
2907  template <typename VectorType>
2908  DEAL_II_DEPRECATED_EARLY void
2909  gather_evaluate(const VectorType &input_vector,
2910  const bool evaluate_values,
2911  const bool evaluate_gradients,
2912  const bool evaluate_hessians = false);
2913 
2924  void
2925  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2926 
2927 
2931  DEAL_II_DEPRECATED_EARLY void
2932  integrate(const bool integrate_values, const bool integrate_gradients);
2933 
2945  void
2946  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2947  VectorizedArrayType * values_array);
2948 
2952  DEAL_II_DEPRECATED_EARLY void
2953  integrate(const bool integrate_values,
2954  const bool integrate_gradients,
2955  VectorizedArrayType *values_array);
2956 
2970  template <typename VectorType>
2971  void
2972  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
2973  VectorType & output_vector);
2974 
2978  template <typename VectorType>
2979  DEAL_II_DEPRECATED_EARLY void
2980  integrate_scatter(const bool integrate_values,
2981  const bool integrate_gradients,
2982  VectorType &output_vector);
2983 
2989  quadrature_point(const unsigned int q_point) const;
2990 
2997  const unsigned int dofs_per_component;
2998 
3005  const unsigned int dofs_per_cell;
3006 
3014  const unsigned int n_q_points;
3015 
3016 private:
3021  void
3022  check_template_arguments(const unsigned int fe_no,
3023  const unsigned int first_selected_component);
3024 };
3025 
3026 
3027 
3063 template <int dim,
3064  int fe_degree,
3065  int n_q_points_1d = fe_degree + 1,
3066  int n_components_ = 1,
3067  typename Number = double,
3068  typename VectorizedArrayType = VectorizedArray<Number>>
3070  n_components_,
3071  Number,
3072  true,
3073  VectorizedArrayType>
3074 {
3075  static_assert(
3076  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
3077  "Type of Number and of VectorizedArrayType do not match.");
3078 
3079 public:
3083  using BaseClass =
3085 
3089  using number_type = Number;
3090 
3097 
3104 
3108  static constexpr unsigned int dimension = dim;
3109 
3114  static constexpr unsigned int n_components = n_components_;
3115 
3123  static constexpr unsigned int static_n_q_points =
3124  Utilities::pow(n_q_points_1d, dim - 1);
3125 
3132  static constexpr unsigned int static_n_q_points_cell =
3133  Utilities::pow(n_q_points_1d, dim);
3134 
3141  static constexpr unsigned int static_dofs_per_component =
3142  Utilities::pow(fe_degree + 1, dim);
3143 
3150  static constexpr unsigned int tensor_dofs_per_cell =
3151  static_dofs_per_component * n_components;
3152 
3159  static constexpr unsigned int static_dofs_per_cell =
3160  static_dofs_per_component * n_components;
3161 
3205  const bool is_interior_face = true,
3206  const unsigned int dof_no = 0,
3207  const unsigned int quad_no = 0,
3208  const unsigned int first_selected_component = 0,
3209  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
3210  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
3211  const unsigned int face_type = numbers::invalid_unsigned_int);
3212 
3222  const std::pair<unsigned int, unsigned int> & range,
3223  const bool is_interior_face = true,
3224  const unsigned int dof_no = 0,
3225  const unsigned int quad_no = 0,
3226  const unsigned int first_selected_component = 0);
3227 
3238  void
3239  reinit(const unsigned int face_batch_number);
3240 
3248  void
3249  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
3250 
3254  static bool
3255  fast_evaluation_supported(const unsigned int given_degree,
3256  const unsigned int give_n_q_points_1d);
3257 
3268  void
3269  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3270 
3274  DEAL_II_DEPRECATED_EARLY void
3275  evaluate(const bool evaluate_values, const bool evaluate_gradients);
3276 
3289  void
3290  evaluate(const VectorizedArrayType * values_array,
3291  const EvaluationFlags::EvaluationFlags evaluation_flag);
3292 
3296  DEAL_II_DEPRECATED_EARLY void
3297  evaluate(const VectorizedArrayType *values_array,
3298  const bool evaluate_values,
3299  const bool evaluate_gradients);
3300 
3312  template <typename VectorType>
3313  void
3314  gather_evaluate(const VectorType & input_vector,
3315  const EvaluationFlags::EvaluationFlags evaluation_flag);
3316 
3320  template <typename VectorType>
3321  DEAL_II_DEPRECATED_EARLY void
3322  gather_evaluate(const VectorType &input_vector,
3323  const bool evaluate_values,
3324  const bool evaluate_gradients);
3325 
3335  void
3336  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3337 
3341  DEAL_II_DEPRECATED_EARLY void
3342  integrate(const bool integrate_values, const bool integrate_gradients);
3343 
3352  void
3353  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
3354  VectorizedArrayType * values_array);
3355 
3359  DEAL_II_DEPRECATED_EARLY void
3360  integrate(const bool integrate_values,
3361  const bool integrate_gradients,
3362  VectorizedArrayType *values_array);
3363 
3375  template <typename VectorType>
3376  void
3377  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
3378  VectorType & output_vector);
3379 
3383  template <typename VectorType>
3384  void
3385  integrate_scatter(const bool integrate_values,
3386  const bool integrate_gradients,
3387  VectorType &output_vector);
3388 
3394  quadrature_point(const unsigned int q_point) const;
3395 
3402  const unsigned int dofs_per_component;
3403 
3410  const unsigned int dofs_per_cell;
3411 
3419  const unsigned int n_q_points;
3420 
3421 
3422 private:
3426  std::array<unsigned int, VectorizedArrayType::size()>
3427  compute_face_no_data();
3428 
3432  std::array<unsigned int, VectorizedArrayType::size()>
3433  compute_face_orientations();
3434 };
3435 
3436 
3437 
3438 namespace internal
3439 {
3440  namespace MatrixFreeFunctions
3441  {
3442  // a helper function to compute the number of DoFs of a DGP element at
3443  // compile time, depending on the degree
3444  template <int dim, int degree>
3446  {
3447  // this division is always without remainder
3448  static constexpr unsigned int value =
3449  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
3450  };
3451 
3452  // base specialization: 1d elements have 'degree+1' degrees of freedom
3453  template <int degree>
3454  struct DGP_dofs_per_component<1, degree>
3455  {
3456  static constexpr unsigned int value = degree + 1;
3457  };
3458  } // namespace MatrixFreeFunctions
3459 } // namespace internal
3460 
3461 
3462 /*----------------------- Inline functions ----------------------------------*/
3463 
3464 #ifndef DOXYGEN
3465 
3466 
3467 /*----------------------- FEEvaluationBaseData ------------------------*/
3468 
3469 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3473  const unsigned int dof_no,
3474  const unsigned int first_selected_component,
3475  const unsigned int quad_no_in,
3476  const unsigned int fe_degree,
3477  const unsigned int n_q_points,
3478  const bool is_interior_face,
3479  const unsigned int active_fe_index_in,
3480  const unsigned int active_quad_index_in,
3481  const unsigned int face_type)
3482  : scratch_data_array(data_in.acquire_scratch_data())
3483  , quad_no(quad_no_in)
3484  , matrix_info(&data_in)
3485  , dof_info(&data_in.get_dof_info(dof_no))
3486  , mapping_data(
3487  internal::MatrixFreeFunctions::
3488  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3489  data_in.get_mapping_info(),
3490  quad_no))
3491  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
3492  data_in.get_dof_info(dof_no).fe_index_from_degree(
3493  first_selected_component,
3494  fe_degree) :
3495  (active_fe_index_in != numbers::invalid_unsigned_int ?
3496  active_fe_index_in :
3497  0))
3498  , active_quad_index(
3499  fe_degree != numbers::invalid_unsigned_int ?
3500  (mapping_data->quad_index_from_n_q_points(n_q_points)) :
3501  (active_quad_index_in != numbers::invalid_unsigned_int ?
3502  active_quad_index_in :
3503  std::min<unsigned int>(active_fe_index,
3504  mapping_data->descriptor.size() - 1)))
3505  , descriptor(
3506  &mapping_data->descriptor
3507  [is_face ?
3508  (active_quad_index * std::max<unsigned int>(1, dim - 1) +
3509  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
3510  active_quad_index])
3511  , n_quadrature_points(descriptor->n_q_points)
3512  , data(&data_in.get_shape_info(
3513  dof_no,
3514  quad_no_in,
3515  dof_info->component_to_base_index[first_selected_component],
3516  active_fe_index,
3517  active_quad_index))
3518  , jacobian(nullptr)
3519  , J_value(nullptr)
3520  , normal_vectors(nullptr)
3521  , normal_x_jacobian(nullptr)
3522  , quadrature_weights(descriptor->quadrature_weights.begin())
3523  , cell(numbers::invalid_unsigned_int)
3524  , is_interior_face(is_interior_face)
3525  , dof_access_index(
3526  is_face ?
3527  (is_interior_face ?
3528  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3529  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3530  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3531  , cell_type(internal::MatrixFreeFunctions::general)
3532 {
3533  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3534  AssertDimension(matrix_info->get_task_info().vectorization_length,
3535  VectorizedArrayType::size());
3536  AssertDimension(n_quadrature_points, descriptor->n_q_points);
3537 }
3538 
3539 
3540 
3541 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3544  const Mapping<dim> & mapping,
3545  const FiniteElement<dim> &fe,
3546  const Quadrature<1> & quadrature,
3547  const UpdateFlags update_flags,
3548  const unsigned int first_selected_component,
3550  *other)
3551  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3553  , active_fe_index(numbers::invalid_unsigned_int)
3554  , active_quad_index(numbers::invalid_unsigned_int)
3555  , descriptor(nullptr)
3556  , n_quadrature_points(
3557  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3558  , matrix_info(nullptr)
3559  , dof_info(nullptr)
3560  , mapping_data(nullptr)
3561  ,
3562  // select the correct base element from the given FE component
3564  Quadrature<dim - is_face>(quadrature),
3565  fe,
3566  fe.component_to_base_index(first_selected_component).first))
3567  , jacobian(nullptr)
3568  , J_value(nullptr)
3569  , normal_vectors(nullptr)
3570  , normal_x_jacobian(nullptr)
3571  , quadrature_weights(nullptr)
3572  , cell(0)
3574  , is_interior_face(true)
3576 {
3577  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3578  ExcInternalError());
3579  if (other != nullptr &&
3580  other->mapped_geometry->get_quadrature() == quadrature)
3581  mapped_geometry = other->mapped_geometry;
3582  else
3583  mapped_geometry =
3584  std::make_shared<internal::MatrixFreeFunctions::
3586  mapping, quadrature, update_flags);
3587  cell = 0;
3588 
3589  mapping_data = &mapped_geometry->get_data_storage();
3590  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3591  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3592 }
3593 
3594 
3595 
3596 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3600  &other)
3601  : scratch_data_array(other.matrix_info == nullptr ?
3603  other.matrix_info->acquire_scratch_data())
3604  , quad_no(other.quad_no)
3605  , active_fe_index(other.active_fe_index)
3606  , active_quad_index(other.active_quad_index)
3607  , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
3608  , n_quadrature_points(other.n_quadrature_points)
3609  , matrix_info(other.matrix_info)
3610  , dof_info(other.dof_info)
3611  , mapping_data(other.mapping_data)
3612  , data(other.matrix_info == nullptr ?
3614  *other.data) :
3615  other.data)
3616  , jacobian(nullptr)
3617  , J_value(nullptr)
3618  , normal_vectors(nullptr)
3619  , normal_x_jacobian(nullptr)
3620  , quadrature_weights(other.descriptor == nullptr ?
3621  nullptr :
3622  descriptor->quadrature_weights.begin())
3625  , is_interior_face(other.is_interior_face)
3626  , dof_access_index(other.dof_access_index)
3627 {
3628  // Create deep copy of mapped geometry for use in parallel...
3629  if (other.mapped_geometry.get() != nullptr)
3630  {
3631  mapped_geometry = std::make_shared<
3632  internal::MatrixFreeFunctions::
3633  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3634  other.mapped_geometry->get_fe_values().get_mapping(),
3635  other.mapped_geometry->get_quadrature(),
3636  other.mapped_geometry->get_fe_values().get_update_flags());
3637  mapping_data = &mapped_geometry->get_data_storage();
3638  cell = 0;
3639 
3640  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3641  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3642  }
3643 }
3644 
3645 
3646 
3647 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3651 {
3652  AssertDimension(quad_no, other.quad_no);
3653  AssertDimension(active_fe_index, other.active_fe_index);
3654  AssertDimension(active_quad_index, other.active_quad_index);
3655 
3656  // release old memory
3657  if (matrix_info == nullptr)
3658  {
3659  delete data;
3660  delete scratch_data_array;
3661  }
3662  else
3663  {
3664  matrix_info->release_scratch_data(scratch_data_array);
3665  }
3666 
3667  matrix_info = other.matrix_info;
3668  dof_info = other.dof_info;
3669  descriptor = other.descriptor;
3670  mapping_data = other.mapping_data;
3671  if (other.matrix_info == nullptr)
3672  {
3674  *other.data);
3675  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3676  }
3677  else
3678  {
3679  data = other.data;
3680  scratch_data_array = matrix_info->acquire_scratch_data();
3681  }
3682 
3683  quadrature_weights =
3684  (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr);
3687  is_interior_face = other.is_interior_face;
3688  dof_access_index = other.dof_access_index;
3689 
3690  // Create deep copy of mapped geometry for use in parallel...
3691  if (other.mapped_geometry.get() != nullptr)
3692  {
3693  mapped_geometry = std::make_shared<
3694  internal::MatrixFreeFunctions::
3695  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3696  other.mapped_geometry->get_fe_values().get_mapping(),
3697  other.mapped_geometry->get_quadrature(),
3698  other.mapped_geometry->get_fe_values().get_update_flags());
3699  cell = 0;
3700  mapping_data = &mapped_geometry->get_data_storage();
3701  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3702  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3703  }
3704 
3705  return *this;
3706 }
3707 
3708 
3709 
3710 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3713 {
3714  if (matrix_info != nullptr)
3715  {
3716  try
3717  {
3718  matrix_info->release_scratch_data(scratch_data_array);
3719  }
3720  catch (...)
3721  {}
3722  }
3723  else
3724  {
3725  delete scratch_data_array;
3726  delete data;
3727  data = nullptr;
3728  }
3729  scratch_data_array = nullptr;
3730 }
3731 
3732 
3733 
3734 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3735 inline unsigned int
3738 {
3739  if (matrix_info == nullptr)
3740  return 0;
3741  else
3742  {
3743  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3744  return this->mapping_data->data_index_offsets[cell];
3745  }
3746 }
3747 
3748 
3749 
3750 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3753  const
3754 {
3756  return cell_type;
3757 }
3758 
3759 
3760 
3761 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3764  get_shape_info() const
3765 {
3766  Assert(data != nullptr, ExcInternalError());
3767  return *data;
3768 }
3769 
3770 
3771 
3772 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3775  const
3776 {
3777  Assert(dof_info != nullptr,
3778  ExcMessage(
3779  "FEEvaluation was not initialized with a MatrixFree object!"));
3780  return *dof_info;
3781 }
3782 
3783 
3784 
3785 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3788  get_normal_vector(const unsigned int q_point) const
3789 {
3790  AssertIndexRange(q_point, n_quadrature_points);
3791  Assert(normal_vectors != nullptr,
3793  "update_normal_vectors"));
3794  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3795  return normal_vectors[0];
3796  else
3797  return normal_vectors[q_point];
3798 }
3799 
3800 
3801 
3802 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3803 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3805  const unsigned int q_point) const
3806 {
3807  AssertIndexRange(q_point, n_quadrature_points);
3808  Assert(J_value != nullptr,
3810  "update_values|update_gradients"));
3811  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3812  {
3813  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3814  return J_value[0] * this->quadrature_weights[q_point];
3815  }
3816  else
3817  return J_value[q_point];
3818 }
3819 
3820 
3821 
3822 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3825  inverse_jacobian(const unsigned int q_point) const
3826 {
3827  AssertIndexRange(q_point, n_quadrature_points);
3828  Assert(this->jacobian != nullptr,
3830  "update_gradients"));
3831  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3832  return jacobian[0];
3833  else
3834  return jacobian[q_point];
3835 }
3836 
3837 
3838 
3839 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3840 inline std::array<unsigned int, VectorizedArrayType::size()>
3842  const
3843 {
3844  Assert(this->matrix_info != nullptr, ExcNotInitialized());
3845 
3846  const unsigned int n_lanes = VectorizedArrayType::size();
3847  std::array<unsigned int, n_lanes> cells;
3848 
3849  // initialize array
3850  for (unsigned int i = 0; i < n_lanes; ++i)
3851  cells[i] = numbers::invalid_unsigned_int;
3852 
3853  if ((is_face == false) ||
3854  (is_face &&
3855  this->dof_access_index ==
3857  this->is_interior_face))
3858  {
3859  // cell or interior face face (element-centric loop)
3860  for (unsigned int i = 0; i < n_lanes; ++i)
3861  cells[i] = cell * n_lanes + i;
3862  }
3863  else if (is_face &&
3864  this->dof_access_index ==
3866  this->is_interior_face == false)
3867  {
3868  // exterior face (element-centric loop): for this case, we need to
3869  // look into the FaceInfo field that collects information from both
3870  // sides of a face once for the global mesh, and pick the face id that
3871  // is not the local one (cell_this).
3872  for (unsigned int i = 0; i < n_lanes; ++i)
3873  {
3874  // compute actual (non vectorized) cell ID
3875  const unsigned int cell_this = this->cell * n_lanes + i;
3876  // compute face ID
3877  unsigned int face_index =
3878  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3879  this->face_no,
3880  i);
3881 
3882  if (face_index == numbers::invalid_unsigned_int)
3883  continue; // invalid face ID: no neighbor on boundary
3884 
3885  // get cell ID on both sides of face
3886  auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3887  .cells_interior[face_index % n_lanes];
3888  auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3889  .cells_exterior[face_index % n_lanes];
3890 
3891  // compare the IDs with the given cell ID
3892  if (cell_m == cell_this)
3893  cells[i] = cell_p; // neighbor has the other ID
3894  else if (cell_p == cell_this)
3895  cells[i] = cell_m;
3896  }
3897  }
3898  else if (is_face)
3899  {
3900  // face-centric faces
3901  const unsigned int *cells_ =
3902  is_interior_face ?
3903  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3904  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3905  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3906  if (cells_[i] != numbers::invalid_unsigned_int)
3907  cells[i] = cells_[i];
3908  }
3909 
3910  return cells;
3911 }
3912 
3913 
3914 namespace internal
3915 {
3916  template <int dim,
3917  typename Number,
3918  bool is_face,
3919  typename VectorizedArrayType,
3920  typename VectorizedArrayType2,
3921  typename GlobalVectorType,
3922  typename FU>
3923  inline void
3924  process_cell_data(
3927  GlobalVectorType & array,
3928  VectorizedArrayType2 & out,
3929  const FU & fu)
3930  {
3931  (void)matrix_info;
3932  Assert(matrix_info != nullptr, ExcNotImplemented());
3933  AssertDimension(array.size(),
3934  matrix_info->get_task_info().cell_partition_data.back());
3935 
3936  // 1) collect ids of cell
3937  const auto cells = phi.get_cell_ids();
3938 
3939  // 2) actually gather values
3940  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3941  if (cells[i] != numbers::invalid_unsigned_int)
3942  fu(out[i],
3943  array[cells[i] / VectorizedArrayType::size()]
3944  [cells[i] % VectorizedArrayType::size()]);
3945  }
3946 } // namespace internal
3947 
3948 
3949 
3950 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3951 std::array<unsigned int, VectorizedArrayType::size()>
3953  get_cell_or_face_ids() const
3954 {
3955  const unsigned int v_len = VectorizedArrayType::size();
3956  std::array<unsigned int, VectorizedArrayType::size()> cells;
3957 
3958  // initialize array
3959  for (unsigned int i = 0; i < v_len; ++i)
3960  cells[i] = numbers::invalid_unsigned_int;
3961 
3962  if (is_face &&
3963  this->dof_access_index ==
3965  this->is_interior_face == false)
3966  {
3967  // cell-based face-loop: plus face
3968  for (unsigned int i = 0; i < v_len; ++i)
3969  {
3970  // compute actual (non vectorized) cell ID
3971  const unsigned int cell_this = this->cell * v_len + i;
3972  // compute face ID
3973  unsigned int fn =
3974  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3975  this->face_no,
3976  i);
3977 
3979  continue; // invalid face ID: no neighbor on boundary
3980 
3981  // get cell ID on both sides of face
3982  auto cell_m = this->matrix_info->get_face_info(fn / v_len)
3983  .cells_interior[fn % v_len];
3984  auto cell_p = this->matrix_info->get_face_info(fn / v_len)
3985  .cells_exterior[fn % v_len];
3986 
3987  // compare the IDs with the given cell ID
3988  if (cell_m == cell_this)
3989  cells[i] = cell_p; // neighbor has the other ID
3990  else if (cell_p == cell_this)
3991  cells[i] = cell_m;
3992  }
3993  }
3994  else
3995  {
3996  for (unsigned int i = 0; i < v_len; ++i)
3997  cells[i] = cell * v_len + i;
3998  }
3999 
4000  return cells;
4001 }
4002 
4003 
4004 
4005 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4006 inline VectorizedArrayType
4008  const AlignedVector<VectorizedArrayType> &array) const
4009 {
4010  VectorizedArrayType out = Number(1.);
4011  internal::process_cell_data(
4012  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
4013  local = global;
4014  });
4015  return out;
4016 }
4017 
4018 
4019 
4020 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4021 inline void
4024  const VectorizedArrayType & in) const
4025 {
4026  internal::process_cell_data(
4027  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
4028  global = local;
4029  });
4030 }
4031 
4032 
4033 
4034 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4035 template <typename T>
4036 inline std::array<T, VectorizedArrayType::size()>
4038  const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
4039 {
4040  std::array<T, VectorizedArrayType::size()> out;
4041  internal::process_cell_data(
4042  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
4043  local = global;
4044  });
4045  return out;
4046 }
4047 
4048 
4049 
4050 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4051 template <typename T>
4052 inline void
4054  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
4055  const std::array<T, VectorizedArrayType::size()> & in) const
4056 {
4057  internal::process_cell_data(
4058  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
4059  global = local;
4060  });
4061 }
4062 
4063 
4064 
4065 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4066 inline const std::vector<unsigned int> &
4069 {
4070  return data->lexicographic_numbering;
4071 }
4072 
4073 
4074 
4075 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4078  get_scratch_data() const
4079 {
4081  const_cast<VectorizedArrayType *>(scratch_data),
4082  scratch_data_array->end() - scratch_data);
4083 }
4084 
4085 
4086 
4087 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4088 inline unsigned int
4090  get_quadrature_index() const
4091 {
4092  return this->quad_no;
4093 }
4094 
4095 
4096 
4097 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4098 inline unsigned int
4101 {
4102  if (is_face && this->dof_access_index ==
4104  return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
4105  else
4106  return this->cell;
4107 }
4108 
4109 
4110 
4111 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4112 inline unsigned int
4114  get_active_fe_index() const
4115 {
4116  return active_fe_index;
4117 }
4118 
4119 
4120 
4121 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4122 inline unsigned int
4125 {
4126  return active_quad_index;
4127 }
4128 
4129 
4130 
4131 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4134  get_matrix_free() const
4135 {
4136  Assert(matrix_info != nullptr,
4137  ExcMessage(
4138  "FEEvaluation was not initialized with a MatrixFree object!"));
4139  return *matrix_info;
4140 }
4141 
4142 
4143 /*----------------------- FEEvaluationBase ----------------------------------*/
4144 
4145 template <int dim,
4146  int n_components_,
4147  typename Number,
4148  bool is_face,
4149  typename VectorizedArrayType>
4150 inline FEEvaluationBase<dim,
4151  n_components_,
4152  Number,
4153  is_face,
4154  VectorizedArrayType>::
4155  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
4156  const unsigned int dof_no,
4157  const unsigned int first_selected_component,
4158  const unsigned int quad_no_in,
4159  const unsigned int fe_degree,
4160  const unsigned int n_q_points,
4161  const bool is_interior_face,
4162  const unsigned int active_fe_index,
4163  const unsigned int active_quad_index,
4164  const unsigned int face_type)
4166  data_in,
4167  dof_no,
4168  first_selected_component,
4169  quad_no_in,
4170  fe_degree,
4171  n_q_points,
4172  is_interior_face,
4173  active_fe_index,
4174  active_quad_index,
4175  face_type)
4176  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
4177  , dof_values_initialized(false)
4178  , values_quad_initialized(false)
4179  , gradients_quad_initialized(false)
4180  , hessians_quad_initialized(false)
4181  , values_quad_submitted(false)
4182  , gradients_quad_submitted(false)
4183  , first_selected_component(first_selected_component)
4184 {
4185  set_data_pointers();
4186  Assert(
4187  this->dof_info->start_components.back() == 1 ||
4188  static_cast<int>(n_components_) <=
4189  static_cast<int>(
4190  this->dof_info->start_components
4191  [this->dof_info->component_to_base_index[first_selected_component] +
4192  1]) -
4193  first_selected_component,
4194  ExcMessage(
4195  "You tried to construct a vector-valued evaluator with " +
4196  std::to_string(n_components) +
4197  " components. However, "
4198  "the current base element has only " +
4200  this->dof_info->start_components
4201  [this->dof_info->component_to_base_index[first_selected_component] +
4202  1] -
4203  first_selected_component) +
4204  " components left when starting from local element index " +
4206  first_selected_component -
4207  this->dof_info->start_components
4208  [this->dof_info->component_to_base_index[first_selected_component]]) +
4209  " (global index " + std::to_string(first_selected_component) + ")"));
4210 
4211  // do not check for correct dimensions of data fields here, should be done
4212  // in derived classes
4213 }
4214 
4215 
4216 
4217 template <int dim,
4218  int n_components_,
4219  typename Number,
4220  bool is_face,
4221  typename VectorizedArrayType>
4222 inline FEEvaluationBase<dim,
4223  n_components_,
4224  Number,
4225  is_face,
4226  VectorizedArrayType>::
4227  FEEvaluationBase(
4228  const Mapping<dim> & mapping,
4229  const FiniteElement<dim> &fe,
4230  const Quadrature<1> & quadrature,
4231  const UpdateFlags update_flags,
4232  const unsigned int first_selected_component,
4234  *other)
4236  mapping,
4237  fe,
4238  quadrature,
4239  update_flags,
4240  first_selected_component,
4241  other)
4242  , n_fe_components(n_components_)
4243  , dof_values_initialized(false)
4244  , values_quad_initialized(false)
4245  , gradients_quad_initialized(false)
4246  , hessians_quad_initialized(false)
4247  , values_quad_submitted(false)
4248  , gradients_quad_submitted(false)
4249  // keep the number of the selected component within the current base element
4250  // for reading dof values
4251  , first_selected_component(first_selected_component)
4252 {
4253  set_data_pointers();
4254 
4255  const unsigned int base_element_number =
4256  fe.component_to_base_index(first_selected_component).first;
4257  Assert(fe.element_multiplicity(base_element_number) == 1 ||
4258  fe.element_multiplicity(base_element_number) -
4259  first_selected_component >=
4260  n_components_,
4261  ExcMessage("The underlying element must at least contain as many "
4262  "components as requested by this class"));
4263  (void)base_element_number;
4264 }
4265 
4266 
4267 
4268 template <int dim,
4269  int n_components_,
4270  typename Number,
4271  bool is_face,
4272  typename VectorizedArrayType>
4273 inline FEEvaluationBase<dim,
4274  n_components_,
4275  Number,
4276  is_face,
4277  VectorizedArrayType>::
4278  FEEvaluationBase(const FEEvaluationBase<dim,
4279  n_components_,
4280  Number,
4281  is_face,
4282  VectorizedArrayType> &other)
4284  , n_fe_components(other.n_fe_components)
4285  , dof_values_initialized(false)
4286  , values_quad_initialized(false)
4287  , gradients_quad_initialized(false)
4288  , hessians_quad_initialized(false)
4289  , values_quad_submitted(false)
4290  , gradients_quad_submitted(false)
4291  , first_selected_component(other.first_selected_component)
4292 {
4293  set_data_pointers();
4294 }
4295 
4296 
4297 
4298 template <int dim,
4299  int n_components_,
4300  typename Number,
4301  bool is_face,
4302  typename VectorizedArrayType>
4303 inline FEEvaluationBase<dim,
4304  n_components_,
4305  Number,
4306  is_face,
4307  VectorizedArrayType> &
4309 operator=(const FEEvaluationBase<dim,
4310  n_components_,
4311  Number,
4312  is_face,
4313  VectorizedArrayType> &other)
4314 {
4316  operator=(other);
4317  AssertDimension(n_fe_components, other.n_fe_components);
4318  AssertDimension(first_selected_component, other.first_selected_component);
4319 
4320  return *this;
4321 }
4322 
4323 
4324 
4325 template <int dim,
4326  int n_components_,
4327  typename Number,
4328  bool is_face,
4329  typename VectorizedArrayType>
4330 inline void
4333 {
4334  Assert(this->scratch_data_array != nullptr, ExcInternalError());
4335 
4336  const unsigned int tensor_dofs_per_component =
4337  Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4338  const unsigned int dofs_per_component =
4339  this->data->dofs_per_component_on_cell;
4340  const unsigned int n_quadrature_points = this->n_quadrature_points;
4341 
4342  const unsigned int shift =
4343  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4344  n_components_ * 3 +
4345  2 * n_quadrature_points;
4346  const unsigned int allocated_size =
4347  shift + n_components_ * dofs_per_component +
4348  (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
4349  n_quadrature_points);
4350  this->scratch_data_array->resize_fast(allocated_size);
4351 
4352  // set the pointers to the correct position in the data array
4353  for (unsigned int c = 0; c < n_components_; ++c)
4354  {
4355  values_dofs[c] =
4356  this->scratch_data_array->begin() + c * dofs_per_component;
4357  }
4358  values_quad =
4359  this->scratch_data_array->begin() + n_components * dofs_per_component;
4360  gradients_quad = this->scratch_data_array->begin() +
4361  n_components * (dofs_per_component + n_quadrature_points);
4362  gradients_from_hessians_quad =
4363  this->scratch_data_array->begin() +
4364  n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4365  hessians_quad =
4366  this->scratch_data_array->begin() +
4367  n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
4368  this->scratch_data = this->scratch_data_array->begin() +
4369  n_components_ * dofs_per_component +
4370  (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
4371  n_quadrature_points);
4372 }
4373 
4374 
4375 
4376 namespace internal
4377 {
4378  // allows to select between block vectors and non-block vectors, which
4379  // allows to use a unified interface for extracting blocks on block vectors
4380  // and doing nothing on usual vectors
4381  template <typename VectorType, bool>
4382  struct BlockVectorSelector
4383  {};
4384 
4385  template <typename VectorType>
4386  struct BlockVectorSelector<VectorType, true>
4387  {
4388  using BaseVectorType = typename VectorType::BlockType;
4389 
4390  static BaseVectorType *
4391  get_vector_component(VectorType &vec, const unsigned int component)
4392  {
4393  AssertIndexRange(component, vec.n_blocks());
4394  return &vec.block(component);
4395  }
4396  };
4397 
4398  template <typename VectorType>
4399  struct BlockVectorSelector<VectorType, false>
4400  {
4401  using BaseVectorType = VectorType;
4402 
4403  static BaseVectorType *
4404  get_vector_component(VectorType &vec, const unsigned int component)
4405  {
4406  // FEEvaluation allows to combine several vectors from a scalar
4407  // FiniteElement into a "vector-valued" FEEvaluation object with
4408  // multiple components. These components can be extracted with the other
4409  // get_vector_component functions. If we do not get a vector of vectors
4410  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
4411  // must make sure that we do not duplicate the components in input
4412  // and/or duplicate the resulting integrals. In such a case, we should
4413  // only get the zeroth component in the vector contained set nullptr for
4414  // the others which allows us to catch unintended use in
4415  // read_write_operation.
4416  if (component == 0)
4417  return &vec;
4418  else
4419  return nullptr;
4420  }
4421  };
4422 
4423  template <typename VectorType>
4424  struct BlockVectorSelector<std::vector<VectorType>, false>
4425  {
4426  using BaseVectorType = VectorType;
4427 
4428  static BaseVectorType *
4429  get_vector_component(std::vector<VectorType> &vec,
4430  const unsigned int component)
4431  {
4432  AssertIndexRange(component, vec.size());
4433  return &vec[component];
4434  }
4435  };
4436 
4437  template <typename VectorType>
4438  struct BlockVectorSelector<std::vector<VectorType *>, false>
4439  {
4440  using BaseVectorType = VectorType;
4441 
4442  static BaseVectorType *
4443  get_vector_component(std::vector<VectorType *> &vec,
4444  const unsigned int component)
4445  {
4446  AssertIndexRange(component, vec.size());
4447  return vec[component];
4448  }
4449  };
4450 } // namespace internal
4451 
4452 
4453 
4454 template <int dim,
4455  int n_components_,
4456  typename Number,
4457  bool is_face,
4458  typename VectorizedArrayType>
4459 template <typename VectorType, typename VectorOperation>
4460 inline void
4463  const VectorOperation & operation,
4464  const std::array<VectorType *, n_components_> &src,
4465  const std::array<
4467  n_components_> & src_sm,
4468  const std::bitset<VectorizedArrayType::size()> &mask,
4469  const bool apply_constraints) const
4470 {
4471  // Case 1: No MatrixFree object given, simple case because we do not need to
4472  // process constraints and need not care about vectorization -> go to
4473  // separate function
4474  if (this->matrix_info == nullptr)
4475  {
4476  read_write_operation_global(operation, src);
4477  return;
4478  }
4479 
4480  Assert(this->dof_info != nullptr, ExcNotInitialized());
4481  Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
4482  if (n_fe_components == 1)
4483  for (unsigned int comp = 0; comp < n_components; ++comp)
4484  {
4485  Assert(src[comp] != nullptr,
4486  ExcMessage("The finite element underlying this FEEvaluation "
4487  "object is scalar, but you requested " +
4488  std::to_string(n_components) +
4489  " components via the template argument in "
4490  "FEEvaluation. In that case, you must pass an "
4491  "std::vector<VectorType> or a BlockVector to " +
4492  "read_dof_values and distribute_local_to_global."));
4493  internal::check_vector_compatibility(*src[comp], *this->dof_info);
4494  }
4495  else
4496  {
4497  internal::check_vector_compatibility(*src[0], *this->dof_info);
4498  }
4499 
4500  // Case 2: contiguous indices which use reduced storage of indices and can
4501  // use vectorized load/store operations -> go to separate function
4503  this->cell,
4504  this->dof_info->index_storage_variants[this->dof_access_index].size());
4505  if (this->dof_info->index_storage_variants
4506  [is_face ? this->dof_access_index :
4508  [this->cell] >=
4510  {
4511  read_write_operation_contiguous(operation, src, src_sm, mask);
4512  return;
4513  }
4514 
4515  // Case 3: standard operation with one index per degree of freedom -> go on
4516  // here
4517  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4518  Assert(mask.count() == n_lanes,
4519  ExcNotImplemented("Masking currently not implemented for "
4520  "non-contiguous DoF storage"));
4521 
4522  unsigned int n_vectorization_actual =
4523  this->dof_info
4524  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4525 
4526  bool has_hn_constraints = false;
4527 
4528  if (is_face == false)
4529  {
4530  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4531  if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
4532  this->dof_info
4533  ->hanging_node_constraint_masks[(this->cell * n_lanes + v) *
4534  n_fe_components +
4535  first_selected_component] !=
4537  has_hn_constraints = true;
4538  }
4539 
4540  std::integral_constant<bool,
4542  vector_selector;
4543 
4544  const unsigned int dofs_per_component =
4545  this->data->dofs_per_component_on_cell;
4546  if (this->dof_info->index_storage_variants
4547  [is_face ? this->dof_access_index :
4549  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
4550  IndexStorageVariants::interleaved &&
4551  (has_hn_constraints == false))
4552  {
4553  const unsigned int *dof_indices =
4554  this->dof_info->dof_indices_interleaved.data() +
4555  this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4556  .first +
4557  this->dof_info
4558  ->component_dof_indices_offset[this->active_fe_index]
4559  [this->first_selected_component] *
4560  n_lanes;
4561  if (n_components == 1 || n_fe_components == 1)
4562  for (unsigned int i = 0; i < dofs_per_component;
4563  ++i, dof_indices += n_lanes)
4564  for (unsigned int comp = 0; comp < n_components; ++comp)
4565  operation.process_dof_gather(dof_indices,
4566  *src[comp],
4567  0,
4568  values_dofs[comp][i],
4569  vector_selector);
4570  else
4571  for (unsigned int comp = 0; comp < n_components; ++comp)
4572  for (unsigned int i = 0; i < dofs_per_component;
4573  ++i, dof_indices += n_lanes)
4574  operation.process_dof_gather(
4575  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4576  return;
4577  }
4578 
4579  const unsigned int * dof_indices[n_lanes];
4580  VectorizedArrayType **values_dofs =
4581  const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
4582 
4583  // Assign the appropriate cell ids for face/cell case and get the pointers
4584  // to the dof indices of the cells on all lanes
4585  unsigned int cells_copied[n_lanes];
4586  const unsigned int *cells;
4587  bool has_constraints = false;
4588  const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4589 
4590  if (is_face)
4591  {
4592  if (this->dof_access_index ==
4594  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4595  cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4596  cells =
4597  this->dof_access_index ==
4599  &cells_copied[0] :
4600  (this->is_interior_face ?
4601  &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
4602  &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
4603  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4604  {
4605  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4606  ExcInternalError());
4607  const std::pair<unsigned int, unsigned int> *my_index_start =
4608  &this->dof_info->row_starts[cells[v] * n_fe_components +
4609  first_selected_component];
4610 
4611  // check whether any of the SIMD lanes has constraints, i.e., the
4612  // constraint indicator which is the second entry of row_starts
4613  // increments on this cell
4614  if (my_index_start[n_components_read].second !=
4615  my_index_start[0].second)
4616  has_constraints = true;
4617 
4618  dof_indices[v] =
4619  this->dof_info->dof_indices.data() + my_index_start[0].first;
4620  }
4621  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4622  dof_indices[v] = nullptr;
4623  }
4624  else
4625  {
4626  AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
4627  this->dof_info->row_starts.size());
4628  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4629  {
4630  const std::pair<unsigned int, unsigned int> *my_index_start =
4631  &this->dof_info
4632  ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4633  first_selected_component];
4634  if (my_index_start[n_components_read].second !=
4635  my_index_start[0].second)
4636  has_constraints = true;
4637 
4638  if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
4639  this->dof_info
4640  ->hanging_node_constraint_masks[(this->cell * n_lanes + v) *
4641  n_fe_components +
4642  first_selected_component] !=
4644  has_hn_constraints = true;
4645 
4646  Assert(my_index_start[n_components_read].first ==
4647  my_index_start[0].first ||
4648  my_index_start[0].first < this->dof_info->dof_indices.size(),
4649  ExcIndexRange(0,
4650  my_index_start[0].first,
4651  this->dof_info->dof_indices.size()));
4652  dof_indices[v] =
4653  this->dof_info->dof_indices.data() + my_index_start[0].first;
4654  }
4655  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4656  dof_indices[v] = nullptr;
4657  }
4658 
4659  // Case where we have no constraints throughout the whole cell: Can go
4660  // through the list of DoFs directly
4661  if (!has_constraints && apply_constraints)
4662  {
4663  if (n_vectorization_actual < n_lanes)
4664  for (unsigned int comp = 0; comp < n_components; ++comp)
4665  for (unsigned int i = 0; i < dofs_per_component; ++i)
4666  operation.process_empty(values_dofs[comp][i]);
4667  if (n_components == 1 || n_fe_components == 1)
4668  {
4669  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4670  for (unsigned int i = 0; i < dofs_per_component; ++i)
4671  for (unsigned int comp = 0; comp < n_components; ++comp)
4672  operation.process_dof(dof_indices[v][i],
4673  *src[comp],
4674  values_dofs[comp][i][v]);
4675  }
4676  else
4677  {
4678  for (unsigned int comp = 0; comp < n_components; ++comp)
4679  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4680  for (unsigned int i = 0; i < dofs_per_component; ++i)
4681  operation.process_dof(
4682  dof_indices[v][comp * dofs_per_component + i],
4683  *src[0],
4684  values_dofs[comp][i][v]);
4685  }
4686  return;
4687  }
4688 
4689  // In the case where there are some constraints to be resolved, loop over
4690  // all vector components that are filled and then over local dofs. ind_local
4691  // holds local number on cell, index iterates over the elements of
4692  // index_local_to_global and dof_indices points to the global indices stored
4693  // in index_local_to_global
4694  if (n_vectorization_actual < n_lanes)
4695  for (unsigned int comp = 0; comp < n_components; ++comp)
4696  for (unsigned int i = 0; i < dofs_per_component; ++i)
4697  operation.process_empty(values_dofs[comp][i]);
4698  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4699  {
4700  const unsigned int cell_index =
4701  is_face ? cells[v] : this->cell * n_lanes + v;
4702  const unsigned int cell_dof_index =
4703  cell_index * n_fe_components + first_selected_component;
4704  const unsigned int n_components_read =
4705  n_fe_components > 1 ? n_components : 1;
4706  unsigned int index_indicators =
4707  this->dof_info->row_starts[cell_dof_index].second;
4708  unsigned int next_index_indicators =
4709  this->dof_info->row_starts[cell_dof_index + 1].second;
4710 
4711  // For read_dof_values_plain, redirect the dof_indices field to the
4712  // unconstrained indices
4713  if (apply_constraints == false &&
4714  (this->dof_info->row_starts[cell_dof_index].second !=
4715  this->dof_info->row_starts[cell_dof_index + n_components_read]
4716  .second ||
4717  (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
4718  this->dof_info->hanging_node_constraint_masks[cell_dof_index] !=
4720  {
4721  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
4723  ExcNotInitialized());
4724  dof_indices[v] =
4725  this->dof_info->plain_dof_indices.data() +
4726  this->dof_info
4727  ->component_dof_indices_offset[this->active_fe_index]
4728  [this->first_selected_component] +
4729  this->dof_info->row_starts_plain_indices[cell_index];
4730  next_index_indicators = index_indicators;
4731  }
4732 
4733  if (n_components == 1 || n_fe_components == 1)
4734  {
4735  unsigned int ind_local = 0;
4736  for (; index_indicators != next_index_indicators; ++index_indicators)
4737  {
4738  const std::pair<unsigned short, unsigned short> indicator =
4739  this->dof_info->constraint_indicator[index_indicators];
4740  // run through values up to next constraint
4741  for (unsigned int j = 0; j < indicator.first; ++j)
4742  for (unsigned int comp = 0; comp < n_components; ++comp)
4743  operation.process_dof(dof_indices[v][j],
4744  *src[comp],
4745  values_dofs[comp][ind_local + j][v]);
4746 
4747  ind_local += indicator.first;
4748  dof_indices[v] += indicator.first;
4749 
4750  // constrained case: build the local value as a linear
4751  // combination of the global value according to constraints
4752  Number value[n_components];
4753  for (unsigned int comp = 0; comp < n_components; ++comp)
4754  operation.pre_constraints(values_dofs[comp][ind_local][v],
4755  value[comp]);
4756 
4757  const Number *data_val =
4758  this->matrix_info->constraint_pool_begin(indicator.second);
4759  const Number *end_pool =
4760  this->matrix_info->constraint_pool_end(indicator.second);
4761  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4762  for (unsigned int comp = 0; comp < n_components; ++comp)
4763  operation.process_constraint(*dof_indices[v],
4764  *data_val,
4765  *src[comp],
4766  value[comp]);
4767 
4768  for (unsigned int comp = 0; comp < n_components; ++comp)
4769  operation.post_constraints(value[comp],
4770  values_dofs[comp][ind_local][v]);
4771  ind_local++;
4772  }
4773 
4774  AssertIndexRange(ind_local, dofs_per_component + 1);
4775 
4776  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4777  for (unsigned int comp = 0; comp < n_components; ++comp)
4778  operation.process_dof(*dof_indices[v],
4779  *src[comp],
4780  values_dofs[comp][ind_local][v]);
4781  }
4782  else
4783  {
4784  // case with vector-valued finite elements where all components are
4785  // included in one single vector. Assumption: first come all entries
4786  // to the first component, then all entries to the second one, and
4787  // so on. This is ensured by the way MatrixFree reads out the
4788  // indices.
4789  for (unsigned int comp = 0; comp < n_components; ++comp)
4790  {
4791  unsigned int ind_local = 0;
4792 
4793  // check whether there is any constraint on the current cell
4794  for (; index_indicators != next_index_indicators;
4795  ++index_indicators)
4796  {
4797  const std::pair<unsigned short, unsigned short> indicator =
4798  this->dof_info->constraint_indicator[index_indicators];
4799 
4800  // run through values up to next constraint
4801  for (unsigned int j = 0; j < indicator.first; ++j)
4802  operation.process_dof(dof_indices[v][j],
4803  *src[0],
4804  values_dofs[comp][ind_local + j][v]);
4805  ind_local += indicator.first;
4806  dof_indices[v] += indicator.first;
4807 
4808  // constrained case: build the local value as a linear
4809  // combination of the global value according to constraints
4810  Number value;
4811  operation.pre_constraints(values_dofs[comp][ind_local][v],
4812  value);
4813 
4814  const Number *data_val =
4815  this->matrix_info->constraint_pool_begin(indicator.second);
4816  const Number *end_pool =
4817  this->matrix_info->constraint_pool_end(indicator.second);
4818 
4819  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4820  operation.process_constraint(*dof_indices[v],
4821  *data_val,
4822  *src[0],
4823  value);
4824 
4825  operation.post_constraints(value,
4826  values_dofs[comp][ind_local][v]);
4827  ind_local++;
4828  }
4829 
4830  AssertIndexRange(ind_local, dofs_per_component + 1);
4831 
4832  // get the dof values past the last constraint
4833  for (; ind_local < dofs_per_component;
4834  ++dof_indices[v], ++ind_local)
4835  {
4836  AssertIndexRange(*dof_indices[v], src[0]->size());
4837  operation.process_dof(*dof_indices[v],
4838  *src[0],
4839  values_dofs[comp][ind_local][v]);
4840  }
4841 
4842  if (apply_constraints == true && comp + 1 < n_components)
4843  next_index_indicators =
4844  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4845  }
4846  }
4847  }
4848 }
4849 
4850 
4851 
4852 template <int dim,
4853  int n_components_,
4854  typename Number,
4855  bool is_face,
4856  typename VectorizedArrayType>
4857 template <typename VectorType, typename VectorOperation>
4858 inline void
4861  const VectorOperation & operation,
4862  const std::array<VectorType *, n_components_> &src) const
4863 {
4864  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4865 
4866  unsigned int index =
4867  first_selected_component * this->data->dofs_per_component_on_cell;
4868  for (unsigned int comp = 0; comp < n_components; ++comp)
4869  {
4870  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4871  ++i, ++index)
4872  {
4873  operation.process_empty(values_dofs[comp][i]);
4874  operation.process_dof_global(
4875  local_dof_indices[this->data->lexicographic_numbering[index]],
4876  *src[0],
4877  values_dofs[comp][i][0]);
4878  }
4879  }
4880 }
4881 
4882 
4883 
4884 template <int dim,
4885  int n_components_,
4886  typename Number,
4887  bool is_face,
4888  typename VectorizedArrayType>
4889 template <typename VectorType, typename VectorOperation>
4890 inline void
4893  const VectorOperation & operation,
4894  const std::array<VectorType *, n_components_> &src,
4895  const std::array<
4897  n_components_> & vectors_sm,
4898  const std::bitset<VectorizedArrayType::size()> &mask) const
4899 {
4900  // This functions processes the functions read_dof_values,
4901  // distribute_local_to_global, and set_dof_values with the same code for
4902  // contiguous cell indices (DG case). The distinction between these three
4903  // cases is made by the input VectorOperation that either reads values from
4904  // a vector and puts the data into the local data field or write local data
4905  // into the vector. Certain operations are no-ops for the given use case.
4906 
4907  std::integral_constant<bool,
4909  vector_selector;
4911  is_face ? this->dof_access_index :
4913  const unsigned int n_lanes = mask.count();
4914 
4915  const std::vector<unsigned int> &dof_indices_cont =
4916  this->dof_info->dof_indices_contiguous[ind];
4917 
4918  // Simple case: We have contiguous storage, so we can simply copy out the
4919  // data
4920  if ((this->dof_info->index_storage_variants[ind][this->cell] ==
4922  interleaved_contiguous &&
4923  n_lanes == VectorizedArrayType::size()) &&
4924  !(is_face &&
4925  this->dof_access_index ==
4927  this->is_interior_face == false))
4928  {
4929  const unsigned int dof_index =
4930  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4931  this->dof_info->component_dof_indices_offset[this->active_fe_index]
4932  [first_selected_component] *
4933  VectorizedArrayType::size();
4934  if (n_components == 1 || n_fe_components == 1)
4935  for (unsigned int comp = 0; comp < n_components; ++comp)
4936  operation.process_dofs_vectorized(
4937  this->data->dofs_per_component_on_cell,
4938  dof_index,
4939  *src[comp],
4940  values_dofs[comp],
4941  vector_selector);
4942  else
4943  operation.process_dofs_vectorized(
4944  this->data->dofs_per_component_on_cell * n_components,
4945  dof_index,
4946  *src[0],
4947  values_dofs[0],
4948  vector_selector);
4949  return;
4950  }
4951 
4952  std::array<unsigned int, VectorizedArrayType::size()> cells =
4953  this->get_cell_or_face_ids();
4954 
4955  // More general case: Must go through the components one by one and apply
4956  // some transformations
4957  const unsigned int n_filled_lanes =
4958  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4959 
4960  const bool is_ecl =
4961  this->dof_access_index ==
4963  this->is_interior_face == false;
4964 
4965  if (vectors_sm[0] != nullptr)
4966  {
4967  const auto compute_vector_ptrs = [&](const unsigned int comp) {
4968  std::array<typename VectorType::value_type *,
4969  VectorizedArrayType::size()>
4970  vector_ptrs = {};
4971 
4972  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4973  {
4975  ExcNotImplemented());
4976  Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
4977  ExcIndexRange(
4978  ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
4979  Assert(cells[v] <
4980  this->dof_info->dof_indices_contiguous_sm[ind].size(),
4981  ExcIndexRange(
4982  cells[v],
4983  0,
4984  this->dof_info->dof_indices_contiguous_sm[ind].size()));
4985 
4986  const auto &temp =
4987  this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
4988 
4989  if (temp.first != numbers::invalid_unsigned_int)
4990  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
4991  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
4992  this->dof_info->component_dof_indices_offset
4993  [this->active_fe_index][this->first_selected_component]);
4994  else
4995  vector_ptrs[v] = nullptr;
4996  }
4997  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
4998  ++v)
4999  vector_ptrs[v] = nullptr;
5000 
5001  return vector_ptrs;
5002  };
5003 
5004  if (n_filled_lanes == VectorizedArrayType::size() &&
5005  n_lanes == VectorizedArrayType::size() && !is_ecl)
5006  {
5007  if (n_components == 1 || n_fe_components == 1)
5008  {
5009  for (unsigned int comp = 0; comp < n_components; ++comp)
5010  {
5011  auto vector_ptrs = compute_vector_ptrs(comp);
5012  operation.process_dofs_vectorized_transpose(
5013  this->data->dofs_per_component_on_cell,
5014  vector_ptrs,
5015  values_dofs[comp],
5016  vector_selector);
5017  }
5018  }
5019  else
5020  {
5021  auto vector_ptrs = compute_vector_ptrs(0);
5022  operation.process_dofs_vectorized_transpose(
5023  this->data->dofs_per_component_on_cell * n_components,
5024  vector_ptrs,
5025  &values_dofs[0][0],
5026  vector_selector);
5027  }
5028  }
5029  else
5030  for (unsigned int comp = 0; comp < n_components; ++comp)
5031  {
5032  auto vector_ptrs = compute_vector_ptrs(
5033  (n_components == 1 || n_fe_components == 1) ? comp : 0);
5034 
5035  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5036  ++i)
5037  operation.process_empty(values_dofs[comp][i]);
5038 
5039  if (n_components == 1 || n_fe_components == 1)
5040  {
5041  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5042  if (mask[v] == true)
5043  for (unsigned int i = 0;
5044  i < this->data->dofs_per_component_on_cell;
5045  ++i)
5046  operation.process_dof(vector_ptrs[v][i],
5047  values_dofs[comp][i][v]);
5048  }
5049  else
5050  {
5051  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5052  if (mask[v] == true)
5053  for (unsigned int i = 0;
5054  i < this->data->dofs_per_component_on_cell;
5055  ++i)
5056  operation.process_dof(
5057  vector_ptrs[v]
5058  [i + comp * this->data
5059  ->dofs_per_component_on_cell],
5060  values_dofs[comp][i][v]);
5061  }
5062  }
5063  return;
5064  }
5065 
5066  unsigned int dof_indices[VectorizedArrayType::size()];
5067 
5068  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5069  {
5071  dof_indices[v] =
5072  dof_indices_cont[cells[v]] +
5073  this->dof_info
5074  ->component_dof_indices_offset[this->active_fe_index]
5075  [this->first_selected_component] *
5076  this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
5077  }
5078 
5079  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
5080  dof_indices[v] = numbers::invalid_unsigned_int;
5081 
5082  // In the case with contiguous cell indices, we know that there are no
5083  // constraints and that the indices within each element are contiguous
5084  if (n_filled_lanes == VectorizedArrayType::size() &&
5085  n_lanes == VectorizedArrayType::size() && !is_ecl)
5086  {
5087  if (this->dof_info->index_storage_variants[ind][this->cell] ==
5089  contiguous)
5090  {
5091  if (n_components == 1 || n_fe_components == 1)
5092  for (unsigned int comp = 0; comp < n_components; ++comp)
5093  operation.process_dofs_vectorized_transpose(
5094  this->data->dofs_per_component_on_cell,
5095  dof_indices,
5096  *src[comp],
5097  values_dofs[comp],
5098  vector_selector);
5099  else
5100  operation.process_dofs_vectorized_transpose(
5101  this->data->dofs_per_component_on_cell * n_components,
5102  dof_indices,
5103  *src[0],
5104  &values_dofs[0][0],
5105  vector_selector);
5106  }
5107  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
5109  interleaved_contiguous_strided)
5110  {
5111  if (n_components == 1 || n_fe_components == 1)
5112  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5113  ++i)
5114  {
5115  for (unsigned int comp = 0; comp < n_components; ++comp)
5116  operation.process_dof_gather(dof_indices,
5117  *src[comp],
5118  i * VectorizedArrayType::size(),
5119  values_dofs[comp][i],
5120  vector_selector);
5121  }
5122  else
5123  for (unsigned int comp = 0; comp < n_components; ++comp)
5124  for (unsigned int i = 0;
5125  i < this->data->dofs_per_component_on_cell;
5126  ++i)
5127  {
5128  operation.process_dof_gather(
5129  dof_indices,
5130  *src[0],
5131  (comp * this->data->dofs_per_component_on_cell + i) *
5132  VectorizedArrayType::size(),
5133  values_dofs[comp][i],
5134  vector_selector);
5135  }
5136  }
5137  else
5138  {
5139  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
5141  IndexStorageVariants::interleaved_contiguous_mixed_strides,
5142  ExcNotImplemented());
5143  const unsigned int *offsets =
5144  &this->dof_info->dof_indices_interleave_strides
5145  [ind][VectorizedArrayType::size() * this->cell];
5146  if (n_components == 1 || n_fe_components == 1)
5147  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5148  ++i)
5149  {
5150  for (unsigned int comp = 0; comp < n_components; ++comp)
5151  operation.process_dof_gather(dof_indices,
5152  *src[comp],
5153  0,
5154  values_dofs[comp][i],
5155  vector_selector);
5157  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5158  dof_indices[v] += offsets[v];
5159  }
5160  else
5161  for (unsigned int comp = 0; comp < n_components; ++comp)
5162  for (unsigned int i = 0;
5163  i < this->data->dofs_per_component_on_cell;
5164  ++i)
5165  {
5166  operation.process_dof_gather(dof_indices,
5167  *src[0],
5168  0,
5169  values_dofs[comp][i],
5170  vector_selector);
5172  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5173  dof_indices[v] += offsets[v];
5174  }
5175  }
5176  }
5177  else
5178  for (unsigned int comp = 0; comp < n_components; ++comp)
5179  {
5180  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5181  ++i)
5182  operation.process_empty(values_dofs[comp][i]);
5183  if (this->dof_info->index_storage_variants[ind][this->cell] ==
5185  contiguous)
5186  {
5187  if (n_components == 1 || n_fe_components == 1)
5188  {
5189  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5190  if (mask[v] == true)
5191  for (unsigned int i = 0;
5192  i < this->data->dofs_per_component_on_cell;
5193  ++i)
5194  operation.process_dof(dof_indices[v] + i,
5195  *src[comp],
5196  values_dofs[comp][i][v]);
5197  }
5198  else
5199  {
5200  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5201  if (mask[v] == true)
5202  for (unsigned int i = 0;
5203  i < this->data->dofs_per_component_on_cell;
5204  ++i)
5205  operation.process_dof(
5206  dof_indices[v] + i +
5207  comp * this->data->dofs_per_component_on_cell,
5208  *src[0],
5209  values_dofs[comp][i][v]);
5210  }
5211  }
5212  else
5213  {
5214  const unsigned int *offsets =
5215  &this->dof_info->dof_indices_interleave_strides
5216  [ind][VectorizedArrayType::size() * this->cell];
5217  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5218  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
5219  if (n_components == 1 || n_fe_components == 1)
5220  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5221  {
5222  if (mask[v] == true)
5223  for (unsigned int i = 0;
5224  i < this->data->dofs_per_component_on_cell;
5225  ++i)
5226  operation.process_dof(dof_indices[v] + i * offsets[v],
5227  *src[comp],
5228  values_dofs[comp][i][v]);
5229  }
5230  else
5231  {
5232  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5233  if (mask[v] == true)
5234  for (unsigned int i = 0;
5235  i < this->data->dofs_per_component_on_cell;
5236  ++i)
5237  operation.process_dof(
5238  dof_indices[v] +
5239  (i + comp * this->data->dofs_per_component_on_cell) *
5240  offsets[v],
5241  *src[0],
5242  values_dofs[comp][i][v]);
5243  }
5244  }
5245  }
5246 }
5247 
5248 namespace internal
5249 {
5250  template <typename Number,
5251  typename VectorType,
5252  typename std::enable_if<!IsBlockVector<VectorType>::value,
5253  VectorType>::type * = nullptr>
5254  decltype(std::declval<VectorType>().begin())
5255  get_beginning(VectorType &vec)
5256  {
5257  return vec.begin();
5258  }
5259 
5260  template <typename Number,
5261  typename VectorType,
5262  typename std::enable_if<IsBlockVector<VectorType>::value,
5263  VectorType>::type * = nullptr>
5264  typename VectorType::value_type *
5265  get_beginning(VectorType &)
5266  {
5267  return nullptr;
5268  }
5269 
5270  template <typename VectorType,
5271  typename std::enable_if<has_shared_vector_data<VectorType>::value,
5272  VectorType>::type * = nullptr>
5273  const std::vector<ArrayView<const typename VectorType::value_type>> *
5274  get_shared_vector_data(VectorType & vec,
5275  const bool is_valid_mode_for_sm,
5276  const unsigned int active_fe_index,
5278  {
5279  // note: no hp is supported
5280  if (is_valid_mode_for_sm &&
5281  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
5282  .size() > 0 &&
5283  active_fe_index == 0)
5284  return &vec.shared_vector_data();
5285  else
5286  return nullptr;
5287  }
5288 
5289  template <typename VectorType,
5290  typename std::enable_if<!has_shared_vector_data<VectorType>::value,
5291  VectorType>::type * = nullptr>
5292  const std::vector<ArrayView<const typename VectorType::value_type>> *
5293  get_shared_vector_data(VectorType &,
5294  const bool,
5295  const unsigned int,
5297  {
5298  return nullptr;
5299  }
5300 
5301  template <int n_components, typename VectorType>
5302  std::pair<
5303  std::array<typename internal::BlockVectorSelector<
5304  typename std::remove_const<VectorType>::type,
5306  value>::BaseVectorType *,
5307  n_components>,
5308  std::array<
5309  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
5310  typename std::remove_const<VectorType>::type,
5312  BaseVectorType::value_type>> *,
5313  n_components>>
5314  get_vector_data(VectorType & src,
5315  const unsigned int first_index,
5316  const bool is_valid_mode_for_sm,
5317  const unsigned int active_fe_index,
5319  {
5320  // select between block vectors and non-block vectors. Note that the number
5321  // of components is checked in the internal data
5322  std::pair<
5323  std::array<typename internal::BlockVectorSelector<
5324  typename std::remove_const<VectorType>::type,
5326  value>::BaseVectorType *,
5327  n_components>,
5328  std::array<
5329  const std::vector<
5330  ArrayView<const typename internal::BlockVectorSelector<
5331  typename std::remove_const<VectorType>::type,
5333  value>::BaseVectorType::value_type>> *,
5334  n_components>>
5335  src_data;
5336 
5337  for (unsigned int d = 0; d < n_components; ++d)
5338  src_data.first[d] = internal::BlockVectorSelector<
5339  typename std::remove_const<VectorType>::type,
5340  IsBlockVector<typename std::remove_const<VectorType>::type>::value>::
5341  get_vector_component(
5342  const_cast<typename std::remove_const<VectorType>::type &>(src),
5343  d + first_index);
5344 
5345  for (unsigned int d = 0; d < n_components; ++d)
5346  src_data.second[d] = get_shared_vector_data(*src_data.first[d],
5347  is_valid_mode_for_sm,
5348  active_fe_index,
5349  dof_info);
5350 
5351  return src_data;
5352  }
5353 } // namespace internal
5354 
5355 
5356 
5357 template <int dim,
5358  int n_components_,
5359  typename Number,
5360  bool is_face,
5361  typename VectorizedArrayType>
5362 template <bool transpose>
5363 inline void
5366 {
5367  if (this->dof_info == nullptr ||
5368  this->dof_info->hanging_node_constraint_masks.size() == 0)
5369  return; // nothing to do with faces
5370 
5371  unsigned int n_vectorization_actual =
5372  this->dof_info
5373  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
5374 
5375  constexpr unsigned int n_lanes = VectorizedArrayType::size();
5376  std::array<internal::MatrixFreeFunctions::ConstraintKinds, n_lanes>
5377  constraint_mask;
5378 
5379  bool hn_available = false;
5380 
5381  unsigned int cells_copied[n_lanes];
5382  const unsigned int *cells;
5383 
5384  if (is_face)
5385  {
5386  if (this->dof_access_index ==
5388  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
5389  cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
5390  cells =
5391  this->dof_access_index ==
5393  &cells_copied[0] :
5394  (this->is_interior_face ?
5395  &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
5396  &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
5397  }
5398 
5399  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
5400  {
5401  const unsigned int cell_index =
5402  is_face ? cells[v] : this->cell * n_lanes + v;
5403  const unsigned int cell_dof_index =
5404  cell_index * n_fe_components + first_selected_component;
5405 
5406  const auto mask =
5407  this->dof_info->hanging_node_constraint_masks[cell_dof_index];
5408  constraint_mask[v] = mask;
5409 
5410  hn_available |=
5412  }
5413 
5414  if (hn_available == false)
5415  return; // no hanging node on cell batch -> nothing to do
5416 
5417  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
5418  constraint_mask[v] =
5420 
5422  apply(n_components,
5423  this->data->data.front().fe_degree,
5424  *this,
5425  transpose,
5426  constraint_mask,
5427  values_dofs[0]);
5428 }
5429 
5430 
5431 
5432 template <int dim,
5433  int n_components_,
5434  typename Number,
5435  bool is_face,
5436  typename VectorizedArrayType>
5437 template <typename VectorType>
5438 inline void
5440  read_dof_values(const VectorType &src, const unsigned int first_index)
5441 {
5442  const auto src_data = internal::get_vector_data<n_components_>(
5443  src,
5444  first_index,
5445  this->dof_access_index ==
5447  this->active_fe_index,
5448  this->dof_info);
5449 
5451  read_write_operation(reader,
5452  src_data.first,
5453  src_data.second,
5454  std::bitset<VectorizedArrayType::size()>().flip(),
5455  true);
5456 
5457  apply_hanging_node_constraints<false>();
5458 
5459 # ifdef DEBUG
5460  dof_values_initialized = true;
5461 # endif
5462 }
5463 
5464 
5465 
5466 template <int dim,
5467  int n_components_,
5468  typename Number,
5469  bool is_face,
5470  typename VectorizedArrayType>
5471 template <typename VectorType>
5472 inline void
5474  read_dof_values_plain(const VectorType &src, const unsigned int first_index)
5475 {
5476  const auto src_data = internal::get_vector_data<n_components_>(
5477  src,
5478  first_index,
5479  this->dof_access_index ==
5481  this->active_fe_index,
5482  this->dof_info);
5483 
5485  read_write_operation(reader,
5486  src_data.first,
5487  src_data.second,
5488  std::bitset<VectorizedArrayType::size()>().flip(),
5489  false);
5490 
5491 # ifdef DEBUG
5492  dof_values_initialized = true;
5493 # endif
5494 }
5495 
5496 
5497 
5498 template <int dim,
5499  int n_components_,
5500  typename Number,
5501  bool is_face,
5502  typename VectorizedArrayType>
5503 template <typename VectorType>
5504 inline void
5507  VectorType & dst,
5508  const unsigned int first_index,
5509  const std::bitset<VectorizedArrayType::size()> &mask) const
5510 {
5511 # ifdef DEBUG
5512  Assert(dof_values_initialized == true,
5514 # endif
5515 
5516  apply_hanging_node_constraints<true>();
5517 
5518  const auto dst_data = internal::get_vector_data<n_components_>(
5519  dst,
5520  first_index,
5521  this->dof_access_index ==
5523  this->active_fe_index,
5524  this->dof_info);
5525 
5527  distributor;
5528  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
5529 }
5530 
5531 
5532 
5533 template <int dim,
5534  int n_components_,
5535  typename Number,
5536  bool is_face,
5537  typename VectorizedArrayType>
5538 template <typename VectorType>
5539 inline void
5542  const unsigned int first_index,
5543  const std::bitset<VectorizedArrayType::size()> &mask) const
5544 {
5545 # ifdef DEBUG
5546  Assert(dof_values_initialized == true,
5548 # endif
5549 
5550  const auto dst_data = internal::get_vector_data<n_components_>(
5551  dst,
5552  first_index,
5553  this->dof_access_index ==
5555  this->active_fe_index,
5556  this->dof_info);
5557 
5559  read_write_operation(setter, dst_data.first, dst_data.second, mask);
5560 }
5561 
5562 
5563 
5564 template <int dim,
5565  int n_components_,
5566  typename Number,
5567  bool is_face,
5568  typename VectorizedArrayType>
5569 template <typename VectorType>
5570 inline void
5573  VectorType & dst,
5574  const unsigned int first_index,
5575  const std::bitset<VectorizedArrayType::size()> &mask) const
5576 {
5577 # ifdef DEBUG
5578  Assert(dof_values_initialized == true,
5580 # endif
5581 
5582  const auto dst_data = internal::get_vector_data<n_components_>(
5583  dst,
5584  first_index,
5585  this->dof_access_index ==
5587  this->active_fe_index,
5588  this->dof_info);
5589 
5591  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
5592 }
5593 
5594 
5595 
5596 /*------------------------------ access to data fields ----------------------*/
5597 
5598 
5599 
5600 template <int dim,
5601  int n_components,
5602  typename Number,
5603  bool is_face,
5604  typename VectorizedArrayType>
5605 inline const VectorizedArrayType *
5607  begin_dof_values() const
5608 {
5609  return &values_dofs[0][0];
5610 }
5611 
5612 
5613 
5614 template <int dim,
5615  int n_components,
5616  typename Number,
5617  bool is_face,
5618  typename VectorizedArrayType>
5619 inline VectorizedArrayType *
5622 {
5623 # ifdef DEBUG
5624  dof_values_initialized = true;
5625 # endif
5626  return &values_dofs[0][0];
5627 }
5628 
5629 
5630 
5631 template <int dim,
5632  int n_components,
5633  typename Number,
5634  bool is_face,
5635  typename VectorizedArrayType>
5636 inline const VectorizedArrayType *
5638  begin_values() const
5639 {
5640 # ifdef DEBUG
5641  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
5642 # endif
5643  return values_quad;
5644 }
5645 
5646 
5647 
5648 template <int dim,
5649  int n_components,
5650  typename Number,
5651  bool is_face,
5652  typename VectorizedArrayType>
5653 inline VectorizedArrayType *
5655  begin_values()
5656 {
5657 # ifdef DEBUG
5658  values_quad_initialized = true;
5659  values_quad_submitted = true;
5660 # endif
5661  return values_quad;
5662 }
5663 
5664 
5665 
5666 template <int dim,
5667  int n_components,
5668  typename Number,
5669  bool is_face,
5670  typename VectorizedArrayType>
5671 inline const VectorizedArrayType *
5673  begin_gradients() const
5674 {
5675 # ifdef DEBUG
5676  Assert(gradients_quad_initialized || gradients_quad_submitted,
5677  ExcNotInitialized());
5678 # endif
5679  return gradients_quad;
5680 }
5681 
5682 
5683 
5684 template <int dim,
5685  int n_components,
5686  typename Number,
5687  bool is_face,
5688  typename VectorizedArrayType>
5689 inline VectorizedArrayType *
5692 {
5693 # ifdef DEBUG
5694  gradients_quad_submitted = true;
5695  gradients_quad_initialized = true;
5696 # endif
5697  return gradients_quad;
5698 }
5699 
5700 
5701 
5702 template <int dim,
5703  int n_components,
5704  typename Number,
5705  bool is_face,
5706  typename VectorizedArrayType>
5707 inline const VectorizedArrayType *
5709  begin_hessians() const
5710 {
5711 # ifdef DEBUG
5712  Assert(hessians_quad_initialized, ExcNotInitialized());
5713 # endif
5714  return hessians_quad;
5715 }
5716 
5717 
5718 
5719 template <int dim,
5720  int n_components,
5721  typename Number,
5722  bool is_face,
5723  typename VectorizedArrayType>
5724 inline VectorizedArrayType *
5727 {
5728 # ifdef DEBUG
5729  hessians_quad_initialized = true;
5730 # endif
5731  return hessians_quad;
5732 }
5733 
5734 
5735 
5736 template <int dim,
5737  int n_components,
5738  typename Number,
5739  bool is_face,
5740  typename VectorizedArrayType>
5741 inline unsigned int
5744 {
5745  return first_selected_component;
5746 }
5747 
5748 
5749 
5750 template <int dim,
5751  int n_components_,
5752  typename Number,
5753  bool is_face,
5754  typename VectorizedArrayType>
5757  get_dof_value(const unsigned int dof) const
5758 {
5759  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5761  for (unsigned int comp = 0; comp < n_components; ++comp)
5762  return_value[comp] = this->values_dofs[comp][dof];
5763  return return_value;
5764 }
5765 
5766 
5767 
5768 template <int dim,
5769  int n_components_,
5770  typename Number,
5771  bool is_face,
5772  typename VectorizedArrayType>
5775  get_value(const unsigned int q_point) const
5776 {
5777 # ifdef DEBUG
5778  Assert(this->values_quad_initialized == true,
5780 # endif
5781 
5782  AssertIndexRange(q_point, this->n_quadrature_points);
5783  const std::size_t nqp = this->n_quadrature_points;
5785  for (unsigned int comp = 0; comp < n_components; ++comp)
5786  return_value[comp] = values_quad[comp * nqp + q_point];
5787  return return_value;
5788 }
5789 
5790 
5791 
5792 template <int dim,
5793  int n_components_,
5794  typename Number,
5795  bool is_face,
5796  typename VectorizedArrayType>
5797 inline DEAL_II_ALWAYS_INLINE
5800  get_gradient(const unsigned int q_point) const
5801 {
5802 # ifdef DEBUG
5803  Assert(this->gradients_quad_initialized == true,
5805 # endif
5806 
5807  AssertIndexRange(q_point, this->n_quadrature_points);
5808  Assert(this->jacobian != nullptr,
5810  "update_gradients"));
5811  const std::size_t nqp = this->n_quadrature_points;
5813 
5814  // Cartesian cell
5815  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5816  {
5817  for (unsigned int d = 0; d < dim; ++d)
5818  for (unsigned int comp = 0; comp < n_components; ++comp)
5819  grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] *
5820  this->jacobian[0][d][d];
5821  }
5822  // cell with general/affine Jacobian
5823  else
5824  {
5826  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5827  q_point :
5828  0];
5829  for (unsigned int comp = 0; comp < n_components; ++comp)
5830  for (unsigned int d = 0; d < dim; ++d)
5831  {
5832  grad_out[comp][d] =
5833  jac[d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5834  for (unsigned int e = 1; e < dim; ++e)
5835  grad_out[comp][d] +=
5836  jac[d][e] * gradients_quad[(comp * dim + e) * nqp + q_point];
5837  }
5838  }
5839  return grad_out;
5840 }
5841 
5842 
5843 
5844 template <int dim,
5845  int n_components_,
5846  typename Number,
5847  bool is_face,
5848  typename VectorizedArrayType>
5851  get_normal_derivative(const unsigned int q_point) const
5852 {
5853  AssertIndexRange(q_point, this->n_quadrature_points);
5854 # ifdef DEBUG
5855  Assert(this->gradients_quad_initialized == true,
5857 # endif
5858 
5859  Assert(this->normal_x_jacobian != nullptr,
5861  "update_gradients"));
5862 
5863  const std::size_t nqp = this->n_quadrature_points;
5865 
5866  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5867  for (unsigned int comp = 0; comp < n_components; ++comp)
5868  grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5869  (this->normal_x_jacobian[0][dim - 1]);
5870  else
5871  {
5872  const std::size_t index =
5873  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5874  for (unsigned int comp = 0; comp < n_components; ++comp)
5875  {
5876  grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5877  this->normal_x_jacobian[index][0];
5878  for (unsigned int d = 1; d < dim; ++d)
5879  grad_out[comp] += gradients_quad[(comp * dim + d) * nqp + q_point] *
5880  this->normal_x_jacobian[index][d];
5881  }
5882  }
5883  return grad_out;
5884 }
5885 
5886 
5887 
5888 namespace internal
5889 {
5890  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5891  // store the lower diagonal because of symmetry
5892  template <typename VectorizedArrayType>
5893  inline void
5894  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
5895  const VectorizedArrayType *const hessians,
5896  const unsigned int,
5897  VectorizedArrayType (&tmp)[1][1])
5898  {
5899  tmp[0][0] = jac[0][0] * hessians[0];
5900  }
5901 
5902  template <typename VectorizedArrayType>
5903  inline void
5904  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
5905  const VectorizedArrayType *const hessians,
5906  const unsigned int nqp,
5907  VectorizedArrayType (&tmp)[2][2])
5908  {
5909  for (unsigned int d = 0; d < 2; ++d)
5910  {
5911  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
5912  tmp[1][d] =
5913  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
5914  }
5915  }
5916 
5917  template <typename VectorizedArrayType>
5918  inline void
5919  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
5920  const VectorizedArrayType *const hessians,
5921  const unsigned int nqp,
5922  VectorizedArrayType (&tmp)[3][3])
5923  {
5924  for (unsigned int d = 0; d < 3; ++d)
5925  {
5926  tmp[0][d] =
5927  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
5928  jac[d][2] * hessians[4 * nqp]);
5929  tmp[1][d] =
5930  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
5931  jac[d][2] * hessians[5 * nqp]);
5932  tmp[2][d] =
5933  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
5934  jac[d][2] * hessians[2 * nqp]);
5935  }
5936  }
5937 } // namespace internal
5938 
5939 
5940 
5941 template <int dim,
5942  int n_components_,
5943  typename Number,
5944  bool is_face,
5945  typename VectorizedArrayType>
5948  get_hessian(const unsigned int q_point) const
5949 {
5950 # ifdef DEBUG
5951  Assert(this->hessians_quad_initialized == true,
5953 # endif
5954  AssertIndexRange(q_point, this->n_quadrature_points);
5955 
5956  Assert(this->jacobian != nullptr,
5958  "update_hessian"));
5960  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5961  0 :
5962  q_point];
5963 
5965 
5966  const std::size_t nqp = this->n_quadrature_points;
5967  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5968 
5969  // Cartesian cell
5970  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5971  {
5972  for (unsigned int comp = 0; comp < n_components; ++comp)
5973  {
5974  for (unsigned int d = 0; d < dim; ++d)
5975  hessian_out[comp][d][d] =
5976  hessians_quad[(comp * hdim + d) * nqp + q_point] *
5977  (jac[d][d] * jac[d][d]);
5978  switch (dim)
5979  {
5980  case 1:
5981  break;
5982  case 2:
5983  hessian_out[comp][0][1] =
5984  hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5985  (jac[0][0] * jac[1][1]);
5986  break;
5987  case 3:
5988  hessian_out[comp][0][1] =
5989  hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5990  (jac[0][0] * jac[1][1]);
5991  hessian_out[comp][0][2] =
5992  hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5993  (jac[0][0] * jac[2][2]);
5994  hessian_out[comp][1][2] =
5995  hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5996  (jac[1][1] * jac[2][2]);
5997  break;
5998  default:
5999  Assert(false, ExcNotImplemented());
6000  }
6001  for (unsigned int d = 0; d < dim; ++d)
6002  for (unsigned int e = d + 1; e < dim; ++e)
6003  hessian_out[comp][e][d] = hessian_out[comp][d][e];
6004  }
6005  }
6006  // cell with general Jacobian, but constant within the cell
6007  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6008  {
6009  for (unsigned int comp = 0; comp < n_components; ++comp)
6010  {
6011  VectorizedArrayType tmp[dim][dim];
6012  internal::hessian_unit_times_jac(
6013  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
6014 
6015  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
6016  for (unsigned int d = 0; d < dim; ++d)
6017  for (unsigned int e = d; e < dim; ++e)
6018  {
6019  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
6020  for (unsigned int f = 1; f < dim; ++f)
6021  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
6022  }
6023 
6024  // no J' * grad(u) part here because the Jacobian is constant
6025  // throughout the cell and hence, its derivative is zero
6026 
6027  // take symmetric part
6028  for (unsigned int d = 0; d < dim; ++d)
6029  for (unsigned int e = d + 1; e < dim; ++e)
6030  hessian_out[comp][e][d] = hessian_out[comp][d][e];
6031  }
6032  }
6033  // cell with general Jacobian
6034  else
6035  {
6036  const auto &jac_grad =
6037  this->mapping_data->jacobian_gradients
6038  [1 - this->is_interior_face]
6039  [this->mapping_data->data_index_offsets[this->cell] + q_point];
6040  for (unsigned int comp = 0; comp < n_components; ++comp)
6041  {
6042  VectorizedArrayType tmp[dim][dim];
6043  internal::hessian_unit_times_jac(
6044  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
6045 
6046  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
6047  for (unsigned int d = 0; d < dim; ++d)
6048  for (unsigned int e = d; e < dim; ++e)
6049  {
6050  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
6051  for (unsigned int f = 1; f < dim; ++f)
6052  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
6053  }
6054 
6055  // add diagonal part of J' * grad(u)
6056  for (unsigned int d = 0; d < dim; ++d)
6057  for (unsigned int e = 0; e < dim; ++e)
6058  hessian_out[comp][d][d] +=
6059  jac_grad[d][e] *
6060  gradients_from_hessians_quad[(comp * dim + e) * nqp + q_point];
6061 
6062  // add off-diagonal part of J' * grad(u)
6063  for (unsigned int d = 0, count = dim; d < dim; ++d)
6064  for (unsigned int e = d + 1; e < dim; ++e, ++count)
6065  for (unsigned int f = 0; f < dim; ++f)
6066  hessian_out[comp][d][e] +=
6067  jac_grad[count][f] *
6068  gradients_from_hessians_quad[(comp * dim + f) * nqp +
6069  q_point];
6070 
6071  // take symmetric part
6072  for (unsigned int d = 0; d < dim; ++d)
6073  for (unsigned int e = d + 1; e < dim; ++e)
6074  hessian_out[comp][e][d] = hessian_out[comp][d][e];
6075  }
6076  }
6077  return hessian_out;
6078 }
6079 
6080 
6081 
6082 template <int dim,
6083  int n_components_,
6084  typename Number,
6085  bool is_face,
6086  typename VectorizedArrayType>
6089  get_hessian_diagonal(const unsigned int q_point) const
6090 {
6091  Assert(!is_face, ExcNotImplemented());
6092 # ifdef DEBUG
6093  Assert(this->hessians_quad_initialized == true,
6095 # endif
6096  AssertIndexRange(q_point, this->n_quadrature_points);
6097 
6098  Assert(this->jacobian != nullptr, ExcNotImplemented());
6100  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
6101  0 :
6102  q_point];
6103 
6104  const std::size_t nqp = this->n_quadrature_points;
6105  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
6107 
6108  // Cartesian cell
6109  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
6110  {
6111  for (unsigned int comp = 0; comp < n_components; ++comp)
6112  for (unsigned int d = 0; d < dim; ++d)
6113  hessian_out[comp][d] =
6114  hessians_quad[(comp * hdim + d) * nqp + q_point] *
6115  (jac[d][d] * jac[d][d]);
6116  }
6117  // cell with general Jacobian, but constant within the cell
6118  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
6119  {
6120  for (unsigned int comp = 0; comp < n_components; ++comp)
6121  {
6122  // compute laplacian before the gradient because it needs to access
6123  // unscaled gradient data
6124  VectorizedArrayType tmp[dim][dim];
6125  internal::hessian_unit_times_jac(
6126  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
6127 
6128  // compute only the trace part of hessian, J * tmp = J *
6129  // hess_unit(u) * J^T
6130  for (unsigned int d = 0; d < dim; ++d)
6131  {
6132  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
6133  for (unsigned int f = 1; f < dim; ++f)
6134  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
6135  }
6136  }
6137  }
6138  // cell with general Jacobian
6139  else
6140  {
6141  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
6142  &jac_grad =
6143  this->mapping_data->jacobian_gradients
6144  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
6145  for (unsigned int comp = 0; comp < n_components; ++comp)
6146  {
6147  // compute laplacian before the gradient because it needs to access
6148  // unscaled gradient data
6149  VectorizedArrayType tmp[dim][dim];
6150  internal::hessian_unit_times_jac(
6151  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
6152 
6153  // compute only the trace part of hessian, J * tmp = J *
6154  // hess_unit(u) * J^T
6155  for (unsigned int d = 0; d < dim; ++d)
6156  {
6157  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
6158  for (unsigned int f = 1; f < dim; ++f)
6159  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
6160  }
6161 
6162  for (unsigned int d = 0; d < dim; ++d)
6163  for (unsigned int e = 0; e < dim; ++e)
6164  hessian_out[comp][d] +=
6165  jac_grad[d][e] *
6166  gradients_from_hessians_quad[(comp * dim + e) * nqp + q_point];
6167  }
6168  }
6169  return hessian_out;
6170 }
6171 
6172 
6173 
6174 template <int dim,
6175  int n_components_,
6176  typename Number,
6177  bool is_face,
6178  typename VectorizedArrayType>
6181  get_laplacian(const unsigned int q_point) const
6182 {
6183  Assert(is_face == false, ExcNotImplemented());
6184 # ifdef DEBUG
6185  Assert(this->hessians_quad_initialized == true,
6187 # endif
6188  AssertIndexRange(q_point, this->n_quadrature_points);
6189 
6191  const auto hess_diag = get_hessian_diagonal(q_point);
6192  for (unsigned int comp = 0; comp < n_components; ++comp)
6193  {
6194  laplacian_out[comp] = hess_diag[comp][0];
6195  for (unsigned int d = 1; d < dim; ++d)
6196  laplacian_out[comp] += hess_diag[comp][d];
6197  }
6198  return laplacian_out;
6199 }
6200 
6201 
6202 
6203 template <int dim,
6204  int n_components_,
6205  typename Number,
6206  bool is_face,
6207  typename VectorizedArrayType>
6208 inline DEAL_II_ALWAYS_INLINE void
6211  const unsigned int dof)
6212 {
6213 # ifdef DEBUG
6214  this->dof_values_initialized = true;
6215 # endif
6216  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6217  for (unsigned int comp = 0; comp < n_components; ++comp)
6218  this->values_dofs[comp][dof] = val_in[comp];
6219 }
6220 
6221 
6222 
6223 template <int dim,
6224  int n_components_,
6225  typename Number,
6226  bool is_face,
6227  typename VectorizedArrayType>
6228 inline DEAL_II_ALWAYS_INLINE void
6231  const unsigned int q_point)
6232 {
6234  AssertIndexRange(q_point, this->n_quadrature_points);
6235  Assert(this->J_value != nullptr,
6237  "update_values"));
6238 # ifdef DEBUG
6239  this->values_quad_submitted = true;
6240 # endif
6241 
6242  const std::size_t nqp = this->n_quadrature_points;
6243  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6244  {
6245  const VectorizedArrayType JxW =
6246  this->J_value[0] * this->quadrature_weights[q_point];
6247  for (unsigned int comp = 0; comp < n_components; ++comp)
6248  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6249  }
6250  else
6251  {
6252  const VectorizedArrayType JxW = this->J_value[q_point];
6253  for (unsigned int comp = 0; comp < n_components; ++comp)
6254  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6255  }
6256 }
6257 
6258 
6259 
6260 template <int dim,
6261  int n_components_,
6262  typename Number,
6263  bool is_face,
6264  typename VectorizedArrayType>
6265 inline DEAL_II_ALWAYS_INLINE void
6268  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
6269  const unsigned int q_point)
6270 {
6272  AssertIndexRange(q_point, this->n_quadrature_points);
6273  Assert(this->J_value != nullptr,
6275  "update_gradients"));
6276  Assert(this->jacobian != nullptr,
6278  "update_gradients"));
6279 # ifdef DEBUG
6280  this->gradients_quad_submitted = true;
6281 # endif
6282 
6283  const std::size_t nqp = this->n_quadrature_points;
6284  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6285  {
6286  const VectorizedArrayType JxW =
6287  this->J_value[0] * this->quadrature_weights[q_point];
6288  for (unsigned int d = 0; d < dim; ++d)
6289  {
6290  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
6291  for (unsigned int comp = 0; comp < n_components; ++comp)
6292  gradients_quad[(comp * dim + d) * nqp + q_point] =
6293  grad_in[comp][d] * factor;
6294  }
6295  }
6296  else
6297  {
6299  this->cell_type > internal::MatrixFreeFunctions::affine ?
6300  this->jacobian[q_point] :
6301  this->jacobian[0];
6302  const VectorizedArrayType JxW =
6303  this->cell_type > internal::MatrixFreeFunctions::affine ?
6304  this->J_value[q_point] :
6305  this->J_value[0] * this->quadrature_weights[q_point];
6306  for (unsigned int comp = 0; comp < n_components; ++comp)
6307  for (unsigned int d = 0; d < dim; ++d)
6308  {
6309  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
6310  for (unsigned int e = 1; e < dim; ++e)
6311  new_val += (jac[e][d] * grad_in[comp][e]);
6312  gradients_quad[(comp * dim + d) * nqp + q_point] = new_val * JxW;
6313  }
6314  }
6315 }
6316 
6317 
6318 
6319 template <int dim,
6320  int n_components_,
6321  typename Number,
6322  bool is_face,
6323  typename VectorizedArrayType>
6324 inline DEAL_II_ALWAYS_INLINE void
6328  const unsigned int q_point)
6329 {
6330  AssertIndexRange(q_point, this->n_quadrature_points);
6331  Assert(this->normal_x_jacobian != nullptr,
6333  "update_gradients"));
6334 # ifdef DEBUG
6335  this->gradients_quad_submitted = true;
6336 # endif
6337 
6338  const std::size_t nqp = this->n_quadrature_points;
6339  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
6340  for (unsigned int comp = 0; comp < n_components; ++comp)
6341  {
6342  for (unsigned int d = 0; d < dim - 1; ++d)
6343  gradients_quad[(comp * dim + d) * nqp + q_point] =
6344  VectorizedArrayType();
6345  gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
6346  grad_in[comp] *
6347  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
6348  this->quadrature_weights[q_point]);
6349  }
6350  else
6351  {
6352  const unsigned int index =
6353  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
6355  this->normal_x_jacobian[index];
6356  for (unsigned int comp = 0; comp < n_components; ++comp)
6357  {
6358  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
6359  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6360  factor = factor * this->quadrature_weights[q_point];
6361  for (unsigned int d = 0; d < dim; ++d)
6362  gradients_quad[(comp * dim + d) * nqp + q_point] = factor * jac[d];
6363  }
6364  }
6365 }
6366 
6367 
6368 
6369 template <int dim,
6370  int n_components_,
6371  typename Number,
6372  bool is_face,
6373  typename VectorizedArrayType>
6374 inline DEAL_II_ALWAYS_INLINE void
6377  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
6378  hessian_in,
6379  const unsigned int q_point)
6380 {
6382  AssertIndexRange(q_point, this->n_quadrature_points);
6383  Assert(this->J_value != nullptr,
6385  "update_hessians"));
6386  Assert(this->jacobian != nullptr,
6388  "update_hessians"));
6389 # ifdef DEBUG
6390  this->hessians_quad_submitted = true;
6391 # endif
6392 
6393  // compute hessian_unit = J^T * hessian_in(u) * J
6394  const std::size_t nqp = this->n_quadrature_points;
6395  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
6396  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6397  {
6398  const VectorizedArrayType JxW =
6399  this->J_value[0] * this->quadrature_weights[q_point];
6400 
6401  // diagonal part
6402  for (unsigned int d = 0; d < dim; ++d)
6403  {
6404  const auto jac_d = this->jacobian[0][d][d];
6405  const VectorizedArrayType factor = jac_d * jac_d * JxW;
6406  for (unsigned int comp = 0; comp < n_components; ++comp)
6407  hessians_quad[(comp * hdim + d) * nqp + q_point] =
6408  hessian_in[comp][d][d] * factor;
6409  }
6410 
6411  // off diagonal part
6412  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
6413  for (unsigned int e = 0; e < d; ++e, ++off_dia)
6414  {
6415  const auto jac_d = this->jacobian[0][d][d];
6416  const auto jac_e = this->jacobian[0][e][e];
6417  const VectorizedArrayType factor = jac_d * jac_e * JxW;
6418  for (unsigned int comp = 0; comp < n_components; ++comp)
6419  hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
6420  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
6421  }
6422  }
6423  // cell with general Jacobian, but constant within the cell
6424  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6425  {
6426  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
6427  const VectorizedArrayType JxW =
6428  this->J_value[0] * this->quadrature_weights[q_point];
6429  for (unsigned int comp = 0; comp < n_components; ++comp)
6430  {
6431  // 1. tmp = hessian_in(u) * J
6432  VectorizedArrayType tmp[dim][dim];
6433  for (unsigned int i = 0; i < dim; ++i)
6434  for (unsigned int j = 0; j < dim; ++j)
6435  {
6436  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
6437  for (unsigned int k = 1; k < dim; ++k)
6438  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
6439  }
6440 
6441  // 2. hessian_unit = J^T * tmp
6442  VectorizedArrayType tmp2[dim][dim];
6443  for (unsigned int i = 0; i < dim; ++i)
6444  for (unsigned int j = 0; j < dim; ++j)
6445  {
6446  tmp2[i][j] = jac[0][i] * tmp[0][j];
6447  for (unsigned int k = 1; k < dim; ++k)
6448  tmp2[i][j] += jac[k][i] * tmp[k][j];
6449  }
6450 
6451  // diagonal part
6452  for (unsigned int d = 0; d < dim; ++d)
6453  hessians_quad[(comp * hdim + d) * nqp + q_point] = tmp2[d][d] * JxW;
6454 
6455  // off diagonal part
6456  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
6457  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
6458  hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
6459  (tmp2[d][e] + tmp2[e][d]) * JxW;
6460  }
6461  }
6462  else
6463  {
6464  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
6465  const VectorizedArrayType JxW = this->J_value[q_point];
6466  const auto & jac_grad =
6467  this->mapping_data->jacobian_gradients
6468  [1 - this->is_interior_face]
6469  [this->mapping_data->data_index_offsets[this->cell] + q_point];
6470  for (unsigned int comp = 0; comp < n_components; ++comp)
6471  {
6472  // 1. tmp = hessian_in(u) * J
6473  VectorizedArrayType tmp[dim][dim];
6474  for (unsigned int i = 0; i < dim; ++i)
6475  for (unsigned int j = 0; j < dim; ++j)
6476  {
6477  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
6478  for (unsigned int k = 1; k < dim; ++k)
6479  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
6480  }
6481 
6482  // 2. hessian_unit = J^T * tmp
6483  VectorizedArrayType tmp2[dim][dim];
6484  for (unsigned int i = 0; i < dim; ++i)
6485  for (unsigned int j = 0; j < dim; ++j)
6486  {
6487  tmp2[i][j] = jac[0][i] * tmp[0][j];
6488  for (unsigned int k = 1; k < dim; ++k)
6489  tmp2[i][j] += jac[k][i] * tmp[k][j];
6490  }
6491 
6492  // diagonal part
6493  for (unsigned int d = 0; d < dim; ++d)
6494  hessians_quad[(comp * hdim + d) * nqp + q_point] = tmp2[d][d] * JxW;
6495 
6496  // off diagonal part
6497  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
6498  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
6499  hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
6500  (tmp2[d][e] + tmp2[e][d]) * JxW;
6501 
6502  // 3. gradient_unit = J' ** hessian_in
6503  for (unsigned int d = 0; d < dim; ++d)
6504  {
6505  VectorizedArrayType sum = 0;
6506  for (unsigned int e = 0; e < dim; ++e)
6507  sum += hessian_in[comp][e][e] * jac_grad[e][d];
6508  for (unsigned int e = 0, count = dim; e < dim; ++e)
6509  for (unsigned int f = e + 1; f < dim; ++f, ++count)
6510  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
6511  jac_grad[count][d];
6512  gradients_from_hessians_quad[(comp * dim + d) * nqp + q_point] =
6513  sum * JxW;
6514  }
6515  }
6516  }
6517 }
6518 
6519 
6520 
6521 template <int dim,
6522  int n_components_,
6523  typename Number,
6524  bool is_face,
6525  typename VectorizedArrayType>
6528  integrate_value() const
6529 {
6531 # ifdef DEBUG
6532  Assert(this->values_quad_submitted == true,
6534 # endif
6535 
6537  const std::size_t nqp = this->n_quadrature_points;
6538  for (unsigned int q = 0; q < nqp; ++q)
6539  for (unsigned int comp = 0; comp < n_components; ++comp)
6540  return_value[comp] += this->values_quad[comp * nqp + q];
6541  return (return_value);
6542 }
6543 
6544 
6545 
6546 /*----------------------- FEEvaluationAccess --------------------------------*/
6547 
6548 
6549 template <int dim,
6550  int n_components_,
6551  typename Number,
6552  bool is_face,
6553  typename VectorizedArrayType>
6554 inline FEEvaluationAccess<dim,
6555  n_components_,
6556  Number,
6557  is_face,
6558  VectorizedArrayType>::
6559  FEEvaluationAccess(
6561  const unsigned int dof_no,
6562  const unsigned int first_selected_component,
6563  const unsigned int quad_no_in,
6564  const unsigned int fe_degree,
6565  const unsigned int n_q_points,
6566  const bool is_interior_face,
6567  const unsigned int active_fe_index,
6568  const unsigned int active_quad_index,
6569  const unsigned int face_type)
6571  data_in,
6572  dof_no,
6573  first_selected_component,
6574  quad_no_in,
6575  fe_degree,
6576  n_q_points,
6577  is_interior_face,
6578  active_fe_index,
6579  active_quad_index,
6580  face_type)
6581 {}
6582 
6583 
6584 
6585 template <int dim,
6586  int n_components_,
6587  typename Number,
6588  bool is_face,
6589  typename VectorizedArrayType>
6590 inline FEEvaluationAccess<dim,
6591  n_components_,
6592  Number,
6593  is_face,
6594  VectorizedArrayType>::
6595  FEEvaluationAccess(
6596  const Mapping<dim> & mapping,
6597  const FiniteElement<dim> &fe,
6598  const Quadrature<1> & quadrature,
6599  const UpdateFlags update_flags,
6600  const unsigned int first_selected_component,
6602  *other)
6604  mapping,
6605  fe,
6606  quadrature,
6607  update_flags,
6608  first_selected_component,
6609  other)
6610 {}
6611 
6612 
6613 
6614 template <int dim,
6615  int n_components_,
6616  typename Number,
6617  bool is_face,
6618  typename VectorizedArrayType>
6619 inline FEEvaluationAccess<dim,
6620  n_components_,
6621  Number,
6622  is_face,
6623  VectorizedArrayType>::
6624  FEEvaluationAccess(const FEEvaluationAccess<dim,
6625  n_components_,
6626  Number,
6627  is_face,
6628  VectorizedArrayType> &other)
6630  other)
6631 {}
6632 
6633 
6634 
6635 template <int dim,
6636  int n_components_,
6637  typename Number,
6638  bool is_face,
6639  typename VectorizedArrayType>
6640 inline FEEvaluationAccess<dim,
6641  n_components_,
6642  Number,
6643  is_face,
6644  VectorizedArrayType> &
6646 operator=(const FEEvaluationAccess<dim,
6647  n_components_,
6648  Number,
6649  is_face,
6650  VectorizedArrayType> &other)
6651 {
6652  this->FEEvaluationBase<dim,
6653  n_components_,
6654  Number,
6655  is_face,
6656  VectorizedArrayType>::operator=(other);
6657  return *this;
6658 }
6659 
6660 
6661 
6662 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
6663 
6664 
6665 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6669  const unsigned int dof_no,
6670  const unsigned int first_selected_component,
6671  const unsigned int quad_no_in,
6672  const unsigned int fe_degree,
6673  const unsigned int n_q_points,
6674  const bool is_interior_face,
6675  const unsigned int active_fe_index,
6676  const unsigned int active_quad_index,
6677  const unsigned int face_type)
6679  data_in,
6680  dof_no,
6681  first_selected_component,
6682  quad_no_in,
6683  fe_degree,
6684  n_q_points,
6685  is_interior_face,
6686  active_fe_index,
6687  active_quad_index,
6688  face_type)
6689 {}
6690 
6691 
6692 
6693 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6696  const Mapping<dim> & mapping,
6697  const FiniteElement<dim> &fe,
6698  const Quadrature<1> & quadrature,
6699  const UpdateFlags update_flags,
6700  const unsigned int first_selected_component,
6702  *other)
6704  mapping,
6705  fe,
6706  quadrature,
6707  update_flags,
6708  first_selected_component,
6709  other)
6710 {}
6711 
6712 
6713 
6714 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6718  &other)
6720 {}
6721 
6722 
6723 
6724 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6728 {
6729  this
6730  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
6731  other);
6732  return *this;
6733 }
6734 
6735 
6736 
6737 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6738 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6740  const unsigned int dof) const
6741 {
6742  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6743  return this->values_dofs[0][dof];
6744 }
6745 
6746 
6747 
6748 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6749 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6751  const unsigned int q_point) const
6752 {
6753 # ifdef DEBUG
6754  Assert(this->values_quad_initialized == true,
6756 # endif
6757  AssertIndexRange(q_point, this->n_quadrature_points);
6758  return this->values_quad[q_point];
6759 }
6760 
6761 
6762 
6763 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6764 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6766  get_normal_derivative(const unsigned int q_point) const
6767 {
6768  return BaseClass::get_normal_derivative(q_point)[0];
6769 }
6770 
6771 
6772 
6773 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6776  const unsigned int q_point) const
6777 {
6778  // could use the base class gradient, but that involves too many expensive
6779  // initialization operations on tensors
6780 
6781 # ifdef DEBUG
6782  Assert(this->gradients_quad_initialized == true,
6784 # endif
6785  AssertIndexRange(q_point, this->n_quadrature_points);
6786 
6787  Assert(this->jacobian != nullptr,
6789  "update_gradients"));
6790 
6792 
6793  const std::size_t nqp = this->n_quadrature_points;
6794  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6795  {
6796  for (unsigned int d = 0; d < dim; ++d)
6797  grad_out[d] =
6798  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
6799  }
6800  // cell with general/affine Jacobian
6801  else
6802  {
6804  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
6805  q_point :
6806  0];
6807  for (unsigned int d = 0; d < dim; ++d)
6808  {
6809  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
6810  for (unsigned int e = 1; e < dim; ++e)
6811  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
6812  }
6813  }
6814  return grad_out;
6815 }
6816 
6817 
6818 
6819 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6822  const unsigned int q_point) const
6823 {
6824  return BaseClass::get_hessian(q_point)[0];
6825 }
6826 
6827 
6828 
6829 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6832  get_hessian_diagonal(const unsigned int q_point) const
6833 {
6834  return BaseClass::get_hessian_diagonal(q_point)[0];
6835 }
6836 
6837 
6838 
6839 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6840 inline VectorizedArrayType
6842  const unsigned int q_point) const
6843 {
6844  return BaseClass::get_laplacian(q_point)[0];
6845 }
6846 
6847 
6848 
6849 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6850 inline void DEAL_II_ALWAYS_INLINE
6852  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6853 {
6854 # ifdef DEBUG
6855  this->dof_values_initialized = true;
6856  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6857 # endif
6858  this->values_dofs[0][dof] = val_in;
6859 }
6860 
6861 
6862 
6863 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6864 inline void DEAL_II_ALWAYS_INLINE
6866  const VectorizedArrayType val_in,
6867  const unsigned int q_point)
6868 {
6870  AssertIndexRange(q_point, this->n_quadrature_points);
6871  Assert(this->J_value != nullptr,
6873  "update_value"));
6874 # ifdef DEBUG
6875  this->values_quad_submitted = true;
6876 # endif
6877 
6878  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6879  {
6880  const VectorizedArrayType JxW =
6881  this->J_value[0] * this->quadrature_weights[q_point];
6882  this->values_quad[q_point] = val_in * JxW;
6883  }
6884  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
6885  {
6886  this->values_quad[q_point] = val_in * this->J_value[q_point];
6887  }
6888 }
6889 
6890 
6891 
6892 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6893 inline DEAL_II_ALWAYS_INLINE void
6895  const Tensor<1, 1, VectorizedArrayType> val_in,
6896  const unsigned int q_point)
6897 {
6898  submit_value(val_in[0], q_point);
6899 }
6900 
6901 
6902 
6903 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6904 inline DEAL_II_ALWAYS_INLINE void
6906  submit_normal_derivative(const VectorizedArrayType grad_in,
6907  const unsigned int q_point)
6908 {
6910  grad[0] = grad_in;
6911  BaseClass::submit_normal_derivative(grad, q_point);
6912 }
6913 
6914 
6915 
6916 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6917 inline DEAL_II_ALWAYS_INLINE void
6920  const unsigned int q_point)
6921 {
6923  AssertIndexRange(q_point, this->n_quadrature_points);
6924  Assert(this->J_value != nullptr,
6926  "update_gradients"));
6927  Assert(this->jacobian != nullptr,
6929  "update_gradients"));
6930 # ifdef DEBUG
6931  this->gradients_quad_submitted = true;
6932 # endif
6933 
6934  const std::size_t nqp = this->n_quadrature_points;
6935  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6936  {
6937  const VectorizedArrayType JxW =
6938  this->J_value[0] * this->quadrature_weights[q_point];
6939  for (unsigned int d = 0; d < dim; ++d)
6940  this->gradients_quad[d * nqp + q_point] =
6941  (grad_in[d] * this->jacobian[0][d][d] * JxW);
6942  }
6943  // general/affine cell type
6944  else
6945  {
6947  this->cell_type > internal::MatrixFreeFunctions::affine ?
6948  this->jacobian[q_point] :
6949  this->jacobian[0];
6950  const VectorizedArrayType JxW =
6951  this->cell_type > internal::MatrixFreeFunctions::affine ?
6952  this->J_value[q_point] :
6953  this->J_value[0] * this->quadrature_weights[q_point];
6954  for (unsigned int d = 0; d < dim; ++d)
6955  {
6956  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
6957  for (unsigned int e = 1; e < dim; ++e)
6958  new_val += jac[e][d] * grad_in[e];
6959  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
6960  }
6961  }
6962 }
6963 
6964 
6965 
6966 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6967 inline DEAL_II_ALWAYS_INLINE void
6970  const unsigned int q_point)
6971 {
6973  hessian[0] = hessian_in;
6974  BaseClass::submit_hessian(hessian, q_point);
6975 }
6976 
6977 
6978 
6979 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6980 inline VectorizedArrayType
6982  integrate_value() const
6983 {
6984  return BaseClass::integrate_value()[0];
6985 }
6986 
6987 
6988 
6989 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
6990 
6991 
6992 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6996  const unsigned int dof_no,
6997  const unsigned int first_selected_component,
6998  const unsigned int quad_no_in,
6999  const unsigned int fe_degree,
7000  const unsigned int n_q_points,
7001  const bool is_interior_face,
7002  const unsigned int active_fe_index,
7003  const unsigned int active_quad_index,
7004  const unsigned int face_type)
7006  data_in,
7007  dof_no,
7008  first_selected_component,
7009  quad_no_in,
7010  fe_degree,
7011  n_q_points,
7012  is_interior_face,
7013  active_fe_index,
7014  active_quad_index,
7015  face_type)
7016 {}
7017 
7018 
7019 
7020 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7023  const Mapping<dim> & mapping,
7024  const FiniteElement<dim> &fe,
7025  const Quadrature<1> & quadrature,
7026  const UpdateFlags update_flags,
7027  const unsigned int first_selected_component,
7029  *other)
7031  mapping,
7032  fe,
7033  quadrature,
7034  update_flags,
7035  first_selected_component,
7036  other)
7037 {}
7038 
7039 
7040 
7041 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7045  &other)
7047 {}
7048 
7049 
7050 
7051 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7055  &other)
7056 {
7058  operator=(other);
7059  return *this;
7060 }
7061 
7062 
7063 
7064 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7067  get_gradient(const unsigned int q_point) const
7068 {
7069  return BaseClass::get_gradient(q_point);
7070 }
7071 
7072 
7073 
7074 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7075 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7077  get_divergence(const unsigned int q_point) const
7078 {
7079 # ifdef DEBUG
7080  Assert(this->gradients_quad_initialized == true,
7082 # endif
7083  AssertIndexRange(q_point, this->n_quadrature_points);
7084  Assert(this->jacobian != nullptr,
7086  "update_gradients"));
7087 
7088  VectorizedArrayType divergence;
7089  const std::size_t nqp = this->n_quadrature_points;
7090 
7091  // Cartesian cell
7092  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
7093  {
7094  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
7095  for (unsigned int d = 1; d < dim; ++d)
7096  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
7097  this->jacobian[0][d][d];
7098  }
7099  // cell with general/constant Jacobian
7100  else
7101  {
7103  this->cell_type == internal::MatrixFreeFunctions::general ?
7104  this->jacobian[q_point] :
7105  this->jacobian[0];
7106  divergence = jac[0][0] * this->gradients_quad[q_point];
7107  for (unsigned int e = 1; e < dim; ++e)
7108  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
7109  for (unsigned int d = 1; d < dim; ++d)
7110  for (unsigned int e = 0; e < dim; ++e)
7111  divergence +=
7112  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
7113  }
7114  return divergence;
7115 }
7116 
7117 
7118 
7119 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7122  get_symmetric_gradient(const unsigned int q_point) const
7123 {
7124  // copy from generic function into dim-specialization function
7125  const auto grad = get_gradient(q_point);
7126  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
7127  VectorizedArrayType half = Number(0.5);
7128  for (unsigned int d = 0; d < dim; ++d)
7129  symmetrized[d] = grad[d][d];
7130  switch (dim)
7131  {
7132  case 1:
7133  break;
7134  case 2:
7135  symmetrized[2] = grad[0][1] + grad[1][0];
7136  symmetrized[2] *= half;
7137  break;
7138  case 3:
7139  symmetrized[3] = grad[0][1] + grad[1][0];
7140  symmetrized[3] *= half;
7141  symmetrized[4] = grad[0][2] + grad[2][0];
7142  symmetrized[4] *= half;
7143  symmetrized[5] = grad[1][2] + grad[2][1];
7144  symmetrized[5] *= half;
7145  break;
7146  default:
7147  Assert(false, ExcNotImplemented());
7148  }
7150 }
7151 
7152 
7153 
7154 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7155 inline DEAL_II_ALWAYS_INLINE
7158  const unsigned int q_point) const
7159 {
7160  // copy from generic function into dim-specialization function
7161  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
7163  switch (dim)
7164  {
7165  case 1:
7166  Assert(false,
7167  ExcMessage(
7168  "Computing the curl in 1d is not a useful operation"));
7169  break;
7170  case 2:
7171  curl[0] = grad[1][0] - grad[0][1];
7172  break;
7173  case 3:
7174  curl[0] = grad[2][1] - grad[1][2];
7175  curl[1] = grad[0][2] - grad[2][0];
7176  curl[2] = grad[1][0] - grad[0][1];
7177  break;
7178  default:
7179  Assert(false, ExcNotImplemented());
7180  }
7181  return curl;
7182 }
7183 
7184 
7185 
7186 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7189  get_hessian_diagonal(const unsigned int q_point) const
7190 {
7191  return BaseClass::get_hessian_diagonal(q_point);
7192 }
7193 
7194 
7195 
7196 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7199  const unsigned int q_point) const
7200 {
7201 # ifdef DEBUG
7202  Assert(this->hessians_quad_initialized == true,
7204 # endif
7205  AssertIndexRange(q_point, this->n_quadrature_points);
7206  return BaseClass::get_hessian(q_point);
7207 }
7208 
7209 
7210 
7211 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7212 inline DEAL_II_ALWAYS_INLINE void
7215  const unsigned int q_point)
7216 {
7217  BaseClass::submit_gradient(grad_in, q_point);
7218 }
7219 
7220 
7221 
7222 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7223 inline DEAL_II_ALWAYS_INLINE void
7226  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
7227  const unsigned int q_point)
7228 {
7229  BaseClass::submit_gradient(grad_in, q_point);
7230 }
7231 
7232 
7233 
7234 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7235 inline DEAL_II_ALWAYS_INLINE void
7237  submit_divergence(const VectorizedArrayType div_in,
7238  const unsigned int q_point)
7239 {
7241  AssertIndexRange(q_point, this->n_quadrature_points);
7242  Assert(this->J_value != nullptr,
7244  "update_gradients"));
7245  Assert(this->jacobian != nullptr,
7247  "update_gradients"));
7248 # ifdef DEBUG
7249  this->gradients_quad_submitted = true;
7250 # endif
7251 
7252  const std::size_t nqp = this->n_quadrature_points;
7253  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
7254  {
7255  const VectorizedArrayType fac =
7256  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
7257  for (unsigned int d = 0; d < dim; ++d)
7258  {
7259  this->gradients_quad[(d * dim + d) * nqp + q_point] =
7260  (fac * this->jacobian[0][d][d]);
7261  for (unsigned int e = d + 1; e < dim; ++e)
7262  {
7263  this->gradients_quad[(d * dim + e) * nqp + q_point] =
7264  VectorizedArrayType();
7265  this->gradients_quad[(e * dim + d) * nqp + q_point] =
7266  VectorizedArrayType();
7267  }
7268  }
7269  }
7270  else
7271  {
7273  this->cell_type == internal::MatrixFreeFunctions::general ?
7274  this->jacobian[q_point] :
7275  this->jacobian[0];
7276  const VectorizedArrayType fac =
7277  (this->cell_type == internal::MatrixFreeFunctions::general ?
7278  this->J_value[q_point] :
7279  this->J_value[0] * this->quadrature_weights[q_point]) *
7280  div_in;
7281  for (unsigned int d = 0; d < dim; ++d)
7282  {
7283  for (unsigned int e = 0; e < dim; ++e)
7284  this->gradients_quad[(d * dim + e) * nqp + q_point] =
7285  jac[d][e] * fac;
7286  }
7287  }
7288 }
7289 
7290 
7291 
7292 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7293 inline DEAL_II_ALWAYS_INLINE void
7297  const unsigned int q_point)
7298 {
7299  // could have used base class operator, but that involves some overhead
7300  // which is inefficient. it is nice to have the symmetric tensor because
7301  // that saves some operations
7303  AssertIndexRange(q_point, this->n_quadrature_points);
7304  Assert(this->J_value != nullptr,
7306  "update_gradients"));
7307  Assert(this->jacobian != nullptr,
7309  "update_gradients"));
7310 # ifdef DEBUG
7311  this->gradients_quad_submitted = true;
7312 # endif
7313 
7314  const std::size_t nqp = this->n_quadrature_points;
7315  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
7316  {
7317  const VectorizedArrayType JxW =
7318  this->J_value[0] * this->quadrature_weights[q_point];
7319  for (unsigned int d = 0; d < dim; ++d)
7320  this->gradients_quad[(d * dim + d) * nqp + q_point] =
7321  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
7322  for (unsigned int e = 0, counter = dim; e < dim; ++e)
7323  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
7324  {
7325  const VectorizedArrayType value =
7326  sym_grad.access_raw_entry(counter) * JxW;
7327  this->gradients_quad[(e * dim + d) * nqp + q_point] =
7328  value * this->jacobian[0][d][d];
7329  this->gradients_quad[(d * dim + e) * nqp + q_point] =
7330  value * this->jacobian[0][e][e];
7331  }
7332  }
7333  // general/affine cell type
7334  else
7335  {
7336  const VectorizedArrayType JxW =
7337  this->cell_type == internal::MatrixFreeFunctions::general ?
7338  this->J_value[q_point] :
7339  this->J_value[0] * this->quadrature_weights[q_point];
7341  this->cell_type == internal::MatrixFreeFunctions::general ?
7342  this->jacobian[q_point] :
7343  this->jacobian[0];
7344  VectorizedArrayType weighted[dim][dim];
7345  for (unsigned int i = 0; i < dim; ++i)
7346  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
7347  for (unsigned int i = 0, counter = dim; i < dim; ++i)
7348  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
7349  {
7350  const VectorizedArrayType value =
7351  sym_grad.access_raw_entry(counter) * JxW;
7352  weighted[i][j] = value;
7353  weighted[j][i] = value;
7354  }
7355  for (unsigned int comp = 0; comp < dim; ++comp)
7356  for (unsigned int d = 0; d < dim; ++d)
7357  {
7358  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
7359  for (unsigned int e = 1; e < dim; ++e)
7360  new_val += jac[e][d] * weighted[comp][e];
7361  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
7362  }
7363  }
7364 }
7365 
7366 
7367 
7368 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
7369 inline DEAL_II_ALWAYS_INLINE void
7372  const unsigned int q_point)
7373 {
7375  switch (dim)
7376  {
7377  case 1:
7378  Assert(false,
7379  ExcMessage(
7380  "Testing by the curl in 1d is not a useful operation"));
7381  break;
7382  case 2:
7383  grad[1][0] = curl[0];
7384  grad[0][1] = -curl[0];
7385  break;
7386  case 3:
7387  grad[2][1] = curl[0];
7388  grad[1][2] = -curl[0];
7389  grad[0][2] = curl[1];
7390  grad[2][0] = -curl[1];
7391  grad[1][0] = curl[2];
7392  grad[0][1] = -curl[2];
7393  break;
7394  default:
7395  Assert(false, ExcNotImplemented());
7396  }
7397  submit_gradient(grad, q_point);
7398 }
7399 
7400 
7401 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
7402 
7403 
7404 template <typename Number, bool is_face, typename VectorizedArrayType>
7407  const unsigned int dof_no,
7408  const unsigned int first_selected_component,
7409  const unsigned int quad_no_in,
7410  const unsigned int fe_degree,
7411  const unsigned int n_q_points,
7412  const bool is_interior_face,
7413  const unsigned int active_fe_index,
7414  const unsigned int active_quad_index,
7415  const unsigned int face_type)
7417  data_in,
7418  dof_no,
7419  first_selected_component,
7420  quad_no_in,
7421  fe_degree,
7422  n_q_points,
7423  is_interior_face,
7424  active_fe_index,
7425  active_quad_index,
7426  face_type)
7427 {}
7428 
7429 
7430 
7431 template <typename Number, bool is_face, typename VectorizedArrayType>
7434  const Mapping<1> & mapping,
7435  const FiniteElement<1> &fe,
7436  const Quadrature<1> & quadrature,
7437  const UpdateFlags update_flags,
7438  const unsigned int first_selected_component,
7441  mapping,
7442  fe,
7443  quadrature,
7444  update_flags,
7445  first_selected_component,
7446  other)
7447 {}
7448 
7449 
7450 
7451 template <typename Number, bool is_face, typename VectorizedArrayType>
7456 {}
7457 
7458 
7459 
7460 template <typename Number, bool is_face, typename VectorizedArrayType>
7464 {
7466  other);
7467  return *this;
7468 }
7469 
7470 
7471 
7472 template <typename Number, bool is_face, typename VectorizedArrayType>
7473 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7475  const unsigned int dof) const
7476 {
7477  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7478  return this->values_dofs[0][dof];
7479 }
7480 
7481 
7482 
7483 template <typename Number, bool is_face, typename VectorizedArrayType>
7484 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7486  const unsigned int q_point) const
7487 {
7488 # ifdef DEBUG
7489  Assert(this->values_quad_initialized == true,
7491 # endif
7492  AssertIndexRange(q_point, this->n_quadrature_points);
7493  return this->values_quad[q_point];
7494 }
7495 
7496 
7497 
7498 template <typename Number, bool is_face, typename VectorizedArrayType>
7501  const unsigned int q_point) const
7502 {
7503  // could use the base class gradient, but that involves too many inefficient
7504  // initialization operations on tensors
7505 
7506 # ifdef DEBUG
7507  Assert(this->gradients_quad_initialized == true,
7509 # endif
7510  AssertIndexRange(q_point, this->n_quadrature_points);
7511 
7513  this->cell_type == internal::MatrixFreeFunctions::general ?
7514  this->jacobian[q_point] :
7515  this->jacobian[0];
7516 
7518  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
7519 
7520  return grad_out;
7521 }
7522 
7523 
7524 
7525 template <typename Number, bool is_face, typename VectorizedArrayType>
7526 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7528  const unsigned int q_point) const
7529 {
7530  return get_gradient(q_point)[0];
7531 }
7532 
7533 
7534 
7535 template <typename Number, bool is_face, typename VectorizedArrayType>
7536 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7538  get_normal_derivative(const unsigned int q_point) const
7539 {
7540  return BaseClass::get_normal_derivative(q_point)[0];
7541 }
7542 
7543 
7544 
7545 template <typename Number, bool is_face, typename VectorizedArrayType>
7548  const unsigned int q_point) const
7549 {
7550  return BaseClass::get_hessian(q_point)[0];
7551 }
7552 
7553 
7554 
7555 template <typename Number, bool is_face, typename VectorizedArrayType>
7558  get_hessian_diagonal(const unsigned int q_point) const
7559 {
7560  return BaseClass::get_hessian_diagonal(q_point)[0];
7561 }
7562 
7563 
7564 
7565 template <typename Number, bool is_face, typename VectorizedArrayType>
7566 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7568  const unsigned int q_point) const
7569 {
7570  return BaseClass::get_laplacian(q_point)[0];
7571 }
7572 
7573 
7574 
7575 template <typename Number, bool is_face, typename VectorizedArrayType>
7578  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
7579 {
7580 # ifdef DEBUG
7581  this->dof_values_initialized = true;
7582  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7583 # endif
7584  this->values_dofs[0][dof] = val_in;
7585 }
7586 
7587 
7588 
7589 template <typename Number, bool is_face, typename VectorizedArrayType>
7590 inline DEAL_II_ALWAYS_INLINE void
7592  const VectorizedArrayType val_in,
7593  const unsigned int q_point)
7594 {
7596  AssertIndexRange(q_point, this->n_quadrature_points);
7597 # ifdef DEBUG
7598  this->values_quad_submitted = true;
7599 # endif
7600 
7601  if (this->cell_type == internal::MatrixFreeFunctions::general)
7602  {
7603  const VectorizedArrayType JxW = this->J_value[q_point];
7604  this->values_quad[q_point] = val_in * JxW;
7605  }
7606  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
7607  {
7608  const VectorizedArrayType JxW =
7609  this->J_value[0] * this->quadrature_weights[q_point];
7610  this->values_quad[q_point] = val_in * JxW;
7611  }
7612 }
7613 
7614 
7615 
7616 template <typename Number, bool is_face, typename VectorizedArrayType>
7617 inline DEAL_II_ALWAYS_INLINE void
7619  const Tensor<1, 1, VectorizedArrayType> val_in,
7620  const unsigned int q_point)
7621 {
7622  submit_value(val_in[0], q_point);
7623 }
7624 
7625 
7626 
7627 template <typename Number, bool is_face, typename VectorizedArrayType>
7628 inline DEAL_II_ALWAYS_INLINE void
7630  const Tensor<1, 1, VectorizedArrayType> grad_in,
7631  const unsigned int q_point)
7632 {
7633  submit_gradient(grad_in[0], q_point);
7634 }
7635 
7636 
7637 
7638 template <typename Number, bool is_face, typename VectorizedArrayType>
7639 inline DEAL_II_ALWAYS_INLINE void
7641  const VectorizedArrayType grad_in,
7642  const unsigned int q_point)
7643 {
7645  AssertIndexRange(q_point, this->n_quadrature_points);
7646 # ifdef DEBUG
7647  this->gradients_quad_submitted = true;
7648 # endif
7649 
7651  this->cell_type == internal::MatrixFreeFunctions::general ?
7652  this->jacobian[q_point] :
7653  this->jacobian[0];
7654  const VectorizedArrayType JxW =
7655  this->cell_type == internal::MatrixFreeFunctions::general ?
7656  this->J_value[q_point] :
7657  this->J_value[0] * this->quadrature_weights[q_point];
7658 
7659  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7660 }
7661 
7662 
7663 
7664 template <typename Number, bool is_face, typename VectorizedArrayType>
7665 inline DEAL_II_ALWAYS_INLINE void
7667  const Tensor<2, 1, VectorizedArrayType> grad_in,
7668  const unsigned int q_point)
7669 {
7670  submit_gradient(grad_in[0][0], q_point);
7671 }
7672 
7673 
7674 
7675 template <typename Number, bool is_face, typename VectorizedArrayType>
7676 inline DEAL_II_ALWAYS_INLINE void
7678  submit_normal_derivative(const VectorizedArrayType grad_in,
7679  const unsigned int q_point)
7680 {
7682  grad[0] = grad_in;
7683  BaseClass::submit_normal_derivative(grad, q_point);
7684 }
7685 
7686 
7687 
7688 template <typename Number, bool is_face, typename VectorizedArrayType>
7689 inline DEAL_II_ALWAYS_INLINE void
7692  const unsigned int q_point)
7693 {
7694  BaseClass::submit_normal_derivative(grad_in, q_point);
7695 }
7696 
7697 
7698 template <typename Number, bool is_face, typename VectorizedArrayType>
7699 inline DEAL_II_ALWAYS_INLINE void
7701  const Tensor<2, 1, VectorizedArrayType> hessian_in,
7702  const unsigned int q_point)
7703 {
7705  hessian[0] = hessian_in;
7706  BaseClass::submit_hessian(hessian, q_point);
7707 }
7708 
7709 
7710 template <typename Number, bool is_face, typename VectorizedArrayType>
7711 inline VectorizedArrayType
7713  integrate_value() const
7714 {
7715  return BaseClass::integrate_value()[0];
7716 }
7717 
7718 
7719 
7720 /*-------------------------- FEEvaluation -----------------------------------*/
7721 
7722 
7723 template <int dim,
7724  int fe_degree,
7725  int n_q_points_1d,
7726  int n_components_,
7727  typename Number,
7728  typename VectorizedArrayType>
7729 inline FEEvaluation<dim,
7730  fe_degree,
7731  n_q_points_1d,
7732  n_components_,
7733  Number,
7734  VectorizedArrayType>::
7735  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
7736  const unsigned int fe_no,
7737  const unsigned int quad_no,
7738  const unsigned int first_selected_component,
7739  const unsigned int active_fe_index,
7740  const unsigned int active_quad_index)
7741  : BaseClass(data_in,
7742  fe_no,
7743  first_selected_component,
7744  quad_no,
7745  fe_degree,
7746  static_n_q_points,
7747  true /*note: this is not a face*/,
7748  active_fe_index,
7749  active_quad_index)
7750  , dofs_per_component(this->data->dofs_per_component_on_cell)
7751  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7752  , n_q_points(this->data->n_q_points)
7753 {
7754  check_template_arguments(fe_no, 0);
7755 }
7756 
7757 
7758 
7759 template <int dim,
7760  int fe_degree,
7761  int n_q_points_1d,
7762  int n_components_,
7763  typename Number,
7764  typename VectorizedArrayType>
7765 inline FEEvaluation<dim,
7766  fe_degree,
7767  n_q_points_1d,
7768  n_components_,
7769  Number,
7770  VectorizedArrayType>::
7771  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7772  const std::pair<unsigned int, unsigned int> & range,
7773  const unsigned int dof_no,
7774  const unsigned int quad_no,
7775  const unsigned int first_selected_component)
7776  : FEEvaluation(matrix_free,
7777  dof_no,
7778  quad_no,
7779  first_selected_component,
7780  matrix_free.get_cell_active_fe_index(range))
7781 {}
7782 
7783 
7784 
7785 template <int dim,
7786  int fe_degree,
7787  int n_q_points_1d,
7788  int n_components_,
7789  typename Number,
7790  typename VectorizedArrayType>
7791 inline FEEvaluation<dim,
7792  fe_degree,
7793  n_q_points_1d,
7794  n_components_,
7795  Number,
7796  VectorizedArrayType>::
7797  FEEvaluation(const Mapping<dim> & mapping,
7798  const FiniteElement<dim> &fe,
7799  const Quadrature<1> & quadrature,
7800  const UpdateFlags update_flags,
7801  const unsigned int first_selected_component)
7802  : BaseClass(mapping,
7803  fe,
7804  quadrature,
7805  update_flags,
7806  first_selected_component,
7807  nullptr)
7808  , dofs_per_component(this->data->dofs_per_component_on_cell)
7809  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7810  , n_q_points(this->data->n_q_points)
7811 {
7812  check_template_arguments(numbers::invalid_unsigned_int, 0);
7813 }
7814 
7815 
7816 
7817 template <int dim,
7818  int fe_degree,
7819  int n_q_points_1d,
7820  int n_components_,
7821  typename Number,
7822  typename VectorizedArrayType>
7823 inline FEEvaluation<dim,
7824  fe_degree,
7825  n_q_points_1d,
7826  n_components_,
7827  Number,
7828  VectorizedArrayType>::
7829  FEEvaluation(const FiniteElement<dim> &fe,
7830  const Quadrature<1> & quadrature,
7831  const UpdateFlags update_flags,
7832  const unsigned int first_selected_component)
7833  : BaseClass(StaticMappingQ1<dim>::mapping,
7834  fe,
7835  quadrature,
7836  update_flags,
7837  first_selected_component,
7838  nullptr)
7839  , dofs_per_component(this->data->dofs_per_component_on_cell)
7840  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7841  , n_q_points(this->data->n_q_points)
7842 {
7843  check_template_arguments(numbers::invalid_unsigned_int, 0);
7844 }
7845 
7846 
7847 
7848 template <int dim,
7849  int fe_degree,
7850  int n_q_points_1d,
7851  int n_components_,
7852  typename Number,
7853  typename VectorizedArrayType>
7854 inline FEEvaluation<dim,
7855  fe_degree,
7856  n_q_points_1d,
7857  n_components_,
7858  Number,
7859  VectorizedArrayType>::
7860  FEEvaluation(
7861  const FiniteElement<dim> & fe,
7863  const unsigned int first_selected_component)
7864  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7865  fe,
7866  other.mapped_geometry->get_quadrature(),
7867  other.mapped_geometry->get_fe_values().get_update_flags(),
7868  first_selected_component,
7869  &other)
7870  , dofs_per_component(this->data->dofs_per_component_on_cell)
7871  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7872  , n_q_points(this->data->n_q_points)
7873 {
7874  check_template_arguments(numbers::invalid_unsigned_int, 0);
7875 }
7876 
7877 
7878 
7879 template <int dim,
7880  int fe_degree,
7881  int n_q_points_1d,
7882  int n_components_,
7883  typename Number,
7884  typename VectorizedArrayType>
7885 inline FEEvaluation<dim,
7886  fe_degree,
7887  n_q_points_1d,
7888  n_components_,
7889  Number,
7890  VectorizedArrayType>::FEEvaluation(const FEEvaluation
7891  &other)
7892  : BaseClass(other)
7893  , dofs_per_component(this->data->dofs_per_component_on_cell)
7894  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7895  , n_q_points(this->data->n_q_points)
7896 {
7897  check_template_arguments(numbers::invalid_unsigned_int, 0);
7898 }
7899 
7900 
7901 
7902 template <int dim,
7903  int fe_degree,
7904  int n_q_points_1d,
7905  int n_components_,
7906  typename Number,
7907  typename VectorizedArrayType>
7908 inline FEEvaluation<dim,
7909  fe_degree,
7910  n_q_points_1d,
7911  n_components_,
7912  Number,
7913  VectorizedArrayType> &
7914 FEEvaluation<dim,
7915  fe_degree,
7916  n_q_points_1d,
7917  n_components_,
7918  Number,
7919  VectorizedArrayType>::operator=(const FEEvaluation &other)
7920 {
7921  BaseClass::operator=(other);
7922  check_template_arguments(numbers::invalid_unsigned_int, 0);
7923  return *this;
7924 }
7925 
7926 
7927 
7928 template <int dim,
7929  int fe_degree,
7930  int n_q_points_1d,
7931  int n_components_,
7932  typename Number,
7933  typename VectorizedArrayType>
7934 inline void
7935 FEEvaluation<dim,
7936  fe_degree,
7937  n_q_points_1d,
7938  n_components_,
7939  Number,
7940  VectorizedArrayType>::
7941  check_template_arguments(const unsigned int dof_no,
7942  const unsigned int first_selected_component)
7943 {
7944  (void)dof_no;
7945  (void)first_selected_component;
7946 
7947 # ifdef DEBUG
7948  // print error message when the dimensions do not match. Propose a possible
7949  // fix
7950  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7951  static_cast<unsigned int>(fe_degree) !=
7952  this->data->data.front().fe_degree) ||
7953  n_q_points != this->n_quadrature_points)
7954  {
7955  std::string message =
7956  "-------------------------------------------------------\n";
7957  message += "Illegal arguments in constructor/wrong template arguments!\n";
7958  message += " Called --> FEEvaluation<dim,";
7959  message += Utilities::int_to_string(fe_degree) + ",";
7960  message += Utilities::int_to_string(n_q_points_1d);
7961  message += "," + Utilities::int_to_string(n_components);
7962  message += ",Number>(data";
7963  if (first_selected_component != numbers::invalid_unsigned_int)
7964  {
7965  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7966  message += Utilities::int_to_string(this->quad_no) + ", ";
7967  message += Utilities::int_to_string(first_selected_component);
7968  }
7969  message += ")\n";
7970 
7971  // check whether some other vector component has the correct number of
7972  // points
7973  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7974  proposed_fe_comp = numbers::invalid_unsigned_int,
7975  proposed_quad_comp = numbers::invalid_unsigned_int;
7976  if (dof_no != numbers::invalid_unsigned_int)
7977  {
7978  if (static_cast<unsigned int>(fe_degree) ==
7979  this->data->data.front().fe_degree)
7980  {
7981  proposed_dof_comp = dof_no;
7982  proposed_fe_comp = first_selected_component;
7983  }
7984  else
7985  for (unsigned int no = 0; no < this->matrix_info->n_components();
7986  ++no)
7987  for (unsigned int nf = 0;
7988  nf < this->matrix_info->n_base_elements(no);
7989  ++nf)
7990  if (this->matrix_info
7991  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7992  .data.front()
7993  .fe_degree == static_cast<unsigned int>(fe_degree))
7994  {
7995  proposed_dof_comp = no;
7996  proposed_fe_comp = nf;
7997  break;
7998  }
7999  if (n_q_points ==
8000  this->mapping_data->descriptor[this->active_quad_index]
8001  .n_q_points)
8002  proposed_quad_comp = this->quad_no;
8003  else
8004  for (unsigned int no = 0;
8005  no < this->matrix_info->get_mapping_info().cell_data.size();
8006  ++no)
8007  if (this->matrix_info->get_mapping_info()
8008  .cell_data[no]
8009  .descriptor[this->active_quad_index]
8010  .n_q_points == n_q_points)
8011  {
8012  proposed_quad_comp = no;
8013  break;
8014  }
8015  }
8016  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
8017  proposed_quad_comp != numbers::invalid_unsigned_int)
8018  {
8019  if (proposed_dof_comp != first_selected_component)
8020  message += "Wrong vector component selection:\n";
8021  else
8022  message += "Wrong quadrature formula selection:\n";
8023  message += " Did you mean FEEvaluation<dim,";
8024  message += Utilities::int_to_string(fe_degree) + ",";
8025  message += Utilities::int_to_string(n_q_points_1d);
8026  message += "," + Utilities::int_to_string(n_components);
8027  message += ",Number>(data";
8028  if (dof_no != numbers::invalid_unsigned_int)
8029  {
8030  message +=
8031  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
8032  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
8033  message += Utilities::int_to_string(proposed_fe_comp);
8034  }
8035  message += ")?\n";
8036  std::string correct_pos;
8037  if (proposed_dof_comp != dof_no)
8038  correct_pos = " ^ ";
8039  else
8040  correct_pos = " ";
8041  if (proposed_quad_comp != this->quad_no)
8042  correct_pos += " ^ ";
8043  else
8044  correct_pos += " ";
8045  if (proposed_fe_comp != first_selected_component)
8046  correct_pos += " ^\n";
8047  else
8048  correct_pos += " \n";
8049  message += " " +
8050  correct_pos;
8051  }
8052  // ok, did not find the numbers specified by the template arguments in
8053  // the given list. Suggest correct template arguments
8054  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
8055  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
8056  message += "Wrong template arguments:\n";
8057  message += " Did you mean FEEvaluation<dim,";
8058  message +=
8059  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
8060  message += Utilities::int_to_string(proposed_n_q_points_1d);
8061  message += "," + Utilities::int_to_string(n_components);
8062  message += ",Number>(data";
8063  if (dof_no != numbers::invalid_unsigned_int)
8064  {
8065  message += ", " + Utilities::int_to_string(dof_no) + ", ";
8066  message += Utilities::int_to_string(this->quad_no);
8067  message += ", " + Utilities::int_to_string(first_selected_component);
8068  }
8069  message += ")?\n";
8070  std::string correct_pos;
8071  if (this->data->data.front().fe_degree !=
8072  static_cast<unsigned int>(fe_degree))
8073  correct_pos = " ^";
8074  else
8075  correct_pos = " ";
8076  if (proposed_n_q_points_1d != n_q_points_1d)
8077  correct_pos += " ^\n";
8078  else
8079  correct_pos += " \n";
8080  message += " " + correct_pos;
8081 
8082  Assert(static_cast<unsigned int>(fe_degree) ==
8083  this->data->data.front().fe_degree &&
8084  n_q_points == this->n_quadrature_points,
8085  ExcMessage(message));
8086  }
8087  if (dof_no != numbers::invalid_unsigned_int)
8089  n_q_points,
8090  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
8091 # endif
8092 }
8093 
8094 
8095 
8096 template <int dim,
8097  int fe_degree,
8098  int n_q_points_1d,
8099  int n_components_,
8100  typename Number,
8101  typename VectorizedArrayType>
8102 inline void
8103 FEEvaluation<dim,
8104  fe_degree,
8105  n_q_points_1d,
8106  n_components_,
8107  Number,
8108  VectorizedArrayType>::reinit(const unsigned int cell_index)
8109 {
8110  Assert(this->mapped_geometry == nullptr,
8111  ExcMessage("FEEvaluation was initialized without a matrix-free object."
8112  " Integer indexing is not possible"));
8113  if (this->mapped_geometry != nullptr)
8114  return;
8115 
8116  Assert(this->dof_info != nullptr, ExcNotInitialized());
8117  Assert(this->mapping_data != nullptr, ExcNotInitialized());
8118  this->cell = cell_index;
8119  this->cell_type =
8120  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
8121 
8122  const unsigned int offsets =
8123  this->mapping_data->data_index_offsets[cell_index];
8124  this->jacobian = &this->mapping_data->jacobians[0][offsets];
8125  this->J_value = &this->mapping_data->JxW_values[offsets];
8126 
8127 # ifdef DEBUG
8128  this->dof_values_initialized = false;
8129  this->values_quad_initialized = false;
8130  this->gradients_quad_initialized = false;
8131  this->hessians_quad_initialized = false;
8132 # endif
8133 }
8134 
8135 
8136 
8137 template <int dim,
8138  int fe_degree,
8139  int n_q_points_1d,
8140  int n_components_,
8141  typename Number,
8142  typename VectorizedArrayType>
8143 template <bool level_dof_access>
8144 inline void
8145 FEEvaluation<dim,
8146  fe_degree,
8147  n_q_points_1d,
8148  n_components_,
8149  Number,
8150  VectorizedArrayType>::
8152 {
8153  Assert(this->matrix_info == nullptr,
8154  ExcMessage("Cannot use initialization from cell iterator if "
8155  "initialized from MatrixFree object. Use variant for "
8156  "on the fly computation with arguments as for FEValues "
8157  "instead"));
8158  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
8159  this->mapped_geometry->reinit(
8160  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
8161  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
8162  if (level_dof_access)
8163  cell->get_mg_dof_indices(this->local_dof_indices);
8164  else
8165  cell->get_dof_indices(this->local_dof_indices);
8166 }
8167 
8168 
8169 
8170 template <int dim,
8171  int fe_degree,
8172  int n_q_points_1d,
8173  int n_components_,
8174  typename Number,
8175  typename VectorizedArrayType>
8176 inline void
8177 FEEvaluation<dim,
8178  fe_degree,
8179  n_q_points_1d,
8180  n_components_,
8181  Number,
8182  VectorizedArrayType>::
8183  reinit(const typename Triangulation<dim>::cell_iterator &cell)
8184 {
8185  Assert(this->matrix_info == 0,
8186  ExcMessage("Cannot use initialization from cell iterator if "
8187  "initialized from MatrixFree object. Use variant for "
8188  "on the fly computation with arguments as for FEValues "
8189  "instead"));
8190  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
8191  this->mapped_geometry->reinit(cell);
8192 }
8193 
8194 
8195 
8196 template <int dim,
8197  int fe_degree,
8198  int n_q_points_1d,
8199  int n_components_,
8200  typename Number,
8201  typename VectorizedArrayType>
8203 FEEvaluation<dim,
8204  fe_degree,
8205  n_q_points_1d,
8206  n_components_,
8207  Number,
8208  VectorizedArrayType>::quadrature_point(const unsigned int q) const
8209 {
8210  if (this->matrix_info == nullptr)
8211  {
8212  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
8215  "update_quadrature_points"));
8216  }
8217  else
8218  {
8219  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
8221  "update_quadrature_points"));
8222  }
8223 
8224  AssertIndexRange(q, n_q_points);
8225 
8227  &this->mapping_data->quadrature_points
8228  [this->mapping_data->quadrature_point_offsets[this->cell]];
8229 
8230  // Cartesian/affine mesh: only first vertex of cell is stored, we must
8231  // compute it through the Jacobian (which is stored in non-inverted and
8232  // non-transposed form as index '1' in the jacobian field)
8233  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
8234  {
8235  Assert(this->jacobian != nullptr, ExcNotInitialized());
8236  Point<dim, VectorizedArrayType> point = quadrature_points[0];
8237 
8238  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
8239  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
8240  for (unsigned int d = 0; d < dim; ++d)
8241  point[d] += jac[d][d] * static_cast<Number>(
8242  this->descriptor->quadrature.point(q)[d]);
8243  else
8244  for (unsigned int d = 0; d < dim; ++d)
8245  for (unsigned int e = 0; e < dim; ++e)
8246  point[d] += jac[d][e] * static_cast<Number>(
8247  this->descriptor->quadrature.point(q)[e]);
8248  return point;
8249  }
8250  else
8251  return quadrature_points[q];
8252 }
8253 
8254 
8255 
8256 template <int dim,
8257  int fe_degree,
8258  int n_q_points_1d,
8259  int n_components_,
8260  typename Number,
8261  typename VectorizedArrayType>
8262 inline void
8263 FEEvaluation<dim,
8264  fe_degree,
8265  n_q_points_1d,
8266  n_components_,
8267  Number,
8268  VectorizedArrayType>::evaluate(const bool evaluate_values,
8269  const bool evaluate_gradients,
8270  const bool evaluate_hessians)
8271 {
8272 # ifdef DEBUG
8273  Assert(this->dof_values_initialized == true,
8275 # endif
8276  evaluate(this->values_dofs[0],
8277  evaluate_values,
8278  evaluate_gradients,
8279  evaluate_hessians);
8280 }
8281 
8282 
8283 template <int dim,
8284  int fe_degree,
8285  int n_q_points_1d,
8286  int n_components_,
8287  typename Number,
8288  typename VectorizedArrayType>
8289 inline void
8290 FEEvaluation<dim,
8291  fe_degree,
8292  n_q_points_1d,
8293  n_components_,
8294  Number,
8295  VectorizedArrayType>::
8296  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
8297 {
8298 # ifdef DEBUG
8299  Assert(this->dof_values_initialized == true,
8301 # endif
8302  evaluate(this->values_dofs[0], evaluation_flags);
8303 }
8304 
8305 
8306 
8307 template <int dim,
8308  int fe_degree,
8309  int n_q_points_1d,
8310  int n_components_,
8311  typename Number,
8312  typename VectorizedArrayType>
8313 inline void
8314 FEEvaluation<dim,
8315  fe_degree,
8316  n_q_points_1d,
8317  n_components_,
8318  Number,
8319  VectorizedArrayType>::evaluate(const VectorizedArrayType
8320  * values_array,
8321  const bool evaluate_values,
8322  const bool evaluate_gradients,
8323  const bool evaluate_hessians)
8324 {
8326  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8327  ((evaluate_gradients) ? EvaluationFlags::gradients :
8328  EvaluationFlags::nothing) |
8329  ((evaluate_hessians) ? EvaluationFlags::hessians :
8330  EvaluationFlags::nothing);
8331 
8332  evaluate(values_array, flag);
8333 }
8334 
8335 
8336 
8337 template <int dim,
8338  int fe_degree,
8339  int n_q_points_1d,
8340  int n_components_,
8341  typename Number,
8342  typename VectorizedArrayType>
8343 inline void
8344 FEEvaluation<dim,
8345  fe_degree,
8346  n_q_points_1d,
8347  n_components_,
8348  Number,
8349  VectorizedArrayType>::
8350  evaluate(const VectorizedArrayType * values_array,
8351  const EvaluationFlags::EvaluationFlags evaluation_flag)
8352 {
8353  const bool hessians_on_general_cells =
8354  evaluation_flag & EvaluationFlags::hessians &&
8355  (this->cell_type > internal::MatrixFreeFunctions::affine);
8356  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
8357  if (hessians_on_general_cells)
8358  evaluation_flag_actual |= EvaluationFlags::gradients;
8359 
8360  if (fe_degree > -1)
8361