Reference documentation for deal.II version Git ac854ef673 2021-03-06 10:09:45 +0100
\(\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 - 2020 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  static constexpr unsigned int dimension = dim;
553  static constexpr unsigned int n_components = n_components_;
554 
591  template <typename VectorType>
592  void
593  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
594 
623  template <typename VectorType>
624  void
625  read_dof_values_plain(const VectorType & src,
626  const unsigned int first_index = 0);
627 
659  template <typename VectorType>
660  void
661  distribute_local_to_global(
662  VectorType & dst,
663  const unsigned int first_index = 0,
664  const std::bitset<VectorizedArrayType::size()> &mask =
665  std::bitset<VectorizedArrayType::size()>().flip()) const;
666 
705  template <typename VectorType>
706  void
707  set_dof_values(VectorType & dst,
708  const unsigned int first_index = 0,
709  const std::bitset<VectorizedArrayType::size()> &mask =
710  std::bitset<VectorizedArrayType::size()>().flip()) const;
711 
715  template <typename VectorType>
716  void
717  set_dof_values_plain(
718  VectorType & dst,
719  const unsigned int first_index = 0,
720  const std::bitset<VectorizedArrayType::size()> &mask =
721  std::bitset<VectorizedArrayType::size()>().flip()) const;
722 
724 
745  value_type
746  get_dof_value(const unsigned int dof) const;
747 
758  void
759  submit_dof_value(const value_type val_in, const unsigned int dof);
760 
773  value_type
774  get_value(const unsigned int q_point) const;
775 
788  void
789  submit_value(const value_type val_in, const unsigned int q_point);
790 
802  get_gradient(const unsigned int q_point) const;
803 
818  value_type
819  get_normal_derivative(const unsigned int q_point) const;
820 
833  void
834  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
835 
854  void
855  submit_normal_derivative(const value_type grad_in,
856  const unsigned int q_point);
857 
870  get_hessian(const unsigned int q_point) const;
871 
882  get_hessian_diagonal(const unsigned int q_point) const;
883 
895  value_type
896  get_laplacian(const unsigned int q_point) const;
897 
898 #ifdef DOXYGEN
899  // doxygen does not anyhow mention functions coming from partial template
900  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
901  // For now, hack in those functions manually only to fix documentation:
902 
909  VectorizedArrayType
910  get_divergence(const unsigned int q_point) const;
911 
921  get_symmetric_gradient(const unsigned int q_point) const;
922 
930  get_curl(const unsigned int q_point) const;
931 
947  void
948  submit_divergence(const VectorizedArrayType div_in,
949  const unsigned int q_point);
950 
967  void
968  submit_symmetric_gradient(
970  const unsigned int q_point);
971 
984  void
985  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
986  const unsigned int q_point);
987 
988 #endif
989 
1006  value_type
1007  integrate_value() const;
1008 
1010 
1023  const VectorizedArrayType *
1024  begin_dof_values() const;
1025 
1034  VectorizedArrayType *
1035  begin_dof_values();
1036 
1047  const VectorizedArrayType *
1048  begin_values() const;
1049 
1060  VectorizedArrayType *
1061  begin_values();
1062 
1074  const VectorizedArrayType *
1075  begin_gradients() const;
1076 
1088  VectorizedArrayType *
1089  begin_gradients();
1090 
1103  const VectorizedArrayType *
1104  begin_hessians() const;
1105 
1118  VectorizedArrayType *
1119  begin_hessians();
1120 
1122 
1126  unsigned int
1127  get_first_selected_component() const;
1128 
1129 protected:
1140  const unsigned int dof_no,
1141  const unsigned int first_selected_component,
1142  const unsigned int quad_no,
1143  const unsigned int fe_degree,
1144  const unsigned int n_q_points,
1145  const bool is_interior_face,
1146  const unsigned int active_fe_index,
1147  const unsigned int active_quad_index,
1148  const unsigned int face_type);
1149 
1186  const Mapping<dim> & mapping,
1187  const FiniteElement<dim> &fe,
1188  const Quadrature<1> & quadrature,
1189  const UpdateFlags update_flags,
1190  const unsigned int first_selected_component,
1192  *other);
1193 
1200  FEEvaluationBase(const FEEvaluationBase &other);
1201 
1209  operator=(const FEEvaluationBase &other);
1210 
1217  template <typename VectorType, typename VectorOperation>
1218  void
1219  read_write_operation(
1220  const VectorOperation & operation,
1221  const std::array<VectorType *, n_components_> &vectors,
1222  const std::array<
1224  n_components_> & vectors_sm,
1225  const std::bitset<VectorizedArrayType::size()> &mask,
1226  const bool apply_constraints = true) const;
1227 
1235  template <typename VectorType, typename VectorOperation>
1236  void
1237  read_write_operation_contiguous(
1238  const VectorOperation & operation,
1239  const std::array<VectorType *, n_components_> &vectors,
1240  const std::array<
1242  n_components_> & vectors_sm,
1243  const std::bitset<VectorizedArrayType::size()> &mask) const;
1244 
1252  template <typename VectorType, typename VectorOperation>
1253  void
1254  read_write_operation_global(
1255  const VectorOperation & operation,
1256  const std::array<VectorType *, n_components_> &vectors) const;
1257 
1270  VectorizedArrayType *values_dofs[n_components];
1271 
1283  VectorizedArrayType *values_quad;
1284 
1298  VectorizedArrayType *gradients_quad;
1299 
1311  VectorizedArrayType *hessians_quad;
1312 
1317  const unsigned int n_fe_components;
1318 
1325 
1332 
1339 
1346 
1353 
1360 
1365  const unsigned int first_selected_component;
1366 
1371  mutable std::vector<types::global_dof_index> local_dof_indices;
1372 
1373 private:
1378  void
1379  set_data_pointers();
1380 };
1381 
1382 
1383 
1391 template <int dim,
1392  int n_components_,
1393  typename Number,
1394  bool is_face,
1395  typename VectorizedArrayType = VectorizedArray<Number>>
1397  n_components_,
1398  Number,
1399  is_face,
1400  VectorizedArrayType>
1401 {
1402  static_assert(
1403  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1404  "Type of Number and of VectorizedArrayType do not match.");
1405 
1406 public:
1407  using number_type = Number;
1409  using gradient_type =
1411  static constexpr unsigned int dimension = dim;
1412  static constexpr unsigned int n_components = n_components_;
1413  using BaseClass =
1415 
1416 protected:
1426  const unsigned int dof_no,
1427  const unsigned int first_selected_component,
1428  const unsigned int quad_no,
1429  const unsigned int fe_degree,
1430  const unsigned int n_q_points,
1431  const bool is_interior_face = true,
1432  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1433  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1434  const unsigned int face_type = numbers::invalid_unsigned_int);
1435 
1441  const Mapping<dim> & mapping,
1442  const FiniteElement<dim> &fe,
1443  const Quadrature<1> & quadrature,
1444  const UpdateFlags update_flags,
1445  const unsigned int first_selected_component,
1447  *other);
1448 
1452  FEEvaluationAccess(const FEEvaluationAccess &other);
1453 
1458  operator=(const FEEvaluationAccess &other);
1459 };
1460 
1461 
1462 
1471 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1472 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1473  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1474 {
1475  static_assert(
1476  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1477  "Type of Number and of VectorizedArrayType do not match.");
1478 
1479 public:
1480  using number_type = Number;
1481  using value_type = VectorizedArrayType;
1483  static constexpr unsigned int dimension = dim;
1484  using BaseClass =
1486 
1489  value_type
1490  get_dof_value(const unsigned int dof) const;
1491 
1494  void
1495  submit_dof_value(const value_type val_in, const unsigned int dof);
1496 
1499  value_type
1500  get_value(const unsigned int q_point) const;
1501 
1504  void
1505  submit_value(const value_type val_in, const unsigned int q_point);
1506 
1509  void
1510  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
1511  const unsigned int q_point);
1512 
1516  get_gradient(const unsigned int q_point) const;
1517 
1520  value_type
1521  get_normal_derivative(const unsigned int q_point) const;
1522 
1525  void
1526  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1527 
1530  void
1531  submit_normal_derivative(const value_type grad_in,
1532  const unsigned int q_point);
1533 
1537  get_hessian(unsigned int q_point) const;
1538 
1542  get_hessian_diagonal(const unsigned int q_point) const;
1543 
1546  value_type
1547  get_laplacian(const unsigned int q_point) const;
1548 
1551  value_type
1552  integrate_value() const;
1553 
1554 protected:
1564  const unsigned int dof_no,
1565  const unsigned int first_selected_component,
1566  const unsigned int quad_no,
1567  const unsigned int fe_degree,
1568  const unsigned int n_q_points,
1569  const bool is_interior_face = true,
1570  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1571  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1572  const unsigned int face_type = numbers::invalid_unsigned_int);
1573 
1579  const Mapping<dim> & mapping,
1580  const FiniteElement<dim> &fe,
1581  const Quadrature<1> & quadrature,
1582  const UpdateFlags update_flags,
1583  const unsigned int first_selected_component,
1585  *other);
1586 
1590  FEEvaluationAccess(const FEEvaluationAccess &other);
1591 
1596  operator=(const FEEvaluationAccess &other);
1597 };
1598 
1599 
1600 
1610 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1611 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1612  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1613 {
1614  static_assert(
1615  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1616  "Type of Number and of VectorizedArrayType do not match.");
1617 
1618 public:
1619  using number_type = Number;
1622  static constexpr unsigned int dimension = dim;
1623  static constexpr unsigned int n_components = dim;
1624  using BaseClass =
1626 
1630  get_gradient(const unsigned int q_point) const;
1631 
1636  VectorizedArrayType
1637  get_divergence(const unsigned int q_point) const;
1638 
1646  get_symmetric_gradient(const unsigned int q_point) const;
1647 
1653  get_curl(const unsigned int q_point) const;
1654 
1658  get_hessian(const unsigned int q_point) const;
1659 
1663  get_hessian_diagonal(const unsigned int q_point) const;
1664 
1667  void
1668  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1669 
1678  void
1679  submit_gradient(
1680  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1681  const unsigned int q_point);
1682 
1691  void
1692  submit_divergence(const VectorizedArrayType div_in,
1693  const unsigned int q_point);
1694 
1703  void
1704  submit_symmetric_gradient(
1706  const unsigned int q_point);
1707 
1712  void
1713  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1714  const unsigned int q_point);
1715 
1716 protected:
1726  const unsigned int dof_no,
1727  const unsigned int first_selected_component,
1728  const unsigned int quad_no,
1729  const unsigned int dofs_per_cell,
1730  const unsigned int n_q_points,
1731  const bool is_interior_face = true,
1732  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1733  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1734  const unsigned int face_type = numbers::invalid_unsigned_int);
1735 
1741  const Mapping<dim> & mapping,
1742  const FiniteElement<dim> &fe,
1743  const Quadrature<1> & quadrature,
1744  const UpdateFlags update_flags,
1745  const unsigned int first_selected_component,
1747  *other);
1748 
1752  FEEvaluationAccess(const FEEvaluationAccess &other);
1753 
1758  operator=(const FEEvaluationAccess &other);
1759 };
1760 
1761 
1770 template <typename Number, bool is_face, typename VectorizedArrayType>
1771 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1772  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1773 {
1774  static_assert(
1775  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1776  "Type of Number and of VectorizedArrayType do not match.");
1777 
1778 public:
1779  using number_type = Number;
1780  using value_type = VectorizedArrayType;
1782  static constexpr unsigned int dimension = 1;
1783  using BaseClass =
1785 
1788  value_type
1789  get_dof_value(const unsigned int dof) const;
1790 
1793  void
1794  submit_dof_value(const value_type val_in, const unsigned int dof);
1795 
1798  value_type
1799  get_value(const unsigned int q_point) const;
1800 
1803  void
1804  submit_value(const value_type val_in, const unsigned int q_point);
1805 
1808  void
1809  submit_value(const gradient_type val_in, const unsigned int q_point);
1810 
1814  get_gradient(const unsigned int q_point) const;
1815 
1818  value_type
1819  get_normal_derivative(const unsigned int q_point) const;
1820 
1823  void
1824  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1825 
1828  void
1829  submit_gradient(const value_type grad_in, const unsigned int q_point);
1830 
1833  void
1834  submit_normal_derivative(const value_type grad_in,
1835  const unsigned int q_point);
1836 
1839  void
1840  submit_normal_derivative(const gradient_type grad_in,
1841  const unsigned int q_point);
1842 
1846  get_hessian(unsigned int q_point) const;
1847 
1851  get_hessian_diagonal(const unsigned int q_point) const;
1852 
1855  value_type
1856  get_laplacian(const unsigned int q_point) const;
1857 
1860  value_type
1861  integrate_value() const;
1862 
1863 protected:
1872  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
1873  const unsigned int dof_no,
1874  const unsigned int first_selected_component,
1875  const unsigned int quad_no,
1876  const unsigned int fe_degree,
1877  const unsigned int n_q_points,
1878  const bool is_interior_face = true,
1879  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1880  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1881  const unsigned int face_type = numbers::invalid_unsigned_int);
1882 
1888  const Mapping<1> & mapping,
1889  const FiniteElement<1> &fe,
1890  const Quadrature<1> & quadrature,
1891  const UpdateFlags update_flags,
1892  const unsigned int first_selected_component,
1894 
1898  FEEvaluationAccess(const FEEvaluationAccess &other);
1899 
1904  operator=(const FEEvaluationAccess &other);
1905 };
1906 
1907 
1908 
2460 template <int dim,
2461  int fe_degree,
2462  int n_q_points_1d,
2463  int n_components_,
2464  typename Number,
2465  typename VectorizedArrayType>
2466 class FEEvaluation : public FEEvaluationAccess<dim,
2467  n_components_,
2468  Number,
2469  false,
2470  VectorizedArrayType>
2471 {
2472  static_assert(
2473  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2474  "Type of Number and of VectorizedArrayType do not match.");
2475 
2476 public:
2480  using BaseClass =
2482 
2486  using number_type = Number;
2487 
2494 
2501 
2505  static constexpr unsigned int dimension = dim;
2506 
2511  static constexpr unsigned int n_components = n_components_;
2512 
2519  static constexpr unsigned int static_n_q_points =
2520  Utilities::pow(n_q_points_1d, dim);
2521 
2529  static constexpr unsigned int static_dofs_per_component =
2530  Utilities::pow(fe_degree + 1, dim);
2531 
2539  static constexpr unsigned int tensor_dofs_per_cell =
2540  static_dofs_per_component * n_components;
2541 
2549  static constexpr unsigned int static_dofs_per_cell =
2550  static_dofs_per_component * n_components;
2551 
2586  FEEvaluation(
2588  const unsigned int dof_no = 0,
2589  const unsigned int quad_no = 0,
2590  const unsigned int first_selected_component = 0,
2591  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2592  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2593 
2602  const std::pair<unsigned int, unsigned int> & range,
2603  const unsigned int dof_no = 0,
2604  const unsigned int quad_no = 0,
2605  const unsigned int first_selected_component = 0);
2606 
2633  FEEvaluation(const Mapping<dim> & mapping,
2634  const FiniteElement<dim> &fe,
2635  const Quadrature<1> & quadrature,
2636  const UpdateFlags update_flags,
2637  const unsigned int first_selected_component = 0);
2638 
2644  FEEvaluation(const FiniteElement<dim> &fe,
2645  const Quadrature<1> & quadrature,
2646  const UpdateFlags update_flags,
2647  const unsigned int first_selected_component = 0);
2648 
2659  FEEvaluation(
2660  const FiniteElement<dim> & fe,
2662  const unsigned int first_selected_component = 0);
2663 
2670  FEEvaluation(const FEEvaluation &other);
2671 
2678  FEEvaluation &
2679  operator=(const FEEvaluation &other);
2680 
2689  void
2690  reinit(const unsigned int cell_batch_index);
2691 
2704  template <bool level_dof_access>
2705  void
2707 
2718  void
2719  reinit(const typename Triangulation<dim>::cell_iterator &cell);
2720 
2730  void
2731  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2732 
2737  void
2738  evaluate(const bool evaluate_values,
2739  const bool evaluate_gradients,
2740  const bool evaluate_hessians = false);
2741 
2754  void
2755  evaluate(const VectorizedArrayType * values_array,
2756  const EvaluationFlags::EvaluationFlags evaluation_flag);
2757 
2762  void
2763  evaluate(const VectorizedArrayType *values_array,
2764  const bool evaluate_values,
2765  const bool evaluate_gradients,
2766  const bool evaluate_hessians = false);
2767 
2781  template <typename VectorType>
2782  void
2783  gather_evaluate(const VectorType & input_vector,
2784  const EvaluationFlags::EvaluationFlags evaluation_flag);
2785 
2789  template <typename VectorType>
2790  void
2791  gather_evaluate(const VectorType &input_vector,
2792  const bool evaluate_values,
2793  const bool evaluate_gradients,
2794  const bool evaluate_hessians = false);
2795 
2806  void
2807  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2808 
2809 
2813  void
2814  integrate(const bool integrate_values, const bool integrate_gradients);
2815 
2827  void
2828  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2829  VectorizedArrayType * values_array);
2830 
2834  void
2835  integrate(const bool integrate_values,
2836  const bool integrate_gradients,
2837  VectorizedArrayType *values_array);
2838 
2852  template <typename VectorType>
2853  void
2854  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
2855  VectorType & output_vector);
2856 
2860  template <typename VectorType>
2861  void
2862  integrate_scatter(const bool integrate_values,
2863  const bool integrate_gradients,
2864  VectorType &output_vector);
2865 
2871  quadrature_point(const unsigned int q_point) const;
2872 
2879  const unsigned int dofs_per_component;
2880 
2887  const unsigned int dofs_per_cell;
2888 
2896  const unsigned int n_q_points;
2897 
2898 private:
2903  void
2904  check_template_arguments(const unsigned int fe_no,
2905  const unsigned int first_selected_component);
2906 };
2907 
2908 
2909 
2945 template <int dim,
2946  int fe_degree,
2947  int n_q_points_1d = fe_degree + 1,
2948  int n_components_ = 1,
2949  typename Number = double,
2950  typename VectorizedArrayType = VectorizedArray<Number>>
2952  n_components_,
2953  Number,
2954  true,
2955  VectorizedArrayType>
2956 {
2957  static_assert(
2958  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2959  "Type of Number and of VectorizedArrayType do not match.");
2960 
2961 public:
2965  using BaseClass =
2967 
2971  using number_type = Number;
2972 
2979 
2986 
2990  static constexpr unsigned int dimension = dim;
2991 
2996  static constexpr unsigned int n_components = n_components_;
2997 
3005  static constexpr unsigned int static_n_q_points =
3006  Utilities::pow(n_q_points_1d, dim - 1);
3007 
3014  static constexpr unsigned int static_n_q_points_cell =
3015  Utilities::pow(n_q_points_1d, dim);
3016 
3023  static constexpr unsigned int static_dofs_per_component =
3024  Utilities::pow(fe_degree + 1, dim);
3025 
3032  static constexpr unsigned int tensor_dofs_per_cell =
3033  static_dofs_per_component * n_components;
3034 
3041  static constexpr unsigned int static_dofs_per_cell =
3042  static_dofs_per_component * n_components;
3043 
3087  const bool is_interior_face = true,
3088  const unsigned int dof_no = 0,
3089  const unsigned int quad_no = 0,
3090  const unsigned int first_selected_component = 0,
3091  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
3092  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
3093  const unsigned int face_type = numbers::invalid_unsigned_int);
3094 
3104  const std::pair<unsigned int, unsigned int> & range,
3105  const bool is_interior_face = true,
3106  const unsigned int dof_no = 0,
3107  const unsigned int quad_no = 0,
3108  const unsigned int first_selected_component = 0);
3109 
3120  void
3121  reinit(const unsigned int face_batch_number);
3122 
3130  void
3131  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
3132 
3143  void
3144  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3145 
3149  void
3150  evaluate(const bool evaluate_values, const bool evaluate_gradients);
3151 
3164  void
3165  evaluate(const VectorizedArrayType * values_array,
3166  const EvaluationFlags::EvaluationFlags evaluation_flag);
3167 
3171  void
3172  evaluate(const VectorizedArrayType *values_array,
3173  const bool evaluate_values,
3174  const bool evaluate_gradients);
3175 
3187  template <typename VectorType>
3188  void
3189  gather_evaluate(const VectorType & input_vector,
3190  const EvaluationFlags::EvaluationFlags evaluation_flag);
3191 
3195  template <typename VectorType>
3196  void
3197  gather_evaluate(const VectorType &input_vector,
3198  const bool evaluate_values,
3199  const bool evaluate_gradients);
3200 
3210  void
3211  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3212 
3216  void
3217  integrate(const bool integrate_values, const bool integrate_gradients);
3218 
3227  void
3228  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
3229  VectorizedArrayType * values_array);
3230 
3234  void
3235  integrate(const bool integrate_values,
3236  const bool integrate_gradients,
3237  VectorizedArrayType *values_array);
3238 
3250  template <typename VectorType>
3251  void
3252  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
3253  VectorType & output_vector);
3254 
3258  template <typename VectorType>
3259  void
3260  integrate_scatter(const bool integrate_values,
3261  const bool integrate_gradients,
3262  VectorType &output_vector);
3263 
3269  quadrature_point(const unsigned int q_point) const;
3270 
3277  const unsigned int dofs_per_component;
3278 
3285  const unsigned int dofs_per_cell;
3286 
3294  const unsigned int n_q_points;
3295 
3296 
3297 private:
3301  std::array<unsigned int, VectorizedArrayType::size()>
3302  compute_face_no_data();
3303 
3307  std::array<unsigned int, VectorizedArrayType::size()>
3308  compute_face_orientations();
3309 };
3310 
3311 
3312 
3313 namespace internal
3314 {
3315  namespace MatrixFreeFunctions
3316  {
3317  // a helper function to compute the number of DoFs of a DGP element at
3318  // compile time, depending on the degree
3319  template <int dim, int degree>
3321  {
3322  // this division is always without remainder
3323  static constexpr unsigned int value =
3324  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
3325  };
3326 
3327  // base specialization: 1d elements have 'degree+1' degrees of freedom
3328  template <int degree>
3329  struct DGP_dofs_per_component<1, degree>
3330  {
3331  static constexpr unsigned int value = degree + 1;
3332  };
3333  } // namespace MatrixFreeFunctions
3334 } // namespace internal
3335 
3336 
3337 /*----------------------- Inline functions ----------------------------------*/
3338 
3339 #ifndef DOXYGEN
3340 
3341 
3342 /*----------------------- FEEvaluationBaseData ------------------------*/
3343 
3344 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3348  const unsigned int dof_no,
3349  const unsigned int first_selected_component,
3350  const unsigned int quad_no_in,
3351  const unsigned int fe_degree,
3352  const unsigned int n_q_points,
3353  const bool is_interior_face,
3354  const unsigned int active_fe_index_in,
3355  const unsigned int active_quad_index_in,
3356  const unsigned int face_type)
3357  : scratch_data_array(data_in.acquire_scratch_data())
3358  , quad_no(quad_no_in)
3359  , matrix_info(&data_in)
3360  , dof_info(&data_in.get_dof_info(dof_no))
3361  , mapping_data(
3362  internal::MatrixFreeFunctions::
3363  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3364  data_in.get_mapping_info(),
3365  quad_no))
3366  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
3367  data_in.get_dof_info(dof_no).fe_index_from_degree(
3368  first_selected_component,
3369  fe_degree) :
3370  (active_fe_index_in != numbers::invalid_unsigned_int ?
3371  active_fe_index_in :
3372  0))
3373  , active_quad_index(
3374  fe_degree != numbers::invalid_unsigned_int ?
3375  (mapping_data->quad_index_from_n_q_points(n_q_points)) :
3376  (active_quad_index_in != numbers::invalid_unsigned_int ?
3377  active_quad_index_in :
3378  std::min<unsigned int>(active_fe_index,
3379  mapping_data->descriptor.size() - 1)))
3380  , descriptor(
3381  &mapping_data->descriptor
3382  [is_face ?
3383  (active_quad_index * std::max<unsigned int>(1, dim - 1) +
3384  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
3385  active_quad_index])
3386  , n_quadrature_points(descriptor->n_q_points)
3387  , data(&data_in.get_shape_info(
3388  dof_no,
3389  quad_no_in,
3390  dof_info->component_to_base_index[first_selected_component],
3391  active_fe_index,
3392  active_quad_index))
3393  , jacobian(nullptr)
3394  , J_value(nullptr)
3395  , normal_vectors(nullptr)
3396  , normal_x_jacobian(nullptr)
3397  , quadrature_weights(descriptor->quadrature_weights.begin())
3398  , cell(numbers::invalid_unsigned_int)
3399  , is_interior_face(is_interior_face)
3400  , dof_access_index(
3401  is_face ?
3402  (is_interior_face ?
3403  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3404  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3405  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3406  , cell_type(internal::MatrixFreeFunctions::general)
3407 {
3408  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3409  AssertDimension(matrix_info->get_task_info().vectorization_length,
3410  VectorizedArrayType::size());
3411  AssertDimension(n_quadrature_points, descriptor->n_q_points);
3412 }
3413 
3414 
3415 
3416 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3419  const Mapping<dim> & mapping,
3420  const FiniteElement<dim> &fe,
3421  const Quadrature<1> & quadrature,
3422  const UpdateFlags update_flags,
3423  const unsigned int first_selected_component,
3425  *other)
3426  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3428  , active_fe_index(numbers::invalid_unsigned_int)
3429  , active_quad_index(numbers::invalid_unsigned_int)
3430  , descriptor(nullptr)
3431  , n_quadrature_points(
3432  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3433  , matrix_info(nullptr)
3434  , dof_info(nullptr)
3435  , mapping_data(nullptr)
3436  ,
3437  // select the correct base element from the given FE component
3439  Quadrature<dim - is_face>(quadrature),
3440  fe,
3441  fe.component_to_base_index(first_selected_component).first))
3442  , jacobian(nullptr)
3443  , J_value(nullptr)
3444  , normal_vectors(nullptr)
3445  , normal_x_jacobian(nullptr)
3446  , quadrature_weights(nullptr)
3447  , cell(0)
3449  , is_interior_face(true)
3451 {
3452  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3453  ExcInternalError());
3454  if (other != nullptr &&
3455  other->mapped_geometry->get_quadrature() == quadrature)
3456  mapped_geometry = other->mapped_geometry;
3457  else
3458  mapped_geometry =
3459  std::make_shared<internal::MatrixFreeFunctions::
3461  mapping, quadrature, update_flags);
3462  cell = 0;
3463 
3464  mapping_data = &mapped_geometry->get_data_storage();
3465  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3466  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3467 }
3468 
3469 
3470 
3471 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3475  &other)
3476  : scratch_data_array(other.matrix_info == nullptr ?
3478  other.matrix_info->acquire_scratch_data())
3479  , quad_no(other.quad_no)
3480  , active_fe_index(other.active_fe_index)
3481  , active_quad_index(other.active_quad_index)
3482  , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
3483  , n_quadrature_points(other.n_quadrature_points)
3484  , matrix_info(other.matrix_info)
3485  , dof_info(other.dof_info)
3486  , mapping_data(other.mapping_data)
3487  , data(other.matrix_info == nullptr ?
3489  *other.data) :
3490  other.data)
3491  , jacobian(nullptr)
3492  , J_value(nullptr)
3493  , normal_vectors(nullptr)
3494  , normal_x_jacobian(nullptr)
3495  , quadrature_weights(other.descriptor == nullptr ?
3496  nullptr :
3497  descriptor->quadrature_weights.begin())
3500  , is_interior_face(other.is_interior_face)
3501  , dof_access_index(other.dof_access_index)
3502 {
3503  // Create deep copy of mapped geometry for use in parallel...
3504  if (other.mapped_geometry.get() != nullptr)
3505  {
3506  mapped_geometry = std::make_shared<
3507  internal::MatrixFreeFunctions::
3508  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3509  other.mapped_geometry->get_fe_values().get_mapping(),
3510  other.mapped_geometry->get_quadrature(),
3511  other.mapped_geometry->get_fe_values().get_update_flags());
3512  mapping_data = &mapped_geometry->get_data_storage();
3513  cell = 0;
3514 
3515  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3516  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3517  }
3518 }
3519 
3520 
3521 
3522 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3526 {
3527  AssertDimension(quad_no, other.quad_no);
3528  AssertDimension(active_fe_index, other.active_fe_index);
3529  AssertDimension(active_quad_index, other.active_quad_index);
3530 
3531  // release old memory
3532  if (matrix_info == nullptr)
3533  {
3534  delete data;
3535  delete scratch_data_array;
3536  }
3537  else
3538  {
3539  matrix_info->release_scratch_data(scratch_data_array);
3540  }
3541 
3542  matrix_info = other.matrix_info;
3543  dof_info = other.dof_info;
3544  descriptor = other.descriptor;
3545  mapping_data = other.mapping_data;
3546  if (other.matrix_info == nullptr)
3547  {
3549  *other.data);
3550  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3551  }
3552  else
3553  {
3554  data = other.data;
3555  scratch_data_array = matrix_info->acquire_scratch_data();
3556  }
3557 
3558  quadrature_weights =
3559  (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr);
3562  is_interior_face = other.is_interior_face;
3563  dof_access_index = other.dof_access_index;
3564 
3565  // Create deep copy of mapped geometry for use in parallel...
3566  if (other.mapped_geometry.get() != nullptr)
3567  {
3568  mapped_geometry = std::make_shared<
3569  internal::MatrixFreeFunctions::
3570  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3571  other.mapped_geometry->get_fe_values().get_mapping(),
3572  other.mapped_geometry->get_quadrature(),
3573  other.mapped_geometry->get_fe_values().get_update_flags());
3574  cell = 0;
3575  mapping_data = &mapped_geometry->get_data_storage();
3576  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3577  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3578  }
3579 
3580  return *this;
3581 }
3582 
3583 
3584 
3585 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3588 {
3589  if (matrix_info != nullptr)
3590  {
3591  try
3592  {
3593  matrix_info->release_scratch_data(scratch_data_array);
3594  }
3595  catch (...)
3596  {}
3597  }
3598  else
3599  {
3600  delete scratch_data_array;
3601  delete data;
3602  data = nullptr;
3603  }
3604  scratch_data_array = nullptr;
3605 }
3606 
3607 
3608 
3609 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3610 inline unsigned int
3613 {
3614  if (matrix_info == nullptr)
3615  return 0;
3616  else
3617  {
3618  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3619  return this->mapping_data->data_index_offsets[cell];
3620  }
3621 }
3622 
3623 
3624 
3625 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3628  const
3629 {
3631  return cell_type;
3632 }
3633 
3634 
3635 
3636 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3639  get_shape_info() const
3640 {
3641  Assert(data != nullptr, ExcInternalError());
3642  return *data;
3643 }
3644 
3645 
3646 
3647 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3650  const
3651 {
3652  Assert(dof_info != nullptr,
3653  ExcMessage(
3654  "FEEvaluation was not initialized with a MatrixFree object!"));
3655  return *dof_info;
3656 }
3657 
3658 
3659 
3660 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3663  get_normal_vector(const unsigned int q_point) const
3664 {
3665  AssertIndexRange(q_point, n_quadrature_points);
3666  Assert(normal_vectors != nullptr,
3668  "update_normal_vectors"));
3669  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3670  return normal_vectors[0];
3671  else
3672  return normal_vectors[q_point];
3673 }
3674 
3675 
3676 
3677 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3678 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3680  const unsigned int q_point) const
3681 {
3682  AssertIndexRange(q_point, n_quadrature_points);
3683  Assert(J_value != nullptr,
3685  "update_values|update_gradients"));
3686  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3687  {
3688  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3689  return J_value[0] * this->quadrature_weights[q_point];
3690  }
3691  else
3692  return J_value[q_point];
3693 }
3694 
3695 
3696 
3697 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3700  inverse_jacobian(const unsigned int q_point) const
3701 {
3702  AssertIndexRange(q_point, n_quadrature_points);
3703  Assert(this->jacobian != nullptr,
3705  "update_gradients"));
3706  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3707  return jacobian[0];
3708  else
3709  return jacobian[q_point];
3710 }
3711 
3712 
3713 
3714 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3715 inline std::array<unsigned int, VectorizedArrayType::size()>
3717  const
3718 {
3719  Assert(this->matrix_info != nullptr, ExcNotInitialized());
3720 
3721  const unsigned int n_lanes = VectorizedArrayType::size();
3722  std::array<unsigned int, n_lanes> cells;
3723 
3724  // initialize array
3725  for (unsigned int i = 0; i < n_lanes; ++i)
3726  cells[i] = numbers::invalid_unsigned_int;
3727 
3728  if ((is_face == false) ||
3729  (is_face &&
3730  this->dof_access_index ==
3732  this->is_interior_face))
3733  {
3734  // cell or interior face face (element-centric loop)
3735  for (unsigned int i = 0; i < n_lanes; ++i)
3736  cells[i] = cell * n_lanes + i;
3737  }
3738  else if (is_face &&
3739  this->dof_access_index ==
3741  this->is_interior_face == false)
3742  {
3743  // exterior face (element-centric loop): for this case, we need to
3744  // look into the FaceInfo field that collects information from both
3745  // sides of a face once for the global mesh, and pick the face id that
3746  // is not the local one (cell_this).
3747  for (unsigned int i = 0; i < n_lanes; i++)
3748  {
3749  // compute actual (non vectorized) cell ID
3750  const unsigned int cell_this = this->cell * n_lanes + i;
3751  // compute face ID
3752  unsigned int face_index =
3753  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3754  this->face_no,
3755  i);
3756 
3757  if (face_index == numbers::invalid_unsigned_int)
3758  continue; // invalid face ID: no neighbor on boundary
3759 
3760  // get cell ID on both sides of face
3761  auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3762  .cells_interior[face_index % n_lanes];
3763  auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3764  .cells_exterior[face_index % n_lanes];
3765 
3766  // compare the IDs with the given cell ID
3767  if (cell_m == cell_this)
3768  cells[i] = cell_p; // neighbor has the other ID
3769  else if (cell_p == cell_this)
3770  cells[i] = cell_m;
3771  }
3772  }
3773  else if (is_face)
3774  {
3775  // face-centric faces
3776  const unsigned int *cells_ =
3777  is_interior_face ?
3778  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3779  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3780  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3781  if (cells_[i] != numbers::invalid_unsigned_int)
3782  cells[i] = cells_[i];
3783  }
3784 
3785  return cells;
3786 }
3787 
3788 
3789 namespace internal
3790 {
3791  template <int dim,
3792  typename Number,
3793  bool is_face,
3794  typename VectorizedArrayType,
3795  typename VectorizedArrayType2,
3796  typename GlobalVectorType,
3797  typename FU>
3798  inline void
3799  process_cell_data(
3802  GlobalVectorType & array,
3803  VectorizedArrayType2 & out,
3804  const FU & fu)
3805  {
3806  (void)matrix_info;
3807  Assert(matrix_info != nullptr, ExcNotImplemented());
3808  AssertDimension(array.size(),
3809  matrix_info->get_task_info().cell_partition_data.back());
3810 
3811  // 1) collect ids of cell
3812  const auto cells = phi.get_cell_ids();
3813 
3814  // 2) actually gather values
3815  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3816  if (cells[i] != numbers::invalid_unsigned_int)
3817  fu(out[i],
3818  array[cells[i] / VectorizedArrayType::size()]
3819  [cells[i] % VectorizedArrayType::size()]);
3820  }
3821 } // namespace internal
3822 
3823 
3824 
3825 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3826 std::array<unsigned int, VectorizedArrayType::size()>
3828  get_cell_or_face_ids() const
3829 {
3830  const unsigned int v_len = VectorizedArrayType::size();
3831  std::array<unsigned int, VectorizedArrayType::size()> cells;
3832 
3833  // initialize array
3834  for (unsigned int i = 0; i < v_len; ++i)
3835  cells[i] = numbers::invalid_unsigned_int;
3836 
3837  if (is_face &&
3838  this->dof_access_index ==
3840  this->is_interior_face == false)
3841  {
3842  // cell-based face-loop: plus face
3843  for (unsigned int i = 0; i < v_len; i++)
3844  {
3845  // compute actual (non vectorized) cell ID
3846  const unsigned int cell_this = this->cell * v_len + i;
3847  // compute face ID
3848  unsigned int fn =
3849  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3850  this->face_no,
3851  i);
3852 
3854  continue; // invalid face ID: no neighbor on boundary
3855 
3856  // get cell ID on both sides of face
3857  auto cell_m = this->matrix_info->get_face_info(fn / v_len)
3858  .cells_interior[fn % v_len];
3859  auto cell_p = this->matrix_info->get_face_info(fn / v_len)
3860  .cells_exterior[fn % v_len];
3861 
3862  // compare the IDs with the given cell ID
3863  if (cell_m == cell_this)
3864  cells[i] = cell_p; // neighbor has the other ID
3865  else if (cell_p == cell_this)
3866  cells[i] = cell_m;
3867  }
3868  }
3869  else
3870  {
3871  for (unsigned int i = 0; i < v_len; ++i)
3872  cells[i] = cell * v_len + i;
3873  }
3874 
3875  return cells;
3876 }
3877 
3878 
3879 
3880 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3881 inline VectorizedArrayType
3883  const AlignedVector<VectorizedArrayType> &array) const
3884 {
3885  VectorizedArrayType out = Number(1.);
3886  internal::process_cell_data(
3887  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3888  local = global;
3889  });
3890  return out;
3891 }
3892 
3893 
3894 
3895 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3896 inline void
3899  const VectorizedArrayType & in) const
3900 {
3901  internal::process_cell_data(
3902  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3903  global = local;
3904  });
3905 }
3906 
3907 
3908 
3909 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3910 template <typename T>
3911 inline std::array<T, VectorizedArrayType::size()>
3913  const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
3914 {
3915  std::array<T, VectorizedArrayType::size()> out;
3916  internal::process_cell_data(
3917  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3918  local = global;
3919  });
3920  return out;
3921 }
3922 
3923 
3924 
3925 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3926 template <typename T>
3927 inline void
3929  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
3930  const std::array<T, VectorizedArrayType::size()> & in) const
3931 {
3932  internal::process_cell_data(
3933  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3934  global = local;
3935  });
3936 }
3937 
3938 
3939 
3940 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3941 inline const std::vector<unsigned int> &
3944 {
3945  return data->lexicographic_numbering;
3946 }
3947 
3948 
3949 
3950 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3953  get_scratch_data() const
3954 {
3956  const_cast<VectorizedArrayType *>(scratch_data),
3957  scratch_data_array->end() - scratch_data);
3958 }
3959 
3960 
3961 
3962 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3963 inline unsigned int
3965  get_quadrature_index() const
3966 {
3967  return this->quad_no;
3968 }
3969 
3970 
3971 
3972 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3973 inline unsigned int
3976 {
3977  if (is_face && this->dof_access_index ==
3979  return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
3980  else
3981  return this->cell;
3982 }
3983 
3984 
3985 
3986 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3987 inline unsigned int
3989  get_active_fe_index() const
3990 {
3991  return active_fe_index;
3992 }
3993 
3994 
3995 
3996 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3997 inline unsigned int
4000 {
4001  return active_quad_index;
4002 }
4003 
4004 
4005 
4006 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
4009  get_matrix_free() const
4010 {
4011  Assert(matrix_info != nullptr,
4012  ExcMessage(
4013  "FEEvaluation was not initialized with a MatrixFree object!"));
4014  return *matrix_info;
4015 }
4016 
4017 
4018 /*----------------------- FEEvaluationBase ----------------------------------*/
4019 
4020 template <int dim,
4021  int n_components_,
4022  typename Number,
4023  bool is_face,
4024  typename VectorizedArrayType>
4025 inline FEEvaluationBase<dim,
4026  n_components_,
4027  Number,
4028  is_face,
4029  VectorizedArrayType>::
4030  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
4031  const unsigned int dof_no,
4032  const unsigned int first_selected_component,
4033  const unsigned int quad_no_in,
4034  const unsigned int fe_degree,
4035  const unsigned int n_q_points,
4036  const bool is_interior_face,
4037  const unsigned int active_fe_index,
4038  const unsigned int active_quad_index,
4039  const unsigned int face_type)
4041  data_in,
4042  dof_no,
4043  first_selected_component,
4044  quad_no_in,
4045  fe_degree,
4046  n_q_points,
4047  is_interior_face,
4048  active_fe_index,
4049  active_quad_index,
4050  face_type)
4051  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
4052  , dof_values_initialized(false)
4053  , values_quad_initialized(false)
4054  , gradients_quad_initialized(false)
4055  , hessians_quad_initialized(false)
4056  , values_quad_submitted(false)
4057  , gradients_quad_submitted(false)
4058  , first_selected_component(first_selected_component)
4059 {
4060  set_data_pointers();
4061  Assert(
4062  this->dof_info->start_components.back() == 1 ||
4063  static_cast<int>(n_components_) <=
4064  static_cast<int>(
4065  this->dof_info->start_components
4066  [this->dof_info->component_to_base_index[first_selected_component] +
4067  1]) -
4068  first_selected_component,
4069  ExcMessage(
4070  "You tried to construct a vector-valued evaluator with " +
4071  std::to_string(n_components) +
4072  " components. However, "
4073  "the current base element has only " +
4075  this->dof_info->start_components
4076  [this->dof_info->component_to_base_index[first_selected_component] +
4077  1] -
4078  first_selected_component) +
4079  " components left when starting from local element index " +
4081  first_selected_component -
4082  this->dof_info->start_components
4083  [this->dof_info->component_to_base_index[first_selected_component]]) +
4084  " (global index " + std::to_string(first_selected_component) + ")"));
4085 
4086  // do not check for correct dimensions of data fields here, should be done
4087  // in derived classes
4088 }
4089 
4090 
4091 
4092 template <int dim,
4093  int n_components_,
4094  typename Number,
4095  bool is_face,
4096  typename VectorizedArrayType>
4097 inline FEEvaluationBase<dim,
4098  n_components_,
4099  Number,
4100  is_face,
4101  VectorizedArrayType>::
4102  FEEvaluationBase(
4103  const Mapping<dim> & mapping,
4104  const FiniteElement<dim> &fe,
4105  const Quadrature<1> & quadrature,
4106  const UpdateFlags update_flags,
4107  const unsigned int first_selected_component,
4109  *other)
4111  mapping,
4112  fe,
4113  quadrature,
4114  update_flags,
4115  first_selected_component,
4116  other)
4117  , n_fe_components(n_components_)
4118  , dof_values_initialized(false)
4119  , values_quad_initialized(false)
4120  , gradients_quad_initialized(false)
4121  , hessians_quad_initialized(false)
4122  , values_quad_submitted(false)
4123  , gradients_quad_submitted(false)
4124  // keep the number of the selected component within the current base element
4125  // for reading dof values
4126  , first_selected_component(first_selected_component)
4127 {
4128  set_data_pointers();
4129 
4130  const unsigned int base_element_number =
4131  fe.component_to_base_index(first_selected_component).first;
4132  Assert(fe.element_multiplicity(base_element_number) == 1 ||
4133  fe.element_multiplicity(base_element_number) -
4134  first_selected_component >=
4135  n_components_,
4136  ExcMessage("The underlying element must at least contain as many "
4137  "components as requested by this class"));
4138  (void)base_element_number;
4139 }
4140 
4141 
4142 
4143 template <int dim,
4144  int n_components_,
4145  typename Number,
4146  bool is_face,
4147  typename VectorizedArrayType>
4148 inline FEEvaluationBase<dim,
4149  n_components_,
4150  Number,
4151  is_face,
4152  VectorizedArrayType>::
4153  FEEvaluationBase(const FEEvaluationBase<dim,
4154  n_components_,
4155  Number,
4156  is_face,
4157  VectorizedArrayType> &other)
4159  , n_fe_components(other.n_fe_components)
4160  , dof_values_initialized(false)
4161  , values_quad_initialized(false)
4162  , gradients_quad_initialized(false)
4163  , hessians_quad_initialized(false)
4164  , values_quad_submitted(false)
4165  , gradients_quad_submitted(false)
4166  , first_selected_component(other.first_selected_component)
4167 {
4168  set_data_pointers();
4169 }
4170 
4171 
4172 
4173 template <int dim,
4174  int n_components_,
4175  typename Number,
4176  bool is_face,
4177  typename VectorizedArrayType>
4178 inline FEEvaluationBase<dim,
4179  n_components_,
4180  Number,
4181  is_face,
4182  VectorizedArrayType> &
4184 operator=(const FEEvaluationBase<dim,
4185  n_components_,
4186  Number,
4187  is_face,
4188  VectorizedArrayType> &other)
4189 {
4191  operator=(other);
4192  AssertDimension(n_fe_components, other.n_fe_components);
4193  AssertDimension(first_selected_component, other.first_selected_component);
4194 
4195  return *this;
4196 }
4197 
4198 
4199 
4200 template <int dim,
4201  int n_components_,
4202  typename Number,
4203  bool is_face,
4204  typename VectorizedArrayType>
4205 inline void
4208 {
4209  Assert(this->scratch_data_array != nullptr, ExcInternalError());
4210 
4211  const unsigned int tensor_dofs_per_component =
4212  Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4213  const unsigned int dofs_per_component =
4214  this->data->dofs_per_component_on_cell;
4215  const unsigned int n_quadrature_points = this->n_quadrature_points;
4216 
4217  const unsigned int shift =
4218  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4219  n_components_ * 3 +
4220  2 * n_quadrature_points;
4221  const unsigned int allocated_size =
4222  shift + n_components_ * dofs_per_component +
4223  (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4224  this->scratch_data_array->resize_fast(allocated_size);
4225 
4226  // set the pointers to the correct position in the data array
4227  for (unsigned int c = 0; c < n_components_; ++c)
4228  {
4229  values_dofs[c] =
4230  this->scratch_data_array->begin() + c * dofs_per_component;
4231  }
4232  values_quad =
4233  this->scratch_data_array->begin() + n_components * dofs_per_component;
4234  gradients_quad = this->scratch_data_array->begin() +
4235  n_components * (dofs_per_component + n_quadrature_points);
4236  hessians_quad =
4237  this->scratch_data_array->begin() +
4238  n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4239  this->scratch_data =
4240  this->scratch_data_array->begin() + n_components_ * dofs_per_component +
4241  (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4242 }
4243 
4244 
4245 
4246 namespace internal
4247 {
4248  // allows to select between block vectors and non-block vectors, which
4249  // allows to use a unified interface for extracting blocks on block vectors
4250  // and doing nothing on usual vectors
4251  template <typename VectorType, bool>
4252  struct BlockVectorSelector
4253  {};
4254 
4255  template <typename VectorType>
4256  struct BlockVectorSelector<VectorType, true>
4257  {
4258  using BaseVectorType = typename VectorType::BlockType;
4259 
4260  static BaseVectorType *
4261  get_vector_component(VectorType &vec, const unsigned int component)
4262  {
4263  AssertIndexRange(component, vec.n_blocks());
4264  return &vec.block(component);
4265  }
4266  };
4267 
4268  template <typename VectorType>
4269  struct BlockVectorSelector<VectorType, false>
4270  {
4271  using BaseVectorType = VectorType;
4272 
4273  static BaseVectorType *
4274  get_vector_component(VectorType &vec, const unsigned int component)
4275  {
4276  // FEEvaluation allows to combine several vectors from a scalar
4277  // FiniteElement into a "vector-valued" FEEvaluation object with
4278  // multiple components. These components can be extracted with the other
4279  // get_vector_component functions. If we do not get a vector of vectors
4280  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
4281  // must make sure that we do not duplicate the components in input
4282  // and/or duplicate the resulting integrals. In such a case, we should
4283  // only get the zeroth component in the vector contained set nullptr for
4284  // the others which allows us to catch unintended use in
4285  // read_write_operation.
4286  if (component == 0)
4287  return &vec;
4288  else
4289  return nullptr;
4290  }
4291  };
4292 
4293  template <typename VectorType>
4294  struct BlockVectorSelector<std::vector<VectorType>, false>
4295  {
4296  using BaseVectorType = VectorType;
4297 
4298  static BaseVectorType *
4299  get_vector_component(std::vector<VectorType> &vec,
4300  const unsigned int component)
4301  {
4302  AssertIndexRange(component, vec.size());
4303  return &vec[component];
4304  }
4305  };
4306 
4307  template <typename VectorType>
4308  struct BlockVectorSelector<std::vector<VectorType *>, false>
4309  {
4310  using BaseVectorType = VectorType;
4311 
4312  static BaseVectorType *
4313  get_vector_component(std::vector<VectorType *> &vec,
4314  const unsigned int component)
4315  {
4316  AssertIndexRange(component, vec.size());
4317  return vec[component];
4318  }
4319  };
4320 } // namespace internal
4321 
4322 
4323 
4324 template <int dim,
4325  int n_components_,
4326  typename Number,
4327  bool is_face,
4328  typename VectorizedArrayType>
4329 template <typename VectorType, typename VectorOperation>
4330 inline void
4333  const VectorOperation & operation,
4334  const std::array<VectorType *, n_components_> &src,
4335  const std::array<
4337  n_components_> & src_sm,
4338  const std::bitset<VectorizedArrayType::size()> &mask,
4339  const bool apply_constraints) const
4340 {
4341  // Case 1: No MatrixFree object given, simple case because we do not need to
4342  // process constraints and need not care about vectorization -> go to
4343  // separate function
4344  if (this->matrix_info == nullptr)
4345  {
4346  read_write_operation_global(operation, src);
4347  return;
4348  }
4349 
4350  Assert(this->dof_info != nullptr, ExcNotInitialized());
4351  Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
4352  if (n_fe_components == 1)
4353  for (unsigned int comp = 0; comp < n_components; ++comp)
4354  {
4355  Assert(src[comp] != nullptr,
4356  ExcMessage("The finite element underlying this FEEvaluation "
4357  "object is scalar, but you requested " +
4358  std::to_string(n_components) +
4359  " components via the template argument in "
4360  "FEEvaluation. In that case, you must pass an "
4361  "std::vector<VectorType> or a BlockVector to " +
4362  "read_dof_values and distribute_local_to_global."));
4363  internal::check_vector_compatibility(*src[comp], *this->dof_info);
4364  }
4365  else
4366  {
4367  internal::check_vector_compatibility(*src[0], *this->dof_info);
4368  }
4369 
4370  // Case 2: contiguous indices which use reduced storage of indices and can
4371  // use vectorized load/store operations -> go to separate function
4373  this->cell,
4374  this->dof_info->index_storage_variants[this->dof_access_index].size());
4375  if (this->dof_info->index_storage_variants
4376  [is_face ? this->dof_access_index :
4378  [this->cell] >=
4380  {
4381  read_write_operation_contiguous(operation, src, src_sm, mask);
4382  return;
4383  }
4384 
4385  // Case 3: standard operation with one index per degree of freedom -> go on
4386  // here
4387  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4388  Assert(mask.count() == n_lanes,
4389  ExcNotImplemented("Masking currently not implemented for "
4390  "non-contiguous DoF storage"));
4391 
4392  std::integral_constant<bool,
4394  vector_selector;
4395 
4396  const unsigned int dofs_per_component =
4397  this->data->dofs_per_component_on_cell;
4398  if (this->dof_info->index_storage_variants
4399  [is_face ? this->dof_access_index :
4401  [this->cell] ==
4403  {
4404  const unsigned int *dof_indices =
4405  this->dof_info->dof_indices_interleaved.data() +
4406  this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4407  .first +
4408  this->dof_info
4409  ->component_dof_indices_offset[this->active_fe_index]
4410  [this->first_selected_component] *
4411  n_lanes;
4412  if (n_components == 1 || n_fe_components == 1)
4413  for (unsigned int i = 0; i < dofs_per_component;
4414  ++i, dof_indices += n_lanes)
4415  for (unsigned int comp = 0; comp < n_components; ++comp)
4416  operation.process_dof_gather(dof_indices,
4417  *src[comp],
4418  0,
4419  values_dofs[comp][i],
4420  vector_selector);
4421  else
4422  for (unsigned int comp = 0; comp < n_components; ++comp)
4423  for (unsigned int i = 0; i < dofs_per_component;
4424  ++i, dof_indices += n_lanes)
4425  operation.process_dof_gather(
4426  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4427  return;
4428  }
4429 
4430  const unsigned int * dof_indices[n_lanes];
4431  VectorizedArrayType **values_dofs =
4432  const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
4433 
4434  // Assign the appropriate cell ids for face/cell case and get the pointers
4435  // to the dof indices of the cells on all lanes
4436  unsigned int cells_copied[n_lanes];
4437  const unsigned int *cells;
4438  unsigned int n_vectorization_actual =
4439  this->dof_info
4440  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4441  bool has_constraints = false;
4442  const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4443  if (is_face)
4444  {
4445  if (this->dof_access_index ==
4447  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4448  cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4449  cells =
4450  this->dof_access_index ==
4452  &cells_copied[0] :
4453  (this->is_interior_face ?
4454  &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
4455  &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
4456  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4457  {
4458  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4459  ExcInternalError());
4460  const std::pair<unsigned int, unsigned int> *my_index_start =
4461  &this->dof_info->row_starts[cells[v] * n_fe_components +
4462  first_selected_component];
4463 
4464  // check whether any of the SIMD lanes has constraints, i.e., the
4465  // constraint indicator which is the second entry of row_starts
4466  // increments on this cell
4467  if (my_index_start[n_components_read].second !=
4468  my_index_start[0].second)
4469  has_constraints = true;
4470 
4471  dof_indices[v] =
4472  this->dof_info->dof_indices.data() + my_index_start[0].first;
4473  }
4474  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4475  dof_indices[v] = nullptr;
4476  }
4477  else
4478  {
4479  AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
4480  this->dof_info->row_starts.size());
4481  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4482  {
4483  const std::pair<unsigned int, unsigned int> *my_index_start =
4484  &this->dof_info
4485  ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4486  first_selected_component];
4487  if (my_index_start[n_components_read].second !=
4488  my_index_start[0].second)
4489  has_constraints = true;
4490  Assert(my_index_start[n_components_read].first ==
4491  my_index_start[0].first ||
4492  my_index_start[0].first < this->dof_info->dof_indices.size(),
4493  ExcIndexRange(0,
4494  my_index_start[0].first,
4495  this->dof_info->dof_indices.size()));
4496  dof_indices[v] =
4497  this->dof_info->dof_indices.data() + my_index_start[0].first;
4498  }
4499  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4500  dof_indices[v] = nullptr;
4501  }
4502 
4503  // Case where we have no constraints throughout the whole cell: Can go
4504  // through the list of DoFs directly
4505  if (!has_constraints)
4506  {
4507  if (n_vectorization_actual < n_lanes)
4508  for (unsigned int comp = 0; comp < n_components; ++comp)
4509  for (unsigned int i = 0; i < dofs_per_component; ++i)
4510  operation.process_empty(values_dofs[comp][i]);
4511  if (n_components == 1 || n_fe_components == 1)
4512  {
4513  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4514  for (unsigned int i = 0; i < dofs_per_component; ++i)
4515  for (unsigned int comp = 0; comp < n_components; ++comp)
4516  operation.process_dof(dof_indices[v][i],
4517  *src[comp],
4518  values_dofs[comp][i][v]);
4519  }
4520  else
4521  {
4522  for (unsigned int comp = 0; comp < n_components; ++comp)
4523  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4524  for (unsigned int i = 0; i < dofs_per_component; ++i)
4525  operation.process_dof(
4526  dof_indices[v][comp * dofs_per_component + i],
4527  *src[0],
4528  values_dofs[comp][i][v]);
4529  }
4530  return;
4531  }
4532 
4533  // In the case where there are some constraints to be resolved, loop over
4534  // all vector components that are filled and then over local dofs. ind_local
4535  // holds local number on cell, index iterates over the elements of
4536  // index_local_to_global and dof_indices points to the global indices stored
4537  // in index_local_to_global
4538  if (n_vectorization_actual < n_lanes)
4539  for (unsigned int comp = 0; comp < n_components; ++comp)
4540  for (unsigned int i = 0; i < dofs_per_component; ++i)
4541  operation.process_empty(values_dofs[comp][i]);
4542  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4543  {
4544  const unsigned int cell_index =
4545  is_face ? cells[v] : this->cell * n_lanes + v;
4546  const unsigned int cell_dof_index =
4547  cell_index * n_fe_components + first_selected_component;
4548  const unsigned int n_components_read =
4549  n_fe_components > 1 ? n_components : 1;
4550  unsigned int index_indicators =
4551  this->dof_info->row_starts[cell_dof_index].second;
4552  unsigned int next_index_indicators =
4553  this->dof_info->row_starts[cell_dof_index + 1].second;
4554 
4555  // For read_dof_values_plain, redirect the dof_indices field to the
4556  // unconstrained indices
4557  if (apply_constraints == false &&
4558  this->dof_info->row_starts[cell_dof_index].second !=
4559  this->dof_info->row_starts[cell_dof_index + n_components_read]
4560  .second)
4561  {
4562  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
4564  ExcNotInitialized());
4565  dof_indices[v] =
4566  this->dof_info->plain_dof_indices.data() +
4567  this->dof_info
4568  ->component_dof_indices_offset[this->active_fe_index]
4569  [this->first_selected_component] +
4570  this->dof_info->row_starts_plain_indices[cell_index];
4571  next_index_indicators = index_indicators;
4572  }
4573 
4574  if (n_components == 1 || n_fe_components == 1)
4575  {
4576  unsigned int ind_local = 0;
4577  for (; index_indicators != next_index_indicators; ++index_indicators)
4578  {
4579  const std::pair<unsigned short, unsigned short> indicator =
4580  this->dof_info->constraint_indicator[index_indicators];
4581  // run through values up to next constraint
4582  for (unsigned int j = 0; j < indicator.first; ++j)
4583  for (unsigned int comp = 0; comp < n_components; ++comp)
4584  operation.process_dof(dof_indices[v][j],
4585  *src[comp],
4586  values_dofs[comp][ind_local + j][v]);
4587 
4588  ind_local += indicator.first;
4589  dof_indices[v] += indicator.first;
4590 
4591  // constrained case: build the local value as a linear
4592  // combination of the global value according to constraints
4593  Number value[n_components];
4594  for (unsigned int comp = 0; comp < n_components; ++comp)
4595  operation.pre_constraints(values_dofs[comp][ind_local][v],
4596  value[comp]);
4597 
4598  const Number *data_val =
4599  this->matrix_info->constraint_pool_begin(indicator.second);
4600  const Number *end_pool =
4601  this->matrix_info->constraint_pool_end(indicator.second);
4602  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4603  for (unsigned int comp = 0; comp < n_components; ++comp)
4604  operation.process_constraint(*dof_indices[v],
4605  *data_val,
4606  *src[comp],
4607  value[comp]);
4608 
4609  for (unsigned int comp = 0; comp < n_components; ++comp)
4610  operation.post_constraints(value[comp],
4611  values_dofs[comp][ind_local][v]);
4612  ind_local++;
4613  }
4614 
4615  AssertIndexRange(ind_local, dofs_per_component + 1);
4616 
4617  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4618  for (unsigned int comp = 0; comp < n_components; ++comp)
4619  operation.process_dof(*dof_indices[v],
4620  *src[comp],
4621  values_dofs[comp][ind_local][v]);
4622  }
4623  else
4624  {
4625  // case with vector-valued finite elements where all components are
4626  // included in one single vector. Assumption: first come all entries
4627  // to the first component, then all entries to the second one, and
4628  // so on. This is ensured by the way MatrixFree reads out the
4629  // indices.
4630  for (unsigned int comp = 0; comp < n_components; ++comp)
4631  {
4632  unsigned int ind_local = 0;
4633 
4634  // check whether there is any constraint on the current cell
4635  for (; index_indicators != next_index_indicators;
4636  ++index_indicators)
4637  {
4638  const std::pair<unsigned short, unsigned short> indicator =
4639  this->dof_info->constraint_indicator[index_indicators];
4640 
4641  // run through values up to next constraint
4642  for (unsigned int j = 0; j < indicator.first; ++j)
4643  operation.process_dof(dof_indices[v][j],
4644  *src[0],
4645  values_dofs[comp][ind_local + j][v]);
4646  ind_local += indicator.first;
4647  dof_indices[v] += indicator.first;
4648 
4649  // constrained case: build the local value as a linear
4650  // combination of the global value according to constraints
4651  Number value;
4652  operation.pre_constraints(values_dofs[comp][ind_local][v],
4653  value);
4654 
4655  const Number *data_val =
4656  this->matrix_info->constraint_pool_begin(indicator.second);
4657  const Number *end_pool =
4658  this->matrix_info->constraint_pool_end(indicator.second);
4659 
4660  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4661  operation.process_constraint(*dof_indices[v],
4662  *data_val,
4663  *src[0],
4664  value);
4665 
4666  operation.post_constraints(value,
4667  values_dofs[comp][ind_local][v]);
4668  ind_local++;
4669  }
4670 
4671  AssertIndexRange(ind_local, dofs_per_component + 1);
4672 
4673  // get the dof values past the last constraint
4674  for (; ind_local < dofs_per_component;
4675  ++dof_indices[v], ++ind_local)
4676  {
4677  AssertIndexRange(*dof_indices[v], src[0]->size());
4678  operation.process_dof(*dof_indices[v],
4679  *src[0],
4680  values_dofs[comp][ind_local][v]);
4681  }
4682 
4683  if (apply_constraints == true && comp + 1 < n_components)
4684  next_index_indicators =
4685  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4686  }
4687  }
4688  }
4689 }
4690 
4691 
4692 
4693 template <int dim,
4694  int n_components_,
4695  typename Number,
4696  bool is_face,
4697  typename VectorizedArrayType>
4698 template <typename VectorType, typename VectorOperation>
4699 inline void
4702  const VectorOperation & operation,
4703  const std::array<VectorType *, n_components_> &src) const
4704 {
4705  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4706 
4707  unsigned int index =
4708  first_selected_component * this->data->dofs_per_component_on_cell;
4709  for (unsigned int comp = 0; comp < n_components; ++comp)
4710  {
4711  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4712  ++i, ++index)
4713  {
4714  operation.process_empty(values_dofs[comp][i]);
4715  operation.process_dof_global(
4716  local_dof_indices[this->data->lexicographic_numbering[index]],
4717  *src[0],
4718  values_dofs[comp][i][0]);
4719  }
4720  }
4721 }
4722 
4723 
4724 
4725 template <int dim,
4726  int n_components_,
4727  typename Number,
4728  bool is_face,
4729  typename VectorizedArrayType>
4730 template <typename VectorType, typename VectorOperation>
4731 inline void
4734  const VectorOperation & operation,
4735  const std::array<VectorType *, n_components_> &src,
4736  const std::array<
4738  n_components_> & vectors_sm,
4739  const std::bitset<VectorizedArrayType::size()> &mask) const
4740 {
4741  // This functions processes the functions read_dof_values,
4742  // distribute_local_to_global, and set_dof_values with the same code for
4743  // contiguous cell indices (DG case). The distinction between these three
4744  // cases is made by the input VectorOperation that either reads values from
4745  // a vector and puts the data into the local data field or write local data
4746  // into the vector. Certain operations are no-ops for the given use case.
4747 
4748  std::integral_constant<bool,
4750  vector_selector;
4752  is_face ? this->dof_access_index :
4754  const unsigned int n_lanes = mask.count();
4755 
4756  const std::vector<unsigned int> &dof_indices_cont =
4757  this->dof_info->dof_indices_contiguous[ind];
4758 
4759  // Simple case: We have contiguous storage, so we can simply copy out the
4760  // data
4761  if ((this->dof_info->index_storage_variants[ind][this->cell] ==
4763  interleaved_contiguous &&
4764  n_lanes == VectorizedArrayType::size()) &&
4765  !(is_face &&
4766  this->dof_access_index ==
4768  this->is_interior_face == false))
4769  {
4770  const unsigned int dof_index =
4771  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4772  this->dof_info->component_dof_indices_offset[this->active_fe_index]
4773  [first_selected_component] *
4774  VectorizedArrayType::size();
4775  if (n_components == 1 || n_fe_components == 1)
4776  for (unsigned int comp = 0; comp < n_components; ++comp)
4777  operation.process_dofs_vectorized(
4778  this->data->dofs_per_component_on_cell,
4779  dof_index,
4780  *src[comp],
4781  values_dofs[comp],
4782  vector_selector);
4783  else
4784  operation.process_dofs_vectorized(
4785  this->data->dofs_per_component_on_cell * n_components,
4786  dof_index,
4787  *src[0],
4788  values_dofs[0],
4789  vector_selector);
4790  return;
4791  }
4792 
4793  std::array<unsigned int, VectorizedArrayType::size()> cells =
4794  this->get_cell_or_face_ids();
4795 
4796  // More general case: Must go through the components one by one and apply
4797  // some transformations
4798  const unsigned int n_filled_lanes =
4799  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4800 
4801  const bool is_ecl =
4802  this->dof_access_index ==
4804  this->is_interior_face == false;
4805 
4806  if (vectors_sm[0] != nullptr)
4807  {
4808  const auto compute_vector_ptrs = [&](const unsigned int comp) {
4809  std::array<typename VectorType::value_type *,
4810  VectorizedArrayType::size()>
4811  vector_ptrs = {};
4812 
4813  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4814  {
4816  ExcNotImplemented());
4817  Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
4818  ExcIndexRange(
4819  ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
4820  Assert(cells[v] <
4821  this->dof_info->dof_indices_contiguous_sm[ind].size(),
4822  ExcIndexRange(
4823  cells[v],
4824  0,
4825  this->dof_info->dof_indices_contiguous_sm[ind].size()));
4826 
4827  const auto &temp =
4828  this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
4829 
4830  if (temp.first != numbers::invalid_unsigned_int)
4831  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
4832  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
4833  this->dof_info->component_dof_indices_offset
4834  [this->active_fe_index][this->first_selected_component]);
4835  else
4836  vector_ptrs[v] = nullptr;
4837  }
4838  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
4839  ++v)
4840  vector_ptrs[v] = nullptr;
4841 
4842  return vector_ptrs;
4843  };
4844 
4845  if (n_filled_lanes == VectorizedArrayType::size() &&
4846  n_lanes == VectorizedArrayType::size() && !is_ecl)
4847  {
4848  if (n_components == 1 || n_fe_components == 1)
4849  {
4850  for (unsigned int comp = 0; comp < n_components; ++comp)
4851  {
4852  auto vector_ptrs = compute_vector_ptrs(comp);
4853  operation.process_dofs_vectorized_transpose(
4854  this->data->dofs_per_component_on_cell,
4855  vector_ptrs,
4856  values_dofs[comp],
4857  vector_selector);
4858  }
4859  }
4860  else
4861  {
4862  auto vector_ptrs = compute_vector_ptrs(0);
4863  operation.process_dofs_vectorized_transpose(
4864  this->data->dofs_per_component_on_cell * n_components,
4865  vector_ptrs,
4866  &values_dofs[0][0],
4867  vector_selector);
4868  }
4869  }
4870  else
4871  for (unsigned int comp = 0; comp < n_components; ++comp)
4872  {
4873  auto vector_ptrs = compute_vector_ptrs(
4874  (n_components == 1 || n_fe_components == 1) ? comp : 0);
4875 
4876  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4877  ++i)
4878  operation.process_empty(values_dofs[comp][i]);
4879 
4880  if (n_components == 1 || n_fe_components == 1)
4881  {
4882  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4883  if (mask[v] == true)
4884  for (unsigned int i = 0;
4885  i < this->data->dofs_per_component_on_cell;
4886  ++i)
4887  operation.process_dof(vector_ptrs[v][i],
4888  values_dofs[comp][i][v]);
4889  }
4890  else
4891  {
4892  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4893  if (mask[v] == true)
4894  for (unsigned int i = 0;
4895  i < this->data->dofs_per_component_on_cell;
4896  ++i)
4897  operation.process_dof(
4898  vector_ptrs[v]
4899  [i + comp * this->data
4900  ->dofs_per_component_on_cell],
4901  values_dofs[comp][i][v]);
4902  }
4903  }
4904  return;
4905  }
4906 
4907  unsigned int dof_indices[VectorizedArrayType::size()];
4908 
4909  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4910  {
4912  dof_indices[v] =
4913  dof_indices_cont[cells[v]] +
4914  this->dof_info
4915  ->component_dof_indices_offset[this->active_fe_index]
4916  [this->first_selected_component] *
4917  this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
4918  }
4919 
4920  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4921  dof_indices[v] = numbers::invalid_unsigned_int;
4922 
4923  // In the case with contiguous cell indices, we know that there are no
4924  // constraints and that the indices within each element are contiguous
4925  if (n_filled_lanes == VectorizedArrayType::size() &&
4926  n_lanes == VectorizedArrayType::size() && !is_ecl)
4927  {
4928  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4930  contiguous)
4931  {
4932  if (n_components == 1 || n_fe_components == 1)
4933  for (unsigned int comp = 0; comp < n_components; ++comp)
4934  operation.process_dofs_vectorized_transpose(
4935  this->data->dofs_per_component_on_cell,
4936  dof_indices,
4937  *src[comp],
4938  values_dofs[comp],
4939  vector_selector);
4940  else
4941  operation.process_dofs_vectorized_transpose(
4942  this->data->dofs_per_component_on_cell * n_components,
4943  dof_indices,
4944  *src[0],
4945  &values_dofs[0][0],
4946  vector_selector);
4947  }
4948  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
4950  interleaved_contiguous_strided)
4951  {
4952  if (n_components == 1 || n_fe_components == 1)
4953  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4954  ++i)
4955  {
4956  for (unsigned int comp = 0; comp < n_components; ++comp)
4957  operation.process_dof_gather(dof_indices,
4958  *src[comp],
4959  i * VectorizedArrayType::size(),
4960  values_dofs[comp][i],
4961  vector_selector);
4962  }
4963  else
4964  for (unsigned int comp = 0; comp < n_components; ++comp)
4965  for (unsigned int i = 0;
4966  i < this->data->dofs_per_component_on_cell;
4967  ++i)
4968  {
4969  operation.process_dof_gather(
4970  dof_indices,
4971  *src[0],
4972  (comp * this->data->dofs_per_component_on_cell + i) *
4973  VectorizedArrayType::size(),
4974  values_dofs[comp][i],
4975  vector_selector);
4976  }
4977  }
4978  else
4979  {
4980  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
4982  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4983  ExcNotImplemented());
4984  const unsigned int *offsets =
4985  &this->dof_info->dof_indices_interleave_strides
4986  [ind][VectorizedArrayType::size() * this->cell];
4987  if (n_components == 1 || n_fe_components == 1)
4988  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4989  ++i)
4990  {
4991  for (unsigned int comp = 0; comp < n_components; ++comp)
4992  operation.process_dof_gather(dof_indices,
4993  *src[comp],
4994  0,
4995  values_dofs[comp][i],
4996  vector_selector);
4998  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4999  dof_indices[v] += offsets[v];
5000  }
5001  else
5002  for (unsigned int comp = 0; comp < n_components; ++comp)
5003  for (unsigned int i = 0;
5004  i < this->data->dofs_per_component_on_cell;
5005  ++i)
5006  {
5007  operation.process_dof_gather(dof_indices,
5008  *src[0],
5009  0,
5010  values_dofs[comp][i],
5011  vector_selector);
5013  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
5014  dof_indices[v] += offsets[v];
5015  }
5016  }
5017  }
5018  else
5019  for (unsigned int comp = 0; comp < n_components; ++comp)
5020  {
5021  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
5022  ++i)
5023  operation.process_empty(values_dofs[comp][i]);
5024  if (this->dof_info->index_storage_variants[ind][this->cell] ==
5026  contiguous)
5027  {
5028  if (n_components == 1 || n_fe_components == 1)
5029  {
5030  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5031  if (mask[v] == true)
5032  for (unsigned int i = 0;
5033  i < this->data->dofs_per_component_on_cell;
5034  ++i)
5035  operation.process_dof(dof_indices[v] + i,
5036  *src[comp],
5037  values_dofs[comp][i][v]);
5038  }
5039  else
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(
5047  dof_indices[v] + i +
5048  comp * this->data->dofs_per_component_on_cell,
5049  *src[0],
5050  values_dofs[comp][i][v]);
5051  }
5052  }
5053  else
5054  {
5055  const unsigned int *offsets =
5056  &this->dof_info->dof_indices_interleave_strides
5057  [ind][VectorizedArrayType::size() * this->cell];
5058  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5059  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
5060  if (n_components == 1 || n_fe_components == 1)
5061  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5062  {
5063  if (mask[v] == true)
5064  for (unsigned int i = 0;
5065  i < this->data->dofs_per_component_on_cell;
5066  ++i)
5067  operation.process_dof(dof_indices[v] + i * offsets[v],
5068  *src[comp],
5069  values_dofs[comp][i][v]);
5070  }
5071  else
5072  {
5073  for (unsigned int v = 0; v < n_filled_lanes; ++v)
5074  if (mask[v] == true)
5075  for (unsigned int i = 0;
5076  i < this->data->dofs_per_component_on_cell;
5077  ++i)
5078  operation.process_dof(
5079  dof_indices[v] +
5080  (i + comp * this->data->dofs_per_component_on_cell) *
5081  offsets[v],
5082  *src[0],
5083  values_dofs[comp][i][v]);
5084  }
5085  }
5086  }
5087 }
5088 
5089 namespace internal
5090 {
5091  template <typename Number,
5092  typename VectorType,
5093  typename std::enable_if<!IsBlockVector<VectorType>::value,
5094  VectorType>::type * = nullptr>
5095  decltype(std::declval<VectorType>().begin())
5096  get_beginning(VectorType &vec)
5097  {
5098  return vec.begin();
5099  }
5100 
5101  template <typename Number,
5102  typename VectorType,
5103  typename std::enable_if<IsBlockVector<VectorType>::value,
5104  VectorType>::type * = nullptr>
5105  typename VectorType::value_type *
5106  get_beginning(VectorType &)
5107  {
5108  return nullptr;
5109  }
5110 
5111  template <typename VectorType,
5112  typename std::enable_if<has_shared_vector_data<VectorType>::value,
5113  VectorType>::type * = nullptr>
5114  const std::vector<ArrayView<const typename VectorType::value_type>> *
5115  get_shared_vector_data(VectorType & vec,
5116  const bool is_valid_mode_for_sm,
5117  const unsigned int active_fe_index,
5119  {
5120  // note: no hp is supported
5121  if (is_valid_mode_for_sm &&
5122  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
5123  .size() > 0 &&
5124  active_fe_index == 0)
5125  return &vec.shared_vector_data();
5126  else
5127  return nullptr;
5128  }
5129 
5130  template <typename VectorType,
5131  typename std::enable_if<!has_shared_vector_data<VectorType>::value,
5132  VectorType>::type * = nullptr>
5133  const std::vector<ArrayView<const typename VectorType::value_type>> *
5134  get_shared_vector_data(VectorType &,
5135  const bool,
5136  const unsigned int,
5138  {
5139  return nullptr;
5140  }
5141 
5142  template <int n_components, typename VectorType>
5143  std::pair<
5144  std::array<typename internal::BlockVectorSelector<
5145  typename std::remove_const<VectorType>::type,
5147  value>::BaseVectorType *,
5148  n_components>,
5149  std::array<
5150  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
5151  typename std::remove_const<VectorType>::type,
5153  BaseVectorType::value_type>> *,
5154  n_components>>
5155  get_vector_data(VectorType & src,
5156  const unsigned int first_index,
5157  const bool is_valid_mode_for_sm,
5158  const unsigned int active_fe_index,
5160  {
5161  // select between block vectors and non-block vectors. Note that the number
5162  // of components is checked in the internal data
5163  std::pair<
5164  std::array<typename internal::BlockVectorSelector<
5165  typename std::remove_const<VectorType>::type,
5167  value>::BaseVectorType *,
5168  n_components>,
5169  std::array<
5170  const std::vector<
5171  ArrayView<const typename internal::BlockVectorSelector<
5172  typename std::remove_const<VectorType>::type,
5174  value>::BaseVectorType::value_type>> *,
5175  n_components>>
5176  src_data;
5177 
5178  for (unsigned int d = 0; d < n_components; ++d)
5179  src_data.first[d] = internal::BlockVectorSelector<
5180  typename std::remove_const<VectorType>::type,
5181  IsBlockVector<typename std::remove_const<VectorType>::type>::value>::
5182  get_vector_component(
5183  const_cast<typename std::remove_const<VectorType>::type &>(src),
5184  d + first_index);
5185 
5186  for (unsigned int d = 0; d < n_components; ++d)
5187  src_data.second[d] = get_shared_vector_data(*src_data.first[d],
5188  is_valid_mode_for_sm,
5189  active_fe_index,
5190  dof_info);
5191 
5192  return src_data;
5193  }
5194 } // namespace internal
5195 
5196 
5197 
5198 template <int dim,
5199  int n_components_,
5200  typename Number,
5201  bool is_face,
5202  typename VectorizedArrayType>
5203 template <typename VectorType>
5204 inline void
5206  read_dof_values(const VectorType &src, const unsigned int first_index)
5207 {
5208  const auto src_data = internal::get_vector_data<n_components_>(
5209  src,
5210  first_index,
5211  this->dof_access_index ==
5213  this->active_fe_index,
5214  this->dof_info);
5215 
5217  read_write_operation(reader,
5218  src_data.first,
5219  src_data.second,
5220  std::bitset<VectorizedArrayType::size()>().flip(),
5221  true);
5222 
5223 # ifdef DEBUG
5224  dof_values_initialized = true;
5225 # endif
5226 }
5227 
5228 
5229 
5230 template <int dim,
5231  int n_components_,
5232  typename Number,
5233  bool is_face,
5234  typename VectorizedArrayType>
5235 template <typename VectorType>
5236 inline void
5238  read_dof_values_plain(const VectorType &src, const unsigned int first_index)
5239 {
5240  const auto src_data = internal::get_vector_data<n_components_>(
5241  src,
5242  first_index,
5243  this->dof_access_index ==
5245  this->active_fe_index,
5246  this->dof_info);
5247 
5249  read_write_operation(reader,
5250  src_data.first,
5251  src_data.second,
5252  std::bitset<VectorizedArrayType::size()>().flip(),
5253  false);
5254 
5255 # ifdef DEBUG
5256  dof_values_initialized = true;
5257 # endif
5258 }
5259 
5260 
5261 
5262 template <int dim,
5263  int n_components_,
5264  typename Number,
5265  bool is_face,
5266  typename VectorizedArrayType>
5267 template <typename VectorType>
5268 inline void
5271  VectorType & dst,
5272  const unsigned int first_index,
5273  const std::bitset<VectorizedArrayType::size()> &mask) const
5274 {
5275 # ifdef DEBUG
5276  Assert(dof_values_initialized == true,
5278 # endif
5279 
5280  const auto dst_data = internal::get_vector_data<n_components_>(
5281  dst,
5282  first_index,
5283  this->dof_access_index ==
5285  this->active_fe_index,
5286  this->dof_info);
5287 
5289  distributor;
5290  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
5291 }
5292 
5293 
5294 
5295 template <int dim,
5296  int n_components_,
5297  typename Number,
5298  bool is_face,
5299  typename VectorizedArrayType>
5300 template <typename VectorType>
5301 inline void
5304  const unsigned int first_index,
5305  const std::bitset<VectorizedArrayType::size()> &mask) const
5306 {
5307 # ifdef DEBUG
5308  Assert(dof_values_initialized == true,
5310 # endif
5311 
5312  const auto dst_data = internal::get_vector_data<n_components_>(
5313  dst,
5314  first_index,
5315  this->dof_access_index ==
5317  this->active_fe_index,
5318  this->dof_info);
5319 
5321  read_write_operation(setter, dst_data.first, dst_data.second, mask);
5322 }
5323 
5324 
5325 
5326 template <int dim,
5327  int n_components_,
5328  typename Number,
5329  bool is_face,
5330  typename VectorizedArrayType>
5331 template <typename VectorType>
5332 inline void
5335  VectorType & dst,
5336  const unsigned int first_index,
5337  const std::bitset<VectorizedArrayType::size()> &mask) const
5338 {
5339 # ifdef DEBUG
5340  Assert(dof_values_initialized == true,
5342 # endif
5343 
5344  const auto dst_data = internal::get_vector_data<n_components_>(
5345  dst,
5346  first_index,
5347  this->dof_access_index ==
5349  this->active_fe_index,
5350  this->dof_info);
5351 
5353  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
5354 }
5355 
5356 
5357 
5358 /*------------------------------ access to data fields ----------------------*/
5359 
5360 
5361 
5362 template <int dim,
5363  int n_components,
5364  typename Number,
5365  bool is_face,
5366  typename VectorizedArrayType>
5367 inline const VectorizedArrayType *
5369  begin_dof_values() const
5370 {
5371  return &values_dofs[0][0];
5372 }
5373 
5374 
5375 
5376 template <int dim,
5377  int n_components,
5378  typename Number,
5379  bool is_face,
5380  typename VectorizedArrayType>
5381 inline VectorizedArrayType *
5384 {
5385 # ifdef DEBUG
5386  dof_values_initialized = true;
5387 # endif
5388  return &values_dofs[0][0];
5389 }
5390 
5391 
5392 
5393 template <int dim,
5394  int n_components,
5395  typename Number,
5396  bool is_face,
5397  typename VectorizedArrayType>
5398 inline const VectorizedArrayType *
5400  begin_values() const
5401 {
5402 # ifdef DEBUG
5403  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
5404 # endif
5405  return values_quad;
5406 }
5407 
5408 
5409 
5410 template <int dim,
5411  int n_components,
5412  typename Number,
5413  bool is_face,
5414  typename VectorizedArrayType>
5415 inline VectorizedArrayType *
5417  begin_values()
5418 {
5419 # ifdef DEBUG
5420  values_quad_initialized = true;
5421  values_quad_submitted = true;
5422 # endif
5423  return values_quad;
5424 }
5425 
5426 
5427 
5428 template <int dim,
5429  int n_components,
5430  typename Number,
5431  bool is_face,
5432  typename VectorizedArrayType>
5433 inline const VectorizedArrayType *
5435  begin_gradients() const
5436 {
5437 # ifdef DEBUG
5438  Assert(gradients_quad_initialized || gradients_quad_submitted,
5439  ExcNotInitialized());
5440 # endif
5441  return gradients_quad;
5442 }
5443 
5444 
5445 
5446 template <int dim,
5447  int n_components,
5448  typename Number,
5449  bool is_face,
5450  typename VectorizedArrayType>
5451 inline VectorizedArrayType *
5454 {
5455 # ifdef DEBUG
5456  gradients_quad_submitted = true;
5457  gradients_quad_initialized = true;
5458 # endif
5459  return gradients_quad;
5460 }
5461 
5462 
5463 
5464 template <int dim,
5465  int n_components,
5466  typename Number,
5467  bool is_face,
5468  typename VectorizedArrayType>
5469 inline const VectorizedArrayType *
5471  begin_hessians() const
5472 {
5473 # ifdef DEBUG
5474  Assert(hessians_quad_initialized, ExcNotInitialized());
5475 # endif
5476  return hessians_quad;
5477 }
5478 
5479 
5480 
5481 template <int dim,
5482  int n_components,
5483  typename Number,
5484  bool is_face,
5485  typename VectorizedArrayType>
5486 inline VectorizedArrayType *
5489 {
5490 # ifdef DEBUG
5491  hessians_quad_initialized = true;
5492 # endif
5493  return hessians_quad;
5494 }
5495 
5496 
5497 
5498 template <int dim,
5499  int n_components,
5500  typename Number,
5501  bool is_face,
5502  typename VectorizedArrayType>
5503 inline unsigned int
5506 {
5507  return first_selected_component;
5508 }
5509 
5510 
5511 
5512 template <int dim,
5513  int n_components_,
5514  typename Number,
5515  bool is_face,
5516  typename VectorizedArrayType>
5519  get_dof_value(const unsigned int dof) const
5520 {
5521  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5523  for (unsigned int comp = 0; comp < n_components; comp++)
5524  return_value[comp] = this->values_dofs[comp][dof];
5525  return return_value;
5526 }
5527 
5528 
5529 
5530 template <int dim,
5531  int n_components_,
5532  typename Number,
5533  bool is_face,
5534  typename VectorizedArrayType>
5537  get_value(const unsigned int q_point) const
5538 {
5539 # ifdef DEBUG
5540  Assert(this->values_quad_initialized == true,
5542 # endif
5543 
5544  AssertIndexRange(q_point, this->n_quadrature_points);
5545  const std::size_t nqp = this->n_quadrature_points;
5547  for (unsigned int comp = 0; comp < n_components; comp++)
5548  return_value[comp] = values_quad[comp * nqp + q_point];
5549  return return_value;
5550 }
5551 
5552 
5553 
5554 template <int dim,
5555  int n_components_,
5556  typename Number,
5557  bool is_face,
5558  typename VectorizedArrayType>
5559 inline DEAL_II_ALWAYS_INLINE
5562  get_gradient(const unsigned int q_point) const
5563 {
5564 # ifdef DEBUG
5565  Assert(this->gradients_quad_initialized == true,
5567 # endif
5568 
5569  AssertIndexRange(q_point, this->n_quadrature_points);
5570  Assert(this->jacobian != nullptr,
5572  "update_gradients"));
5573  const std::size_t nqp = this->n_quadrature_points;
5575 
5576  // Cartesian cell
5577  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5578  {
5579  for (unsigned int d = 0; d < dim; ++d)
5580  for (unsigned int comp = 0; comp < n_components; comp++)
5581  grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] *
5582  this->jacobian[0][d][d];
5583  }
5584  // cell with general/affine Jacobian
5585  else
5586  {
5588  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5589  q_point :
5590  0];
5591  for (unsigned int comp = 0; comp < n_components; comp++)
5592  for (unsigned int d = 0; d < dim; ++d)
5593  {
5594  grad_out[comp][d] =
5595  jac[d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5596  for (unsigned int e = 1; e < dim; ++e)
5597  grad_out[comp][d] +=
5598  jac[d][e] * gradients_quad[(comp * dim + e) * nqp + q_point];
5599  }
5600  }
5601  return grad_out;
5602 }
5603 
5604 
5605 
5606 template <int dim,
5607  int n_components_,
5608  typename Number,
5609  bool is_face,
5610  typename VectorizedArrayType>
5613  get_normal_derivative(const unsigned int q_point) const
5614 {
5615  AssertIndexRange(q_point, this->n_quadrature_points);
5616 # ifdef DEBUG
5617  Assert(this->gradients_quad_initialized == true,
5619 # endif
5620 
5621  Assert(this->normal_x_jacobian != nullptr,
5623  "update_gradients"));
5624 
5625  const std::size_t nqp = this->n_quadrature_points;
5627 
5628  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5629  for (unsigned int comp = 0; comp < n_components; comp++)
5630  grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5631  (this->normal_x_jacobian[0][dim - 1]);
5632  else
5633  {
5634  const std::size_t index =
5635  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5636  for (unsigned int comp = 0; comp < n_components; comp++)
5637  {
5638  grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5639  this->normal_x_jacobian[index][0];
5640  for (unsigned int d = 1; d < dim; ++d)
5641  grad_out[comp] += gradients_quad[(comp * dim + d) * nqp + q_point] *
5642  this->normal_x_jacobian[index][d];
5643  }
5644  }
5645  return grad_out;
5646 }
5647 
5648 
5649 
5650 namespace internal
5651 {
5652  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5653  // store the lower diagonal because of symmetry
5654  template <typename VectorizedArrayType>
5655  inline void
5656  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
5657  const VectorizedArrayType *const hessians,
5658  const unsigned int,
5659  VectorizedArrayType (&tmp)[1][1])
5660  {
5661  tmp[0][0] = jac[0][0] * hessians[0];
5662  }
5663 
5664  template <typename VectorizedArrayType>
5665  inline void
5666  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
5667  const VectorizedArrayType *const hessians,
5668  const unsigned int nqp,
5669  VectorizedArrayType (&tmp)[2][2])
5670  {
5671  for (unsigned int d = 0; d < 2; ++d)
5672  {
5673  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
5674  tmp[1][d] =
5675  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
5676  }
5677  }
5678 
5679  template <typename VectorizedArrayType>
5680  inline void
5681  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
5682  const VectorizedArrayType *const hessians,
5683  const unsigned int nqp,
5684  VectorizedArrayType (&tmp)[3][3])
5685  {
5686  for (unsigned int d = 0; d < 3; ++d)
5687  {
5688  tmp[0][d] =
5689  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
5690  jac[d][2] * hessians[4 * nqp]);
5691  tmp[1][d] =
5692  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
5693  jac[d][2] * hessians[5 * nqp]);
5694  tmp[2][d] =
5695  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
5696  jac[d][2] * hessians[2 * nqp]);
5697  }
5698  }
5699 } // namespace internal
5700 
5701 
5702 
5703 template <int dim,
5704  int n_components_,
5705  typename Number,
5706  bool is_face,
5707  typename VectorizedArrayType>
5710  get_hessian(const unsigned int q_point) const
5711 {
5712  Assert(!is_face, ExcNotImplemented());
5713 # ifdef DEBUG
5714  Assert(this->hessians_quad_initialized == true,
5716 # endif
5717  AssertIndexRange(q_point, this->n_quadrature_points);
5718 
5719  Assert(this->jacobian != nullptr,
5721  "update_hessian"));
5723  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5724  0 :
5725  q_point];
5726 
5728 
5729  const std::size_t nqp = this->n_quadrature_points;
5730  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5731 
5732  // Cartesian cell
5733  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5734  {
5735  for (unsigned int comp = 0; comp < n_components; comp++)
5736  {
5737  for (unsigned int d = 0; d < dim; ++d)
5738  hessian_out[comp][d][d] =
5739  hessians_quad[(comp * hdim + d) * nqp + q_point] *
5740  (jac[d][d] * jac[d][d]);
5741  switch (dim)
5742  {
5743  case 1:
5744  break;
5745  case 2:
5746  hessian_out[comp][0][1] =
5747  hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5748  (jac[0][0] * jac[1][1]);
5749  break;
5750  case 3:
5751  hessian_out[comp][0][1] =
5752  hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5753  (jac[0][0] * jac[1][1]);
5754  hessian_out[comp][0][2] =
5755  hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5756  (jac[0][0] * jac[2][2]);
5757  hessian_out[comp][1][2] =
5758  hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5759  (jac[1][1] * jac[2][2]);
5760  break;
5761  default:
5762  Assert(false, ExcNotImplemented());
5763  }
5764  for (unsigned int d = 0; d < dim; ++d)
5765  for (unsigned int e = d + 1; e < dim; ++e)
5766  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5767  }
5768  }
5769  // cell with general Jacobian, but constant within the cell
5770  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5771  {
5772  for (unsigned int comp = 0; comp < n_components; comp++)
5773  {
5774  VectorizedArrayType tmp[dim][dim];
5775  internal::hessian_unit_times_jac(
5776  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5777 
5778  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5779  for (unsigned int d = 0; d < dim; ++d)
5780  for (unsigned int e = d; e < dim; ++e)
5781  {
5782  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5783  for (unsigned int f = 1; f < dim; ++f)
5784  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5785  }
5786 
5787  // no J' * grad(u) part here because the Jacobian is constant
5788  // throughout the cell and hence, its derivative is zero
5789 
5790  // take symmetric part
5791  for (unsigned int d = 0; d < dim; ++d)
5792  for (unsigned int e = d + 1; e < dim; ++e)
5793  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5794  }
5795  }
5796  // cell with general Jacobian
5797  else
5798  {
5799  const auto &jac_grad =
5800  this->mapping_data->jacobian_gradients
5801  [1 - this->is_interior_face]
5802  [this->mapping_data->data_index_offsets[this->cell] + q_point];
5803  for (unsigned int comp = 0; comp < n_components; comp++)
5804  {
5805  // compute laplacian before the gradient because it needs to access
5806  // unscaled gradient data
5807  VectorizedArrayType tmp[dim][dim];
5808  internal::hessian_unit_times_jac(
5809  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5810 
5811  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5812  for (unsigned int d = 0; d < dim; ++d)
5813  for (unsigned int e = d; e < dim; ++e)
5814  {
5815  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5816  for (unsigned int f = 1; f < dim; ++f)
5817  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5818  }
5819 
5820  // add diagonal part of J' * grad(u)
5821  for (unsigned int d = 0; d < dim; ++d)
5822  for (unsigned int e = 0; e < dim; ++e)
5823  hessian_out[comp][d][d] +=
5824  jac_grad[d][e] *
5825  gradients_quad[(comp * dim + e) * nqp + q_point];
5826 
5827  // add off-diagonal part of J' * grad(u)
5828  for (unsigned int d = 0, count = dim; d < dim; ++d)
5829  for (unsigned int e = d + 1; e < dim; ++e, ++count)
5830  for (unsigned int f = 0; f < dim; ++f)
5831  hessian_out[comp][d][e] +=
5832  jac_grad[count][f] *
5833  gradients_quad[(comp * dim + f) * nqp + q_point];
5834 
5835  // take symmetric part
5836  for (unsigned int d = 0; d < dim; ++d)
5837  for (unsigned int e = d + 1; e < dim; ++e)
5838  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5839  }
5840  }
5841  return hessian_out;
5842 }
5843 
5844 
5845 
5846 template <int dim,
5847  int n_components_,
5848  typename Number,
5849  bool is_face,
5850  typename VectorizedArrayType>
5853  get_hessian_diagonal(const unsigned int q_point) const
5854 {
5855  Assert(!is_face, ExcNotImplemented());
5856 # ifdef DEBUG
5857  Assert(this->hessians_quad_initialized == true,
5859 # endif
5860  AssertIndexRange(q_point, this->n_quadrature_points);
5861 
5862  Assert(this->jacobian != nullptr, ExcNotImplemented());
5864  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5865  0 :
5866  q_point];
5867 
5868  const std::size_t nqp = this->n_quadrature_points;
5869  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5871 
5872  // Cartesian cell
5873  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5874  {
5875  for (unsigned int comp = 0; comp < n_components; comp++)
5876  for (unsigned int d = 0; d < dim; ++d)
5877  hessian_out[comp][d] =
5878  hessians_quad[(comp * hdim + d) * nqp + q_point] *
5879  (jac[d][d] * jac[d][d]);
5880  }
5881  // cell with general Jacobian, but constant within the cell
5882  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5883  {
5884  for (unsigned int comp = 0; comp < n_components; comp++)
5885  {
5886  // compute laplacian before the gradient because it needs to access
5887  // unscaled gradient data
5888  VectorizedArrayType tmp[dim][dim];
5889  internal::hessian_unit_times_jac(
5890  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5891 
5892  // compute only the trace part of hessian, J * tmp = J *
5893  // hess_unit(u) * J^T
5894  for (unsigned int d = 0; d < dim; ++d)
5895  {
5896  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5897  for (unsigned int f = 1; f < dim; ++f)
5898  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5899  }
5900  }
5901  }
5902  // cell with general Jacobian
5903  else
5904  {
5905  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5906  &jac_grad =
5907  this->mapping_data->jacobian_gradients
5908  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5909  for (unsigned int comp = 0; comp < n_components; comp++)
5910  {
5911  // compute laplacian before the gradient because it needs to access
5912  // unscaled gradient data
5913  VectorizedArrayType tmp[dim][dim];
5914  internal::hessian_unit_times_jac(
5915  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5916 
5917  // compute only the trace part of hessian, J * tmp = J *
5918  // hess_unit(u) * J^T
5919  for (unsigned int d = 0; d < dim; ++d)
5920  {
5921  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5922  for (unsigned int f = 1; f < dim; ++f)
5923  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5924  }
5925 
5926  for (unsigned int d = 0; d < dim; ++d)
5927  for (unsigned int e = 0; e < dim; ++e)
5928  hessian_out[comp][d] +=
5929  jac_grad[d][e] *
5930  gradients_quad[(comp * dim + e) * nqp + q_point];
5931  }
5932  }
5933  return hessian_out;
5934 }
5935 
5936 
5937 
5938 template <int dim,
5939  int n_components_,
5940  typename Number,
5941  bool is_face,
5942  typename VectorizedArrayType>
5945  get_laplacian(const unsigned int q_point) const
5946 {
5947  Assert(is_face == false, ExcNotImplemented());
5948 # ifdef DEBUG
5949  Assert(this->hessians_quad_initialized == true,
5951 # endif
5952  AssertIndexRange(q_point, this->n_quadrature_points);
5953 
5955  const auto hess_diag = get_hessian_diagonal(q_point);
5956  for (unsigned int comp = 0; comp < n_components; ++comp)
5957  {
5958  laplacian_out[comp] = hess_diag[comp][0];
5959  for (unsigned int d = 1; d < dim; ++d)
5960  laplacian_out[comp] += hess_diag[comp][d];
5961  }
5962  return laplacian_out;
5963 }
5964 
5965 
5966 
5967 template <int dim,
5968  int n_components_,
5969  typename Number,
5970  bool is_face,
5971  typename VectorizedArrayType>
5972 inline DEAL_II_ALWAYS_INLINE void
5975  const unsigned int dof)
5976 {
5977 # ifdef DEBUG
5978  this->dof_values_initialized = true;
5979 # endif
5980  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5981  for (unsigned int comp = 0; comp < n_components; comp++)
5982  this->values_dofs[comp][dof] = val_in[comp];
5983 }
5984 
5985 
5986 
5987 template <int dim,
5988  int n_components_,
5989  typename Number,
5990  bool is_face,
5991  typename VectorizedArrayType>
5992 inline DEAL_II_ALWAYS_INLINE void
5995  const unsigned int q_point)
5996 {
5998  AssertIndexRange(q_point, this->n_quadrature_points);
5999  Assert(this->J_value != nullptr,
6001  "update_values"));
6002 # ifdef DEBUG
6003  this->values_quad_submitted = true;
6004 # endif
6005 
6006  const std::size_t nqp = this->n_quadrature_points;
6007  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6008  {
6009  const VectorizedArrayType JxW =
6010  this->J_value[0] * this->quadrature_weights[q_point];
6011  for (unsigned int comp = 0; comp < n_components; ++comp)
6012  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6013  }
6014  else
6015  {
6016  const VectorizedArrayType JxW = this->J_value[q_point];
6017  for (unsigned int comp = 0; comp < n_components; ++comp)
6018  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
6019  }
6020 }
6021 
6022 
6023 
6024 template <int dim,
6025  int n_components_,
6026  typename Number,
6027  bool is_face,
6028  typename VectorizedArrayType>
6029 inline DEAL_II_ALWAYS_INLINE void
6032  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
6033  const unsigned int q_point)
6034 {
6036  AssertIndexRange(q_point, this->n_quadrature_points);
6037  Assert(this->J_value != nullptr,
6039  "update_gradients"));
6040  Assert(this->jacobian != nullptr,
6042  "update_gradients"));
6043 # ifdef DEBUG
6044  this->gradients_quad_submitted = true;
6045 # endif
6046 
6047  const std::size_t nqp = this->n_quadrature_points;
6048  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6049  {
6050  const VectorizedArrayType JxW =
6051  this->J_value[0] * this->quadrature_weights[q_point];
6052  for (unsigned int d = 0; d < dim; ++d)
6053  {
6054  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
6055  for (unsigned int comp = 0; comp < n_components; comp++)
6056  gradients_quad[(comp * dim + d) * nqp + q_point] =
6057  grad_in[comp][d] * factor;
6058  }
6059  }
6060  else
6061  {
6063  this->cell_type > internal::MatrixFreeFunctions::affine ?
6064  this->jacobian[q_point] :
6065  this->jacobian[0];
6066  const VectorizedArrayType JxW =
6067  this->cell_type > internal::MatrixFreeFunctions::affine ?
6068  this->J_value[q_point] :
6069  this->J_value[0] * this->quadrature_weights[q_point];
6070  for (unsigned int comp = 0; comp < n_components; ++comp)
6071  for (unsigned int d = 0; d < dim; ++d)
6072  {
6073  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
6074  for (unsigned int e = 1; e < dim; ++e)
6075  new_val += (jac[e][d] * grad_in[comp][e]);
6076  gradients_quad[(comp * dim + d) * nqp + q_point] = new_val * JxW;
6077  }
6078  }
6079 }
6080 
6081 
6082 
6083 template <int dim,
6084  int n_components_,
6085  typename Number,
6086  bool is_face,
6087  typename VectorizedArrayType>
6088 inline DEAL_II_ALWAYS_INLINE void
6092  const unsigned int q_point)
6093 {
6094  AssertIndexRange(q_point, this->n_quadrature_points);
6095  Assert(this->normal_x_jacobian != nullptr,
6097  "update_gradients"));
6098 # ifdef DEBUG
6099  this->gradients_quad_submitted = true;
6100 # endif
6101 
6102  const std::size_t nqp = this->n_quadrature_points;
6103  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
6104  for (unsigned int comp = 0; comp < n_components; comp++)
6105  {
6106  for (unsigned int d = 0; d < dim - 1; ++d)
6107  gradients_quad[(comp * dim + d) * nqp + q_point] =
6108  VectorizedArrayType();
6109  gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
6110  grad_in[comp] *
6111  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
6112  this->quadrature_weights[q_point]);
6113  }
6114  else
6115  {
6116  const unsigned int index =
6117  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
6119  this->normal_x_jacobian[index];
6120  for (unsigned int comp = 0; comp < n_components; comp++)
6121  {
6122  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
6123  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6124  factor = factor * this->quadrature_weights[q_point];
6125  for (unsigned int d = 0; d < dim; ++d)
6126  gradients_quad[(comp * dim + d) * nqp + q_point] = factor * jac[d];
6127  }
6128  }
6129 }
6130 
6131 
6132 
6133 template <int dim,
6134  int n_components_,
6135  typename Number,
6136  bool is_face,
6137  typename VectorizedArrayType>
6140  integrate_value() const
6141 {
6143 # ifdef DEBUG
6144  Assert(this->values_quad_submitted == true,
6146 # endif
6147 
6149  const std::size_t nqp = this->n_quadrature_points;
6150  for (unsigned int q = 0; q < nqp; ++q)
6151  for (unsigned int comp = 0; comp < n_components; ++comp)
6152  return_value[comp] += this->values_quad[comp * nqp + q];
6153  return (return_value);
6154 }
6155 
6156 
6157 
6158 /*----------------------- FEEvaluationAccess --------------------------------*/
6159 
6160 
6161 template <int dim,
6162  int n_components_,
6163  typename Number,
6164  bool is_face,
6165  typename VectorizedArrayType>
6166 inline FEEvaluationAccess<dim,
6167  n_components_,
6168  Number,
6169  is_face,
6170  VectorizedArrayType>::
6171  FEEvaluationAccess(
6173  const unsigned int dof_no,
6174  const unsigned int first_selected_component,
6175  const unsigned int quad_no_in,
6176  const unsigned int fe_degree,
6177  const unsigned int n_q_points,
6178  const bool is_interior_face,
6179  const unsigned int active_fe_index,
6180  const unsigned int active_quad_index,
6181  const unsigned int face_type)
6183  data_in,
6184  dof_no,
6185  first_selected_component,
6186  quad_no_in,
6187  fe_degree,
6188  n_q_points,
6189  is_interior_face,
6190  active_fe_index,
6191  active_quad_index,
6192  face_type)
6193 {}
6194 
6195 
6196 
6197 template <int dim,
6198  int n_components_,
6199  typename Number,
6200  bool is_face,
6201  typename VectorizedArrayType>
6202 inline FEEvaluationAccess<dim,
6203  n_components_,
6204  Number,
6205  is_face,
6206  VectorizedArrayType>::
6207  FEEvaluationAccess(
6208  const Mapping<dim> & mapping,
6209  const FiniteElement<dim> &fe,
6210  const Quadrature<1> & quadrature,
6211  const UpdateFlags update_flags,
6212  const unsigned int first_selected_component,
6214  *other)
6216  mapping,
6217  fe,
6218  quadrature,
6219  update_flags,
6220  first_selected_component,
6221  other)
6222 {}
6223 
6224 
6225 
6226 template <int dim,
6227  int n_components_,
6228  typename Number,
6229  bool is_face,
6230  typename VectorizedArrayType>
6231 inline FEEvaluationAccess<dim,
6232  n_components_,
6233  Number,
6234  is_face,
6235  VectorizedArrayType>::
6236  FEEvaluationAccess(const FEEvaluationAccess<dim,
6237  n_components_,
6238  Number,
6239  is_face,
6240  VectorizedArrayType> &other)
6242  other)
6243 {}
6244 
6245 
6246 
6247 template <int dim,
6248  int n_components_,
6249  typename Number,
6250  bool is_face,
6251  typename VectorizedArrayType>
6252 inline FEEvaluationAccess<dim,
6253  n_components_,
6254  Number,
6255  is_face,
6256  VectorizedArrayType> &
6258 operator=(const FEEvaluationAccess<dim,
6259  n_components_,
6260  Number,
6261  is_face,
6262  VectorizedArrayType> &other)
6263 {
6264  this->FEEvaluationBase<dim,
6265  n_components_,
6266  Number,
6267  is_face,
6268  VectorizedArrayType>::operator=(other);
6269  return *this;
6270 }
6271 
6272 
6273 
6274 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
6275 
6276 
6277 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6281  const unsigned int dof_no,
6282  const unsigned int first_selected_component,
6283  const unsigned int quad_no_in,
6284  const unsigned int fe_degree,
6285  const unsigned int n_q_points,
6286  const bool is_interior_face,
6287  const unsigned int active_fe_index,
6288  const unsigned int active_quad_index,
6289  const unsigned int face_type)
6291  data_in,
6292  dof_no,
6293  first_selected_component,
6294  quad_no_in,
6295  fe_degree,
6296  n_q_points,
6297  is_interior_face,
6298  active_fe_index,
6299  active_quad_index,
6300  face_type)
6301 {}
6302 
6303 
6304 
6305 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6308  const Mapping<dim> & mapping,
6309  const FiniteElement<dim> &fe,
6310  const Quadrature<1> & quadrature,
6311  const UpdateFlags update_flags,
6312  const unsigned int first_selected_component,
6314  *other)
6316  mapping,
6317  fe,
6318  quadrature,
6319  update_flags,
6320  first_selected_component,
6321  other)
6322 {}
6323 
6324 
6325 
6326 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6330  &other)
6332 {}
6333 
6334 
6335 
6336 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6340 {
6342  operator=(other);
6343  return *this;
6344 }
6345 
6346 
6347 
6348 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6349 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6351  const unsigned int dof) const
6352 {
6353  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6354  return this->values_dofs[0][dof];
6355 }
6356 
6357 
6358 
6359 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6360 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6362  const unsigned int q_point) const
6363 {
6364 # ifdef DEBUG
6365  Assert(this->values_quad_initialized == true,
6367 # endif
6368  AssertIndexRange(q_point, this->n_quadrature_points);
6369  return this->values_quad[q_point];
6370 }
6371 
6372 
6373 
6374 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6375 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6377  get_normal_derivative(const unsigned int q_point) const
6378 {
6379  return BaseClass::get_normal_derivative(q_point)[0];
6380 }
6381 
6382 
6383 
6384 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6387  const unsigned int q_point) const
6388 {
6389  // could use the base class gradient, but that involves too many expensive
6390  // initialization operations on tensors
6391 
6392 # ifdef DEBUG
6393  Assert(this->gradients_quad_initialized == true,
6395 # endif
6396  AssertIndexRange(q_point, this->n_quadrature_points);
6397 
6398  Assert(this->jacobian != nullptr,
6400  "update_gradients"));
6401 
6403 
6404  const std::size_t nqp = this->n_quadrature_points;
6405  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6406  {
6407  for (unsigned int d = 0; d < dim; ++d)
6408  grad_out[d] =
6409  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
6410  }
6411  // cell with general/affine Jacobian
6412  else
6413  {
6415  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
6416  q_point :
6417  0];
6418  for (unsigned int d = 0; d < dim; ++d)
6419  {
6420  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
6421  for (unsigned int e = 1; e < dim; ++e)
6422  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
6423  }
6424  }
6425  return grad_out;
6426 }
6427 
6428 
6429 
6430 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6433  const unsigned int q_point) const
6434 {
6435  return BaseClass::get_hessian(q_point)[0];
6436 }
6437 
6438 
6439 
6440 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6443  get_hessian_diagonal(const unsigned int q_point) const
6444 {
6445  return BaseClass::get_hessian_diagonal(q_point)[0];
6446 }
6447 
6448 
6449 
6450 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6451 inline VectorizedArrayType
6453  const unsigned int q_point) const
6454 {
6455  return BaseClass::get_laplacian(q_point)[0];
6456 }
6457 
6458 
6459 
6460 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6461 inline void DEAL_II_ALWAYS_INLINE
6463  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6464 {
6465 # ifdef DEBUG
6466  this->dof_values_initialized = true;
6467  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6468 # endif
6469  this->values_dofs[0][dof] = val_in;
6470 }
6471 
6472 
6473 
6474 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6475 inline void DEAL_II_ALWAYS_INLINE
6477  const VectorizedArrayType val_in,
6478  const unsigned int q_point)
6479 {
6481  AssertIndexRange(q_point, this->n_quadrature_points);
6482  Assert(this->J_value != nullptr,
6484  "update_value"));
6485 # ifdef DEBUG
6486  this->values_quad_submitted = true;
6487 # endif
6488 
6489  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6490  {
6491  const VectorizedArrayType JxW =
6492  this->J_value[0] * this->quadrature_weights[q_point];
6493  this->values_quad[q_point] = val_in * JxW;
6494  }
6495  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
6496  {
6497  this->values_quad[q_point] = val_in * this->J_value[q_point];
6498  }
6499 }
6500 
6501 
6502 
6503 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6504 inline DEAL_II_ALWAYS_INLINE void
6506  const Tensor<1, 1, VectorizedArrayType> val_in,
6507  const unsigned int q_point)
6508 {
6509  submit_value(val_in[0], q_point);
6510 }
6511 
6512 
6513 
6514 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6515 inline DEAL_II_ALWAYS_INLINE void
6517  submit_normal_derivative(const VectorizedArrayType grad_in,
6518  const unsigned int q_point)
6519 {
6521  grad[0] = grad_in;
6522  BaseClass::submit_normal_derivative(grad, q_point);
6523 }
6524 
6525 
6526 
6527 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6528 inline DEAL_II_ALWAYS_INLINE void
6531  const unsigned int q_point)
6532 {
6534  AssertIndexRange(q_point, this->n_quadrature_points);
6535  Assert(this->J_value != nullptr,
6537  "update_gradients"));
6538  Assert(this->jacobian != nullptr,
6540  "update_gradients"));
6541 # ifdef DEBUG
6542  this->gradients_quad_submitted = true;
6543 # endif
6544 
6545  const std::size_t nqp = this->n_quadrature_points;
6546  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6547  {
6548  const VectorizedArrayType JxW =
6549  this->J_value[0] * this->quadrature_weights[q_point];
6550  for (unsigned int d = 0; d < dim; ++d)
6551  this->gradients_quad[d * nqp + q_point] =
6552  (grad_in[d] * this->jacobian[0][d][d] * JxW);
6553  }
6554  // general/affine cell type
6555  else
6556  {
6558  this->cell_type > internal::MatrixFreeFunctions::affine ?
6559  this->jacobian[q_point] :
6560  this->jacobian[0];
6561  const VectorizedArrayType JxW =
6562  this->cell_type > internal::MatrixFreeFunctions::affine ?
6563  this->J_value[q_point] :
6564  this->J_value[0] * this->quadrature_weights[q_point];
6565  for (unsigned int d = 0; d < dim; ++d)
6566  {
6567  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
6568  for (unsigned int e = 1; e < dim; ++e)
6569  new_val += jac[e][d] * grad_in[e];
6570  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
6571  }
6572  }
6573 }
6574 
6575 
6576 
6577 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6578 inline VectorizedArrayType
6580  integrate_value() const
6581 {
6582  return BaseClass::integrate_value()[0];
6583 }
6584 
6585 
6586 
6587 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
6588 
6589 
6590 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6594  const unsigned int dof_no,
6595  const unsigned int first_selected_component,
6596  const unsigned int quad_no_in,
6597  const unsigned int fe_degree,
6598  const unsigned int n_q_points,
6599  const bool is_interior_face,
6600  const unsigned int active_fe_index,
6601  const unsigned int active_quad_index,
6602  const unsigned int face_type)
6604  data_in,
6605  dof_no,
6606  first_selected_component,
6607  quad_no_in,
6608  fe_degree,
6609  n_q_points,
6610  is_interior_face,
6611  active_fe_index,
6612  active_quad_index,
6613  face_type)
6614 {}
6615 
6616 
6617 
6618 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6621  const Mapping<dim> & mapping,
6622  const FiniteElement<dim> &fe,
6623  const Quadrature<1> & quadrature,
6624  const UpdateFlags update_flags,
6625  const unsigned int first_selected_component,
6627  *other)
6629  mapping,
6630  fe,
6631  quadrature,
6632  update_flags,
6633  first_selected_component,
6634  other)
6635 {}
6636 
6637 
6638 
6639 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6643  &other)
6645 {}
6646 
6647 
6648 
6649 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6653  &other)
6654 {
6656  operator=(other);
6657  return *this;
6658 }
6659 
6660 
6661 
6662 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6665  get_gradient(const unsigned int q_point) const
6666 {
6667  return BaseClass::get_gradient(q_point);
6668 }
6669 
6670 
6671 
6672 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6673 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6675  get_divergence(const unsigned int q_point) const
6676 {
6677 # ifdef DEBUG
6678  Assert(this->gradients_quad_initialized == true,
6680 # endif
6681  AssertIndexRange(q_point, this->n_quadrature_points);
6682  Assert(this->jacobian != nullptr,
6684  "update_gradients"));
6685 
6686  VectorizedArrayType divergence;
6687  const std::size_t nqp = this->n_quadrature_points;
6688 
6689  // Cartesian cell
6690  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6691  {
6692  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6693  for (unsigned int d = 1; d < dim; ++d)
6694  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
6695  this->jacobian[0][d][d];
6696  }
6697  // cell with general/constant Jacobian
6698  else
6699  {
6701  this->cell_type == internal::MatrixFreeFunctions::general ?
6702  this->jacobian[q_point] :
6703  this->jacobian[0];
6704  divergence = jac[0][0] * this->gradients_quad[q_point];
6705  for (unsigned int e = 1; e < dim; ++e)
6706  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
6707  for (unsigned int d = 1; d < dim; ++d)
6708  for (unsigned int e = 0; e < dim; ++e)
6709  divergence +=
6710  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
6711  }
6712  return divergence;
6713 }
6714 
6715 
6716 
6717 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6720  get_symmetric_gradient(const unsigned int q_point) const
6721 {
6722  // copy from generic function into dim-specialization function
6723  const auto grad = get_gradient(q_point);
6724  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6725  VectorizedArrayType half = Number(0.5);
6726  for (unsigned int d = 0; d < dim; ++d)
6727  symmetrized[d] = grad[d][d];
6728  switch (dim)
6729  {
6730  case 1:
6731  break;
6732  case 2:
6733  symmetrized[2] = grad[0][1] + grad[1][0];
6734  symmetrized[2] *= half;
6735  break;
6736  case 3:
6737  symmetrized[3] = grad[0][1] + grad[1][0];
6738  symmetrized[3] *= half;
6739  symmetrized[4] = grad[0][2] + grad[2][0];
6740  symmetrized[4] *= half;
6741  symmetrized[5] = grad[1][2] + grad[2][1];
6742  symmetrized[5] *= half;
6743  break;
6744  default:
6745  Assert(false, ExcNotImplemented());
6746  }
6748 }
6749 
6750 
6751 
6752 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6753 inline DEAL_II_ALWAYS_INLINE
6756  const unsigned int q_point) const
6757 {
6758  // copy from generic function into dim-specialization function
6759  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6761  switch (dim)
6762  {
6763  case 1:
6764  Assert(false,
6765  ExcMessage(
6766  "Computing the curl in 1d is not a useful operation"));
6767  break;
6768  case 2:
6769  curl[0] = grad[1][0] - grad[0][1];
6770  break;
6771  case 3:
6772  curl[0] = grad[2][1] - grad[1][2];
6773  curl[1] = grad[0][2] - grad[2][0];
6774  curl[2] = grad[1][0] - grad[0][1];
6775  break;
6776  default:
6777  Assert(false, ExcNotImplemented());
6778  }
6779  return curl;
6780 }
6781 
6782 
6783 
6784 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6787  get_hessian_diagonal(const unsigned int q_point) const
6788 {
6789  return BaseClass::get_hessian_diagonal(q_point);
6790 }
6791 
6792 
6793 
6794 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6797  const unsigned int q_point) const
6798 {
6799 # ifdef DEBUG
6800  Assert(this->hessians_quad_initialized == true,
6802 # endif
6803  AssertIndexRange(q_point, this->n_quadrature_points);
6804  return BaseClass::get_hessian(q_point);
6805 }
6806 
6807 
6808 
6809 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6810 inline DEAL_II_ALWAYS_INLINE void
6813  const unsigned int q_point)
6814 {
6815  BaseClass::submit_gradient(grad_in, q_point);
6816 }
6817 
6818 
6819 
6820 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6821 inline DEAL_II_ALWAYS_INLINE void
6824  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6825  const unsigned int q_point)
6826 {
6827  BaseClass::submit_gradient(grad_in, q_point);
6828 }
6829 
6830 
6831 
6832 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6833 inline DEAL_II_ALWAYS_INLINE void
6835  submit_divergence(const VectorizedArrayType div_in,
6836  const unsigned int q_point)
6837 {
6839  AssertIndexRange(q_point, this->n_quadrature_points);
6840  Assert(this->J_value != nullptr,
6842  "update_gradients"));
6843  Assert(this->jacobian != nullptr,
6845  "update_gradients"));
6846 # ifdef DEBUG
6847  this->gradients_quad_submitted = true;
6848 # endif
6849 
6850  const std::size_t nqp = this->n_quadrature_points;
6851  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6852  {
6853  const VectorizedArrayType fac =
6854  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6855  for (unsigned int d = 0; d < dim; ++d)
6856  {
6857  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6858  (fac * this->jacobian[0][d][d]);
6859  for (unsigned int e = d + 1; e < dim; ++e)
6860  {
6861  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6862  VectorizedArrayType();
6863  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6864  VectorizedArrayType();
6865  }
6866  }
6867  }
6868  else
6869  {
6871  this->cell_type == internal::MatrixFreeFunctions::general ?
6872  this->jacobian[q_point] :
6873  this->jacobian[0];
6874  const VectorizedArrayType fac =
6875  (this->cell_type == internal::MatrixFreeFunctions::general ?
6876  this->J_value[q_point] :
6877  this->J_value[0] * this->quadrature_weights[q_point]) *
6878  div_in;
6879  for (unsigned int d = 0; d < dim; ++d)
6880  {
6881  for (unsigned int e = 0; e < dim; ++e)
6882  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6883  jac[d][e] * fac;
6884  }
6885  }
6886 }
6887 
6888 
6889 
6890 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6891 inline DEAL_II_ALWAYS_INLINE void
6895  const unsigned int q_point)
6896 {
6897  // could have used base class operator, but that involves some overhead
6898  // which is inefficient. it is nice to have the symmetric tensor because
6899  // that saves some operations
6901  AssertIndexRange(q_point, this->n_quadrature_points);
6902  Assert(this->J_value != nullptr,
6904  "update_gradients"));
6905  Assert(this->jacobian != nullptr,
6907  "update_gradients"));
6908 # ifdef DEBUG
6909  this->gradients_quad_submitted = true;
6910 # endif
6911 
6912  const std::size_t nqp = this->n_quadrature_points;
6913  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6914  {
6915  const VectorizedArrayType JxW =
6916  this->J_value[0] * this->quadrature_weights[q_point];
6917  for (unsigned int d = 0; d < dim; ++d)
6918  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6919  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6920  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6921  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6922  {
6923  const VectorizedArrayType value =
6924  sym_grad.access_raw_entry(counter) * JxW;
6925  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6926  value * this->jacobian[0][d][d];
6927  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6928  value * this->jacobian[0][e][e];
6929  }
6930  }
6931  // general/affine cell type
6932  else
6933  {
6934  const VectorizedArrayType JxW =
6935  this->cell_type == internal::MatrixFreeFunctions::general ?
6936  this->J_value[q_point] :
6937  this->J_value[0] * this->quadrature_weights[q_point];
6939  this->cell_type == internal::MatrixFreeFunctions::general ?
6940  this->jacobian[q_point] :
6941  this->jacobian[0];
6942  VectorizedArrayType weighted[dim][dim];
6943  for (unsigned int i = 0; i < dim; ++i)
6944  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6945  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6946  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6947  {
6948  const VectorizedArrayType value =
6949  sym_grad.access_raw_entry(counter) * JxW;
6950  weighted[i][j] = value;
6951  weighted[j][i] = value;
6952  }
6953  for (unsigned int comp = 0; comp < dim; ++comp)
6954  for (unsigned int d = 0; d < dim; ++d)
6955  {
6956  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6957  for (unsigned int e = 1; e < dim; ++e)
6958  new_val += jac[e][d] * weighted[comp][e];
6959  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
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  switch (dim)
6974  {
6975  case 1:
6976  Assert(false,
6977  ExcMessage(
6978  "Testing by the curl in 1d is not a useful operation"));
6979  break;
6980  case 2:
6981  grad[1][0] = curl[0];
6982  grad[0][1] = -curl[0];
6983  break;
6984  case 3:
6985  grad[2][1] = curl[0];
6986  grad[1][2] = -curl[0];
6987  grad[0][2] = curl[1];
6988  grad[2][0] = -curl[1];
6989  grad[1][0] = curl[2];
6990  grad[0][1] = -curl[2];
6991  break;
6992  default:
6993  Assert(false, ExcNotImplemented());
6994  }
6995  submit_gradient(grad, q_point);
6996 }
6997 
6998 
6999 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
7000 
7001 
7002 template <typename Number, bool is_face, typename VectorizedArrayType>
7005  const unsigned int dof_no,
7006  const unsigned int first_selected_component,
7007  const unsigned int quad_no_in,
7008  const unsigned int fe_degree,
7009  const unsigned int n_q_points,
7010  const bool is_interior_face,
7011  const unsigned int active_fe_index,
7012  const unsigned int active_quad_index,
7013  const unsigned int face_type)
7015  data_in,
7016  dof_no,
7017  first_selected_component,
7018  quad_no_in,
7019  fe_degree,
7020  n_q_points,
7021  is_interior_face,
7022  active_fe_index,
7023  active_quad_index,
7024  face_type)
7025 {}
7026 
7027 
7028 
7029 template <typename Number, bool is_face, typename VectorizedArrayType>
7032  const Mapping<1> & mapping,
7033  const FiniteElement<1> &fe,
7034  const Quadrature<1> & quadrature,
7035  const UpdateFlags update_flags,
7036  const unsigned int first_selected_component,
7039  mapping,
7040  fe,
7041  quadrature,
7042  update_flags,
7043  first_selected_component,
7044  other)
7045 {}
7046 
7047 
7048 
7049 template <typename Number, bool is_face, typename VectorizedArrayType>
7054 {}
7055 
7056 
7057 
7058 template <typename Number, bool is_face, typename VectorizedArrayType>
7062 {
7064  other);
7065  return *this;
7066 }
7067 
7068 
7069 
7070 template <typename Number, bool is_face, typename VectorizedArrayType>
7071 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7073  const unsigned int dof) const
7074 {
7075  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7076  return this->values_dofs[0][dof];
7077 }
7078 
7079 
7080 
7081 template <typename Number, bool is_face, typename VectorizedArrayType>
7082 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7084  const unsigned int q_point) const
7085 {
7086 # ifdef DEBUG
7087  Assert(this->values_quad_initialized == true,
7089 # endif
7090  AssertIndexRange(q_point, this->n_quadrature_points);
7091  return this->values_quad[q_point];
7092 }
7093 
7094 
7095 
7096 template <typename Number, bool is_face, typename VectorizedArrayType>
7099  const unsigned int q_point) const
7100 {
7101  // could use the base class gradient, but that involves too many inefficient
7102  // initialization operations on tensors
7103 
7104 # ifdef DEBUG
7105  Assert(this->gradients_quad_initialized == true,
7107 # endif
7108  AssertIndexRange(q_point, this->n_quadrature_points);
7109 
7111  this->cell_type == internal::MatrixFreeFunctions::general ?
7112  this->jacobian[q_point] :
7113  this->jacobian[0];
7114 
7116  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
7117 
7118  return grad_out;
7119 }
7120 
7121 
7122 
7123 template <typename Number, bool is_face, typename VectorizedArrayType>
7124 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7126  get_normal_derivative(const unsigned int q_point) const
7127 {
7128  return BaseClass::get_normal_derivative(q_point)[0];
7129 }
7130 
7131 
7132 
7133 template <typename Number, bool is_face, typename VectorizedArrayType>
7136  const unsigned int q_point) const
7137 {
7138  return BaseClass::get_hessian(q_point)[0];
7139 }
7140 
7141 
7142 
7143 template <typename Number, bool is_face, typename VectorizedArrayType>
7146  get_hessian_diagonal(const unsigned int q_point) const
7147 {
7148  return BaseClass::get_hessian_diagonal(q_point)[0];
7149 }
7150 
7151 
7152 
7153 template <typename Number, bool is_face, typename VectorizedArrayType>
7154 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
7156  const unsigned int q_point) const
7157 {
7158  return BaseClass::get_laplacian(q_point)[0];
7159 }
7160 
7161 
7162 
7163 template <typename Number, bool is_face, typename VectorizedArrayType>
7166  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
7167 {
7168 # ifdef DEBUG
7169  this->dof_values_initialized = true;
7170  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
7171 # endif
7172  this->values_dofs[0][dof] = val_in;
7173 }
7174 
7175 
7176 
7177 template <typename Number, bool is_face, typename VectorizedArrayType>
7178 inline DEAL_II_ALWAYS_INLINE void
7180  const VectorizedArrayType val_in,
7181  const unsigned int q_point)
7182 {
7184  AssertIndexRange(q_point, this->n_quadrature_points);
7185 # ifdef DEBUG
7186  this->values_quad_submitted = true;
7187 # endif
7188 
7189  if (this->cell_type == internal::MatrixFreeFunctions::general)
7190  {
7191  const VectorizedArrayType JxW = this->J_value[q_point];
7192  this->values_quad[q_point] = val_in * JxW;
7193  }
7194  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
7195  {
7196  const VectorizedArrayType JxW =
7197  this->J_value[0] * this->quadrature_weights[q_point];
7198  this->values_quad[q_point] = val_in * JxW;
7199  }
7200 }
7201 
7202 
7203 
7204 template <typename Number, bool is_face, typename VectorizedArrayType>
7205 inline DEAL_II_ALWAYS_INLINE void
7207  const Tensor<1, 1, VectorizedArrayType> val_in,
7208  const unsigned int q_point)
7209 {
7210  submit_value(val_in[0], q_point);
7211 }
7212 
7213 
7214 
7215 template <typename Number, bool is_face, typename VectorizedArrayType>
7216 inline DEAL_II_ALWAYS_INLINE void
7218  const Tensor<1, 1, VectorizedArrayType> grad_in,
7219  const unsigned int q_point)
7220 {
7221  submit_gradient(grad_in[0], q_point);
7222 }
7223 
7224 
7225 
7226 template <typename Number, bool is_face, typename VectorizedArrayType>
7227 inline DEAL_II_ALWAYS_INLINE void
7229  const VectorizedArrayType grad_in,
7230  const unsigned int q_point)
7231 {
7233  AssertIndexRange(q_point, this->n_quadrature_points);
7234 # ifdef DEBUG
7235  this->gradients_quad_submitted = true;
7236 # endif
7237 
7239  this->cell_type == internal::MatrixFreeFunctions::general ?
7240  this->jacobian[q_point] :
7241  this->jacobian[0];
7242  const VectorizedArrayType JxW =
7243  this->cell_type == internal::MatrixFreeFunctions::general ?
7244  this->J_value[q_point] :
7245  this->J_value[0] * this->quadrature_weights[q_point];
7246 
7247  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7248 }
7249 
7250 
7251 
7252 template <typename Number, bool is_face, typename VectorizedArrayType>
7253 inline DEAL_II_ALWAYS_INLINE void
7255  submit_normal_derivative(const VectorizedArrayType grad_in,
7256  const unsigned int q_point)
7257 {
7259  grad[0] = grad_in;
7260  BaseClass::submit_normal_derivative(grad, q_point);
7261 }
7262 
7263 
7264 
7265 template <typename Number, bool is_face, typename VectorizedArrayType>
7266 inline DEAL_II_ALWAYS_INLINE void
7269  const unsigned int q_point)
7270 {
7271  BaseClass::submit_normal_derivative(grad_in, q_point);
7272 }
7273 
7274 
7275 
7276 template <typename Number, bool is_face, typename VectorizedArrayType>
7277 inline VectorizedArrayType
7279  integrate_value() const
7280 {
7281  return BaseClass::integrate_value()[0];
7282 }
7283 
7284 
7285 
7286 /*-------------------------- FEEvaluation -----------------------------------*/
7287 
7288 
7289 template <int dim,
7290  int fe_degree,
7291  int n_q_points_1d,
7292  int n_components_,
7293  typename Number,
7294  typename VectorizedArrayType>
7295 inline FEEvaluation<dim,
7296  fe_degree,
7297  n_q_points_1d,
7298  n_components_,
7299  Number,
7300  VectorizedArrayType>::
7301  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
7302  const unsigned int fe_no,
7303  const unsigned int quad_no,
7304  const unsigned int first_selected_component,
7305  const unsigned int active_fe_index,
7306  const unsigned int active_quad_index)
7307  : BaseClass(data_in,
7308  fe_no,
7309  first_selected_component,
7310  quad_no,
7311  fe_degree,
7312  static_n_q_points,
7313  true /*note: this is not a face*/,
7314  active_fe_index,
7315  active_quad_index)
7316  , dofs_per_component(this->data->dofs_per_component_on_cell)
7317  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7318  , n_q_points(this->data->n_q_points)
7319 {
7320  check_template_arguments(fe_no, 0);
7321 }
7322 
7323 
7324 
7325 template <int dim,
7326  int fe_degree,
7327  int n_q_points_1d,
7328  int n_components_,
7329  typename Number,
7330  typename VectorizedArrayType>
7331 inline FEEvaluation<dim,
7332  fe_degree,
7333  n_q_points_1d,
7334  n_components_,
7335  Number,
7336  VectorizedArrayType>::
7337  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7338  const std::pair<unsigned int, unsigned int> & range,
7339  const unsigned int dof_no,
7340  const unsigned int quad_no,
7341  const unsigned int first_selected_component)
7342  : FEEvaluation(matrix_free,
7343  dof_no,
7344  quad_no,
7345  first_selected_component,
7346  matrix_free.get_cell_active_fe_index(range))
7347 {}
7348 
7349 
7350 
7351 template <int dim,
7352  int fe_degree,
7353  int n_q_points_1d,
7354  int n_components_,
7355  typename Number,
7356  typename VectorizedArrayType>
7357 inline FEEvaluation<dim,
7358  fe_degree,
7359  n_q_points_1d,
7360  n_components_,
7361  Number,
7362  VectorizedArrayType>::
7363  FEEvaluation(const Mapping<dim> & mapping,
7364  const FiniteElement<dim> &fe,
7365  const Quadrature<1> & quadrature,
7366  const UpdateFlags update_flags,
7367  const unsigned int first_selected_component)
7368  : BaseClass(mapping,
7369  fe,
7370  quadrature,
7371  update_flags,
7372  first_selected_component,
7373  nullptr)
7374  , dofs_per_component(this->data->dofs_per_component_on_cell)
7375  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7376  , n_q_points(this->data->n_q_points)
7377 {
7378  check_template_arguments(numbers::invalid_unsigned_int, 0);
7379 }
7380 
7381 
7382 
7383 template <int dim,
7384  int fe_degree,
7385  int n_q_points_1d,
7386  int n_components_,
7387  typename Number,
7388  typename VectorizedArrayType>
7389 inline FEEvaluation<dim,
7390  fe_degree,
7391  n_q_points_1d,
7392  n_components_,
7393  Number,
7394  VectorizedArrayType>::
7395  FEEvaluation(const FiniteElement<dim> &fe,
7396  const Quadrature<1> & quadrature,
7397  const UpdateFlags update_flags,
7398  const unsigned int first_selected_component)
7399  : BaseClass(StaticMappingQ1<dim>::mapping,
7400  fe,
7401  quadrature,
7402  update_flags,
7403  first_selected_component,
7404  nullptr)
7405  , dofs_per_component(this->data->dofs_per_component_on_cell)
7406  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7407  , n_q_points(this->data->n_q_points)
7408 {
7409  check_template_arguments(numbers::invalid_unsigned_int, 0);
7410 }
7411 
7412 
7413 
7414 template <int dim,
7415  int fe_degree,
7416  int n_q_points_1d,
7417  int n_components_,
7418  typename Number,
7419  typename VectorizedArrayType>
7420 inline FEEvaluation<dim,
7421  fe_degree,
7422  n_q_points_1d,
7423  n_components_,
7424  Number,
7425  VectorizedArrayType>::
7426  FEEvaluation(
7427  const FiniteElement<dim> & fe,
7429  const unsigned int first_selected_component)
7430  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7431  fe,
7432  other.mapped_geometry->get_quadrature(),
7433  other.mapped_geometry->get_fe_values().get_update_flags(),
7434  first_selected_component,
7435  &other)
7436  , dofs_per_component(this->data->dofs_per_component_on_cell)
7437  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7438  , n_q_points(this->data->n_q_points)
7439 {
7440  check_template_arguments(numbers::invalid_unsigned_int, 0);
7441 }
7442 
7443 
7444 
7445 template <int dim,
7446  int fe_degree,
7447  int n_q_points_1d,
7448  int n_components_,
7449  typename Number,
7450  typename VectorizedArrayType>
7451 inline FEEvaluation<dim,
7452  fe_degree,
7453  n_q_points_1d,
7454  n_components_,
7455  Number,
7456  VectorizedArrayType>::FEEvaluation(const FEEvaluation
7457  &other)
7458  : BaseClass(other)
7459  , dofs_per_component(this->data->dofs_per_component_on_cell)
7460  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7461  , n_q_points(this->data->n_q_points)
7462 {
7463  check_template_arguments(numbers::invalid_unsigned_int, 0);
7464 }
7465 
7466 
7467 
7468 template <int dim,
7469  int fe_degree,
7470  int n_q_points_1d,
7471  int n_components_,
7472  typename Number,
7473  typename VectorizedArrayType>
7474 inline FEEvaluation<dim,
7475  fe_degree,
7476  n_q_points_1d,
7477  n_components_,
7478  Number,
7479  VectorizedArrayType> &
7480 FEEvaluation<dim,
7481  fe_degree,
7482  n_q_points_1d,
7483  n_components_,
7484  Number,
7485  VectorizedArrayType>::operator=(const FEEvaluation &other)
7486 {
7487  BaseClass::operator=(other);
7488  check_template_arguments(numbers::invalid_unsigned_int, 0);
7489  return *this;
7490 }
7491 
7492 
7493 
7494 template <int dim,
7495  int fe_degree,
7496  int n_q_points_1d,
7497  int n_components_,
7498  typename Number,
7499  typename VectorizedArrayType>
7500 inline void
7501 FEEvaluation<dim,
7502  fe_degree,
7503  n_q_points_1d,
7504  n_components_,
7505  Number,
7506  VectorizedArrayType>::
7507  check_template_arguments(const unsigned int dof_no,
7508  const unsigned int first_selected_component)
7509 {
7510  (void)dof_no;
7511  (void)first_selected_component;
7512 
7513 # ifdef DEBUG
7514  // print error message when the dimensions do not match. Propose a possible
7515  // fix
7516  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7517  static_cast<unsigned int>(fe_degree) !=
7518  this->data->data.front().fe_degree) ||
7519  n_q_points != this->n_quadrature_points)
7520  {
7521  std::string message =
7522  "-------------------------------------------------------\n";
7523  message += "Illegal arguments in constructor/wrong template arguments!\n";
7524  message += " Called --> FEEvaluation<dim,";
7525  message += Utilities::int_to_string(fe_degree) + ",";
7526  message += Utilities::int_to_string(n_q_points_1d);
7527  message += "," + Utilities::int_to_string(n_components);
7528  message += ",Number>(data";
7529  if (first_selected_component != numbers::invalid_unsigned_int)
7530  {
7531  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7532  message += Utilities::int_to_string(this->quad_no) + ", ";
7533  message += Utilities::int_to_string(first_selected_component);
7534  }
7535  message += ")\n";
7536 
7537  // check whether some other vector component has the correct number of
7538  // points
7539  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7540  proposed_fe_comp = numbers::invalid_unsigned_int,
7541  proposed_quad_comp = numbers::invalid_unsigned_int;
7542  if (dof_no != numbers::invalid_unsigned_int)
7543  {
7544  if (static_cast<unsigned int>(fe_degree) ==
7545  this->data->data.front().fe_degree)
7546  {
7547  proposed_dof_comp = dof_no;
7548  proposed_fe_comp = first_selected_component;
7549  }
7550  else
7551  for (unsigned int no = 0; no < this->matrix_info->n_components();
7552  ++no)
7553  for (unsigned int nf = 0;
7554  nf < this->matrix_info->n_base_elements(no);
7555  ++nf)
7556  if (this->matrix_info
7557  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7558  .data.front()
7559  .fe_degree == static_cast<unsigned int>(fe_degree))
7560  {
7561  proposed_dof_comp = no;
7562  proposed_fe_comp = nf;
7563  break;
7564  }
7565  if (n_q_points ==
7566  this->mapping_data->descriptor[this->active_quad_index]
7567  .n_q_points)
7568  proposed_quad_comp = this->quad_no;
7569  else
7570  for (unsigned int no = 0;
7571  no < this->matrix_info->get_mapping_info().cell_data.size();
7572  ++no)
7573  if (this->matrix_info->get_mapping_info()
7574  .cell_data[no]
7575  .descriptor[this->active_quad_index]
7576  .n_q_points == n_q_points)
7577  {
7578  proposed_quad_comp = no;
7579  break;
7580  }
7581  }
7582  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7583  proposed_quad_comp != numbers::invalid_unsigned_int)
7584  {
7585  if (proposed_dof_comp != first_selected_component)
7586  message += "Wrong vector component selection:\n";
7587  else
7588  message += "Wrong quadrature formula selection:\n";
7589  message += " Did you mean FEEvaluation<dim,";
7590  message += Utilities::int_to_string(fe_degree) + ",";
7591  message += Utilities::int_to_string(n_q_points_1d);
7592  message += "," + Utilities::int_to_string(n_components);
7593  message += ",Number>(data";
7594  if (dof_no != numbers::invalid_unsigned_int)
7595  {
7596  message +=
7597  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7598  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7599  message += Utilities::int_to_string(proposed_fe_comp);
7600  }
7601  message += ")?\n";
7602  std::string correct_pos;
7603  if (proposed_dof_comp != dof_no)
7604  correct_pos = " ^ ";
7605  else
7606  correct_pos = " ";
7607  if (proposed_quad_comp != this->quad_no)
7608  correct_pos += " ^ ";
7609  else
7610  correct_pos += " ";
7611  if (proposed_fe_comp != first_selected_component)
7612  correct_pos += " ^\n";
7613  else
7614  correct_pos += " \n";
7615  message += " " +
7616  correct_pos;
7617  }
7618  // ok, did not find the numbers specified by the template arguments in
7619  // the given list. Suggest correct template arguments
7620  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7621  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7622  message += "Wrong template arguments:\n";
7623  message += " Did you mean FEEvaluation<dim,";
7624  message +=
7625  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7626  message += Utilities::int_to_string(proposed_n_q_points_1d);
7627  message += "," + Utilities::int_to_string(n_components);
7628  message += ",Number>(data";
7629  if (dof_no != numbers::invalid_unsigned_int)
7630  {
7631  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7632  message += Utilities::int_to_string(this->quad_no);
7633  message += ", " + Utilities::int_to_string(first_selected_component);
7634  }
7635  message += ")?\n";
7636  std::string correct_pos;
7637  if (this->data->data.front().fe_degree !=
7638  static_cast<unsigned int>(fe_degree))
7639  correct_pos = " ^";
7640  else
7641  correct_pos = " ";
7642  if (proposed_n_q_points_1d != n_q_points_1d)
7643  correct_pos += " ^\n";
7644  else
7645  correct_pos += " \n";
7646  message += " " + correct_pos;
7647 
7648  Assert(static_cast<unsigned int>(fe_degree) ==
7649  this->data->data.front().fe_degree &&
7650  n_q_points == this->n_quadrature_points,
7651  ExcMessage(message));
7652  }
7653  if (dof_no != numbers::invalid_unsigned_int)
7655  n_q_points,
7656  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7657 # endif
7658 }
7659 
7660 
7661 
7662 template <int dim,
7663  int fe_degree,
7664  int n_q_points_1d,
7665  int n_components_,
7666  typename Number,
7667  typename VectorizedArrayType>
7668 inline void
7669 FEEvaluation<dim,
7670  fe_degree,
7671  n_q_points_1d,
7672  n_components_,
7673  Number,
7674  VectorizedArrayType>::reinit(const unsigned int cell_index)
7675 {
7676  Assert(this->mapped_geometry == nullptr,
7677  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7678  " Integer indexing is not possible"));
7679  if (this->mapped_geometry != nullptr)
7680  return;
7681 
7682  Assert(this->dof_info != nullptr, ExcNotInitialized());
7683  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7684  this->cell = cell_index;
7685  this->cell_type =
7686  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
7687 
7688  const unsigned int offsets =
7689  this->mapping_data->data_index_offsets[cell_index];
7690  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7691  this->J_value = &this->mapping_data->JxW_values[offsets];
7692 
7693 # ifdef DEBUG
7694  this->dof_values_initialized = false;
7695  this->values_quad_initialized = false;
7696  this->gradients_quad_initialized = false;
7697  this->hessians_quad_initialized = false;
7698 # endif
7699 }
7700 
7701 
7702 
7703 template <int dim,
7704  int fe_degree,
7705  int n_q_points_1d,
7706  int n_components_,
7707  typename Number,
7708  typename VectorizedArrayType>
7709 template <bool level_dof_access>
7710 inline void
7711 FEEvaluation<dim,
7712  fe_degree,
7713  n_q_points_1d,
7714  n_components_,
7715  Number,
7716  VectorizedArrayType>::
7718 {
7719  Assert(this->matrix_info == nullptr,
7720  ExcMessage("Cannot use initialization from cell iterator if "
7721  "initialized from MatrixFree object. Use variant for "
7722  "on the fly computation with arguments as for FEValues "
7723  "instead"));
7724  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7725  this->mapped_geometry->reinit(
7726  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7727  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7728  if (level_dof_access)
7729  cell->get_mg_dof_indices(this->local_dof_indices);
7730  else
7731  cell->get_dof_indices(this->local_dof_indices);
7732 }
7733 
7734 
7735 
7736 template <int dim,
7737  int fe_degree,
7738  int n_q_points_1d,
7739  int n_components_,
7740  typename Number,
7741  typename VectorizedArrayType>
7742 inline void
7743 FEEvaluation<dim,
7744  fe_degree,
7745  n_q_points_1d,
7746  n_components_,
7747  Number,
7748  VectorizedArrayType>::
7749  reinit(const typename Triangulation<dim>::cell_iterator &cell)
7750 {
7751  Assert(this->matrix_info == 0,
7752  ExcMessage("Cannot use initialization from cell iterator if "
7753  "initialized from MatrixFree object. Use variant for "
7754  "on the fly computation with arguments as for FEValues "
7755  "instead"));
7756  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7757  this->mapped_geometry->reinit(cell);
7758 }
7759 
7760 
7761 
7762 template <int dim,
7763  int fe_degree,
7764  int n_q_points_1d,
7765  int n_components_,
7766  typename Number,
7767  typename VectorizedArrayType>
7769 FEEvaluation<dim,
7770  fe_degree,
7771  n_q_points_1d,
7772  n_components_,
7773  Number,
7774  VectorizedArrayType>::quadrature_point(const unsigned int q) const
7775 {
7776  if (this->matrix_info == nullptr)
7777  {
7778  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7781  "update_quadrature_points"));
7782  }
7783  else
7784  {
7785  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
7787  "update_quadrature_points"));
7788  }
7789 
7790  AssertIndexRange(q, n_q_points);
7791 
7793  &this->mapping_data->quadrature_points
7794  [this->mapping_data->quadrature_point_offsets[this->cell]];
7795 
7796  // Cartesian/affine mesh: only first vertex of cell is stored, we must
7797  // compute it through the Jacobian (which is stored in non-inverted and
7798  // non-transposed form as index '1' in the jacobian field)
7799  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
7800  {
7801  Assert(this->jacobian != nullptr, ExcNotInitialized());
7802  Point<dim, VectorizedArrayType> point = quadrature_points[0];
7803 
7804  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
7805  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
7806  for (unsigned int d = 0; d < dim; ++d)
7807  point[d] += jac[d][d] * static_cast<Number>(
7808  this->descriptor->quadrature.point(q)[d]);
7809  else
7810  for (unsigned int d = 0; d < dim; ++d)
7811  for (unsigned int e = 0; e < dim; ++e)
7812  point[d] += jac[d][e] * static_cast<Number>(
7813  this->descriptor->quadrature.point(q)[e]);
7814  return point;
7815  }
7816  else
7817  return quadrature_points[q];
7818 }
7819 
7820 
7821 
7822 template <int dim,
7823  int fe_degree,
7824  int n_q_points_1d,
7825  int n_components_,
7826  typename Number,
7827  typename VectorizedArrayType>
7828 inline void
7829 FEEvaluation<dim,
7830  fe_degree,
7831  n_q_points_1d,
7832  n_components_,
7833  Number,
7834  VectorizedArrayType>::evaluate(const bool evaluate_values,
7835  const bool evaluate_gradients,
7836  const bool evaluate_hessians)
7837 {
7838 # ifdef DEBUG
7839  Assert(this->dof_values_initialized == true,
7841 # endif
7842  evaluate(this->values_dofs[0],
7843  evaluate_values,
7844  evaluate_gradients,
7845  evaluate_hessians);
7846 }
7847 
7848 
7849 template <int dim,
7850  int fe_degree,
7851  int n_q_points_1d,
7852  int n_components_,
7853  typename Number,
7854  typename VectorizedArrayType>
7855 inline void
7856 FEEvaluation<dim,
7857  fe_degree,
7858  n_q_points_1d,
7859  n_components_,
7860  Number,
7861  VectorizedArrayType>::
7862  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7863 {
7864 # ifdef DEBUG
7865  Assert(this->dof_values_initialized == true,
7867 # endif
7868  evaluate(this->values_dofs[0], evaluation_flags);
7869 }
7870 
7871 
7872 
7873 template <int dim,
7874  int fe_degree,
7875  int n_q_points_1d,
7876  int n_components_,
7877  typename Number,
7878  typename VectorizedArrayType>
7879 inline void
7880 FEEvaluation<dim,
7881  fe_degree,
7882  n_q_points_1d,
7883  n_components_,
7884  Number,
7885  VectorizedArrayType>::evaluate(const VectorizedArrayType
7886  * values_array,
7887  const bool evaluate_values,
7888  const bool evaluate_gradients,
7889  const bool evaluate_hessians)
7890 {
7892  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7893  ((evaluate_gradients) ? EvaluationFlags::gradients :
7894  EvaluationFlags::nothing) |
7895  ((evaluate_hessians) ? EvaluationFlags::hessians :
7896  EvaluationFlags::nothing);
7897 
7898  evaluate(values_array, flag);
7899 }
7900 
7901 
7902 
7903 template <int dim,
7904  int fe_degree,
7905  int n_q_points_1d,
7906  int n_components_,
7907  typename Number,
7908  typename VectorizedArrayType>
7909 inline void
7910 FEEvaluation<dim,
7911  fe_degree,
7912  n_q_points_1d,
7913  n_components_,
7914  Number,
7915  VectorizedArrayType>::
7916  evaluate(const VectorizedArrayType * values_array,
7917  const EvaluationFlags::EvaluationFlags evaluation_flags)
7918 {
7919  if (fe_degree > -1)
7921  evaluate(n_components,
7922  evaluation_flags,
7923  *this->data,
7924  const_cast<VectorizedArrayType *>(values_array),
7925  this->values_quad,
7926  this->gradients_quad,
7927  this->hessians_quad,
7928  this->scratch_data);
7929  else
7931  n_components,
7932  evaluation_flags,
7933  *this->data,
7934  const_cast<VectorizedArrayType *>(values_array),
7935  this->values_quad,
7936  this->gradients_quad,
7937  this->hessians_quad,
7938  this->scratch_data);
7939 
7940 # ifdef DEBUG
7941  if (evaluation_flags & EvaluationFlags::values)
7942  this->values_quad_initialized = true;
7943  if (evaluation_flags & EvaluationFlags::gradients)
7944  this->gradients_quad_initialized = true;
7945  if (evaluation_flags & EvaluationFlags::hessians)
7946  this->hessians_quad_initialized = true;
7947 # endif
7948 }
7949 
7950 
7951 
7952 template <int dim,
7953  int fe_degree,
7954  int n_q_points_1d,
7955  int n_components_,
7956  typename Number,
7957  typename VectorizedArrayType>
7958 template <typename VectorType>
7959 inline void
7960 FEEvaluation<
7961  dim,
7962  fe_degree,
7963  n_q_points_1d,
7964  n_components_,
7965  Number,
7966  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7967  const bool evaluate_values,
7968  const bool evaluate_gradients,
7969  const bool evaluate_hessians)
7970 {
7972  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7973  ((evaluate_gradients) ? EvaluationFlags::gradients :
7974  EvaluationFlags::nothing) |
7975  ((evaluate_hessians) ? EvaluationFlags::hessians :
7976  EvaluationFlags::nothing);
7977 
7978  gather_evaluate(input_vector, flag);
7979 }
7980 
7981 
7982 namespace internal
7983 {
7987  template <typename Number,
7988  typename VectorizedArrayType,
7989  typename VectorType,
7990  typename T,
7991  typename std::enable_if<
7993  std::is_same<decltype(std::declval<VectorType>().begin()),
7994  Number *>::value,
7995  VectorType>::type * = nullptr>
7996  bool
7997  try_gather_evaluate_inplace(
7998  T phi,
7999  const VectorType & input_vector,
8000  const unsigned int cell,
8001  const unsigned int active_fe_index,
8002  const unsigned int first_selected_component,
8004  const EvaluationFlags::EvaluationFlags evaluation_flag)
8005  {
8006  // If the index storage is interleaved and contiguous and the vector storage
8007  // has the correct alignment, we can directly pass the pointer into the
8008  // vector to the evaluate() call, without reading the vector entries into a
8009  // separate data field. This saves some operations.
8010  if (std::is_same<typename VectorType::value_type, Number>::value &&
8011  dof_info->index_storage_variants
8014  interleaved_contiguous &&
8015  reinterpret_cast<std::size_t>(
8016  input_vector.begin() +
8017  dof_info->dof_indices_contiguous
8019  [cell * VectorizedArrayType::size()]) %
8020  sizeof(VectorizedArrayType) ==
8021  0)
8022  {
8023  const VectorizedArrayType *vec_values =
8024  reinterpret_cast<const VectorizedArrayType *>(
8025  input_vector.begin() +
8026  dof_info->dof_indices_contiguous
8028  [cell * VectorizedArrayType::size()] +
8029  dof_info->component_dof_indices_offset[active_fe_index]
8030  [first_selected_component] *
8031  VectorizedArrayType::size());
8032 
8033  phi->evaluate(vec_values, evaluation_flag);
8034 
8035  return true;
8036  }
8037 
8038  return false;
8039  }
8040 
8044  template <typename Number,
8045  typename VectorizedArrayType,
8046  typename VectorType,
8047  typename T,
8048  typename std::enable_if<
8050  !std::is_same<decltype(std::declval<VectorType>().begin()),
8051  Number *>::value,
8052  VectorType>::type * = nullptr>
8053  bool
8054  try_gather_evaluate_inplace(T,
8055  const VectorType &,
8056  const unsigned int,
8057  const unsigned int,
8058  const unsigned int,
8061  {
8062  return false;
8063  }
8064 
8068  template <int dim,
8069  int fe_degree,
8070  int n_q_points_1d,
8071  typename Number,
8072  typename VectorizedArrayType,
8073  typename VectorType,
8074  typename std::enable_if<
8076  std::is_same<decltype(std::declval<VectorType>().begin()),
8077  Number *>::value,
8078  VectorType>::type * = nullptr>
8079  bool
8080  try_integrate_scatter_inplace(
8081  VectorType & destination,
8082  const unsigned int cell,
8083  const unsigned int n_components,
8084  const unsigned int active_fe_index,
8085  const unsigned int first_selected_component,
8087  VectorizedArrayType * values_quad,
8088  VectorizedArrayType * gradients_quad,
8089  VectorizedArrayType * scratch_data,
8091  const EvaluationFlags::EvaluationFlags integration_flag)
8092  {
8093  // If the index storage is interleaved and contiguous and the vector storage
8094  // has the correct alignment, we can directly pass the pointer into the
8095  // vector to the integrate() call, without writing temporary results into a
8096  // separate data field that will later be added into the vector. This saves
8097  // some operations.
8098  if (std::is_same<typename VectorType::value_type, Number>::value &&
8099  dof_info->index_storage_variants
8102  interleaved_contiguous &&
8103  reinterpret_cast<std::size_t>(
8104  destination.begin() +
8105  dof_info->dof_indices_contiguous
8107  [cell * VectorizedArrayType::size()]) %
8108  sizeof(VectorizedArrayType) ==
8109  0)
8110  {
8111  VectorizedArrayType *vec_values =
8112  reinterpret_cast<VectorizedArrayType *>(
8113  destination.begin() +
8114  dof_info->dof_indices_contiguous
8116  [cell * VectorizedArrayType::size()] +
8117  dof_info->component_dof_indices_offset[active_fe_index]
8118  [first_selected_component] *
8119  VectorizedArrayType::size());
8120  if (fe_degree > -1)
8122  integrate(n_components,
8123  integration_flag,
8124  *data,
8125  vec_values,
8126  values_quad,
8127  gradients_quad,
8128  scratch_data,
8129  true);
8130  else
8132  n_components,
8133  integration_flag,
8134  *data,
8135  vec_values,
8136  values_quad,
8137  gradients_quad,
8138  scratch_data,
8139  true);
8140 
8141  return true;
8142  }
8143 
8144  return false;
8145  }
8146 
8150  template <int dim,
8151  int fe_degree,
8152  int n_q_points_1d,
8153  typename Number,
8154  typename VectorizedArrayType,
8155  typename VectorType,
8156  typename std::enable_if<
8158  !std::is_same<decltype(std::declval<VectorType>().begin()),
8159  Number *>::value,
8160  VectorType>::type * = nullptr>
8161  bool
8162  try_integrate_scatter_inplace(
8163  VectorType &,
8164  const unsigned int,
8165  const unsigned int,
8166  const unsigned int,
8167  const unsigned int,
8169  const VectorizedArrayType *,
8170  const VectorizedArrayType *,
8171  const VectorizedArrayType *,
8174  {
8175  return false;
8176  }
8177 } // namespace internal
8178 
8179 
8180 
8181 template <int dim,
8182  int fe_degree,
8183  int n_q_points_1d,
8184  int n_components_,
8185  typename Number,
8186  typename VectorizedArrayType>
8187 template <typename VectorType>
8188 inline void
8189 FEEvaluation<dim,
8190  fe_degree,
8191  n_q_points_1d,
8192  n_components_,
8193  Number,
8194  VectorizedArrayType>::
8195  gather_evaluate(const VectorType & input_vector,
8196  const EvaluationFlags::EvaluationFlags evaluation_flag)
8197 {
8198  if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
8199  this,
8200  input_vector,
8201  this->cell,
8202  this->active_fe_index,
8203  this->first_selected_component,
8204  this->dof_info,
8205  evaluation_flag) == false)
8206  {
8207  this->read_dof_values(input_vector);
8208  evaluate(this->begin_dof_values(), evaluation_flag);
8209  }
8210 }
8211 
8212 
8213 
8214 template <int dim,
8215  int fe_degree,
8216  int n_q_points_1d,
8217  int n_components_,
8218  typename Number,
8219  typename VectorizedArrayType>
8220 inline void
8221 FEEvaluation<dim,
8222  fe_degree,
8223  n_q_points_1d,
8224  n_components_,
8225  Number,
8226  VectorizedArrayType>::integrate(const bool integrate_values,
8227  const bool integrate_gradients)
8228 {
8229  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8230 
8231 # ifdef DEBUG
8232  this->dof_values_initialized = true;
8233 # endif
8234 }
8235 
8236 
8237 
8238 template <int dim,
8239  int fe_degree,
8240  int n_q_points_1d,
8241  int n_components_,
8242  typename Number,
8243  typename VectorizedArrayType>
8244 inline void
8245 FEEvaluation<dim,
8246  fe_degree,
8247  n_q_points_1d,
8248  n_components_,
8249  Number,
8250  VectorizedArrayType>::
8251  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
8252 {
8253  integrate(integration_flag, this->values_dofs[0]);
8254 
8255 # ifdef DEBUG
8256  this->dof_values_initialized = true;
8257 # endif
8258 }
8259 
8260 
8261 
8262 template <int dim,
8263  int fe_degree,
8264  int n_q_points_1d,
8265  int n_components_,
8266  typename Number,
8267  typename VectorizedArrayType>
8268 inline void
8269 FEEvaluation<dim,
8270  fe_degree,
8271  n_q_points_1d,
8272  n_components_,
8273  Number,
8274  VectorizedArrayType>::integrate(const bool integrate_values,
8275  const bool integrate_gradients,
8276  VectorizedArrayType *values_array)
8277 {
8279  (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
8280  (integrate_gradients ? EvaluationFlags::gradients :
8282  integrate(flag, values_array);
8283 }
8284 
8285 
8286 
8287 template <int dim,
8288  int fe_degree,
8289  int n_q_points_1d,
8290  int n_components_,
8291  typename Number,
8292  typename VectorizedArrayType>
8293 inline void
8294 FEEvaluation<dim,
8295  fe_degree,
8296  n_q_points_1d,
8297  n_components_,
8298  Number,
8299  VectorizedArrayType>::
8300  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
8301  VectorizedArrayType * values_array)
8302 {
8303 # ifdef DEBUG
8304  if (integration_flag & EvaluationFlags::values)
8305  Assert(this->values_quad_submitted == true,
8307  if (integration_flag & EvaluationFlags::gradients)
8308  Assert(this->gradients_quad_submitted == true,
8310 # endif
8311  Assert(this->matrix_info != nullptr ||
8312  this->mapped_geometry->is_initialized(),
8313  ExcNotInitialized());
8314 
8315  Assert(
8316  (integration_flag &
8317  ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
8318  ExcMessage(
8319  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8320 
8321  if (fe_degree > -1)
8323  integrate(n_components,
8324  integration_flag,
8325  *this->data,
8326  values_array,
8327  this->values_quad,
8328  this->gradients_quad,
8329  this->scratch_data,
8330  false);
8331  else
8333  n_components,
8334  integration_flag,
8335  *this->data,
8336  values_array,
8337  this->values_quad,
8338  this->gradients_quad,
8339  this->scratch_data,
8340  false);
8341 
8342 # ifdef DEBUG
8343  this->dof_values_initialized = true;
8344 # endif
8345 }
8346 
8347 
8348 
8349 template <int dim,
8350  int fe_degree,
8351  int n_q_points_1d,
8352  int n_components_,
8353  typename Number,
8354  typename VectorizedArrayType>
8355 template <typename VectorType>
8356 inline void
8357 FEEvaluation<
8358