Reference documentation for deal.II version Git d63a52b21e 2020-07-11 21:58:20 -0400
\(\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 
40 
41 
43 
44 
45 
46 namespace internal
47 {
49 }
50 
51 template <int dim,
52  int fe_degree,
53  int n_q_points_1d = fe_degree + 1,
54  int n_components_ = 1,
55  typename Number = double,
56  typename VectorizedArrayType = VectorizedArray<Number>>
58 
59 
66 namespace EvaluationFlags
67 {
76  {
80  nothing = 0,
84  values = 0x1,
88  gradients = 0x2,
92  hessians = 0x4
93  };
94 
95 
105  inline EvaluationFlags
107  {
108  return static_cast<EvaluationFlags>(static_cast<unsigned int>(f1) |
109  static_cast<unsigned int>(f2));
110  }
111 
112 
113 
120  inline EvaluationFlags &
122  {
123  f1 = f1 | f2;
124  return f1;
125  }
126 
127 
138  const EvaluationFlags f2)
139  {
140  return static_cast<EvaluationFlags>(static_cast<unsigned int>(f1) &
141  static_cast<unsigned int>(f2));
142  }
143 
144 
151  inline EvaluationFlags &
153  {
154  f1 = f1 & f2;
155  return f1;
156  }
157 
158 } // namespace EvaluationFlags
159 
190 template <int dim,
191  int n_components_,
192  typename Number,
193  bool is_face = false,
194  typename VectorizedArrayType = VectorizedArray<Number>>
196 {
197  static_assert(
198  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
199  "Type of Number and of VectorizedArrayType do not match.");
200 
201 public:
202  using number_type = Number;
204  using gradient_type =
206  static constexpr unsigned int dimension = dim;
207  static constexpr unsigned int n_components = n_components_;
208 
216  ~FEEvaluationBase();
217 
226  unsigned int
227  get_mapping_data_index_offset() const;
228 
236  get_cell_type() const;
237 
242  get_shape_info() const;
243 
245 
282  template <typename VectorType>
283  void
284  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
285 
314  template <typename VectorType>
315  void
316  read_dof_values_plain(const VectorType & src,
317  const unsigned int first_index = 0);
318 
350  template <typename VectorType>
351  void
352  distribute_local_to_global(
353  VectorType & dst,
354  const unsigned int first_index = 0,
355  const std::bitset<VectorizedArrayType::size()> &mask =
356  std::bitset<VectorizedArrayType::size()>().flip()) const;
357 
396  template <typename VectorType>
397  void
398  set_dof_values(VectorType & dst,
399  const unsigned int first_index = 0,
400  const std::bitset<VectorizedArrayType::size()> &mask =
401  std::bitset<VectorizedArrayType::size()>().flip()) const;
402 
404 
425  value_type
426  get_dof_value(const unsigned int dof) const;
427 
438  void
439  submit_dof_value(const value_type val_in, const unsigned int dof);
440 
452  value_type
453  get_value(const unsigned int q_point) const;
454 
466  void
467  submit_value(const value_type val_in, const unsigned int q_point);
468 
479  get_gradient(const unsigned int q_point) const;
480 
494  value_type
495  get_normal_derivative(const unsigned int q_point) const;
496 
509  void
510  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
511 
529  void
530  submit_normal_derivative(const value_type grad_in,
531  const unsigned int q_point);
532 
544  get_hessian(const unsigned int q_point) const;
545 
555  get_hessian_diagonal(const unsigned int q_point) const;
556 
567  value_type
568  get_laplacian(const unsigned int q_point) const;
569 
570 #ifdef DOXYGEN
571  // doxygen does not anyhow mention functions coming from partial template
572  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
573  // For now, hack-in those functions manually only to fix documentation:
574 
581  VectorizedArrayType
582  get_divergence(const unsigned int q_point) const;
583 
593  get_symmetric_gradient(const unsigned int q_point) const;
594 
602  get_curl(const unsigned int q_point) const;
603 
619  void
620  submit_divergence(const VectorizedArrayType div_in,
621  const unsigned int q_point);
622 
639  void
640  submit_symmetric_gradient(
642  const unsigned int q_point);
643 
656  void
657  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
658  const unsigned int q_point);
659 
660 #endif
661 
678  value_type
679  integrate_value() const;
680 
685  VectorizedArrayType
686  JxW(const unsigned int q_index) const;
687 
695  inverse_jacobian(const unsigned int q_index) const;
696 
710  get_normal_vector(const unsigned int q_point) const;
711 
718  VectorizedArrayType
719  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
720 
727  void
728  set_cell_data(AlignedVector<VectorizedArrayType> &array,
729  const VectorizedArrayType & value) const;
730 
735  template <typename T>
736  std::array<T, VectorizedArrayType::size()>
737  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
738  &array) const;
739 
744  template <typename T>
745  void
746  set_cell_data(
747  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
748  const std::array<T, VectorizedArrayType::size()> & value) const;
749 
754  std::array<unsigned int, VectorizedArrayType::size()>
755  get_cell_ids() const;
756 
758 
771  const VectorizedArrayType *
772  begin_dof_values() const;
773 
782  VectorizedArrayType *
783  begin_dof_values();
784 
795  const VectorizedArrayType *
796  begin_values() const;
797 
808  VectorizedArrayType *
809  begin_values();
810 
822  const VectorizedArrayType *
823  begin_gradients() const;
824 
836  VectorizedArrayType *
837  begin_gradients();
838 
851  const VectorizedArrayType *
852  begin_hessians() const;
853 
866  VectorizedArrayType *
867  begin_hessians();
868 
874  const std::vector<unsigned int> &
875  get_internal_dof_numbering() const;
876 
884  get_scratch_data() const;
885 
887 
888 protected:
899  const unsigned int dof_no,
900  const unsigned int first_selected_component,
901  const unsigned int quad_no,
902  const unsigned int fe_degree,
903  const unsigned int n_q_points,
904  const bool is_interior_face);
905 
940  template <int n_components_other>
941  FEEvaluationBase(const Mapping<dim> & mapping,
942  const FiniteElement<dim> &fe,
943  const Quadrature<1> & quadrature,
944  const UpdateFlags update_flags,
945  const unsigned int first_selected_component,
946  const FEEvaluationBase<dim,
947  n_components_other,
948  Number,
949  is_face,
950  VectorizedArrayType> *other);
951 
958  FEEvaluationBase(const FEEvaluationBase &other);
959 
967  operator=(const FEEvaluationBase &other);
968 
975  template <typename VectorType, typename VectorOperation>
976  void
977  read_write_operation(const VectorOperation &operation,
978  VectorType * vectors[],
979  const std::bitset<VectorizedArrayType::size()> &mask,
980  const bool apply_constraints = true) const;
981 
989  template <typename VectorType, typename VectorOperation>
990  void
991  read_write_operation_contiguous(
992  const VectorOperation & operation,
993  VectorType * vectors[],
994  const std::bitset<VectorizedArrayType::size()> &mask) const;
995 
1003  template <typename VectorType, typename VectorOperation>
1004  void
1005  read_write_operation_global(const VectorOperation &operation,
1006  VectorType * vectors[]) const;
1007 
1012 
1018  VectorizedArrayType *scratch_data;
1019 
1032  VectorizedArrayType *values_dofs[n_components];
1033 
1045  VectorizedArrayType *values_quad[n_components];
1046 
1060  VectorizedArrayType *gradients_quad[n_components][dim];
1061 
1073  VectorizedArrayType *hessians_quad[n_components][(dim * (dim + 1)) / 2];
1074 
1078  const unsigned int quad_no;
1079 
1084  const unsigned int n_fe_components;
1085 
1090  const unsigned int active_fe_index;
1091 
1096  const unsigned int active_quad_index;
1097 
1101  const unsigned int n_quadrature_points;
1102 
1107 
1114 
1122  (is_face ? dim - 1 : dim),
1123  dim,
1124  Number,
1125  VectorizedArrayType> *mapping_data;
1126 
1134 
1140 
1147  const VectorizedArrayType *J_value;
1148 
1153 
1158 
1162  const Number *quadrature_weights;
1163 
1168  unsigned int cell;
1169 
1175 
1181 
1186  unsigned int face_no;
1187 
1192  unsigned int face_orientation;
1193 
1201  unsigned int subface_index;
1202 
1210 
1217 
1224 
1231 
1238 
1245 
1252 
1257  std::shared_ptr<internal::MatrixFreeFunctions::
1258  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
1260 
1265  const unsigned int first_selected_component;
1266 
1271  mutable std::vector<types::global_dof_index> local_dof_indices;
1272 
1273 private:
1278  void
1279  set_data_pointers();
1280 
1281  // Make other FEEvaluationBase as well as FEEvaluation objects friends.
1282  template <int, int, typename, bool, typename>
1283  friend class FEEvaluationBase;
1284  template <int, int, int, int, typename, typename>
1285  friend class FEEvaluation;
1286 };
1287 
1288 
1289 
1297 template <int dim,
1298  int n_components_,
1299  typename Number,
1300  bool is_face,
1301  typename VectorizedArrayType = VectorizedArray<Number>>
1303  n_components_,
1304  Number,
1305  is_face,
1306  VectorizedArrayType>
1307 {
1308  static_assert(
1309  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1310  "Type of Number and of VectorizedArrayType do not match.");
1311 
1312 public:
1313  using number_type = Number;
1315  using gradient_type =
1317  static constexpr unsigned int dimension = dim;
1318  static constexpr unsigned int n_components = n_components_;
1319  using BaseClass =
1321 
1322 protected:
1332  const unsigned int dof_no,
1333  const unsigned int first_selected_component,
1334  const unsigned int quad_no,
1335  const unsigned int fe_degree,
1336  const unsigned int n_q_points,
1337  const bool is_interior_face = true);
1338 
1343  template <int n_components_other>
1344  FEEvaluationAccess(const Mapping<dim> & mapping,
1345  const FiniteElement<dim> &fe,
1346  const Quadrature<1> & quadrature,
1347  const UpdateFlags update_flags,
1348  const unsigned int first_selected_component,
1349  const FEEvaluationBase<dim,
1350  n_components_other,
1351  Number,
1352  is_face,
1353  VectorizedArrayType> *other);
1354 
1358  FEEvaluationAccess(const FEEvaluationAccess &other);
1359 
1364  operator=(const FEEvaluationAccess &other);
1365 };
1366 
1367 
1368 
1377 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1378 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1379  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1380 {
1381  static_assert(
1382  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1383  "Type of Number and of VectorizedArrayType do not match.");
1384 
1385 public:
1386  using number_type = Number;
1387  using value_type = VectorizedArrayType;
1389  static constexpr unsigned int dimension = dim;
1390  using BaseClass =
1392 
1395  value_type
1396  get_dof_value(const unsigned int dof) const;
1397 
1400  void
1401  submit_dof_value(const value_type val_in, const unsigned int dof);
1402 
1405  value_type
1406  get_value(const unsigned int q_point) const;
1407 
1410  void
1411  submit_value(const value_type val_in, const unsigned int q_point);
1412 
1415  void
1416  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
1417  const unsigned int q_point);
1418 
1422  get_gradient(const unsigned int q_point) const;
1423 
1426  value_type
1427  get_normal_derivative(const unsigned int q_point) const;
1428 
1431  void
1432  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1433 
1436  void
1437  submit_normal_derivative(const value_type grad_in,
1438  const unsigned int q_point);
1439 
1443  get_hessian(unsigned int q_point) const;
1444 
1448  get_hessian_diagonal(const unsigned int q_point) const;
1449 
1452  value_type
1453  get_laplacian(const unsigned int q_point) const;
1454 
1457  value_type
1458  integrate_value() const;
1459 
1460 protected:
1470  const unsigned int dof_no,
1471  const unsigned int first_selected_component,
1472  const unsigned int quad_no,
1473  const unsigned int fe_degree,
1474  const unsigned int n_q_points,
1475  const bool is_interior_face = true);
1476 
1481  template <int n_components_other>
1482  FEEvaluationAccess(const Mapping<dim> & mapping,
1483  const FiniteElement<dim> &fe,
1484  const Quadrature<1> & quadrature,
1485  const UpdateFlags update_flags,
1486  const unsigned int first_selected_component,
1487  const FEEvaluationBase<dim,
1488  n_components_other,
1489  Number,
1490  is_face,
1491  VectorizedArrayType> *other);
1492 
1496  FEEvaluationAccess(const FEEvaluationAccess &other);
1497 
1502  operator=(const FEEvaluationAccess &other);
1503 };
1504 
1505 
1506 
1516 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1517 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1518  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1519 {
1520  static_assert(
1521  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1522  "Type of Number and of VectorizedArrayType do not match.");
1523 
1524 public:
1525  using number_type = Number;
1528  static constexpr unsigned int dimension = dim;
1529  static constexpr unsigned int n_components = dim;
1530  using BaseClass =
1532 
1536  get_gradient(const unsigned int q_point) const;
1537 
1542  VectorizedArrayType
1543  get_divergence(const unsigned int q_point) const;
1544 
1552  get_symmetric_gradient(const unsigned int q_point) const;
1553 
1559  get_curl(const unsigned int q_point) const;
1560 
1564  get_hessian(const unsigned int q_point) const;
1565 
1569  get_hessian_diagonal(const unsigned int q_point) const;
1570 
1573  void
1574  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1575 
1584  void
1585  submit_gradient(
1586  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1587  const unsigned int q_point);
1588 
1597  void
1598  submit_divergence(const VectorizedArrayType div_in,
1599  const unsigned int q_point);
1600 
1609  void
1610  submit_symmetric_gradient(
1612  const unsigned int q_point);
1613 
1618  void
1619  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1620  const unsigned int q_point);
1621 
1622 protected:
1632  const unsigned int dof_no,
1633  const unsigned int first_selected_component,
1634  const unsigned int quad_no,
1635  const unsigned int dofs_per_cell,
1636  const unsigned int n_q_points,
1637  const bool is_interior_face = true);
1638 
1643  template <int n_components_other>
1644  FEEvaluationAccess(const Mapping<dim> & mapping,
1645  const FiniteElement<dim> &fe,
1646  const Quadrature<1> & quadrature,
1647  const UpdateFlags update_flags,
1648  const unsigned int first_selected_component,
1649  const FEEvaluationBase<dim,
1650  n_components_other,
1651  Number,
1652  is_face,
1653  VectorizedArrayType> *other);
1654 
1658  FEEvaluationAccess(const FEEvaluationAccess &other);
1659 
1664  operator=(const FEEvaluationAccess &other);
1665 };
1666 
1667 
1678 template <typename Number, bool is_face, typename VectorizedArrayType>
1679 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1680  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1681 {
1682  static_assert(
1683  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1684  "Type of Number and of VectorizedArrayType do not match.");
1685 
1686 public:
1687  using number_type = Number;
1688  using value_type = VectorizedArrayType;
1690  static constexpr unsigned int dimension = 1;
1691  using BaseClass =
1693 
1696  value_type
1697  get_dof_value(const unsigned int dof) const;
1698 
1701  void
1702  submit_dof_value(const value_type val_in, const unsigned int dof);
1703 
1706  value_type
1707  get_value(const unsigned int q_point) const;
1708 
1711  void
1712  submit_value(const value_type val_in, const unsigned int q_point);
1713 
1716  void
1717  submit_value(const gradient_type val_in, const unsigned int q_point);
1718 
1722  get_gradient(const unsigned int q_point) const;
1723 
1726  value_type
1727  get_normal_derivative(const unsigned int q_point) const;
1728 
1731  void
1732  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1733 
1736  void
1737  submit_gradient(const value_type grad_in, const unsigned int q_point);
1738 
1741  void
1742  submit_normal_derivative(const value_type grad_in,
1743  const unsigned int q_point);
1744 
1747  void
1748  submit_normal_derivative(const gradient_type grad_in,
1749  const unsigned int q_point);
1750 
1754  get_hessian(unsigned int q_point) const;
1755 
1759  get_hessian_diagonal(const unsigned int q_point) const;
1760 
1763  value_type
1764  get_laplacian(const unsigned int q_point) const;
1765 
1768  value_type
1769  integrate_value() const;
1770 
1771 protected:
1780  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
1781  const unsigned int dof_no,
1782  const unsigned int first_selected_component,
1783  const unsigned int quad_no,
1784  const unsigned int fe_degree,
1785  const unsigned int n_q_points,
1786  const bool is_interior_face = true);
1787 
1792  template <int n_components_other>
1793  FEEvaluationAccess(const Mapping<1> & mapping,
1794  const FiniteElement<1> &fe,
1795  const Quadrature<1> & quadrature,
1796  const UpdateFlags update_flags,
1797  const unsigned int first_selected_component,
1798  const FEEvaluationBase<1,
1799  n_components_other,
1800  Number,
1801  is_face,
1802  VectorizedArrayType> *other);
1803 
1807  FEEvaluationAccess(const FEEvaluationAccess &other);
1808 
1813  operator=(const FEEvaluationAccess &other);
1814 };
1815 
1816 
1817 
2361 template <int dim,
2362  int fe_degree,
2363  int n_q_points_1d,
2364  int n_components_,
2365  typename Number,
2366  typename VectorizedArrayType>
2367 class FEEvaluation : public FEEvaluationAccess<dim,
2368  n_components_,
2369  Number,
2370  false,
2371  VectorizedArrayType>
2372 {
2373  static_assert(
2374  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2375  "Type of Number and of VectorizedArrayType do not match.");
2376 
2377 public:
2381  using BaseClass =
2383 
2387  using number_type = Number;
2388 
2395 
2402 
2406  static constexpr unsigned int dimension = dim;
2407 
2412  static constexpr unsigned int n_components = n_components_;
2413 
2420  static constexpr unsigned int static_n_q_points =
2421  Utilities::pow(n_q_points_1d, dim);
2422 
2430  static constexpr unsigned int static_dofs_per_component =
2431  Utilities::pow(fe_degree + 1, dim);
2432 
2440  static constexpr unsigned int tensor_dofs_per_cell =
2441  static_dofs_per_component * n_components;
2442 
2450  static constexpr unsigned int static_dofs_per_cell =
2451  static_dofs_per_component * n_components;
2452 
2479  const unsigned int dof_no = 0,
2480  const unsigned int quad_no = 0,
2481  const unsigned int first_selected_component = 0);
2482 
2509  FEEvaluation(const Mapping<dim> & mapping,
2510  const FiniteElement<dim> &fe,
2511  const Quadrature<1> & quadrature,
2512  const UpdateFlags update_flags,
2513  const unsigned int first_selected_component = 0);
2514 
2520  FEEvaluation(const FiniteElement<dim> &fe,
2521  const Quadrature<1> & quadrature,
2522  const UpdateFlags update_flags,
2523  const unsigned int first_selected_component = 0);
2524 
2535  template <int n_components_other>
2536  FEEvaluation(const FiniteElement<dim> & fe,
2537  const FEEvaluationBase<dim,
2538  n_components_other,
2539  Number,
2540  false,
2541  VectorizedArrayType> &other,
2542  const unsigned int first_selected_component = 0);
2543 
2550  FEEvaluation(const FEEvaluation &other);
2551 
2558  FEEvaluation &
2559  operator=(const FEEvaluation &other);
2560 
2569  void
2570  reinit(const unsigned int cell_batch_index);
2571 
2584  template <bool level_dof_access>
2585  void
2587 
2598  void
2599  reinit(const typename Triangulation<dim>::cell_iterator &cell);
2600 
2610  void
2611  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2612 
2617  void
2618  evaluate(const bool evaluate_values,
2619  const bool evaluate_gradients,
2620  const bool evaluate_hessians = false);
2621 
2634  void
2635  evaluate(const VectorizedArrayType * values_array,
2636  const EvaluationFlags::EvaluationFlags evaluation_flag);
2637 
2642  void
2643  evaluate(const VectorizedArrayType *values_array,
2644  const bool evaluate_values,
2645  const bool evaluate_gradients,
2646  const bool evaluate_hessians = false);
2647 
2661  template <typename VectorType>
2662  void
2663  gather_evaluate(const VectorType & input_vector,
2664  const EvaluationFlags::EvaluationFlags evaluation_flag);
2665 
2669  template <typename VectorType>
2670  void
2671  gather_evaluate(const VectorType &input_vector,
2672  const bool evaluate_values,
2673  const bool evaluate_gradients,
2674  const bool evaluate_hessians = false);
2675 
2686  void
2687  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2688 
2689 
2693  void
2694  integrate(const bool integrate_values, const bool integrate_gradients);
2695 
2707  void
2708  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2709  VectorizedArrayType * values_array);
2710 
2714  void
2715  integrate(const bool integrate_values,
2716  const bool integrate_gradients,
2717  VectorizedArrayType *values_array);
2718 
2732  template <typename VectorType>
2733  void
2734  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
2735  VectorType & output_vector);
2736 
2740  template <typename VectorType>
2741  void
2742  integrate_scatter(const bool integrate_values,
2743  const bool integrate_gradients,
2744  VectorType &output_vector);
2745 
2751  quadrature_point(const unsigned int q_point) const;
2752 
2759  const unsigned int dofs_per_component;
2760 
2767  const unsigned int dofs_per_cell;
2768 
2776  const unsigned int n_q_points;
2777 
2778 private:
2783  void
2784  check_template_arguments(const unsigned int fe_no,
2785  const unsigned int first_selected_component);
2786 };
2787 
2788 
2789 
2825 template <int dim,
2826  int fe_degree,
2827  int n_q_points_1d = fe_degree + 1,
2828  int n_components_ = 1,
2829  typename Number = double,
2830  typename VectorizedArrayType = VectorizedArray<Number>>
2832  n_components_,
2833  Number,
2834  true,
2835  VectorizedArrayType>
2836 {
2837  static_assert(
2838  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2839  "Type of Number and of VectorizedArrayType do not match.");
2840 
2841 public:
2845  using BaseClass =
2847 
2851  using number_type = Number;
2852 
2859 
2866 
2870  static constexpr unsigned int dimension = dim;
2871 
2876  static constexpr unsigned int n_components = n_components_;
2877 
2885  static constexpr unsigned int static_n_q_points =
2886  Utilities::pow(n_q_points_1d, dim - 1);
2887 
2894  static constexpr unsigned int static_n_q_points_cell =
2895  Utilities::pow(n_q_points_1d, dim);
2896 
2903  static constexpr unsigned int static_dofs_per_component =
2904  Utilities::pow(fe_degree + 1, dim);
2905 
2912  static constexpr unsigned int tensor_dofs_per_cell =
2913  static_dofs_per_component * n_components;
2914 
2921  static constexpr unsigned int static_dofs_per_cell =
2922  static_dofs_per_component * n_components;
2923 
2955  const bool is_interior_face = true,
2956  const unsigned int dof_no = 0,
2957  const unsigned int quad_no = 0,
2958  const unsigned int first_selected_component = 0);
2959 
2970  void
2971  reinit(const unsigned int face_batch_number);
2972 
2980  void
2981  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2982 
2993  void
2994  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2995 
2999  void
3000  evaluate(const bool evaluate_values, const bool evaluate_gradients);
3001 
3014  void
3015  evaluate(const VectorizedArrayType * values_array,
3016  const EvaluationFlags::EvaluationFlags evaluation_flag);
3017 
3021  void
3022  evaluate(const VectorizedArrayType *values_array,
3023  const bool evaluate_values,
3024  const bool evaluate_gradients);
3025 
3037  template <typename VectorType>
3038  void
3039  gather_evaluate(const VectorType & input_vector,
3040  const EvaluationFlags::EvaluationFlags evaluation_flag);
3041 
3045  template <typename VectorType>
3046  void
3047  gather_evaluate(const VectorType &input_vector,
3048  const bool evaluate_values,
3049  const bool evaluate_gradients);
3050 
3060  void
3061  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3062 
3066  void
3067  integrate(const bool integrate_values, const bool integrate_gradients);
3068 
3077  void
3078  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
3079  VectorizedArrayType * values_array);
3080 
3084  void
3085  integrate(const bool integrate_values,
3086  const bool integrate_gradients,
3087  VectorizedArrayType *values_array);
3088 
3100  template <typename VectorType>
3101  void
3102  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
3103  VectorType & output_vector);
3104 
3108  template <typename VectorType>
3109  void
3110  integrate_scatter(const bool integrate_values,
3111  const bool integrate_gradients,
3112  VectorType &output_vector);
3113 
3119  quadrature_point(const unsigned int q_point) const;
3120 
3127  const unsigned int dofs_per_component;
3128 
3135  const unsigned int dofs_per_cell;
3136 
3144  const unsigned int n_q_points;
3145 };
3146 
3147 
3148 
3149 namespace internal
3150 {
3151  namespace MatrixFreeFunctions
3152  {
3153  // a helper function to compute the number of DoFs of a DGP element at
3154  // compile time, depending on the degree
3155  template <int dim, int degree>
3157  {
3158  // this division is always without remainder
3159  static constexpr unsigned int value =
3160  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
3161  };
3162 
3163  // base specialization: 1d elements have 'degree+1' degrees of freedom
3164  template <int degree>
3165  struct DGP_dofs_per_component<1, degree>
3166  {
3167  static constexpr unsigned int value = degree + 1;
3168  };
3169  } // namespace MatrixFreeFunctions
3170 } // namespace internal
3171 
3172 
3173 /*----------------------- Inline functions ----------------------------------*/
3174 
3175 #ifndef DOXYGEN
3176 
3177 
3178 
3179 /*----------------------- FEEvaluationBase ----------------------------------*/
3180 
3181 template <int dim,
3182  int n_components_,
3183  typename Number,
3184  bool is_face,
3185  typename VectorizedArrayType>
3186 inline FEEvaluationBase<dim,
3187  n_components_,
3188  Number,
3189  is_face,
3190  VectorizedArrayType>::
3191  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
3192  const unsigned int dof_no,
3193  const unsigned int first_selected_component,
3194  const unsigned int quad_no_in,
3195  const unsigned int fe_degree,
3196  const unsigned int n_q_points,
3197  const bool is_interior_face)
3198  : scratch_data_array(data_in.acquire_scratch_data())
3199  , quad_no(quad_no_in)
3200  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
3201  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
3202  data_in.get_dof_info(dof_no).fe_index_from_degree(
3203  first_selected_component,
3204  fe_degree) :
3205  0)
3206  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
3207  (is_face ? data_in.get_mapping_info()
3208  .face_data[quad_no_in]
3209  .quad_index_from_n_q_points(n_q_points) :
3210  data_in.get_mapping_info()
3211  .cell_data[quad_no_in]
3212  .quad_index_from_n_q_points(n_q_points)) :
3213  0)
3214  , n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ?
3215  n_q_points :
3216  (is_face ? data_in
3217  .get_shape_info(dof_no,
3218  quad_no_in,
3219  active_fe_index,
3220  active_quad_index)
3221  .n_q_points_face :
3222  data_in
3223  .get_shape_info(dof_no,
3224  quad_no_in,
3225  active_fe_index,
3226  active_quad_index)
3227  .n_q_points))
3228  , matrix_info(&data_in)
3229  , dof_info(&data_in.get_dof_info(dof_no))
3230  , mapping_data(
3231  internal::MatrixFreeFunctions::
3232  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3233  data_in.get_mapping_info(),
3234  quad_no))
3235  , data(&data_in.get_shape_info(
3236  dof_no,
3237  quad_no_in,
3238  dof_info->component_to_base_index[first_selected_component],
3239  active_fe_index,
3240  active_quad_index))
3241  , jacobian(nullptr)
3242  , J_value(nullptr)
3243  , normal_vectors(nullptr)
3244  , normal_x_jacobian(nullptr)
3245  , quadrature_weights(
3246  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3247  , cell(numbers::invalid_unsigned_int)
3248  , is_interior_face(is_interior_face)
3249  , dof_access_index(
3250  is_face ?
3251  (is_interior_face ?
3252  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3253  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3254  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3255  , cell_type(internal::MatrixFreeFunctions::general)
3256  , dof_values_initialized(false)
3257  , values_quad_initialized(false)
3258  , gradients_quad_initialized(false)
3259  , hessians_quad_initialized(false)
3260  , values_quad_submitted(false)
3261  , gradients_quad_submitted(false)
3262  , first_selected_component(first_selected_component)
3263 {
3264  set_data_pointers();
3265  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3266  AssertDimension(matrix_info->get_task_info().vectorization_length,
3267  VectorizedArrayType::size());
3268  AssertDimension((is_face ? data->n_q_points_face : data->n_q_points),
3269  n_quadrature_points);
3270  AssertDimension(n_quadrature_points,
3271  mapping_data->descriptor[active_quad_index].n_q_points);
3272  Assert(
3273  dof_info->start_components.back() == 1 ||
3274  static_cast<int>(n_components_) <=
3275  static_cast<int>(
3276  dof_info->start_components
3277  [dof_info->component_to_base_index[first_selected_component] + 1]) -
3278  first_selected_component,
3279  ExcMessage(
3280  "You tried to construct a vector-valued evaluator with " +
3281  std::to_string(n_components) +
3282  " components. However, "
3283  "the current base element has only " +
3285  dof_info->start_components
3286  [dof_info->component_to_base_index[first_selected_component] + 1] -
3287  first_selected_component) +
3288  " components left when starting from local element index " +
3290  first_selected_component -
3291  dof_info->start_components
3292  [dof_info->component_to_base_index[first_selected_component]]) +
3293  " (global index " + std::to_string(first_selected_component) + ")"));
3294 
3295  // do not check for correct dimensions of data fields here, should be done
3296  // in derived classes
3297 }
3298 
3299 
3300 
3301 template <int dim,
3302  int n_components_,
3303  typename Number,
3304  bool is_face,
3305  typename VectorizedArrayType>
3306 template <int n_components_other>
3307 inline FEEvaluationBase<dim,
3308  n_components_,
3309  Number,
3310  is_face,
3311  VectorizedArrayType>::
3312  FEEvaluationBase(const Mapping<dim> & mapping,
3313  const FiniteElement<dim> &fe,
3314  const Quadrature<1> & quadrature,
3315  const UpdateFlags update_flags,
3316  const unsigned int first_selected_component,
3317  const FEEvaluationBase<dim,
3318  n_components_other,
3319  Number,
3320  is_face,
3321  VectorizedArrayType> *other)
3322  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3324  , n_fe_components(n_components_)
3325  , active_fe_index(numbers::invalid_unsigned_int)
3326  , active_quad_index(numbers::invalid_unsigned_int)
3327  , n_quadrature_points(
3328  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3329  , matrix_info(nullptr)
3330  , dof_info(nullptr)
3331  , mapping_data(nullptr)
3332  ,
3333  // select the correct base element from the given FE component
3335  quadrature,
3336  fe,
3337  fe.component_to_base_index(first_selected_component).first))
3338  , jacobian(nullptr)
3339  , J_value(nullptr)
3340  , normal_vectors(nullptr)
3341  , normal_x_jacobian(nullptr)
3342  , quadrature_weights(nullptr)
3343  , cell(0)
3345  , is_interior_face(true)
3347  , dof_values_initialized(false)
3348  , values_quad_initialized(false)
3349  , gradients_quad_initialized(false)
3350  , hessians_quad_initialized(false)
3351  , values_quad_submitted(false)
3352  , gradients_quad_submitted(false)
3353  ,
3354  // keep the number of the selected component within the current base element
3355  // for reading dof values
3356  first_selected_component(first_selected_component)
3357 {
3358  set_data_pointers();
3359 
3360  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3361  ExcInternalError());
3362  if (other != nullptr &&
3363  other->mapped_geometry->get_quadrature() == quadrature)
3364  mapped_geometry = other->mapped_geometry;
3365  else
3366  mapped_geometry =
3367  std::make_shared<internal::MatrixFreeFunctions::
3369  mapping, quadrature, update_flags);
3370  cell = 0;
3371 
3372  mapping_data = &mapped_geometry->get_data_storage();
3373  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3374  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3375 
3376  const unsigned int base_element_number =
3377  fe.component_to_base_index(first_selected_component).first;
3378  Assert(fe.element_multiplicity(base_element_number) == 1 ||
3379  fe.element_multiplicity(base_element_number) -
3380  first_selected_component >=
3381  n_components_,
3382  ExcMessage("The underlying element must at least contain as many "
3383  "components as requested by this class"));
3384  (void)base_element_number;
3385 }
3386 
3387 
3388 
3389 template <int dim,
3390  int n_components_,
3391  typename Number,
3392  bool is_face,
3393  typename VectorizedArrayType>
3394 inline FEEvaluationBase<dim,
3395  n_components_,
3396  Number,
3397  is_face,
3398  VectorizedArrayType>::
3399  FEEvaluationBase(const FEEvaluationBase<dim,
3400  n_components_,
3401  Number,
3402  is_face,
3403  VectorizedArrayType> &other)
3404  : scratch_data_array(other.matrix_info == nullptr ?
3406  other.matrix_info->acquire_scratch_data())
3407  , quad_no(other.quad_no)
3408  , n_fe_components(other.n_fe_components)
3409  , active_fe_index(other.active_fe_index)
3410  , active_quad_index(other.active_quad_index)
3411  , n_quadrature_points(other.n_quadrature_points)
3412  , matrix_info(other.matrix_info)
3413  , dof_info(other.dof_info)
3414  , mapping_data(other.mapping_data)
3415  , data(other.matrix_info == nullptr ?
3417  *other.data) :
3418  other.data)
3419  , jacobian(nullptr)
3420  , J_value(nullptr)
3421  , normal_vectors(nullptr)
3422  , normal_x_jacobian(nullptr)
3423  , quadrature_weights(
3424  other.matrix_info == nullptr ?
3425  nullptr :
3426  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3429  , is_interior_face(other.is_interior_face)
3430  , dof_access_index(other.dof_access_index)
3431  , dof_values_initialized(false)
3432  , values_quad_initialized(false)
3433  , gradients_quad_initialized(false)
3434  , hessians_quad_initialized(false)
3435  , values_quad_submitted(false)
3436  , gradients_quad_submitted(false)
3437  , first_selected_component(other.first_selected_component)
3438 {
3439  set_data_pointers();
3440 
3441  // Create deep copy of mapped geometry for use in parallel...
3442  if (other.mapped_geometry.get() != nullptr)
3443  {
3444  mapped_geometry = std::make_shared<
3445  internal::MatrixFreeFunctions::
3446  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3447  other.mapped_geometry->get_fe_values().get_mapping(),
3448  other.mapped_geometry->get_quadrature(),
3449  other.mapped_geometry->get_fe_values().get_update_flags());
3450  mapping_data = &mapped_geometry->get_data_storage();
3451  cell = 0;
3452 
3453  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3454  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3455  }
3456 }
3457 
3458 
3459 
3460 template <int dim,
3461  int n_components_,
3462  typename Number,
3463  bool is_face,
3464  typename VectorizedArrayType>
3465 inline FEEvaluationBase<dim,
3466  n_components_,
3467  Number,
3468  is_face,
3469  VectorizedArrayType> &
3471 operator=(const FEEvaluationBase<dim,
3472  n_components_,
3473  Number,
3474  is_face,
3475  VectorizedArrayType> &other)
3476 {
3477  AssertDimension(quad_no, other.quad_no);
3478  AssertDimension(n_fe_components, other.n_fe_components);
3479  AssertDimension(active_fe_index, other.active_fe_index);
3480  AssertDimension(active_quad_index, other.active_quad_index);
3481  AssertDimension(first_selected_component, other.first_selected_component);
3482 
3483  // release old memory
3484  if (matrix_info == nullptr)
3485  {
3486  delete data;
3487  delete scratch_data_array;
3488  }
3489  else
3490  {
3491  matrix_info->release_scratch_data(scratch_data_array);
3492  }
3493 
3494  matrix_info = other.matrix_info;
3495  dof_info = other.dof_info;
3496  mapping_data = other.mapping_data;
3497  if (other.matrix_info == nullptr)
3498  {
3500  *other.data);
3501  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3502  }
3503  else
3504  {
3505  data = other.data;
3506  scratch_data_array = matrix_info->acquire_scratch_data();
3507  }
3508  set_data_pointers();
3509 
3510  quadrature_weights =
3511  (mapping_data != nullptr ?
3512  mapping_data->descriptor[active_quad_index].quadrature_weights.begin() :
3513  nullptr);
3516  is_interior_face = other.is_interior_face;
3517  dof_access_index = other.dof_access_index;
3518 
3519  // Create deep copy of mapped geometry for use in parallel...
3520  if (other.mapped_geometry.get() != nullptr)
3521  {
3522  mapped_geometry = std::make_shared<
3523  internal::MatrixFreeFunctions::
3524  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3525  other.mapped_geometry->get_fe_values().get_mapping(),
3526  other.mapped_geometry->get_quadrature(),
3527  other.mapped_geometry->get_fe_values().get_update_flags());
3528  cell = 0;
3529  mapping_data = &mapped_geometry->get_data_storage();
3530  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3531  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3532  }
3533 
3534  return *this;
3535 }
3536 
3537 
3538 
3539 template <int dim,
3540  int n_components_,
3541  typename Number,
3542  bool is_face,
3543  typename VectorizedArrayType>
3544 inline FEEvaluationBase<dim,
3545  n_components_,
3546  Number,
3547  is_face,
3548  VectorizedArrayType>::~FEEvaluationBase()
3549 {
3550  if (matrix_info != nullptr)
3551  {
3552  try
3553  {
3554  matrix_info->release_scratch_data(scratch_data_array);
3555  }
3556  catch (...)
3557  {}
3558  }
3559  else
3560  {
3561  delete scratch_data_array;
3562  delete data;
3563  data = nullptr;
3564  }
3565  scratch_data_array = nullptr;
3566 }
3567 
3568 
3569 
3570 template <int dim,
3571  int n_components_,
3572  typename Number,
3573  bool is_face,
3574  typename VectorizedArrayType>
3575 inline void
3578 {
3579  Assert(scratch_data_array != nullptr, ExcInternalError());
3580 
3581  const unsigned int tensor_dofs_per_component =
3582  Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
3583  const unsigned int dofs_per_component =
3584  this->data->dofs_per_component_on_cell;
3585  const unsigned int n_quadrature_points =
3586  is_face ? this->data->n_q_points_face : this->data->n_q_points;
3587 
3588  const unsigned int shift =
3589  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
3590  n_components_ * 3 +
3591  2 * n_quadrature_points;
3592  const unsigned int allocated_size =
3593  shift + n_components_ * dofs_per_component +
3594  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3595  scratch_data_array->resize_fast(allocated_size);
3596 
3597  // set the pointers to the correct position in the data array
3598  for (unsigned int c = 0; c < n_components_; ++c)
3599  {
3600  this->values_dofs[c] =
3601  scratch_data_array->begin() + c * dofs_per_component;
3602  this->values_quad[c] = scratch_data_array->begin() +
3603  n_components * dofs_per_component +
3604  c * n_quadrature_points;
3605  for (unsigned int d = 0; d < dim; ++d)
3606  this->gradients_quad[c][d] =
3607  scratch_data_array->begin() +
3608  n_components * (dofs_per_component + n_quadrature_points) +
3609  (c * dim + d) * n_quadrature_points;
3610  for (unsigned int d = 0; d < (dim * dim + dim) / 2; ++d)
3611  this->hessians_quad[c][d] =
3612  scratch_data_array->begin() +
3613  n_components *
3614  ((dim + 1) * n_quadrature_points + dofs_per_component) +
3615  (c * (dim * dim + dim) + d) * n_quadrature_points;
3616  }
3617  scratch_data =
3618  scratch_data_array->begin() + n_components_ * dofs_per_component +
3619  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3620 }
3621 
3622 
3623 
3624 template <int dim,
3625  int n_components_,
3626  typename Number,
3627  bool is_face,
3628  typename VectorizedArrayType>
3629 inline unsigned int
3632 {
3633  if (matrix_info == nullptr)
3634  return 0;
3635  else
3636  {
3637  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3638  return this->mapping_data->data_index_offsets[cell];
3639  }
3640 }
3641 
3642 
3643 
3644 template <int dim,
3645  int n_components_,
3646  typename Number,
3647  bool is_face,
3648  typename VectorizedArrayType>
3651  get_cell_type() const
3652 {
3654  return cell_type;
3655 }
3656 
3657 
3658 
3659 template <int dim,
3660  int n_components_,
3661  typename Number,
3662  bool is_face,
3663  typename VectorizedArrayType>
3666  get_shape_info() const
3667 {
3668  Assert(data != nullptr, ExcInternalError());
3669  return *data;
3670 }
3671 
3672 
3673 
3674 template <int dim,
3675  int n_components_,
3676  typename Number,
3677  bool is_face,
3678  typename VectorizedArrayType>
3681  get_normal_vector(const unsigned int q_index) const
3682 {
3683  AssertIndexRange(q_index, n_quadrature_points);
3684  Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!"));
3685  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3686  return normal_vectors[0];
3687  else
3688  return normal_vectors[q_index];
3689 }
3690 
3691 
3692 
3693 template <int dim,
3694  int n_components_,
3695  typename Number,
3696  bool is_face,
3697  typename VectorizedArrayType>
3698 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3700  const unsigned int q_index) const
3701 {
3702  AssertIndexRange(q_index, n_quadrature_points);
3703  Assert(J_value != nullptr, ExcNotInitialized());
3704  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3705  {
3706  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3707  return J_value[0] * this->quadrature_weights[q_index];
3708  }
3709  else
3710  return J_value[q_index];
3711 }
3712 
3713 
3714 
3715 template <int dim,
3716  int n_components_,
3717  typename Number,
3718  bool is_face,
3719  typename VectorizedArrayType>
3722  inverse_jacobian(const unsigned int q_index) const
3723 {
3724  AssertIndexRange(q_index, n_quadrature_points);
3725  Assert(this->jacobian != nullptr, ExcNotImplemented());
3726  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3727  return jacobian[0];
3728  else
3729  return jacobian[q_index];
3730 }
3731 
3732 
3733 
3734 template <int dim,
3735  int n_components_,
3736  typename Number,
3737  bool is_face,
3738  typename VectorizedArrayType>
3739 inline std::array<unsigned int, VectorizedArrayType::size()>
3741  get_cell_ids() const
3742 {
3743  const unsigned int n_lanes = VectorizedArrayType::size();
3744  std::array<unsigned int, n_lanes> cells;
3745 
3746  // initialize array
3747  for (unsigned int i = 0; i < n_lanes; ++i)
3748  cells[i] = numbers::invalid_unsigned_int;
3749 
3750  if ((is_face == false) ||
3751  (is_face &&
3752  this->dof_access_index ==
3754  this->is_interior_face))
3755  {
3756  // cell or interior face face (element-centric loop)
3757  for (unsigned int i = 0; i < n_lanes; ++i)
3758  cells[i] = cell * n_lanes + i;
3759  }
3760  else if (is_face &&
3761  this->dof_access_index ==
3763  this->is_interior_face == false)
3764  {
3765  // exterior face (element-centric loop): for this case, we need to
3766  // look into the FaceInfo field that collects information from both
3767  // sides of a face once for the global mesh, and pick the face id that
3768  // is not the local one (cell_this).
3769  for (unsigned int i = 0; i < n_lanes; i++)
3770  {
3771  // compute actual (non vectorized) cell ID
3772  const unsigned int cell_this = this->cell * n_lanes + i;
3773  // compute face ID
3774  unsigned int face_index =
3775  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3776  this->face_no,
3777  i);
3778 
3779  if (face_index == numbers::invalid_unsigned_int)
3780  continue; // invalid face ID: no neighbor on boundary
3781 
3782  // get cell ID on both sides of face
3783  auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3784  .cells_interior[face_index % n_lanes];
3785  auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3786  .cells_exterior[face_index % n_lanes];
3787 
3788  // compare the IDs with the given cell ID
3789  if (cell_m == cell_this)
3790  cells[i] = cell_p; // neighbor has the other ID
3791  else if (cell_p == cell_this)
3792  cells[i] = cell_m;
3793  }
3794  }
3795  else if (is_face)
3796  {
3797  // face-centric faces
3798  const unsigned int *cells_ =
3799  is_interior_face ?
3800  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3801  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3802  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3803  if (cells_[i] != numbers::invalid_unsigned_int)
3804  cells[i] = cells_[i];
3805  }
3806 
3807  return cells;
3808 }
3809 
3810 
3811 namespace internal
3812 {
3813  template <int dim,
3814  int n_components_,
3815  typename Number,
3816  bool is_face,
3817  typename VectorizedArrayType,
3818  typename VectorizedArrayType2,
3819  typename GlobalVectorType,
3820  typename FU>
3821  inline void
3822  process_cell_data(
3823  const FEEvaluationBase<dim,
3824  n_components_,
3825  Number,
3826  is_face,
3827  VectorizedArrayType> & phi,
3829  GlobalVectorType & array,
3830  VectorizedArrayType2 & out,
3831  const FU & fu)
3832  {
3833  (void)matrix_info;
3834  Assert(matrix_info != nullptr, ExcNotImplemented());
3835  AssertDimension(array.size(),
3836  matrix_info->get_task_info().cell_partition_data.back());
3837 
3838  // 1) collect ids of cell
3839  const auto cells = phi.get_cell_ids();
3840 
3841  // 2) actually gather values
3842  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3843  if (cells[i] != numbers::invalid_unsigned_int)
3844  fu(out[i],
3845  array[cells[i] / VectorizedArrayType::size()]
3846  [cells[i] % VectorizedArrayType::size()]);
3847  }
3848 } // namespace internal
3849 
3850 
3851 
3852 template <int dim,
3853  int n_components_,
3854  typename Number,
3855  bool is_face,
3856  typename VectorizedArrayType>
3857 inline VectorizedArrayType
3860 {
3861  VectorizedArrayType out = Number(1.);
3862  internal::process_cell_data(
3863  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3864  local = global;
3865  });
3866  return out;
3867 }
3868 
3869 
3870 
3871 template <int dim,
3872  int n_components_,
3873  typename Number,
3874  bool is_face,
3875  typename VectorizedArrayType>
3876 inline void
3879  const VectorizedArrayType & in) const
3880 {
3881  internal::process_cell_data(
3882  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3883  global = local;
3884  });
3885 }
3886 
3887 
3888 
3889 template <int dim,
3890  int n_components_,
3891  typename Number,
3892  bool is_face,
3893  typename VectorizedArrayType>
3894 template <typename T>
3895 inline std::array<T, VectorizedArrayType::size()>
3897  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
3898  &array) const
3899 {
3900  std::array<T, VectorizedArrayType::size()> out;
3901  internal::process_cell_data(
3902  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3903  local = global;
3904  });
3905  return out;
3906 }
3907 
3908 
3909 
3910 template <int dim,
3911  int n_components_,
3912  typename Number,
3913  bool is_face,
3914  typename VectorizedArrayType>
3915 template <typename T>
3916 inline void
3919  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
3920  const std::array<T, VectorizedArrayType::size()> & in) const
3921 {
3922  internal::process_cell_data(
3923  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3924  global = local;
3925  });
3926 }
3927 
3928 
3929 
3930 namespace internal
3931 {
3932  // allows to select between block vectors and non-block vectors, which
3933  // allows to use a unified interface for extracting blocks on block vectors
3934  // and doing nothing on usual vectors
3935  template <typename VectorType, bool>
3936  struct BlockVectorSelector
3937  {};
3938 
3939  template <typename VectorType>
3940  struct BlockVectorSelector<VectorType, true>
3941  {
3942  using BaseVectorType = typename VectorType::BlockType;
3943 
3944  static BaseVectorType *
3945  get_vector_component(VectorType &vec, const unsigned int component)
3946  {
3947  AssertIndexRange(component, vec.n_blocks());
3948  return &vec.block(component);
3949  }
3950  };
3951 
3952  template <typename VectorType>
3953  struct BlockVectorSelector<VectorType, false>
3954  {
3955  using BaseVectorType = VectorType;
3956 
3957  static BaseVectorType *
3958  get_vector_component(VectorType &vec, const unsigned int component)
3959  {
3960  // FEEvaluation allows to combine several vectors from a scalar
3961  // FiniteElement into a "vector-valued" FEEvaluation object with
3962  // multiple components. These components can be extracted with the other
3963  // get_vector_component functions. If we do not get a vector of vectors
3964  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3965  // must make sure that we do not duplicate the components in input
3966  // and/or duplicate the resulting integrals. In such a case, we should
3967  // only get the zeroth component in the vector contained set nullptr for
3968  // the others which allows us to catch unintended use in
3969  // read_write_operation.
3970  if (component == 0)
3971  return &vec;
3972  else
3973  return nullptr;
3974  }
3975  };
3976 
3977  template <typename VectorType>
3978  struct BlockVectorSelector<std::vector<VectorType>, false>
3979  {
3980  using BaseVectorType = VectorType;
3981 
3982  static BaseVectorType *
3983  get_vector_component(std::vector<VectorType> &vec,
3984  const unsigned int component)
3985  {
3986  AssertIndexRange(component, vec.size());
3987  return &vec[component];
3988  }
3989  };
3990 
3991  template <typename VectorType>
3992  struct BlockVectorSelector<std::vector<VectorType *>, false>
3993  {
3994  using BaseVectorType = VectorType;
3995 
3996  static BaseVectorType *
3997  get_vector_component(std::vector<VectorType *> &vec,
3998  const unsigned int component)
3999  {
4000  AssertIndexRange(component, vec.size());
4001  return vec[component];
4002  }
4003  };
4004 } // namespace internal
4005 
4006 
4007 
4008 template <int dim,
4009  int n_components_,
4010  typename Number,
4011  bool is_face,
4012  typename VectorizedArrayType>
4013 template <typename VectorType, typename VectorOperation>
4014 inline void
4016  read_write_operation(const VectorOperation &operation,
4017  VectorType * src[],
4018  const std::bitset<VectorizedArrayType::size()> &mask,
4019  const bool apply_constraints) const
4020 {
4021  // Case 1: No MatrixFree object given, simple case because we do not need to
4022  // process constraints and need not care about vectorization -> go to
4023  // separate function
4024  if (matrix_info == nullptr)
4025  {
4026  read_write_operation_global(operation, src);
4027  return;
4028  }
4029 
4030  Assert(dof_info != nullptr, ExcNotInitialized());
4031  Assert(matrix_info->indices_initialized() == true, ExcNotInitialized());
4032  if (n_fe_components == 1)
4033  for (unsigned int comp = 0; comp < n_components; ++comp)
4034  {
4035  Assert(src[comp] != nullptr,
4036  ExcMessage("The finite element underlying this FEEvaluation "
4037  "object is scalar, but you requested " +
4038  std::to_string(n_components) +
4039  " components via the template argument in "
4040  "FEEvaluation. In that case, you must pass an "
4041  "std::vector<VectorType> or a BlockVector to " +
4042  "read_dof_values and distribute_local_to_global."));
4043  internal::check_vector_compatibility(*src[comp], *dof_info);
4044  }
4045  else
4046  {
4047  internal::check_vector_compatibility(*src[0], *dof_info);
4048  }
4049 
4050  // Case 2: contiguous indices which use reduced storage of indices and can
4051  // use vectorized load/store operations -> go to separate function
4052  AssertIndexRange(cell,
4053  dof_info->index_storage_variants[dof_access_index].size());
4054  if (dof_info->index_storage_variants
4055  [is_face ? dof_access_index :
4057  [cell] >=
4059  {
4060  read_write_operation_contiguous(operation, src, mask);
4061  return;
4062  }
4063 
4064  // Case 3: standard operation with one index per degree of freedom -> go on
4065  // here
4066  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4067  Assert(mask.count() == n_lanes,
4068  ExcNotImplemented("Masking currently not implemented for "
4069  "non-contiguous DoF storage"));
4070 
4071  std::integral_constant<bool,
4072  internal::is_vectorizable<VectorType, Number>::value>
4073  vector_selector;
4074 
4075  const unsigned int dofs_per_component =
4076  this->data->dofs_per_component_on_cell;
4077  if (dof_info->index_storage_variants
4078  [is_face ? dof_access_index :
4080  [cell] ==
4082  {
4083  const unsigned int *dof_indices =
4084  dof_info->dof_indices_interleaved.data() +
4085  dof_info->row_starts[cell * n_fe_components * n_lanes].first +
4086  dof_info->component_dof_indices_offset[active_fe_index]
4087  [first_selected_component] *
4088  n_lanes;
4089  if (n_components == 1 || n_fe_components == 1)
4090  for (unsigned int i = 0; i < dofs_per_component;
4091  ++i, dof_indices += n_lanes)
4092  for (unsigned int comp = 0; comp < n_components; ++comp)
4093  operation.process_dof_gather(dof_indices,
4094  *src[comp],
4095  0,
4096  values_dofs[comp][i],
4097  vector_selector);
4098  else
4099  for (unsigned int comp = 0; comp < n_components; ++comp)
4100  for (unsigned int i = 0; i < dofs_per_component;
4101  ++i, dof_indices += n_lanes)
4102  operation.process_dof_gather(
4103  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4104  return;
4105  }
4106 
4107  const unsigned int * dof_indices[n_lanes];
4108  VectorizedArrayType **values_dofs =
4109  const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
4110 
4111  // Assign the appropriate cell ids for face/cell case and get the pointers
4112  // to the dof indices of the cells on all lanes
4113  unsigned int cells_copied[n_lanes];
4114  const unsigned int *cells;
4115  unsigned int n_vectorization_actual =
4116  dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
4117  bool has_constraints = false;
4118  if (is_face)
4119  {
4120  if (dof_access_index ==
4122  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4123  cells_copied[v] = cell * VectorizedArrayType::size() + v;
4124  cells = dof_access_index ==
4126  &cells_copied[0] :
4127  (is_interior_face ?
4128  &this->matrix_info->get_face_info(cell).cells_interior[0] :
4129  &this->matrix_info->get_face_info(cell).cells_exterior[0]);
4130  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4131  {
4132  Assert(cells[v] < dof_info->row_starts.size() - 1,
4133  ExcInternalError());
4134  const std::pair<unsigned int, unsigned int> *my_index_start =
4135  &dof_info->row_starts[cells[v] * n_fe_components +
4136  first_selected_component];
4137 
4138  // check whether any of the SIMD lanes has constraints, i.e., the
4139  // constraint indicator which is the second entry of row_starts
4140  // increments on this cell
4141  if (my_index_start[n_components].second != my_index_start[0].second)
4142  has_constraints = true;
4143 
4144  dof_indices[v] =
4145  dof_info->dof_indices.data() + my_index_start[0].first;
4146  }
4147  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4148  dof_indices[v] = nullptr;
4149  }
4150  else
4151  {
4152  AssertIndexRange((cell + 1) * n_lanes * n_fe_components,
4153  dof_info->row_starts.size());
4154  const unsigned int n_components_read =
4155  n_fe_components > 1 ? n_components : 1;
4156  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4157  {
4158  const std::pair<unsigned int, unsigned int> *my_index_start =
4159  &dof_info->row_starts[(cell * n_lanes + v) * n_fe_components +
4160  first_selected_component];
4161  if (my_index_start[n_components_read].second !=
4162  my_index_start[0].second)
4163  has_constraints = true;
4164  Assert(my_index_start[n_components_read].first ==
4165  my_index_start[0].first ||
4166  my_index_start[0].first < dof_info->dof_indices.size(),
4167  ExcIndexRange(0,
4168  my_index_start[0].first,
4169  dof_info->dof_indices.size()));
4170  dof_indices[v] =
4171  dof_info->dof_indices.data() + my_index_start[0].first;
4172  }
4173  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4174  dof_indices[v] = nullptr;
4175  }
4176 
4177  // Case where we have no constraints throughout the whole cell: Can go
4178  // through the list of DoFs directly
4179  if (!has_constraints)
4180  {
4181  if (n_vectorization_actual < n_lanes)
4182  for (unsigned int comp = 0; comp < n_components; ++comp)
4183  for (unsigned int i = 0; i < dofs_per_component; ++i)
4184  operation.process_empty(values_dofs[comp][i]);
4185  if (n_components == 1 || n_fe_components == 1)
4186  {
4187  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4188  for (unsigned int i = 0; i < dofs_per_component; ++i)
4189  for (unsigned int comp = 0; comp < n_components; ++comp)
4190  operation.process_dof(dof_indices[v][i],
4191  *src[comp],
4192  values_dofs[comp][i][v]);
4193  }
4194  else
4195  {
4196  for (unsigned int comp = 0; comp < n_components; ++comp)
4197  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4198  for (unsigned int i = 0; i < dofs_per_component; ++i)
4199  operation.process_dof(
4200  dof_indices[v][comp * dofs_per_component + i],
4201  *src[0],
4202  values_dofs[comp][i][v]);
4203  }
4204  return;
4205  }
4206 
4207  // In the case where there are some constraints to be resolved, loop over
4208  // all vector components that are filled and then over local dofs. ind_local
4209  // holds local number on cell, index iterates over the elements of
4210  // index_local_to_global and dof_indices points to the global indices stored
4211  // in index_local_to_global
4212  if (n_vectorization_actual < n_lanes)
4213  for (unsigned int comp = 0; comp < n_components; ++comp)
4214  for (unsigned int i = 0; i < dofs_per_component; ++i)
4215  operation.process_empty(values_dofs[comp][i]);
4216  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4217  {
4218  const unsigned int cell_index = is_face ? cells[v] : cell * n_lanes + v;
4219  const unsigned int cell_dof_index =
4220  cell_index * n_fe_components + first_selected_component;
4221  const unsigned int n_components_read =
4222  n_fe_components > 1 ? n_components : 1;
4223  unsigned int index_indicators =
4224  dof_info->row_starts[cell_dof_index].second;
4225  unsigned int next_index_indicators =
4226  dof_info->row_starts[cell_dof_index + 1].second;
4227 
4228  // For read_dof_values_plain, redirect the dof_indices field to the
4229  // unconstrained indices
4230  if (apply_constraints == false &&
4231  dof_info->row_starts[cell_dof_index].second !=
4232  dof_info->row_starts[cell_dof_index + n_components_read].second)
4233  {
4234  Assert(dof_info->row_starts_plain_indices[cell_index] !=
4236  ExcNotInitialized());
4237  dof_indices[v] =
4238  dof_info->plain_dof_indices.data() +
4239  dof_info->component_dof_indices_offset[active_fe_index]
4240  [first_selected_component] +
4241  dof_info->row_starts_plain_indices[cell_index];
4242  next_index_indicators = index_indicators;
4243  }
4244 
4245  if (n_components == 1 || n_fe_components == 1)
4246  {
4247  unsigned int ind_local = 0;
4248  for (; index_indicators != next_index_indicators; ++index_indicators)
4249  {
4250  const std::pair<unsigned short, unsigned short> indicator =
4251  dof_info->constraint_indicator[index_indicators];
4252  // run through values up to next constraint
4253  for (unsigned int j = 0; j < indicator.first; ++j)
4254  for (unsigned int comp = 0; comp < n_components; ++comp)
4255  operation.process_dof(dof_indices[v][j],
4256  *src[comp],
4257  values_dofs[comp][ind_local + j][v]);
4258 
4259  ind_local += indicator.first;
4260  dof_indices[v] += indicator.first;
4261 
4262  // constrained case: build the local value as a linear
4263  // combination of the global value according to constraints
4264  Number value[n_components];
4265  for (unsigned int comp = 0; comp < n_components; ++comp)
4266  operation.pre_constraints(values_dofs[comp][ind_local][v],
4267  value[comp]);
4268 
4269  const Number *data_val =
4270  matrix_info->constraint_pool_begin(indicator.second);
4271  const Number *end_pool =
4272  matrix_info->constraint_pool_end(indicator.second);
4273  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4274  for (unsigned int comp = 0; comp < n_components; ++comp)
4275  operation.process_constraint(*dof_indices[v],
4276  *data_val,
4277  *src[comp],
4278  value[comp]);
4279 
4280  for (unsigned int comp = 0; comp < n_components; ++comp)
4281  operation.post_constraints(value[comp],
4282  values_dofs[comp][ind_local][v]);
4283  ind_local++;
4284  }
4285 
4286  AssertIndexRange(ind_local, dofs_per_component + 1);
4287 
4288  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4289  for (unsigned int comp = 0; comp < n_components; ++comp)
4290  operation.process_dof(*dof_indices[v],
4291  *src[comp],
4292  values_dofs[comp][ind_local][v]);
4293  }
4294  else
4295  {
4296  // case with vector-valued finite elements where all components are
4297  // included in one single vector. Assumption: first come all entries
4298  // to the first component, then all entries to the second one, and
4299  // so on. This is ensured by the way MatrixFree reads out the
4300  // indices.
4301  for (unsigned int comp = 0; comp < n_components; ++comp)
4302  {
4303  unsigned int ind_local = 0;
4304 
4305  // check whether there is any constraint on the current cell
4306  for (; index_indicators != next_index_indicators;
4307  ++index_indicators)
4308  {
4309  const std::pair<unsigned short, unsigned short> indicator =
4310  dof_info->constraint_indicator[index_indicators];
4311 
4312  // run through values up to next constraint
4313  for (unsigned int j = 0; j < indicator.first; ++j)
4314  operation.process_dof(dof_indices[v][j],
4315  *src[0],
4316  values_dofs[comp][ind_local + j][v]);
4317  ind_local += indicator.first;
4318  dof_indices[v] += indicator.first;
4319 
4320  // constrained case: build the local value as a linear
4321  // combination of the global value according to constraints
4322  Number value;
4323  operation.pre_constraints(values_dofs[comp][ind_local][v],
4324  value);
4325 
4326  const Number *data_val =
4327  matrix_info->constraint_pool_begin(indicator.second);
4328  const Number *end_pool =
4329  matrix_info->constraint_pool_end(indicator.second);
4330 
4331  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4332  operation.process_constraint(*dof_indices[v],
4333  *data_val,
4334  *src[0],
4335  value);
4336 
4337  operation.post_constraints(value,
4338  values_dofs[comp][ind_local][v]);
4339  ind_local++;
4340  }
4341 
4342  AssertIndexRange(ind_local, dofs_per_component + 1);
4343 
4344  // get the dof values past the last constraint
4345  for (; ind_local < dofs_per_component;
4346  ++dof_indices[v], ++ind_local)
4347  {
4348  AssertIndexRange(*dof_indices[v], src[0]->size());
4349  operation.process_dof(*dof_indices[v],
4350  *src[0],
4351  values_dofs[comp][ind_local][v]);
4352  }
4353 
4354  if (apply_constraints == true && comp + 1 < n_components)
4355  next_index_indicators =
4356  dof_info->row_starts[cell_dof_index + comp + 2].second;
4357  }
4358  }
4359  }
4360 }
4361 
4362 
4363 
4364 template <int dim,
4365  int n_components_,
4366  typename Number,
4367  bool is_face,
4368  typename VectorizedArrayType>
4369 template <typename VectorType, typename VectorOperation>
4370 inline void
4373  VectorType * src[]) const
4374 {
4375  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4376 
4377  unsigned int index =
4378  first_selected_component * data->dofs_per_component_on_cell;
4379  for (unsigned int comp = 0; comp < n_components; ++comp)
4380  {
4381  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4382  ++i, ++index)
4383  {
4384  operation.process_empty(values_dofs[comp][i]);
4385  operation.process_dof_global(
4386  local_dof_indices[data->lexicographic_numbering[index]],
4387  *src[0],
4388  values_dofs[comp][i][0]);
4389  }
4390  }
4391 }
4392 
4393 
4394 
4395 template <int dim,
4396  int n_components_,
4397  typename Number,
4398  bool is_face,
4399  typename VectorizedArrayType>
4400 template <typename VectorType, typename VectorOperation>
4401 inline void
4404  const VectorOperation & operation,
4405  VectorType * src[],
4406  const std::bitset<VectorizedArrayType::size()> &mask) const
4407 {
4408  // This functions processes the functions read_dof_values,
4409  // distribute_local_to_global, and set_dof_values with the same code for
4410  // contiguous cell indices (DG case). The distinction between these three
4411  // cases is made by the input VectorOperation that either reads values from
4412  // a vector and puts the data into the local data field or write local data
4413  // into the vector. Certain operations are no-ops for the given use case.
4414 
4415  std::integral_constant<bool,
4416  internal::is_vectorizable<VectorType, Number>::value>
4417  vector_selector;
4419  is_face ? dof_access_index :
4421  const unsigned int n_lanes = mask.count();
4422 
4423  const std::vector<unsigned int> &dof_indices_cont =
4424  dof_info->dof_indices_contiguous[ind];
4425 
4426  // Simple case: We have contiguous storage, so we can simply copy out the
4427  // data
4428  if (dof_info->index_storage_variants[ind][cell] ==
4430  interleaved_contiguous &&
4431  n_lanes == VectorizedArrayType::size())
4432  {
4433  const unsigned int dof_index =
4434  dof_indices_cont[cell * VectorizedArrayType::size()] +
4435  dof_info->component_dof_indices_offset[active_fe_index]
4436  [first_selected_component] *
4437  VectorizedArrayType::size();
4438  if (n_components == 1 || n_fe_components == 1)
4439  for (unsigned int comp = 0; comp < n_components; ++comp)
4440  operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
4441  dof_index,
4442  *src[comp],
4443  values_dofs[comp],
4444  vector_selector);
4445  else
4446  operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
4447  n_components,
4448  dof_index,
4449  *src[0],
4450  values_dofs[0],
4451  vector_selector);
4452  return;
4453  }
4454 
4455  // More general case: Must go through the components one by one and apply
4456  // some transformations
4457  const unsigned int n_filled_lanes =
4458  dof_info->n_vectorization_lanes_filled[ind][this->cell];
4459 
4460  unsigned int dof_indices[VectorizedArrayType::size()];
4461  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4462  dof_indices[v] =
4463  dof_indices_cont[cell * VectorizedArrayType::size() + v] +
4464  dof_info->component_dof_indices_offset[active_fe_index]
4465  [first_selected_component] *
4466  dof_info->dof_indices_interleave_strides
4467  [ind][cell * VectorizedArrayType::size() + v];
4468 
4469  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4470  dof_indices[v] = numbers::invalid_unsigned_int;
4471 
4472  // In the case with contiguous cell indices, we know that there are no
4473  // constraints and that the indices within each element are contiguous
4474  if (n_filled_lanes == VectorizedArrayType::size() &&
4475  n_lanes == VectorizedArrayType::size())
4476  {
4477  if (dof_info->index_storage_variants[ind][cell] ==
4479  contiguous)
4480  {
4481  if (n_components == 1 || n_fe_components == 1)
4482  for (unsigned int comp = 0; comp < n_components; ++comp)
4483  operation.process_dofs_vectorized_transpose(
4484  data->dofs_per_component_on_cell,
4485  dof_indices,
4486  *src[comp],
4487  values_dofs[comp],
4488  vector_selector);
4489  else
4490  operation.process_dofs_vectorized_transpose(
4491  data->dofs_per_component_on_cell * n_components,
4492  dof_indices,
4493  *src[0],
4494  &values_dofs[0][0],
4495  vector_selector);
4496  }
4497  else if (dof_info->index_storage_variants[ind][cell] ==
4499  interleaved_contiguous_strided)
4500  {
4501  if (n_components == 1 || n_fe_components == 1)
4502  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4503  {
4504  for (unsigned int comp = 0; comp < n_components; ++comp)
4505  operation.process_dof_gather(dof_indices,
4506  *src[comp],
4507  i * VectorizedArrayType::size(),
4508  values_dofs[comp][i],
4509  vector_selector);
4510  }
4511  else
4512  for (unsigned int comp = 0; comp < n_components; ++comp)
4513  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4514  ++i)
4515  {
4516  operation.process_dof_gather(
4517  dof_indices,
4518  *src[0],
4519  (comp * data->dofs_per_component_on_cell + i) *
4520  VectorizedArrayType::size(),
4521  values_dofs[comp][i],
4522  vector_selector);
4523  }
4524  }
4525  else
4526  {
4527  Assert(dof_info->index_storage_variants[ind][cell] ==
4529  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4530  ExcNotImplemented());
4531  const unsigned int *offsets =
4532  &dof_info->dof_indices_interleave_strides
4533  [ind][VectorizedArrayType::size() * cell];
4534  if (n_components == 1 || n_fe_components == 1)
4535  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4536  {
4537  for (unsigned int comp = 0; comp < n_components; ++comp)
4538  operation.process_dof_gather(dof_indices,
4539  *src[comp],
4540  0,
4541  values_dofs[comp][i],
4542  vector_selector);
4544  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4545  dof_indices[v] += offsets[v];
4546  }
4547  else
4548  for (unsigned int comp = 0; comp < n_components; ++comp)
4549  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4550  ++i)
4551  {
4552  operation.process_dof_gather(dof_indices,
4553  *src[0],
4554  0,
4555  values_dofs[comp][i],
4556  vector_selector);
4558  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4559  dof_indices[v] += offsets[v];
4560  }
4561  }
4562  }
4563  else
4564  for (unsigned int comp = 0; comp < n_components; ++comp)
4565  {
4566  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4567  operation.process_empty(values_dofs[comp][i]);
4568  if (dof_info->index_storage_variants[ind][cell] ==
4570  contiguous)
4571  {
4572  if (n_components == 1 || n_fe_components == 1)
4573  {
4574  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4575  if (mask[v] == true)
4576  for (unsigned int i = 0;
4577  i < data->dofs_per_component_on_cell;
4578  ++i)
4579  operation.process_dof(dof_indices[v] + i,
4580  *src[comp],
4581  values_dofs[comp][i][v]);
4582  }
4583  else
4584  {
4585  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4586  if (mask[v] == true)
4587  for (unsigned int i = 0;
4588  i < data->dofs_per_component_on_cell;
4589  ++i)
4590  operation.process_dof(
4591  dof_indices[v] + i +
4592  comp * data->dofs_per_component_on_cell,
4593  *src[0],
4594  values_dofs[comp][i][v]);
4595  }
4596  }
4597  else
4598  {
4599  const unsigned int *offsets =
4600  &dof_info->dof_indices_interleave_strides
4601  [ind][VectorizedArrayType::size() * cell];
4602  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4603  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4604  if (n_components == 1 || n_fe_components == 1)
4605  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4606  {
4607  if (mask[v] == true)
4608  for (unsigned int i = 0;
4609  i < data->dofs_per_component_on_cell;
4610  ++i)
4611  operation.process_dof(dof_indices[v] + i * offsets[v],
4612  *src[comp],
4613  values_dofs[comp][i][v]);
4614  }
4615  else
4616  {
4617  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4618  if (mask[v] == true)
4619  for (unsigned int i = 0;
4620  i < data->dofs_per_component_on_cell;
4621  ++i)
4622  operation.process_dof(
4623  dof_indices[v] +
4624  (i + comp * data->dofs_per_component_on_cell) *
4625  offsets[v],
4626  *src[0],
4627  values_dofs[comp][i][v]);
4628  }
4629  }
4630  }
4631 }
4632 
4633 
4634 
4635 template <int dim,
4636  int n_components_,
4637  typename Number,
4638  bool is_face,
4639  typename VectorizedArrayType>
4640 template <typename VectorType>
4641 inline void
4643  read_dof_values(const VectorType &src, const unsigned int first_index)
4644 {
4645  // select between block vectors and non-block vectors. Note that the number
4646  // of components is checked in the internal data
4647  typename internal::BlockVectorSelector<
4648  VectorType,
4649  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4650  for (unsigned int d = 0; d < n_components; ++d)
4651  src_data[d] =
4652  internal::BlockVectorSelector<VectorType,
4654  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4655 
4657  read_write_operation(reader,
4658  src_data,
4659  std::bitset<VectorizedArrayType::size()>().flip(),
4660  true);
4661 
4662 # ifdef DEBUG
4663  dof_values_initialized = true;
4664 # endif
4665 }
4666 
4667 
4668 
4669 template <int dim,
4670  int n_components_,
4671  typename Number,
4672  bool is_face,
4673  typename VectorizedArrayType>
4674 template <typename VectorType>
4675 inline void
4677  read_dof_values_plain(const VectorType &src, const unsigned int first_index)
4678 {
4679  // select between block vectors and non-block vectors. Note that the number
4680  // of components is checked in the internal data
4681  typename internal::BlockVectorSelector<
4682  VectorType,
4683  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4684  for (unsigned int d = 0; d < n_components; ++d)
4685  src_data[d] =
4686  internal::BlockVectorSelector<VectorType,
4688  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4689 
4691  read_write_operation(reader,
4692  src_data,
4693  std::bitset<VectorizedArrayType::size()>().flip(),
4694  false);
4695 
4696 # ifdef DEBUG
4697  dof_values_initialized = true;
4698 # endif
4699 }
4700 
4701 
4702 
4703 template <int dim,
4704  int n_components_,
4705  typename Number,
4706  bool is_face,
4707  typename VectorizedArrayType>
4708 template <typename VectorType>
4709 inline void
4712  VectorType & dst,
4713  const unsigned int first_index,
4714  const std::bitset<VectorizedArrayType::size()> &mask) const
4715 {
4716 # ifdef DEBUG
4717  Assert(dof_values_initialized == true,
4719 # endif
4720 
4721  // select between block vectors and non-block vectors. Note that the number
4722  // of components is checked in the internal data
4723  typename internal::BlockVectorSelector<
4724  VectorType,
4725  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4726  for (unsigned int d = 0; d < n_components; ++d)
4727  dst_data[d] = internal::BlockVectorSelector<
4728  VectorType,
4729  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4730  d + first_index);
4731 
4733  distributor;
4734  read_write_operation(distributor, dst_data, mask);
4735 }
4736 
4737 
4738 
4739 template <int dim,
4740  int n_components_,
4741  typename Number,
4742  bool is_face,
4743  typename VectorizedArrayType>
4744 template <typename VectorType>
4745 inline void
4748  const unsigned int first_index,
4749  const std::bitset<VectorizedArrayType::size()> &mask) const
4750 {
4751 # ifdef DEBUG
4752  Assert(dof_values_initialized == true,
4754 # endif
4755 
4756  // select between block vectors and non-block vectors. Note that the number
4757  // of components is checked in the internal data
4758  typename internal::BlockVectorSelector<
4759  VectorType,
4760  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4761  for (unsigned int d = 0; d < n_components; ++d)
4762  dst_data[d] = internal::BlockVectorSelector<
4763  VectorType,
4764  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4765  d + first_index);
4766 
4768  read_write_operation(setter, dst_data, mask);
4769 }
4770 
4771 
4772 
4773 /*------------------------------ access to data fields ----------------------*/
4774 
4775 template <int dim,
4776  int n_components,
4777  typename Number,
4778  bool is_face,
4779  typename VectorizedArrayType>
4780 inline const std::vector<unsigned int> &
4783 {
4784  return data->lexicographic_numbering;
4785 }
4786 
4787 
4788 
4789 template <int dim,
4790  int n_components,
4791  typename Number,
4792  bool is_face,
4793  typename VectorizedArrayType>
4796  get_scratch_data() const
4797 {
4799  const_cast<VectorizedArrayType *>(scratch_data),
4800  scratch_data_array->end() - scratch_data);
4801 }
4802 
4803 
4804 
4805 template <int dim,
4806  int n_components,
4807  typename Number,
4808  bool is_face,
4809  typename VectorizedArrayType>
4810 inline const VectorizedArrayType *
4812  begin_dof_values() const
4813 {
4814  return &values_dofs[0][0];
4815 }
4816 
4817 
4818 
4819 template <int dim,
4820  int n_components,
4821  typename Number,
4822  bool is_face,
4823  typename VectorizedArrayType>
4824 inline VectorizedArrayType *
4827 {
4828 # ifdef DEBUG
4829  dof_values_initialized = true;
4830 # endif
4831  return &values_dofs[0][0];
4832 }
4833 
4834 
4835 
4836 template <int dim,
4837  int n_components,
4838  typename Number,
4839  bool is_face,
4840  typename VectorizedArrayType>
4841 inline const VectorizedArrayType *
4843  begin_values() const
4844 {
4845 # ifdef DEBUG
4846  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
4847 # endif
4848  return &values_quad[0][0];
4849 }
4850 
4851 
4852 
4853 template <int dim,
4854  int n_components,
4855  typename Number,
4856  bool is_face,
4857  typename VectorizedArrayType>
4858 inline VectorizedArrayType *
4860  begin_values()
4861 {
4862 # ifdef DEBUG
4863  values_quad_initialized = true;
4864  values_quad_submitted = true;
4865 # endif
4866  return &values_quad[0][0];
4867 }
4868 
4869 
4870 
4871 template <int dim,
4872  int n_components,
4873  typename Number,
4874  bool is_face,
4875  typename VectorizedArrayType>
4876 inline const VectorizedArrayType *
4878  begin_gradients() const
4879 {
4880 # ifdef DEBUG
4881  Assert(gradients_quad_initialized || gradients_quad_submitted,
4882  ExcNotInitialized());
4883 # endif
4884  return &gradients_quad[0][0][0];
4885 }
4886 
4887 
4888 
4889 template <int dim,
4890  int n_components,
4891  typename Number,
4892  bool is_face,
4893  typename VectorizedArrayType>
4894 inline VectorizedArrayType *
4897 {
4898 # ifdef DEBUG
4899  gradients_quad_submitted = true;
4900  gradients_quad_initialized = true;
4901 # endif
4902  return &gradients_quad[0][0][0];
4903 }
4904 
4905 
4906 
4907 template <int dim,
4908  int n_components,
4909  typename Number,
4910  bool is_face,
4911  typename VectorizedArrayType>
4912 inline const VectorizedArrayType *
4914  begin_hessians() const
4915 {
4916 # ifdef DEBUG
4917  Assert(hessians_quad_initialized, ExcNotInitialized());
4918 # endif
4919  return &hessians_quad[0][0][0];
4920 }
4921 
4922 
4923 
4924 template <int dim,
4925  int n_components,
4926  typename Number,
4927  bool is_face,
4928  typename VectorizedArrayType>
4929 inline VectorizedArrayType *
4932 {
4933 # ifdef DEBUG
4934  hessians_quad_initialized = true;
4935 # endif
4936  return &hessians_quad[0][0][0];
4937 }
4938 
4939 
4940 
4941 template <int dim,
4942  int n_components_,
4943  typename Number,
4944  bool is_face,
4945  typename VectorizedArrayType>
4948  get_dof_value(const unsigned int dof) const
4949 {
4950  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4952  for (unsigned int comp = 0; comp < n_components; comp++)
4953  return_value[comp] = this->values_dofs[comp][dof];
4954  return return_value;
4955 }
4956 
4957 
4958 
4959 template <int dim,
4960  int n_components_,
4961  typename Number,
4962  bool is_face,
4963  typename VectorizedArrayType>
4966  get_value(const unsigned int q_point) const
4967 {
4968 # ifdef DEBUG
4969  Assert(this->values_quad_initialized == true,
4971 # endif
4972  AssertIndexRange(q_point, this->n_quadrature_points);
4974  for (unsigned int comp = 0; comp < n_components; comp++)
4975  return_value[comp] = this->values_quad[comp][q_point];
4976  return return_value;
4977 }
4978 
4979 
4980 
4981 template <int dim,
4982  int n_components_,
4983  typename Number,
4984  bool is_face,
4985  typename VectorizedArrayType>
4986 inline DEAL_II_ALWAYS_INLINE
4989  get_gradient(const unsigned int q_point) const
4990 {
4991 # ifdef DEBUG
4992  Assert(this->gradients_quad_initialized == true,
4994 # endif
4995  AssertIndexRange(q_point, this->n_quadrature_points);
4996 
4997  Assert(jacobian != nullptr, ExcNotInitialized());
4998 
5000 
5001  // Cartesian cell
5002  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5003  {
5004  for (unsigned int comp = 0; comp < n_components; comp++)
5005  for (unsigned int d = 0; d < dim; ++d)
5006  grad_out[comp][d] =
5007  (this->gradients_quad[comp][d][q_point] * jacobian[0][d][d]);
5008  }
5009  // cell with general/affine Jacobian
5010  else
5011  {
5013  jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5014  q_point :
5015  0];
5016  for (unsigned int comp = 0; comp < n_components; comp++)
5017  for (unsigned int d = 0; d < dim; ++d)
5018  {
5019  grad_out[comp][d] =
5020  jac[d][0] * this->gradients_quad[comp][0][q_point];
5021  for (unsigned int e = 1; e < dim; ++e)
5022  grad_out[comp][d] +=
5023  jac[d][e] * this->gradients_quad[comp][e][q_point];
5024  }
5025  }
5026  return grad_out;
5027 }
5028 
5029 
5030 
5031 template <int dim,
5032  int n_components_,
5033  typename Number,
5034  bool is_face,
5035  typename VectorizedArrayType>
5038  get_normal_derivative(const unsigned int q_point) const
5039 {
5040  AssertIndexRange(q_point, this->n_quadrature_points);
5041 # ifdef DEBUG
5042  Assert(this->gradients_quad_initialized == true,
5044 # endif
5045 
5046  Assert(normal_x_jacobian != nullptr, ExcNotInitialized());
5047 
5049  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5050  for (unsigned int comp = 0; comp < n_components; comp++)
5051  grad_out[comp] = this->gradients_quad[comp][dim - 1][q_point] *
5052  (this->normal_x_jacobian[0][dim - 1]);
5053  else
5054  {
5055  const unsigned int index =
5056  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5057  for (unsigned int comp = 0; comp < n_components; comp++)
5058  {
5059  grad_out[comp] = this->gradients_quad[comp][0][q_point] *
5060  this->normal_x_jacobian[index][0];
5061  for (unsigned int d = 1; d < dim; ++d)
5062  grad_out[comp] += this->gradients_quad[comp][d][q_point] *
5063  this->normal_x_jacobian[index][d];
5064  }
5065  }
5066  return grad_out;
5067 }
5068 
5069 
5070 
5071 namespace internal
5072 {
5073  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5074  // store the lower diagonal because of symmetry
5075  template <typename VectorizedArrayType>
5076  inline void
5077  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
5078  const VectorizedArrayType *const hessians_quad[1],
5079  const unsigned int q_point,
5080  VectorizedArrayType (&tmp)[1][1])
5081  {
5082  tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
5083  }
5084 
5085  template <typename VectorizedArrayType>
5086  inline void
5087  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
5088  const VectorizedArrayType *const hessians_quad[3],
5089  const unsigned int q_point,
5090  VectorizedArrayType (&tmp)[2][2])
5091  {
5092  for (unsigned int d = 0; d < 2; ++d)
5093  {
5094  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
5095  jac[d][1] * hessians_quad[2][q_point]);
5096  tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
5097  jac[d][1] * hessians_quad[1][q_point]);
5098  }
5099  }
5100 
5101  template <typename VectorizedArrayType>
5102  inline void
5103  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
5104  const VectorizedArrayType *const hessians_quad[6],
5105  const unsigned int q_point,
5106  VectorizedArrayType (&tmp)[3][3])
5107  {
5108  for (unsigned int d = 0; d < 3; ++d)
5109  {
5110  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
5111  jac[d][1] * hessians_quad[3][q_point] +
5112  jac[d][2] * hessians_quad[4][q_point]);
5113  tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
5114  jac[d][1] * hessians_quad[1][q_point] +
5115  jac[d][2] * hessians_quad[5][q_point]);
5116  tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
5117  jac[d][1] * hessians_quad[5][q_point] +
5118  jac[d][2] * hessians_quad[2][q_point]);
5119  }
5120  }
5121 } // namespace internal
5122 
5123 
5124 
5125 template <int dim,
5126  int n_components_,
5127  typename Number,
5128  bool is_face,
5129  typename VectorizedArrayType>
5132  get_hessian(const unsigned int q_point) const
5133 {
5134  Assert(!is_face, ExcNotImplemented());
5135 # ifdef DEBUG
5136  Assert(this->hessians_quad_initialized == true,
5138 # endif
5139  AssertIndexRange(q_point, this->n_quadrature_points);
5140 
5141  Assert(jacobian != nullptr, ExcNotImplemented());
5143  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5144  0 :
5145  q_point];
5146 
5147  Tensor<2, dim, VectorizedArrayType> hessian_out[n_components];
5148 
5149  // Cartesian cell
5150  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5151  {
5152  for (unsigned int comp = 0; comp < n_components; comp++)
5153  for (unsigned int d = 0; d < dim; ++d)
5154  {
5155  hessian_out[comp][d][d] =
5156  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
5157  switch (dim)
5158  {
5159  case 1:
5160  break;
5161  case 2:
5162  hessian_out[comp][0][1] =
5163  (this->hessians_quad[comp][2][q_point] * jac[0][0] *
5164  jac[1][1]);
5165  break;
5166  case 3:
5167  hessian_out[comp][0][1] =
5168  (this->hessians_quad[comp][3][q_point] * jac[0][0] *
5169  jac[1][1]);
5170  hessian_out[comp][0][2] =
5171  (this->hessians_quad[comp][4][q_point] * jac[0][0] *
5172  jac[2][2]);
5173  hessian_out[comp][1][2] =
5174  (this->hessians_quad[comp][5][q_point] * jac[1][1] *
5175  jac[2][2]);
5176  break;
5177  default:
5178  Assert(false, ExcNotImplemented());
5179  }
5180  for (unsigned int e = d + 1; e < dim; ++e)
5181  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5182  }
5183  }
5184  // cell with general Jacobian, but constant within the cell
5185  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5186  {
5187  for (unsigned int comp = 0; comp < n_components; comp++)
5188  {
5189  // compute laplacian before the gradient because it needs to access
5190  // unscaled gradient data
5191  VectorizedArrayType tmp[dim][dim];
5192  internal::hessian_unit_times_jac(jac,
5193  this->hessians_quad[comp],
5194  q_point,
5195  tmp);
5196 
5197  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5198  for (unsigned int d = 0; d < dim; ++d)
5199  for (unsigned int e = d; e < dim; ++e)
5200  {
5201  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5202  for (unsigned int f = 1; f < dim; ++f)
5203  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5204  }
5205 
5206  // no J' * grad(u) part here because the Jacobian is constant
5207  // throughout the cell and hence, its derivative is zero
5208 
5209  // take symmetric part
5210  for (unsigned int d = 0; d < dim; ++d)
5211  for (unsigned int e = d + 1; e < dim; ++e)
5212  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5213  }
5214  }
5215  // cell with general Jacobian
5216  else
5217  {
5218  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5219  &jac_grad =
5220  mapping_data->jacobian_gradients
5221  [1 - this->is_interior_face]
5222  [this->mapping_data->data_index_offsets[this->cell] + q_point];
5223  for (unsigned int comp = 0; comp < n_components; comp++)
5224  {
5225  // compute laplacian before the gradient because it needs to access
5226  // unscaled gradient data
5227  VectorizedArrayType tmp[dim][dim];
5228  internal::hessian_unit_times_jac(jac,
5229  this->hessians_quad[comp],
5230  q_point,
5231  tmp);
5232 
5233  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5234  for (unsigned int d = 0; d < dim; ++d)
5235  for (unsigned int e = d; e < dim; ++e)
5236  {
5237  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5238  for (unsigned int f = 1; f < dim; ++f)
5239  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5240  }
5241 
5242  // add diagonal part of J' * grad(u)
5243  for (unsigned int d = 0; d < dim; ++d)
5244  for (unsigned int e = 0; e < dim; ++e)
5245  hessian_out[comp][d][d] +=
5246  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5247 
5248  // add off-diagonal part of J' * grad(u)
5249  for (unsigned int d = 0, count = dim; d < dim; ++d)
5250  for (unsigned int e = d + 1; e < dim; ++e, ++count)
5251  for (unsigned int f = 0; f < dim; ++f)
5252  hessian_out[comp][d][e] +=
5253  (jac_grad[count][f] * this->gradients_quad[comp][f][q_point]);
5254 
5255  // take symmetric part
5256  for (unsigned int d = 0; d < dim; ++d)
5257  for (unsigned int e = d + 1; e < dim; ++e)
5258  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5259  }
5260  }
5262  hessian_out);
5263 }
5264 
5265 
5266 
5267 template <int dim,
5268  int n_components_,
5269  typename Number,
5270  bool is_face,
5271  typename VectorizedArrayType>
5274  get_hessian_diagonal(const unsigned int q_point) const
5275 {
5276  Assert(!is_face, ExcNotImplemented());
5277 # ifdef DEBUG
5278  Assert(this->hessians_quad_initialized == true,
5280 # endif
5281  AssertIndexRange(q_point, this->n_quadrature_points);
5282 
5283  Assert(jacobian != nullptr, ExcNotImplemented());
5285  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5286  0 :
5287  q_point];
5288 
5290 
5291  // Cartesian cell
5292  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5293  {
5294  for (unsigned int comp = 0; comp < n_components; comp++)
5295  for (unsigned int d = 0; d < dim; ++d)
5296  hessian_out[comp][d] =
5297  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
5298  }
5299  // cell with general Jacobian, but constant within the cell
5300  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5301  {
5302  for (unsigned int comp = 0; comp < n_components; comp++)
5303  {
5304  // compute laplacian before the gradient because it needs to access
5305  // unscaled gradient data
5306  VectorizedArrayType tmp[dim][dim];
5307  internal::hessian_unit_times_jac(jac,
5308  this->hessians_quad[comp],
5309  q_point,
5310  tmp);
5311 
5312  // compute only the trace part of hessian, J * tmp = J *
5313  // hess_unit(u) * J^T
5314  for (unsigned int d = 0; d < dim; ++d)
5315  {
5316  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5317  for (unsigned int f = 1; f < dim; ++f)
5318  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5319  }
5320  }
5321  }
5322  // cell with general Jacobian
5323  else
5324  {
5325  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5326  &jac_grad =
5327  mapping_data->jacobian_gradients
5328  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5329  for (unsigned int comp = 0; comp < n_components; comp++)
5330  {
5331  // compute laplacian before the gradient because it needs to access
5332  // unscaled gradient data
5333  VectorizedArrayType tmp[dim][dim];
5334  internal::hessian_unit_times_jac(jac,
5335  this->hessians_quad[comp],
5336  q_point,
5337  tmp);
5338 
5339  // compute only the trace part of hessian, J * tmp = J *
5340  // hess_unit(u) * J^T
5341  for (unsigned int d = 0; d < dim; ++d)
5342  {
5343  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5344  for (unsigned int f = 1; f < dim; ++f)
5345  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5346  }
5347 
5348  for (unsigned int d = 0; d < dim; ++d)
5349  for (unsigned int e = 0; e < dim; ++e)
5350  hessian_out[comp][d] +=
5351  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5352  }
5353  }
5354  return hessian_out;
5355 }
5356 
5357 
5358 
5359 template <int dim,
5360  int n_components_,
5361  typename Number,
5362  bool is_face,
5363  typename VectorizedArrayType>
5366  get_laplacian(const unsigned int q_point) const
5367 {
5368  Assert(is_face == false, ExcNotImplemented());
5369 # ifdef DEBUG
5370  Assert(this->hessians_quad_initialized == true,
5372 # endif
5373  AssertIndexRange(q_point, this->n_quadrature_points);
5374 
5377  hess_diag = get_hessian_diagonal(q_point);
5378  for (unsigned int comp = 0; comp < n_components; ++comp)
5379  {
5380  laplacian_out[comp] = hess_diag[comp][0];
5381  for (unsigned int d = 1; d < dim; ++d)
5382  laplacian_out[comp] += hess_diag[comp][d];
5383  }
5384  return laplacian_out;
5385 }
5386 
5387 
5388 
5389 template <int dim,
5390  int n_components_,
5391  typename Number,
5392  bool is_face,
5393  typename VectorizedArrayType>
5394 inline DEAL_II_ALWAYS_INLINE void
5397  const unsigned int dof)
5398 {
5399 # ifdef DEBUG
5400  this->dof_values_initialized = true;
5401 # endif
5402  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5403  for (unsigned int comp = 0; comp < n_components; comp++)
5404  this->values_dofs[comp][dof] = val_in[comp];
5405 }
5406 
5407 
5408 
5409 template <int dim,
5410  int n_components_,
5411  typename Number,
5412  bool is_face,
5413  typename VectorizedArrayType>
5414 inline DEAL_II_ALWAYS_INLINE void
5417  const unsigned int q_point)
5418 {
5420  AssertIndexRange(q_point, this->n_quadrature_points);
5421  Assert(this->J_value != nullptr, ExcNotInitialized());
5422 # ifdef DEBUG
5423  this->values_quad_submitted = true;
5424 # endif
5425 
5426  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5427  {
5428  const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
5429  for (unsigned int comp = 0; comp < n_components; ++comp)
5430  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5431  }
5432  else
5433  {
5434  const VectorizedArrayType JxW = J_value[q_point];
5435  for (unsigned int comp = 0; comp < n_components; ++comp)
5436  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5437  }
5438 }
5439 
5440 
5441 
5442 template <int dim,
5443  int n_components_,
5444  typename Number,
5445  bool is_face,
5446  typename VectorizedArrayType>
5447 inline DEAL_II_ALWAYS_INLINE void
5450  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5451  const unsigned int q_point)
5452 {
5454  AssertIndexRange(q_point, this->n_quadrature_points);
5455  Assert(this->J_value != nullptr, ExcNotInitialized());
5456  Assert(this->jacobian != nullptr, ExcNotInitialized());
5457 # ifdef DEBUG
5458  this->gradients_quad_submitted = true;
5459 # endif
5460 
5461  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5462  {
5463  const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
5464  for (unsigned int comp = 0; comp < n_components; comp++)
5465  for (unsigned int d = 0; d < dim; ++d)
5466  this->gradients_quad[comp][d][q_point] =
5467  (grad_in[comp][d] * jacobian[0][d][d] * JxW);
5468  }
5469  else
5470  {
5472  this->cell_type > internal::MatrixFreeFunctions::affine ?
5473  jacobian[q_point] :
5474  jacobian[0];
5475  const VectorizedArrayType JxW =
5476  this->cell_type > internal::MatrixFreeFunctions::affine ?
5477  J_value[q_point] :
5478  J_value[0] * quadrature_weights[q_point];
5479  for (unsigned int comp = 0; comp < n_components; ++comp)
5480  for (unsigned int d = 0; d < dim; ++d)
5481  {
5482  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5483  for (unsigned int e = 1; e < dim; ++e)
5484  new_val += (jac[e][d] * grad_in[comp][e]);
5485  this->gradients_quad[comp][d][q_point] = new_val * JxW;
5486  }
5487  }
5488 }
5489 
5490 
5491 
5492 template <int dim,
5493  int n_components_,
5494  typename Number,
5495  bool is_face,
5496  typename VectorizedArrayType>
5497 inline DEAL_II_ALWAYS_INLINE void
5501  const unsigned int q_point)
5502 {
5503  AssertIndexRange(q_point, this->n_quadrature_points);
5504  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
5505 # ifdef DEBUG
5506  this->gradients_quad_submitted = true;
5507 # endif
5508 
5509  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5510  for (unsigned int comp = 0; comp < n_components; comp++)
5511  {
5512  for (unsigned int d = 0; d < dim - 1; ++d)
5513  this->gradients_quad[comp][d][q_point] = VectorizedArrayType();
5514  this->gradients_quad[comp][dim - 1][q_point] =
5515  grad_in[comp] *
5516  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5517  this->quadrature_weights[q_point]);
5518  }
5519  else
5520  {
5521  const unsigned int index =
5522  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5523  for (unsigned int comp = 0; comp < n_components; comp++)
5524  {
5525  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5526  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5527  factor = factor * this->quadrature_weights[q_point];
5528  for (unsigned int d = 0; d < dim; ++d)
5529  this->gradients_quad[comp][d][q_point] =
5530  factor * this->normal_x_jacobian[index][d];
5531  }
5532  }
5533 }
5534 
5535 
5536 
5537 template <int dim,
5538  int n_components_,
5539  typename Number,
5540  bool is_face,
5541  typename VectorizedArrayType>
5544  integrate_value() const
5545 {
5547 # ifdef DEBUG
5548  Assert(this->values_quad_submitted == true,
5550 # endif
5552  for (unsigned int comp = 0; comp < n_components; ++comp)
5553  return_value[comp] = this->values_quad[comp][0];
5554  const unsigned int n_q_points = this->n_quadrature_points;
5555  for (unsigned int q = 1; q < n_q_points; ++q)
5556  for (unsigned int comp = 0; comp < n_components; ++comp)
5557  return_value[comp] += this->values_quad[comp][q];
5558  return (return_value);
5559 }
5560 
5561 
5562 
5563 /*----------------------- FEEvaluationAccess --------------------------------*/
5564 
5565 
5566 template <int dim,
5567  int n_components_,
5568  typename Number,
5569  bool is_face,
5570  typename VectorizedArrayType>
5571 inline FEEvaluationAccess<dim,
5572  n_components_,
5573  Number,
5574  is_face,
5575  VectorizedArrayType>::
5576  FEEvaluationAccess(
5578  const unsigned int dof_no,
5579  const unsigned int first_selected_component,
5580  const unsigned int quad_no_in,
5581  const unsigned int fe_degree,
5582  const unsigned int n_q_points,
5583  const bool is_interior_face)
5585  data_in,
5586  dof_no,
5587  first_selected_component,
5588  quad_no_in,
5589  fe_degree,
5590  n_q_points,
5591  is_interior_face)
5592 {}
5593 
5594 
5595 
5596 template <int dim,
5597  int n_components_,
5598  typename Number,
5599  bool is_face,
5600  typename VectorizedArrayType>
5601 template <int n_components_other>
5602 inline FEEvaluationAccess<dim,
5603  n_components_,
5604  Number,
5605  is_face,
5606  VectorizedArrayType>::
5607  FEEvaluationAccess(const Mapping<dim> & mapping,
5608  const FiniteElement<dim> &fe,
5609  const Quadrature<1> & quadrature,
5610  const UpdateFlags update_flags,
5611  const unsigned int first_selected_component,
5612  const FEEvaluationBase<dim,
5613  n_components_other,
5614  Number,
5615  is_face,
5616  VectorizedArrayType> *other)
5618  mapping,
5619  fe,
5620  quadrature,
5621  update_flags,
5622  first_selected_component,
5623  other)
5624 {}
5625 
5626 
5627 
5628 template <int dim,
5629  int n_components_,
5630  typename Number,
5631  bool is_face,
5632  typename VectorizedArrayType>
5633 inline FEEvaluationAccess<dim,
5634  n_components_,
5635  Number,
5636  is_face,
5637  VectorizedArrayType>::
5638  FEEvaluationAccess(const FEEvaluationAccess<dim,
5639  n_components_,
5640  Number,
5641  is_face,
5642  VectorizedArrayType> &other)
5644  other)
5645 {}
5646 
5647 
5648 
5649 template <int dim,
5650  int n_components_,
5651  typename Number,
5652  bool is_face,
5653  typename VectorizedArrayType>
5654 inline FEEvaluationAccess<dim,
5655  n_components_,
5656  Number,
5657  is_face,
5658  VectorizedArrayType> &
5660 operator=(const FEEvaluationAccess<dim,
5661  n_components_,
5662  Number,
5663  is_face,
5664  VectorizedArrayType> &other)
5665 {
5666  this->FEEvaluationBase<dim,
5667  n_components_,
5668  Number,
5669  is_face,
5670  VectorizedArrayType>::operator=(other);
5671  return *this;
5672 }
5673 
5674 
5675 
5676 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5677 
5678 
5679 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5683  const unsigned int dof_no,
5684  const unsigned int first_selected_component,
5685  const unsigned int quad_no_in,
5686  const unsigned int fe_degree,
5687  const unsigned int n_q_points,
5688  const bool is_interior_face)
5690  data_in,
5691  dof_no,
5692  first_selected_component,
5693  quad_no_in,
5694  fe_degree,
5695  n_q_points,
5696  is_interior_face)
5697 {}
5698 
5699 
5700 
5701 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5702 template <int n_components_other>
5704  FEEvaluationAccess(const Mapping<dim> & mapping,
5705  const FiniteElement<dim> &fe,
5706  const Quadrature<1> & quadrature,
5707  const UpdateFlags update_flags,
5708  const unsigned int first_selected_component,
5709  const FEEvaluationBase<dim,
5710  n_components_other,
5711  Number,
5712  is_face,
5713  VectorizedArrayType> *other)
5715  mapping,
5716  fe,
5717  quadrature,
5718  update_flags,
5719  first_selected_component,
5720  other)
5721 {}
5722 
5723 
5724 
5725 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5729  &other)
5731 {}
5732 
5733 
5734 
5735 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5739 {
5741  operator=(other);
5742  return *this;
5743 }
5744 
5745 
5746 
5747 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5748 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5750  const unsigned int dof) const
5751 {
5752  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5753  return this->values_dofs[0][dof];
5754 }
5755 
5756 
5757 
5758 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5759 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5761  const unsigned int q_point) const
5762 {
5763 # ifdef DEBUG
5764  Assert(this->values_quad_initialized == true,
5766 # endif
5767  AssertIndexRange(q_point, this->n_quadrature_points);
5768  return this->values_quad[0][q_point];
5769 }
5770 
5771 
5772 
5773 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5774 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5776  get_normal_derivative(const unsigned int q_point) const
5777 {
5778  return BaseClass::get_normal_derivative(q_point)[0];
5779 }
5780 
5781 
5782 
5783 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5786  const unsigned int q_point) const
5787 {
5788  // could use the base class gradient, but that involves too many expensive
5789  // initialization operations on tensors
5790 
5791 # ifdef DEBUG
5792  Assert(this->gradients_quad_initialized == true,
5794 # endif
5795  AssertIndexRange(q_point, this->n_quadrature_points);
5796 
5797  Assert(this->jacobian != nullptr, ExcNotInitialized());
5798 
5800 
5801  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5802  {
5803  for (unsigned int d = 0; d < dim; ++d)
5804  grad_out[d] =
5805  (this->gradients_quad[0][d][q_point] * this->jacobian[0][d][d]);
5806  }
5807  // cell with general/affine Jacobian
5808  else
5809  {
5811  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5812  q_point :
5813  0];
5814  for (unsigned int d = 0; d < dim; ++d)
5815  {
5816  grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point];
5817  for (unsigned int e = 1; e < dim; ++e)
5818  grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point];
5819  }
5820  }
5821  return grad_out;
5822 }
5823 
5824 
5825 
5826 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5829  const unsigned int q_point) const
5830 {
5831  return BaseClass::get_hessian(q_point)[0];
5832 }
5833 
5834 
5835 
5836 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5839  get_hessian_diagonal(const unsigned int q_point) const
5840 {
5841  return BaseClass::get_hessian_diagonal(q_point)[0];
5842 }
5843 
5844 
5845 
5846 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5847 inline VectorizedArrayType
5849  const unsigned int q_point) const
5850 {
5851  return BaseClass::get_laplacian(q_point)[0];
5852 }
5853 
5854 
5855 
5856 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5857 inline void DEAL_II_ALWAYS_INLINE
5859  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5860 {
5861 # ifdef DEBUG
5862  this->dof_values_initialized = true;
5863  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5864 # endif
5865  this->values_dofs[0][dof] = val_in;
5866 }
5867 
5868 
5869 
5870 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5871 inline void DEAL_II_ALWAYS_INLINE
5873  const VectorizedArrayType val_in,
5874  const unsigned int q_index)
5875 {
5877  AssertIndexRange(q_index, this->n_quadrature_points);
5878  Assert(this->J_value != nullptr, ExcNotInitialized());
5879 # ifdef DEBUG
5880  this->values_quad_submitted = true;
5881 # endif
5882 
5883  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5884  {
5885  const VectorizedArrayType JxW =
5886  this->J_value[0] * this->quadrature_weights[q_index];
5887  this->values_quad[0][q_index] = val_in * JxW;
5888  }
5889  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5890  {
5891  this->values_quad[0][q_index] = val_in * this->J_value[q_index];
5892  }
5893 }
5894 
5895 
5896 
5897 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5898 inline DEAL_II_ALWAYS_INLINE void
5900  const Tensor<1, 1, VectorizedArrayType> val_in,
5901  const unsigned int q_point)
5902 {
5903  submit_value(val_in[0], q_point);
5904 }
5905 
5906 
5907 
5908 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5909 inline DEAL_II_ALWAYS_INLINE void
5911  submit_normal_derivative(const VectorizedArrayType grad_in,
5912  const unsigned int q_point)
5913 {
5915  grad[0] = grad_in;
5916  BaseClass::submit_normal_derivative(grad, q_point);
5917 }
5918 
5919 
5920 
5921 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5922 inline DEAL_II_ALWAYS_INLINE void
5925  const unsigned int q_index)
5926 {
5928  AssertIndexRange(q_index, this->n_quadrature_points);
5929  Assert(this->J_value != nullptr, ExcNotInitialized());
5930  Assert(this->jacobian != nullptr, ExcNotInitialized());
5931 # ifdef DEBUG
5932  this->gradients_quad_submitted = true;
5933 # endif
5934 
5935  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5936  {
5937  const VectorizedArrayType JxW =
5938  this->J_value[0] * this->quadrature_weights[q_index];
5939  for (unsigned int d = 0; d < dim; ++d)
5940  this->gradients_quad[0][d][q_index] =
5941  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5942  }
5943  // general/affine cell type
5944  else
5945  {
5947  this->cell_type > internal::MatrixFreeFunctions::affine ?
5948  this->jacobian[q_index] :
5949  this->jacobian[0];
5950  const VectorizedArrayType JxW =
5951  this->cell_type > internal::MatrixFreeFunctions::affine ?
5952  this->J_value[q_index] :
5953  this->J_value[0] * this->quadrature_weights[q_index];
5954  for (unsigned int d = 0; d < dim; ++d)
5955  {
5956  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5957  for (unsigned int e = 1; e < dim; ++e)
5958  new_val += jac[e][d] * grad_in[e];
5959  this->gradients_quad[0][d][q_index] = new_val * JxW;
5960  }
5961  }
5962 }
5963 
5964 
5965 
5966 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5967 inline VectorizedArrayType
5969  integrate_value() const
5970 {
5971  return BaseClass::integrate_value()[0];
5972 }
5973 
5974 
5975 
5976 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5977 
5978 
5979 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5983  const unsigned int dof_no,
5984  const unsigned int first_selected_component,
5985  const unsigned int quad_no_in,
5986  const unsigned int fe_degree,
5987  const unsigned int n_q_points,
5988  const bool is_interior_face)
5990  data_in,
5991  dof_no,
5992  first_selected_component,
5993  quad_no_in,
5994  fe_degree,
5995  n_q_points,
5996  is_interior_face)
5997 {}
5998 
5999 
6000 
6001 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6002 template <int n_components_other>
6004  FEEvaluationAccess(const Mapping<dim> & mapping,
6005  const FiniteElement<dim> &fe,
6006  const Quadrature<1> & quadrature,
6007  const UpdateFlags update_flags,
6008  const unsigned int first_selected_component,
6009  const FEEvaluationBase<dim,
6010  n_components_other,
6011  Number,
6012  is_face,
6013  VectorizedArrayType> *other)
6015  mapping,
6016  fe,
6017  quadrature,
6018  update_flags,
6019  first_selected_component,
6020  other)
6021 {}
6022 
6023 
6024 
6025 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6029  &other)
6031 {}
6032 
6033 
6034 
6035 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6039  &other)
6040 {
6042  operator=(other);
6043  return *this;
6044 }
6045 
6046 
6047 
6048 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6051  get_gradient(const unsigned int q_point) const
6052 {
6053  return BaseClass::get_gradient(q_point);
6054 }
6055 
6056 
6057 
6058 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6059 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6061  get_divergence(const unsigned int q_point) const
6062 {
6063 # ifdef DEBUG
6064  Assert(this->gradients_quad_initialized == true,
6066 # endif
6067  AssertIndexRange(q_point, this->n_quadrature_points);
6068  Assert(this->jacobian != nullptr, ExcNotInitialized());
6069 
6070  VectorizedArrayType divergence;
6071 
6072  // Cartesian cell
6073  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6074  {
6075  divergence =
6076  (this->gradients_quad[0][0][q_point] * this->jacobian[0][0][0]);
6077  for (unsigned int d = 1; d < dim; ++d)
6078  divergence +=
6079  (this->gradients_quad[d][d][q_point] * this->jacobian[0][d][d]);
6080  }
6081  // cell with general/constant Jacobian
6082  else
6083  {
6085  this->cell_type == internal::MatrixFreeFunctions::general ?
6086  this->jacobian[q_point] :
6087  this->jacobian[0];
6088  divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
6089  for (unsigned int e = 1; e < dim; ++e)
6090  divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
6091  for (unsigned int d = 1; d < dim; ++d)
6092  for (unsigned int e = 0; e < dim; ++e)
6093  divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
6094  }
6095  return divergence;
6096 }
6097 
6098 
6099 
6100 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6103  get_symmetric_gradient(const unsigned int q_point) const
6104 {
6105  // copy from generic function into dim-specialization function
6106  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6107  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6108  VectorizedArrayType half = Number(0.5);
6109  for (unsigned int d = 0; d < dim; ++d)
6110  symmetrized[d] = grad[d][d];
6111  switch (dim)
6112  {
6113  case 1:
6114  break;
6115  case 2:
6116  symmetrized[2] = grad[0][1] + grad[1][0];
6117  symmetrized[2] *= half;
6118  break;
6119  case 3:
6120  symmetrized[3] = grad[0][1] + grad[1][0];
6121  symmetrized[3] *= half;
6122  symmetrized[4] = grad[0][2] + grad[2][0];
6123  symmetrized[4] *= half;
6124  symmetrized[5] = grad[1][2] + grad[2][1];
6125  symmetrized[5] *= half;
6126  break;
6127  default:
6128  Assert(false, ExcNotImplemented());
6129  }
6131 }
6132 
6133 
6134 
6135 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6136 inline DEAL_II_ALWAYS_INLINE
6139  const unsigned int q_point) const
6140 {
6141  // copy from generic function into dim-specialization function
6142  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6144  switch (dim)
6145  {
6146  case 1:
6147  Assert(false,
6148  ExcMessage(
6149  "Computing the curl in 1d is not a useful operation"));
6150  break;
6151  case 2:
6152  curl[0] = grad[1][0] - grad[0][1];
6153  break;
6154  case 3:
6155  curl[0] = grad[2][1] - grad[1][2];
6156  curl[1] = grad[0][2] - grad[2][0];
6157  curl[2] = grad[1][0] - grad[0][1];
6158  break;
6159  default:
6160  Assert(false, ExcNotImplemented());
6161  }
6162  return curl;
6163 }
6164 
6165 
6166 
6167 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6170  get_hessian_diagonal(const unsigned int q_point) const
6171 {
6172  return BaseClass::get_hessian_diagonal(q_point);
6173 }
6174 
6175 
6176 
6177 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6180  const unsigned int q_point) const
6181 {
6182 # ifdef DEBUG
6183  Assert(this->hessians_quad_initialized == true,
6185 # endif
6186  AssertIndexRange(q_point, this->n_quadrature_points);
6187  return BaseClass::get_hessian(q_point);
6188 }
6189 
6190 
6191 
6192 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6193 inline DEAL_II_ALWAYS_INLINE void
6196  const unsigned int q_point)
6197 {
6198  BaseClass::submit_gradient(grad_in, q_point);
6199 }
6200 
6201 
6202 
6203 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6204 inline DEAL_II_ALWAYS_INLINE void
6207  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6208  const unsigned int q_point)
6209 {
6210  BaseClass::submit_gradient(grad_in, q_point);
6211 }
6212 
6213 
6214 
6215 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6216 inline DEAL_II_ALWAYS_INLINE void
6218  submit_divergence(const VectorizedArrayType div_in,
6219  const unsigned int q_point)
6220 {
6222  AssertIndexRange(q_point, this->n_quadrature_points);
6223  Assert(this->J_value != nullptr, ExcNotInitialized());
6224  Assert(this->jacobian != nullptr, ExcNotInitialized());
6225 # ifdef DEBUG
6226  this->gradients_quad_submitted = true;
6227 # endif
6228 
6229  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6230  {
6231  const VectorizedArrayType fac =
6232  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6233  for (unsigned int d = 0; d < dim; ++d)
6234  {
6235  this->gradients_quad[d][d][q_point] = (fac * this->jacobian[0][d][d]);
6236  for (unsigned int e = d + 1; e < dim; ++e)
6237  {
6238  this->gradients_quad[d][e][q_point] = VectorizedArrayType();
6239  this->gradients_quad[e][d][q_point] = VectorizedArrayType();
6240  }
6241  }
6242  }
6243  else
6244  {
6246  this->cell_type == internal::MatrixFreeFunctions::general ?
6247  this->jacobian[q_point] :
6248  this->jacobian[0];
6249  const VectorizedArrayType fac =
6250  (this->cell_type == internal::MatrixFreeFunctions::general ?
6251  this->J_value[q_point] :
6252  this->J_value[0] * this->quadrature_weights[q_point]) *
6253  div_in;
6254  for (unsigned int d = 0; d < dim; ++d)
6255  {
6256  for (unsigned int e = 0; e < dim; ++e)
6257  this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
6258  }
6259  }
6260 }
6261 
6262 
6263 
6264 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6265 inline DEAL_II_ALWAYS_INLINE void
6269  const unsigned int q_point)
6270 {
6271  // could have used base class operator, but that involves some overhead
6272  // which is inefficient. it is nice to have the symmetric tensor because
6273  // that saves some operations
6275  AssertIndexRange(q_point, this->n_quadrature_points);
6276  Assert(this->J_value != nullptr, ExcNotInitialized());
6277  Assert(this->jacobian != nullptr, ExcNotInitialized());
6278 # ifdef DEBUG
6279  this->gradients_quad_submitted = true;
6280 # endif
6281 
6282  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6283  {
6284  const VectorizedArrayType JxW =
6285  this->J_value[0] * this->quadrature_weights[q_point];
6286  for (unsigned int d = 0; d < dim; ++d)
6287  this->gradients_quad[d][d][q_point] =
6288  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6289  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6290  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6291  {
6292  const VectorizedArrayType value =
6293  sym_grad.access_raw_entry(counter) * JxW;
6294  this->gradients_quad[e][d][q_point] =
6295  (value * this->jacobian[0][d][d]);
6296  this->gradients_quad[d][e][q_point] =
6297  (value * this->jacobian[0][e][e]);
6298  }
6299  }
6300  // general/affine cell type
6301  else
6302  {
6303  const VectorizedArrayType JxW =
6304  this->cell_type == internal::MatrixFreeFunctions::general ?
6305  this->J_value[q_point] :
6306  this->J_value[0] * this->quadrature_weights[q_point];
6308  this->cell_type == internal::MatrixFreeFunctions::general ?
6309  this->jacobian[q_point] :
6310  this->jacobian[0];
6311  VectorizedArrayType weighted[dim][dim];
6312  for (unsigned int i = 0; i < dim; ++i)
6313  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6314  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6315  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6316  {
6317  const VectorizedArrayType value =
6318  sym_grad.access_raw_entry(counter) * JxW;
6319  weighted[i][j] = value;
6320  weighted[j][i] = value;
6321  }
6322  for (unsigned int comp = 0; comp < dim; ++comp)
6323  for (unsigned int d = 0; d < dim; ++d)
6324  {
6325  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6326  for (unsigned int e = 1; e < dim; ++e)
6327  new_val += jac[e][d] * weighted[comp][e];
6328  this->gradients_quad[comp][d][q_point] = new_val;
6329  }
6330  }
6331 }
6332 
6333 
6334 
6335 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6336 inline DEAL_II_ALWAYS_INLINE void
6339  const unsigned int q_point)
6340 {
6342  switch (dim)
6343  {
6344  case 1:
6345  Assert(false,
6346  ExcMessage(
6347  "Testing by the curl in 1d is not a useful operation"));
6348  break;
6349  case 2:
6350  grad[1][0] = curl[0];
6351  grad[0][1] = -curl[0];
6352  break;
6353  case 3:
6354  grad[2][1] = curl[0];
6355  grad[1][2] = -curl[0];
6356  grad[0][2] = curl[1];
6357  grad[2][0] = -curl[1];
6358  grad[1][0] = curl[2];
6359  grad[0][1] = -curl[2];
6360  break;
6361  default:
6362  Assert(false, ExcNotImplemented());
6363  }
6364  submit_gradient(grad, q_point);
6365 }
6366 
6367 
6368 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6369 
6370 
6371 template <typename Number, bool is_face, typename VectorizedArrayType>
6374  const unsigned int dof_no,
6375  const unsigned int first_selected_component,
6376  const unsigned int quad_no_in,
6377  const unsigned int fe_degree,
6378  const unsigned int n_q_points,
6379  const bool is_interior_face)
6381  data_in,
6382  dof_no,
6383  first_selected_component,
6384  quad_no_in,
6385  fe_degree,
6386  n_q_points,
6387  is_interior_face)
6388 {}
6389 
6390 
6391 
6392 template <typename Number, bool is_face, typename VectorizedArrayType>
6393 template <int n_components_other>
6395  FEEvaluationAccess(const Mapping<1> & mapping,
6396  const FiniteElement<1> &fe,
6397  const Quadrature<1> & quadrature,
6398  const UpdateFlags update_flags,
6399  const unsigned int first_selected_component,
6400  const FEEvaluationBase<1,
6401  n_components_other,
6402  Number,
6403  is_face,
6404  VectorizedArrayType> *other)
6406  mapping,
6407  fe,
6408  quadrature,
6409  update_flags,
6410  first_selected_component,
6411  other)
6412 {}
6413 
6414 
6415 
6416 template <typename Number, bool is_face, typename VectorizedArrayType>
6421 {}
6422 
6423 
6424 
6425 template <typename Number, bool is_face, typename VectorizedArrayType>
6429 {
6431  other);
6432  return *this;
6433 }
6434 
6435 
6436 
6437 template <typename Number, bool is_face, typename VectorizedArrayType>
6438 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6440  const unsigned int dof) const
6441 {
6442  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6443  return this->values_dofs[0][dof];
6444 }
6445 
6446 
6447 
6448 template <typename Number, bool is_face, typename VectorizedArrayType>
6449 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6451  const unsigned int q_point) const
6452 {
6453 # ifdef DEBUG
6454  Assert(this->values_quad_initialized == true,
6456 # endif
6457  AssertIndexRange(q_point, this->n_quadrature_points);
6458  return this->values_quad[0][q_point];
6459 }
6460 
6461 
6462 
6463 template <typename Number, bool is_face, typename VectorizedArrayType>
6466  const unsigned int q_point) const
6467 {
6468  // could use the base class gradient, but that involves too many inefficient
6469  // initialization operations on tensors
6470 
6471 # ifdef DEBUG
6472  Assert(this->gradients_quad_initialized == true,
6474 # endif
6475  AssertIndexRange(q_point, this->n_quadrature_points);
6476 
6478  this->cell_type == internal::MatrixFreeFunctions::general ?
6479  this->jacobian[q_point] :
6480  this->jacobian[0];
6481 
6483  grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];
6484 
6485  return grad_out;
6486 }
6487 
6488 
6489 
6490 template <typename Number, bool is_face, typename VectorizedArrayType>
6491 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6493  get_normal_derivative(const unsigned int q_point) const
6494 {
6495  return BaseClass::get_normal_derivative(q_point)[0];
6496 }
6497 
6498 
6499 
6500 template <typename Number, bool is_face, typename VectorizedArrayType>
6503  const unsigned int q_point) const
6504 {
6505  return BaseClass::get_hessian(q_point)[0];
6506 }
6507 
6508 
6509 
6510 template <typename Number, bool is_face, typename VectorizedArrayType>
6513  get_hessian_diagonal(const unsigned int q_point) const
6514 {
6515  return BaseClass::get_hessian_diagonal(q_point)[0];
6516 }
6517 
6518 
6519 
6520 template <typename Number, bool is_face, typename VectorizedArrayType>
6521 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6523  const unsigned int q_point) const
6524 {
6525  return BaseClass::get_laplacian(q_point)[0];
6526 }
6527 
6528 
6529 
6530 template <typename Number, bool is_face, typename VectorizedArrayType>
6533  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6534 {
6535 # ifdef DEBUG
6536  this->dof_values_initialized = true;
6537  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6538 # endif
6539  this->values_dofs[0][dof] = val_in;
6540 }
6541 
6542 
6543 
6544 template <typename Number, bool is_face, typename VectorizedArrayType>
6545 inline DEAL_II_ALWAYS_INLINE void
6547  const VectorizedArrayType val_in,
6548  const unsigned int q_point)
6549 {
6551  AssertIndexRange(q_point, this->n_quadrature_points);
6552 # ifdef DEBUG
6553  this->values_quad_submitted = true;
6554 # endif
6555 
6556  if (this->cell_type == internal::MatrixFreeFunctions::general)
6557  {
6558  const VectorizedArrayType JxW = this->J_value[q_point];
6559  this->values_quad[0][q_point] = val_in * JxW;
6560  }
6561  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6562  {
6563  const VectorizedArrayType JxW =
6564  this->J_value[0] * this->quadrature_weights[q_point];
6565  this->values_quad[0][q_point] = val_in * JxW;
6566  }
6567 }
6568 
6569 
6570 
6571 template <typename Number, bool is_face, typename VectorizedArrayType>
6572 inline DEAL_II_ALWAYS_INLINE void
6574  const Tensor<1, 1, VectorizedArrayType> val_in,
6575  const unsigned int q_point)
6576 {
6577  submit_value(val_in[0], q_point);
6578 }
6579 
6580 
6581 
6582 template <typename Number, bool is_face, typename VectorizedArrayType>
6583 inline DEAL_II_ALWAYS_INLINE void
6585  const Tensor<1, 1, VectorizedArrayType> grad_in,
6586  const unsigned int q_point)
6587 {
6588  submit_gradient(grad_in[0], q_point);
6589 }
6590 
6591 
6592 
6593 template <typename Number, bool is_face, typename VectorizedArrayType>
6594 inline DEAL_II_ALWAYS_INLINE void
6596  const VectorizedArrayType grad_in,
6597  const unsigned int q_point)
6598 {
6600  AssertIndexRange(q_point, this->n_quadrature_points);
6601 # ifdef DEBUG
6602  this->gradients_quad_submitted = true;
6603 # endif
6604 
6606  this->cell_type == internal::MatrixFreeFunctions::general ?
6607  this->jacobian[q_point] :
6608  this->jacobian[0];
6609  const VectorizedArrayType JxW =
6610  this->cell_type == internal::MatrixFreeFunctions::general ?
6611  this->J_value[q_point] :
6612  this->J_value[0] * this->quadrature_weights[q_point];
6613 
6614  this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in * JxW;
6615 }
6616 
6617 
6618 
6619 template <typename Number, bool is_face, typename VectorizedArrayType>
6620 inline DEAL_II_ALWAYS_INLINE void
6622  submit_normal_derivative(const VectorizedArrayType grad_in,
6623  const unsigned int q_point)
6624 {
6626  grad[0] = grad_in;
6627  BaseClass::submit_normal_derivative(grad, q_point);
6628 }
6629 
6630 
6631 
6632 template <typename Number, bool is_face, typename VectorizedArrayType>
6633 inline DEAL_II_ALWAYS_INLINE void
6636  const unsigned int q_point)
6637 {
6638  BaseClass::submit_normal_derivative(grad_in, q_point);
6639 }
6640 
6641 
6642 
6643 template <typename Number, bool is_face, typename VectorizedArrayType>
6644 inline VectorizedArrayType
6646  integrate_value() const
6647 {
6648  return BaseClass::integrate_value()[0];
6649 }
6650 
6651 
6652 
6653 /*-------------------------- FEEvaluation -----------------------------------*/
6654 
6655 
6656 template <int dim,
6657  int fe_degree,
6658  int n_q_points_1d,
6659  int n_components_,
6660  typename Number,
6661  typename VectorizedArrayType>
6662 inline FEEvaluation<dim,
6663  fe_degree,
6664  n_q_points_1d,
6665  n_components_,
6666  Number,
6667  VectorizedArrayType>::
6668  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
6669  const unsigned int fe_no,
6670  const unsigned int quad_no,
6671  const unsigned int first_selected_component)
6672  : BaseClass(data_in,
6673  fe_no,
6674  first_selected_component,
6675  quad_no,
6676  fe_degree,
6677  static_n_q_points)
6678  , dofs_per_component(this->data->dofs_per_component_on_cell)
6679  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6680  , n_q_points(this->data->n_q_points)
6681 {
6682  check_template_arguments(fe_no, 0);
6683 }
6684 
6685 
6686 
6687 template <int dim,
6688  int fe_degree,
6689  int n_q_points_1d,
6690  int n_components_,
6691  typename Number,
6692  typename VectorizedArrayType>
6693 inline FEEvaluation<dim,
6694  fe_degree,
6695  n_q_points_1d,
6696  n_components_,
6697  Number,
6698  VectorizedArrayType>::
6699  FEEvaluation(const Mapping<dim> & mapping,
6700  const FiniteElement<dim> &fe,
6701  const Quadrature<1> & quadrature,
6702  const UpdateFlags update_flags,
6703  const unsigned int first_selected_component)
6704  : BaseClass(mapping,
6705  fe,
6706  quadrature,
6707  update_flags,
6708  first_selected_component,
6709  static_cast<
6711  nullptr))
6712  , dofs_per_component(this->data->dofs_per_component_on_cell)
6713  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6714  , n_q_points(this->data->n_q_points)
6715 {
6716  check_template_arguments(numbers::invalid_unsigned_int, 0);
6717 }
6718 
6719 
6720 
6721 template <int dim,
6722  int fe_degree,
6723  int n_q_points_1d,
6724  int n_components_,
6725  typename Number,
6726  typename VectorizedArrayType>
6727 inline FEEvaluation<dim,
6728  fe_degree,
6729  n_q_points_1d,
6730  n_components_,
6731  Number,
6732  VectorizedArrayType>::
6733  FEEvaluation(const FiniteElement<dim> &fe,
6734  const Quadrature<1> & quadrature,
6735  const UpdateFlags update_flags,
6736  const unsigned int first_selected_component)
6737  : BaseClass(StaticMappingQ1<dim>::mapping,
6738  fe,
6739  quadrature,
6740  update_flags,
6741  first_selected_component,
6742  static_cast<
6744  nullptr))
6745  , dofs_per_component(this->data->dofs_per_component_on_cell)
6746  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6747  , n_q_points(this->data->n_q_points)
6748 {
6749  check_template_arguments(numbers::invalid_unsigned_int, 0);
6750 }
6751 
6752 
6753 
6754 template <int dim,
6755  int fe_degree,
6756  int n_q_points_1d,
6757  int n_components_,
6758  typename Number,
6759  typename VectorizedArrayType>
6760 template <int n_components_other>
6761 inline FEEvaluation<dim,
6762  fe_degree,
6763  n_q_points_1d,
6764  n_components_,
6765  Number,
6766  VectorizedArrayType>::
6767  FEEvaluation(const FiniteElement<dim> & fe,
6768  const FEEvaluationBase<dim,
6769  n_components_other,
6770  Number,
6771  false,
6772  VectorizedArrayType> &other,
6773  const unsigned int first_selected_component)
6774  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6775  fe,
6776  other.mapped_geometry->get_quadrature(),
6777  other.mapped_geometry->get_fe_values().get_update_flags(),
6778  first_selected_component,
6779  &other)
6780  , dofs_per_component(this->data->dofs_per_component_on_cell)
6781  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6782  , n_q_points(this->data->n_q_points)
6783 {
6784  check_template_arguments(numbers::invalid_unsigned_int, 0);
6785 }
6786 
6787 
6788 
6789 template <int dim,
6790  int fe_degree,
6791  int n_q_points_1d,
6792  int n_components_,
6793  typename Number,
6794  typename VectorizedArrayType>
6795 inline FEEvaluation<dim,
6796  fe_degree,
6797  n_q_points_1d,
6798  n_components_,
6799  Number,
6800  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6801  &other)
6802  : BaseClass(other)
6803  , dofs_per_component(this->data->dofs_per_component_on_cell)
6804  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6805  , n_q_points(this->data->n_q_points)
6806 {
6807  check_template_arguments(numbers::invalid_unsigned_int, 0);
6808 }
6809 
6810 
6811 
6812 template <int dim,
6813  int fe_degree,
6814  int n_q_points_1d,
6815  int n_components_,
6816  typename Number,
6817  typename VectorizedArrayType>
6818 inline FEEvaluation<dim,
6819  fe_degree,
6820  n_q_points_1d,
6821  n_components_,
6822  Number,
6823  VectorizedArrayType> &
6824 FEEvaluation<dim,
6825  fe_degree,
6826  n_q_points_1d,
6827  n_components_,
6828  Number,
6829  VectorizedArrayType>::operator=(const FEEvaluation &other)
6830 {
6831  BaseClass::operator=(other);
6832  check_template_arguments(numbers::invalid_unsigned_int, 0);
6833  return *this;
6834 }
6835 
6836 
6837 
6838 template <int dim,
6839  int fe_degree,
6840  int n_q_points_1d,
6841  int n_components_,
6842  typename Number,
6843  typename VectorizedArrayType>
6844 inline void
6845 FEEvaluation<dim,
6846  fe_degree,
6847  n_q_points_1d,
6848  n_components_,
6849  Number,
6850  VectorizedArrayType>::
6851  check_template_arguments(const unsigned int dof_no,
6852  const unsigned int first_selected_component)
6853 {
6854  (void)dof_no;
6855  (void)first_selected_component;
6856 
6857 # ifdef DEBUG
6858  // print error message when the dimensions do not match. Propose a possible
6859  // fix
6860  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6861  static_cast<unsigned int>(fe_degree) !=
6862  this->data->data.front().fe_degree) ||
6863  n_q_points != this->n_quadrature_points)
6864  {
6865  std::string message =
6866  "-------------------------------------------------------\n";
6867  message += "Illegal arguments in constructor/wrong template arguments!\n";
6868  message += " Called --> FEEvaluation<dim,";
6869  message += Utilities::int_to_string(fe_degree) + ",";
6870  message += Utilities::int_to_string(n_q_points_1d);
6871  message += "," + Utilities::int_to_string(n_components);
6872  message += ",Number>(data";
6873  if (first_selected_component != numbers::invalid_unsigned_int)
6874  {
6875  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6876  message += Utilities::int_to_string(this->quad_no) + ", ";
6877  message += Utilities::int_to_string(first_selected_component);
6878  }
6879  message += ")\n";
6880 
6881  // check whether some other vector component has the correct number of
6882  // points
6883  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6884  proposed_fe_comp = numbers::invalid_unsigned_int,
6885  proposed_quad_comp = numbers::invalid_unsigned_int;
6886  if (dof_no != numbers::invalid_unsigned_int)
6887  {
6888  if (static_cast<unsigned int>(fe_degree) ==
6889  this->data->data.front().fe_degree)
6890  {
6891  proposed_dof_comp = dof_no;
6892  proposed_fe_comp = first_selected_component;
6893  }
6894  else
6895  for (unsigned int no = 0; no < this->matrix_info->n_components();
6896  ++no)
6897  for (unsigned int nf = 0;
6898  nf < this->matrix_info->n_base_elements(no);
6899  ++nf)
6900  if (this->matrix_info
6901  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6902  .data.front()
6903  .fe_degree == static_cast<unsigned int>(fe_degree))
6904  {
6905  proposed_dof_comp = no;
6906  proposed_fe_comp = nf;
6907  break;
6908  }
6909  if (n_q_points ==
6910  this->mapping_data->descriptor[this->active_quad_index]
6911  .n_q_points)
6912  proposed_quad_comp = this->quad_no;
6913  else
6914  for (unsigned int no = 0;
6915  no < this->matrix_info->get_mapping_info().cell_data.size();
6916  ++no)
6917  if (this->matrix_info->get_mapping_info()
6918  .cell_data[no]
6919  .descriptor[this->active_quad_index]
6920  .n_q_points == n_q_points)
6921  {
6922  proposed_quad_comp = no;
6923  break;
6924  }
6925  }
6926  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6927  proposed_quad_comp != numbers::invalid_unsigned_int)
6928  {
6929  if (proposed_dof_comp != first_selected_component)
6930  message += "Wrong vector component selection:\n";
6931  else
6932  message += "Wrong quadrature formula selection:\n";
6933  message += " Did you mean FEEvaluation<dim,";
6934  message += Utilities::int_to_string(fe_degree) + ",";
6935  message += Utilities::int_to_string(n_q_points_1d);
6936  message += "," + Utilities::int_to_string(n_components);
6937  message += ",Number>(data";
6938  if (dof_no != numbers::invalid_unsigned_int)
6939  {
6940  message +=
6941  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6942  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
6943  message += Utilities::int_to_string(proposed_fe_comp);
6944  }
6945  message += ")?\n";
6946  std::string correct_pos;
6947  if (proposed_dof_comp != dof_no)
6948  correct_pos = " ^ ";
6949  else
6950  correct_pos = " ";
6951  if (proposed_quad_comp != this->quad_no)
6952  correct_pos += " ^ ";
6953  else
6954  correct_pos += " ";
6955  if (proposed_fe_comp != first_selected_component)
6956  correct_pos += " ^\n";
6957  else
6958  correct_pos += " \n";
6959  message += " " +
6960  correct_pos;
6961  }
6962  // ok, did not find the numbers specified by the template arguments in
6963  // the given list. Suggest correct template arguments
6964  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6965  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6966  message += "Wrong template arguments:\n";
6967  message += " Did you mean FEEvaluation<dim,";
6968  message +=
6969  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6970  message += Utilities::int_to_string(proposed_n_q_points_1d);
6971  message += "," + Utilities::int_to_string(n_components);
6972  message += ",Number>(data";
6973  if (dof_no != numbers::invalid_unsigned_int)
6974  {
6975  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6976  message += Utilities::int_to_string(this->quad_no);
6977  message += ", " + Utilities::int_to_string(first_selected_component);
6978  }
6979  message += ")?\n";
6980  std::string correct_pos;
6981  if (this->data->data.front().fe_degree !=
6982  static_cast<unsigned int>(fe_degree))
6983  correct_pos = " ^";
6984  else
6985  correct_pos = " ";
6986  if (proposed_n_q_points_1d != n_q_points_1d)
6987  correct_pos += " ^\n";
6988  else
6989  correct_pos += " \n";
6990  message += " " + correct_pos;
6991 
6992  Assert(static_cast<unsigned int>(fe_degree) ==
6993  this->data->data.front().fe_degree &&
6994  n_q_points == this->n_quadrature_points,
6995  ExcMessage(message));
6996  }
6997  if (dof_no != numbers::invalid_unsigned_int)
6999  n_q_points,
7000  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7001 # endif
7002 }
7003 
7004 
7005 
7006 template <int dim,
7007  int fe_degree,
7008  int n_q_points_1d,
7009  int n_components_,
7010  typename Number,
7011  typename VectorizedArrayType>
7012 inline void
7013 FEEvaluation<dim,
7014  fe_degree,
7015  n_q_points_1d,
7016  n_components_,
7017  Number,
7018  VectorizedArrayType>::reinit(const unsigned int cell_index)
7019 {
7020  Assert(this->mapped_geometry == nullptr,
7021  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7022  " Integer indexing is not possible"));
7023  if (this->mapped_geometry != nullptr)
7024  return;
7025 
7026  Assert(this->dof_info != nullptr, ExcNotInitialized());
7027  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7028  this->cell = cell_index;
7029  this->cell_type =
7030  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
7031 
7032  const unsigned int offsets =
7033  this->mapping_data->data_index_offsets[cell_index];
7034  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7035  this->J_value = &this->mapping_data->JxW_values[offsets];
7036 
7037 # ifdef DEBUG
7038  this->dof_values_initialized = false;
7039  this->values_quad_initialized = false;
7040  this->gradients_quad_initialized = false;
7041  this->hessians_quad_initialized = false;
7042 # endif
7043 }
7044 
7045 
7046 
7047 template <int dim,
7048  int fe_degree,
7049  int n_q_points_1d,
7050  int n_components_,
7051  typename Number,
7052  typename VectorizedArrayType>
7053 template <bool level_dof_access>
7054 inline void
7055 FEEvaluation<dim,
7056  fe_degree,
7057  n_q_points_1d,
7058  n_components_,
7059  Number,
7060  VectorizedArrayType>::
7062 {
7063  Assert(this->matrix_info == nullptr,
7064  ExcMessage("Cannot use initialization from cell iterator if "
7065  "initialized from MatrixFree object. Use variant for "
7066  "on the fly computation with arguments as for FEValues "
7067  "instead"));
7068  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7069  this->mapped_geometry->reinit(
7070  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7071  this->local_dof_indices.resize(cell->get_fe().dofs_per_cell);
7072  if (level_dof_access)
7073  cell->get_mg_dof_indices(this->local_dof_indices);
7074  else
7075  cell->get_dof_indices(this->local_dof_indices);
7076 }
7077 
7078 
7079 
7080 template <int dim,
7081  int fe_degree,
7082  int n_q_points_1d,
7083  int n_components_,
7084  typename Number,
7085  typename VectorizedArrayType>
7086 inline void
7087 FEEvaluation<dim,
7088  fe_degree,
7089  n_q_points_1d,
7090  n_components_,
7091  Number,
7092  VectorizedArrayType>::
7093  reinit(const typename Triangulation<dim>::cell_iterator &cell)
7094 {
7095  Assert(this->matrix_info == 0,
7096  ExcMessage("Cannot use initialization from cell iterator if "
7097  "initialized from MatrixFree object. Use variant for "
7098  "on the fly computation with arguments as for FEValues "
7099  "instead"));
7100  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7101  this->mapped_geometry->reinit(cell);
7102 }
7103 
7104 
7105 
7106 template <int dim,
7107  int fe_degree,
7108  int n_q_points_1d,
7109  int n_components_,
7110  typename Number,
7111  typename VectorizedArrayType>
7113 FEEvaluation<dim,
7114  fe_degree,
7115  n_q_points_1d,
7116  n_components_,
7117  Number,
7118  VectorizedArrayType>::quadrature_point(const unsigned int q) const
7119 {
7120  if (this->matrix_info == nullptr)
7121  {
7122  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7124  ExcNotInitialized());
7125  }
7126  else
7127  {
7128  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
7129  ExcNotInitialized());
7130  }
7131 
7132  AssertIndexRange(q, n_q_points);
7133 
7135  &this->mapping_data->quadrature_points
7136  [this->mapping_data->quadrature_point_offsets[this->cell]];
7137 
7138  // Cartesian/affine mesh: only first vertex of cell is stored, we must
7139  // compute it through the Jacobian (which is stored in non-inverted and
7140  // non-transposed form as index '1' in the jacobian field)
7141  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
7142  {
7143  Assert(this->jacobian != nullptr, ExcNotInitialized());
7144  Point<dim, VectorizedArrayType> point = quadrature_points[0];
7145 
7146  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
7147  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
7148  for (unsigned int d = 0; d < dim; ++d)
7149  point[d] += jac[d][d] *
7150  static_cast<Number>(
7151  this->mapping_data->descriptor[this->active_quad_index]
7152  .quadrature.point(q)[d]);
7153  else
7154  for (unsigned int d = 0; d < dim; ++d)
7155  for (unsigned int e = 0; e < dim; ++e)
7156  point[d] += jac[d][e] * static_cast<Number>(
7157  this->mapping_data
7158  ->descriptor[this->active_quad_index]
7159  .quadrature.point(q)[e]);
7160  return point;
7161  }
7162  else
7163  return quadrature_points[q];
7164 }
7165 
7166 
7167 
7168 template <int dim,
7169  int fe_degree,
7170  int n_q_points_1d,
7171  int n_components_,
7172  typename Number,
7173  typename VectorizedArrayType>
7174 inline void
7175 FEEvaluation<dim,
7176  fe_degree,
7177  n_q_points_1d,
7178  n_components_,
7179  Number,
7180  VectorizedArrayType>::evaluate(const bool evaluate_values,
7181  const bool evaluate_gradients,
7182  const bool evaluate_hessians)
7183 {
7184 # ifdef DEBUG
7185  Assert(this->dof_values_initialized == true,
7187 # endif
7188  evaluate(this->values_dofs[0],
7189  evaluate_values,
7190  evaluate_gradients,
7191  evaluate_hessians);
7192 }
7193 
7194 
7195 template <int dim,
7196  int fe_degree,
7197  int n_q_points_1d,
7198  int n_components_,
7199  typename Number,
7200  typename VectorizedArrayType>
7201 inline void
7202 FEEvaluation<dim,
7203  fe_degree,
7204  n_q_points_1d,
7205  n_components_,
7206  Number,
7207  VectorizedArrayType>::
7208  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7209 {
7210 # ifdef DEBUG
7211  Assert(this->dof_values_initialized == true,
7213 # endif
7214  evaluate(this->values_dofs[0], evaluation_flags);
7215 }
7216 
7217 
7218 
7219 template <int dim,
7220  int fe_degree,
7221  int n_q_points_1d,
7222  int n_components_,
7223  typename Number,
7224  typename VectorizedArrayType>
7225 inline void
7226 FEEvaluation<dim,
7227  fe_degree,
7228  n_q_points_1d,
7229  n_components_,
7230  Number,
7231  VectorizedArrayType>::evaluate(const VectorizedArrayType
7232  * values_array,
7233  const bool evaluate_values,
7234  const bool evaluate_gradients,
7235  const bool evaluate_hessians)
7236 {
7238  dim,
7239  fe_degree,
7240  n_q_points_1d,
7241  n_components,
7242  VectorizedArrayType>::evaluate(*this->data,
7243  const_cast<VectorizedArrayType *>(
7244  values_array),
7245  this->values_quad[0],
7246  this->gradients_quad[0][0],
7247  this->hessians_quad[0][0],
7248  this->scratch_data,
7249  evaluate_values,
7250  evaluate_gradients,
7251  evaluate_hessians);
7252 
7253 # ifdef DEBUG
7254  if (evaluate_values == true)
7255  this->values_quad_initialized = true;
7256  if (evaluate_gradients == true)
7257  this->gradients_quad_initialized = true;
7258  if (evaluate_hessians == true)
7259  this->hessians_quad_initialized = true;
7260 # endif
7261 }
7262 
7263 
7264 
7265 template <int dim,
7266  int fe_degree,
7267  int n_q_points_1d,
7268  int n_components_,
7269  typename Number,
7270  typename VectorizedArrayType>
7271 inline void
7272 FEEvaluation<dim,
7273  fe_degree,
7274  n_q_points_1d,
7275  n_components_,
7276  Number,
7277  VectorizedArrayType>::
7278  evaluate(const VectorizedArrayType * values_array,
7279  const EvaluationFlags::EvaluationFlags evaluation_flags)
7280 {
7281  SelectEvaluator<dim,
7282  fe_degree,
7283  n_q_points_1d,
7284  n_components,
7285  VectorizedArrayType>::
7286  evaluate(*this->data,
7287  const_cast<VectorizedArrayType *>(values_array),
7288  this->values_quad[0],
7289  this->gradients_quad[0][0],
7290  this->hessians_quad[0][0],
7291  this->scratch_data,
7292  evaluation_flags & EvaluationFlags::values,
7293  evaluation_flags & EvaluationFlags::gradients,
7294  evaluation_flags & EvaluationFlags::hessians);
7295 
7296 # ifdef DEBUG
7297  if (evaluation_flags & EvaluationFlags::values)
7298  this->values_quad_initialized = true;
7299  if (evaluation_flags & EvaluationFlags::gradients)
7300  this->gradients_quad_initialized = true;
7301  if (evaluation_flags & EvaluationFlags::hessians)
7302  this->hessians_quad_initialized = true;
7303 # endif
7304 }
7305 
7306 
7307 
7308 template <int dim,
7309  int fe_degree,
7310  int n_q_points_1d,
7311  int n_components_,
7312  typename Number,
7313  typename VectorizedArrayType>
7314 template <typename VectorType>
7315 inline void
7316 FEEvaluation<
7317  dim,
7318  fe_degree,
7319  n_q_points_1d,
7320  n_components_,
7321  Number,
7322  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7323 
7324  const bool evaluate_values,
7325  const bool evaluate_gradients,
7326  const bool evaluate_hessians)
7327 {
7329  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7330  ((evaluate_gradients) ? EvaluationFlags::gradients :
7331  EvaluationFlags::nothing) |
7332  ((evaluate_hessians) ? EvaluationFlags::hessians :
7333  EvaluationFlags::nothing);
7334 
7335  gather_evaluate(input_vector, flag);
7336 }
7337 
7338 template <int dim,
7339  int fe_degree,
7340  int n_q_points_1d,
7341  int n_components_,
7342  typename Number,
7343  typename VectorizedArrayType>
7344 template <typename VectorType>
7345 inline void
7346 FEEvaluation<dim,
7347  fe_degree,
7348  n_q_points_1d,
7349  n_components_,
7350  Number,
7351  VectorizedArrayType>::
7352  gather_evaluate(const VectorType & input_vector,
7353  const EvaluationFlags::EvaluationFlags evaluation_flag)
7354 {
7355  // If the index storage is interleaved and contiguous and the vector storage
7356  // has the correct alignment, we can directly pass the pointer into the
7357  // vector to the evaluate() call, without reading the vector entries into a
7358  // separate data field. This saves some operations.
7359  if (std::is_same<typename VectorType::value_type, Number>::value &&
7360  this->dof_info->index_storage_variants
7362  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7363  IndexStorageVariants::interleaved_contiguous &&
7364  reinterpret_cast<std::size_t>(
7365  input_vector.begin() +
7366  this->dof_info->dof_indices_contiguous
7368  [this->cell * VectorizedArrayType::size()]) %
7369  sizeof(VectorizedArrayType) ==
7370  0)
7371  {
7372  const VectorizedArrayType *vec_values =
7373  reinterpret_cast<const VectorizedArrayType *>(
7374  input_vector.begin() +
7375  this->dof_info->dof_indices_contiguous
7377  [this->cell * VectorizedArrayType::size()] +
7378  this->dof_info
7379  ->component_dof_indices_offset[this->active_fe_index]
7380  [this->first_selected_component] *
7381  VectorizedArrayType::size());
7382 
7383  evaluate(vec_values,
7384  evaluation_flag & EvaluationFlags::values,
7385  evaluation_flag & EvaluationFlags::gradients,
7386  evaluation_flag & EvaluationFlags::hessians);
7387  }
7388  else
7389  {
7390  this->read_dof_values(input_vector);
7391  evaluate(this->begin_dof_values(),
7392  evaluation_flag & EvaluationFlags::values,
7393  evaluation_flag & EvaluationFlags::gradients,
7394  evaluation_flag & EvaluationFlags::hessians);
7395  }
7396 }
7397 
7398 
7399 
7400 template <int dim,
7401  int fe_degree,
7402  int n_q_points_1d,
7403  int n_components_,
7404  typename Number,
7405  typename VectorizedArrayType>
7406 inline void
7407 FEEvaluation<dim,
7408  fe_degree,
7409  n_q_points_1d,
7410  n_components_,
7411  Number,
7412  VectorizedArrayType>::integrate(const bool integrate_values,
7413  const bool integrate_gradients)
7414 {
7415  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7416 
7417 # ifdef DEBUG
7418  this->dof_values_initialized = true;
7419 # endif
7420 }
7421 
7422 
7423 
7424 template <int dim,
7425  int fe_degree,
7426  int n_q_points_1d,
7427  int n_components_,
7428  typename Number,
7429  typename VectorizedArrayType>
7430 inline void
7431 FEEvaluation<dim,
7432  fe_degree,
7433  n_q_points_1d,
7434  n_components_,
7435  Number,
7436  VectorizedArrayType>::
7437  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7438 {
7439  integrate(integration_flag, this->values_dofs[0]);
7440 
7441 # ifdef DEBUG
7442  this->dof_values_initialized = true;
7443 # endif
7444 }
7445 
7446 
7447 
7448 template <int dim,
7449  int fe_degree,
7450  int n_q_points_1d,
7451  int n_components_,
7452  typename Number,
7453  typename VectorizedArrayType>
7454 inline void
7455 FEEvaluation<dim,
7456  fe_degree,
7457  n_q_points_1d,
7458  n_components_,
7459  Number,
7460  VectorizedArrayType>::integrate(const bool integrate_values,
7461  const bool integrate_gradients,
7462  VectorizedArrayType *values_array)
7463 {
7464 # ifdef DEBUG
7465  if (integrate_values == true)
7466  Assert(this->values_quad_submitted == true,
7468  if (integrate_gradients == true)
7469  Assert(this->gradients_quad_submitted == true,
7471 # endif
7472  Assert(this->matrix_info != nullptr ||
7473  this->mapped_geometry->is_initialized(),
7474  ExcNotInitialized());
7475 
7476  SelectEvaluator<dim,
7477  fe_degree,
7478  n_q_points_1d,
7479  n_components,
7480  VectorizedArrayType>::integrate(*this->data,
7481  values_array,
7482  this->values_quad[0],
7483  this->gradients_quad[0][0],
7484  this->scratch_data,
7485  integrate_values,
7486  integrate_gradients,
7487  false);
7488 
7489 # ifdef DEBUG
7490  this->dof_values_initialized = true;
7491 # endif
7492 }
7493 
7494 
7495 
7496 template <int dim,
7497  int fe_degree,
7498  int n_q_points_1d,
7499  int n_components_,
7500  typename Number,
7501  typename VectorizedArrayType>
7502 inline void
7503 FEEvaluation<dim,
7504  fe_degree,
7505  n_q_points_1d,
7506  n_components_,
7507  Number,
7508  VectorizedArrayType>::
7509  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7510  VectorizedArrayType * values_array)
7511 {
7512 # ifdef DEBUG
7513  if (integration_flag & EvaluationFlags::values)
7514  Assert(this->values_quad_submitted == true,
7516  if (integration_flag & EvaluationFlags::gradients)
7517  Assert(this->gradients_quad_submitted == true,
7519 # endif
7520  Assert(this->matrix_info != nullptr ||
7521  this->mapped_geometry->is_initialized(),
7522  ExcNotInitialized());
7523 
7524  Assert(
7525  (integration_flag &
7526  ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
7527  ExcMessage(
7528  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
7529 
7530  SelectEvaluator<dim,
7531  fe_degree,
7532  n_q_points_1d,
7533  n_components,
7534  VectorizedArrayType>::integrate(*this->data,
7535  values_array,
7536  this->values_quad[0],
7537  this->gradients_quad[0][0],
7538  this->scratch_data,
7539  integration_flag &
7540  EvaluationFlags::values,
7541  integration_flag &
7542  EvaluationFlags::gradients,
7543  false);
7544 
7545 # ifdef DEBUG
7546  this->dof_values_initialized = true;
7547 # endif
7548 }
7549 
7550 
7551 
7552 template <int dim,
7553  int fe_degree,
7554  int n_q_points_1d,
7555  int n_components_,
7556  typename Number,
7557  typename VectorizedArrayType>
7558 template <typename VectorType>
7559 inline void
7560 FEEvaluation<
7561  dim,
7562  fe_degree,
7563  n_q_points_1d,
7564  n_components_,
7565  Number,
7566  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7567  const bool integrate_gradients,
7568  VectorType &destination)
7569 {
7571  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7572  ((integrate_gradients) ? EvaluationFlags::gradients :
7573  EvaluationFlags::nothing);
7574 
7575  integrate_scatter(flag, destination);
7576 }
7577 
7578 
7579 
7580 template <int dim,
7581  int fe_degree,
7582  int n_q_points_1d,
7583  int n_components_,
7584  typename Number,
7585  typename VectorizedArrayType>
7586 template <typename VectorType>
7587 inline void
7588 FEEvaluation<dim,
7589  fe_degree,
7590  n_q_points_1d,
7591  n_components_,
7592  Number,
7593  VectorizedArrayType>::
7594  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
7595  VectorType & destination)
7596 {
7597  // If the index storage is interleaved and contiguous and the vector storage
7598  // has the correct alignment, we can directly pass the pointer into the
7599  // vector to the integrate() call, without writing temporary results into a
7600  // separate data field that will later be added into the vector. This saves
7601  // some operations.
7602  if (std::is_same<typename VectorType::value_type, Number>::value &&
7603  this->dof_info->index_storage_variants
7605  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7606  IndexStorageVariants::interleaved_contiguous &&
7607  reinterpret_cast<std::size_t>(
7608  destination.begin() +
7609  this->dof_info->dof_indices_contiguous
7611  [this->cell * VectorizedArrayType::size()]) %
7612  sizeof(VectorizedArrayType) ==
7613  0)
7614  {
7615  VectorizedArrayType *vec_values = reinterpret_cast<VectorizedArrayType *>(
7616  destination.begin() +
7617  this->dof_info->dof_indices_contiguous
7619  [this->cell * VectorizedArrayType::size()] +
7620  this->dof_info
7621  ->component_dof_indices_offset[this->active_fe_index]
7622  [this->first_selected_component] *
7623  VectorizedArrayType::size());
7625  dim,
7626  fe_degree,
7627  n_q_points_1d,
7628  n_components,
7629  VectorizedArrayType>::integrate(*this->data,
7630  vec_values,
7631  this->values_quad[0],
7632  this->gradients_quad[0][0],
7633  this->scratch_data,
7634  evaluation_flag &
7636  evaluation_flag &
7638  true);
7639  }
7640  else
7641  {
7642  integrate(evaluation_flag & EvaluationFlags::values,
7643  evaluation_flag & EvaluationFlags::gradients,
7644  this->begin_dof_values());
7645  this->distribute_local_to_global(destination);
7646  }
7647 }
7648 
7649 
7650 
7651 /*-------------------------- FEFaceEvaluation ---------------------------*/
7652 
7653 
7654 
7655 template <int dim,
7656  int fe_degree,
7657  int n_q_points_1d,
7658  int n_components_,
7659  typename Number,
7660  typename VectorizedArrayType>
7661 inline FEFaceEvaluation<dim,
7662  fe_degree,
7663  n_q_points_1d,
7664  n_components_,
7665  Number,
7666  VectorizedArrayType>::
7667  FEFaceEvaluation(
7669  const bool is_interior_face,
7670  const unsigned int dof_no,
7671  const unsigned int quad_no,
7672  const unsigned int first_selected_component)
7673  : BaseClass(matrix_free,
7674  dof_no,
7675  first_selected_component,
7676  quad_no,
7677  fe_degree,
7678  static_n_q_points,
7679  is_interior_face)
7680  , dofs_per_component(this->data->dofs_per_component_on_cell)
7681  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7682  , n_q_points(this->data->n_q_points_face)
7683 {}
7684 
7685 
7686 
7687 template <int dim,
7688  int fe_degree,
7689  int n_q_points_1d,
7690  int n_components_,
7691  typename Number,
7692  typename VectorizedArrayType>
7693 inline void
7694 FEFaceEvaluation<dim,
7695  fe_degree,
7696  n_q_points_1d,
7697  n_components_,
7698  Number,
7699  VectorizedArrayType>::reinit(const unsigned int face_index)
7700 {
7701  Assert(this->mapped_geometry == nullptr,
7702  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7703  " Integer indexing is not possible"));
7704  if (this->mapped_geometry != nullptr)
7705  return;
7706 
7707  this->cell = face_index;
7708  this->dof_access_index =
7709  this->is_interior_face ?
7712  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7714  VectorizedArrayType::size()> &faces =
7715  this->matrix_info->get_face_info(face_index);
7716  if (face_index >=
7717  this->matrix_info->get_task_info().face_partition_data.back() &&
7718  face_index <
7719  this->matrix_info->get_task_info().boundary_partition_data.back())
7720  Assert(this->is_interior_face,
7721  ExcMessage("Boundary faces do not have a neighbor"));
7722 
7723  this->face_no =
7724  (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
7725  this->subface_index = faces.subface_index;
7726  if (this->is_interior_face == true)
7727  {
7728  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7729  if (faces.face_orientation > 8)
7730  this->face_orientation = faces.face_orientation - 8;
7731  else
7732  this->face_orientation = 0;
7733  }
7734  else
7735  {
7736  if (faces.face_orientation < 8)
7737  this->face_orientation = faces.face_orientation;
7738  else
7739  this->face_orientation = 0;
7740  }
7741 
7742  this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
7743  const unsigned int offsets =
7744  this->mapping_data->data_index_offsets[face_index];
7745  this->J_value = &this->mapping_data->JxW_values[offsets];
7746  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7747  this->jacobian =
7748  &this->mapping_data->jacobians[!this->is_interior_face][offsets];
7749  this->normal_x_jacobian =
7750  &this->mapping_data
7751  ->normals_times_jacobians[!this->is_interior_face][offsets];
7752 
7753 # ifdef DEBUG
7754  this->dof_values_initialized = false;
7755  this->values_quad_initialized = false;
7756  this->gradients_quad_initialized = false;
7757  this->hessians_quad_initialized = false;
7758 # endif
7759 }
7760 
7761 
7762 
7763 template <int dim,
7764  int fe_degree,
7765  int n_q_points_1d,
7766  int n_components_,
7767  typename Number,
7768  typename VectorizedArrayType>
7769 inline void
7770 FEFaceEvaluation<dim,
7771  fe_degree,
7772  n_q_points_1d,
7773  n_components_,
7774  Number,
7775  VectorizedArrayType>::reinit(const unsigned int cell_index,
7776  const unsigned int face_number)
7777 {
7778  Assert(
7779  this->quad_no <
7780  this->matrix_info->get_mapping_info().face_data_by_cells.size(),
7781  ExcMessage(
7782  "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7784  AssertIndexRange(cell_index,
7785  this->matrix_info->get_mapping_info().cell_type.size());
7786  Assert(this->mapped_geometry == nullptr,
7787  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7788  " Integer indexing is not possible"));
7789  if (this->mapped_geometry != nullptr)
7790  return;
7791  Assert(this->matrix_info != nullptr, ExcNotInitialized());
7792 
7793  this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
7794  this->cell = cell_index;
7795  this->face_orientation = 0;
7796  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7797  this->face_no = face_number;
7798  this->dof_access_index =
7800 
7801  const unsigned int offsets =
7802  this->matrix_info->get_mapping_info()
7803  .face_data_by_cells[this->quad_no]
7804  .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
7805  face_number];
7806  AssertIndexRange(offsets,
7807  this->matrix_info->get_mapping_info()
7808  .face_data_by_cells[this->quad_no]
7809  .JxW_values.size());
7810  this->J_value = &this->matrix_info->get_mapping_info()
7811  .face_data_by_cells[this->quad_no]
7812  .JxW_values[offsets];
7813  this->normal_vectors = &this->matrix_info->get_mapping_info()
7814  .face_data_by_cells[this->quad_no]
7815  .normal_vectors[offsets];
7816  this->jacobian = &this->matrix_info->get_mapping_info()
7817  .face_data_by_cells[this->quad_no]
7818  .jacobians[!this->is_interior_face][offsets];
7819  this->normal_x_jacobian =
7820  &this->matrix_info->get_mapping_info()
7821  .face_data_by_cells[this->quad_no]
7822  .normals_times_jacobians[!this->is_interior_face][offsets];
7823 
7824 # ifdef DEBUG
7825  this->dof_values_initialized = false;
7826  this->values_quad_initialized = false;
7827  this->gradients_quad_initialized = false;
7828  this->hessians_quad_initialized = false;
7829 # endif
7830 }
7831 
7832 
7833 
7834 template <int dim,
7835  int fe_degree,
7836  int n_q_points_1d,
7837  int n_components,
7838  typename Number,
7839  typename VectorizedArrayType>
7840 inline void
7841 FEFaceEvaluation<dim,
7842  fe_degree,
7843  n_q_points_1d,
7844  n_components,
7845  Number,
7846  VectorizedArrayType>::evaluate(const bool evaluate_values,
7847  const bool evaluate_gradients)
7848 {
7849 # ifdef DEBUG
7850  Assert(this->dof_values_initialized, ExcNotInitialized());
7851 # endif
7852 
7853  evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients);
7854 }
7855 
7856 
7857 
7858 template <int dim,
7859  int fe_degree,
7860  int n_q_points_1d,
7861  int n_components,
7862  typename Number,
7863  typename VectorizedArrayType>
7864 inline void
7865 FEFaceEvaluation<dim,
7866  fe_degree,
7867  n_q_points_1d,
7868  n_components,
7869  Number,
7870  VectorizedArrayType>::
7871  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
7872 {
7873 # ifdef DEBUG
7874  Assert(this->dof_values_initialized, ExcNotInitialized());
7875 # endif
7876 
7877  evaluate(this->values_dofs[0], evaluation_flag);
7878 }
7879 
7880 
7881 
7882 template <int dim,
7883  int fe_degree,
7884  int n_q_points_1d,
7885  int n_components,
7886  typename Number,
7887  typename VectorizedArrayType>
7888 inline void
7889 FEFaceEvaluation<dim,
7890  fe_degree,
7891  n_q_points_1d,
7892  n_components,
7893  Number,
7894  VectorizedArrayType>::evaluate(const VectorizedArrayType
7895  * values_array,
7896  const bool evaluate_values,
7897  const bool evaluate_gradients)
7898 {
7900  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7901  ((evaluate_gradients) ? EvaluationFlags::gradients :
7902  EvaluationFlags::nothing);
7903 
7904  evaluate(values_array, flag);
7905 }
7906 
7907 
7908 
7909 template <int dim,
7910  int fe_degree,
7911  int n_q_points_1d,
7912  int n_components,
7913  typename Number,
7914  typename VectorizedArrayType>
7915 inline void
7916 FEFaceEvaluation<dim,
7917  fe_degree,
7918  n_q_points_1d,
7919  n_components,
7920  Number,
7921  VectorizedArrayType>::
7922  evaluate(const VectorizedArrayType * values_array,
7923  const EvaluationFlags::EvaluationFlags evaluation_flag)
7924 {
7925  Assert(
7926  (evaluation_flag &
7928  ExcMessage(
7929  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
7930 
7931  if (!(evaluation_flag & EvaluationFlags::values) &&
7932  !(evaluation_flag & EvaluationFlags::gradients))
7933  return;
7934 
7936  dim,
7937  fe_degree,
7938  n_q_points_1d,
7939  n_components,
7940  Number,
7941  VectorizedArrayType>::evaluate(*this->data,
7942  values_array,
7943  this->begin_values(),
7944  this->begin_gradients(),
7945  this->scratch_data,
7946  evaluation_flag & EvaluationFlags::values,
7947  evaluation_flag & EvaluationFlags::gradients,
7948  this->face_no,
7949  this->subface_index,
7950  this->face_orientation,
7951  this->mapping_data
7952  ->descriptor[this->active_fe_index]
7953  .face_orientations);
7954 
7955 # ifdef DEBUG
7956  if (evaluation_flag & EvaluationFlags::values)
7957  this->values_quad_initialized = true;
7958  if (evaluation_flag & EvaluationFlags::gradients)
7959  this->gradients_quad_initialized = true;
7960 # endif
7961 }
7962 
7963 
7964 
7965 template <int dim,
7966  int fe_degree,
7967  int n_q_points_1d,
7968  int n_components,
7969  typename Number,
7970  typename VectorizedArrayType>
7971 inline void
7972 FEFaceEvaluation<dim,
7973  fe_degree,
7974  n_q_points_1d,
7975  n_components,
7976  Number,
7977  VectorizedArrayType>::
7978  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag)
7979 {
7980  integrate(evaluation_flag, this->values_dofs[0]);
7981 
7982 # ifdef DEBUG
7983  this->dof_values_initialized = true;
7984 # endif
7985 }
7986 
7987 
7988 
7989 template <int dim,
7990  int fe_degree,
7991  int n_q_points_1d,
7992  int n_components,
7993  typename Number,
7994  typename VectorizedArrayType>
7995 inline void
7996 FEFaceEvaluation<dim,
7997  fe_degree,
7998  n_q_points_1d,
7999  n_components,
8000  Number,
8001  VectorizedArrayType>::integrate(const bool integrate_values,
8002  const bool integrate_gradients)
8003 {
8004  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8005 
8006 # ifdef DEBUG
8007  this->dof_values_initialized = true;
8008 # endif
8009 }
8010 
8011 
8012 
8013 template <int dim,
8014  int fe_degree,
8015  int n_q_points_1d,
8016  int n_components,
8017  typename Number,
8018  typename VectorizedArrayType>
8019 inline void
8020 FEFaceEvaluation<dim,
8021  fe_degree,
8022  n_q_points_1d,
8023  n_components,
8024  Number,
8025  VectorizedArrayType>::integrate(const bool integrate_values,
8026  const bool integrate_gradients,
8027  VectorizedArrayType
8028  *values_array)
8029 {
8031  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8032  ((integrate_gradients) ? EvaluationFlags::gradients :
8033  EvaluationFlags::nothing);
8034 
8035  integrate(flag, values_array);
8036 }
8037 
8038 
8039 
8040 template <int dim,
8041  int fe_degree,
8042  int n_q_points_1d,
8043  int n_components,
8044  typename Number,
8045  typename VectorizedArrayType>
8046 inline void
8047 FEFaceEvaluation<dim,
8048  fe_degree,
8049  n_q_points_1d,
8050  n_components,
8051  Number,
8052  VectorizedArrayType>::
8053  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
8054  VectorizedArrayType * values_array)
8055 {
8056  Assert(
8057  (evaluation_flag &
8059  ExcMessage(
8060  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8061 
8062  if (!(evaluation_flag & EvaluationFlags::values) &&
8063  !(evaluation_flag & EvaluationFlags::gradients))
8064  return;
8065 
8067  fe_degree,
8068  n_q_points_1d,
8069  n_components,
8070  Number,
8071  VectorizedArrayType>::
8072  integrate(
8073  *this->data,
8074  values_array,
8075  this->begin_values(),
8076  this->begin_gradients(),
8077  this->scratch_data,
8078  evaluation_flag & EvaluationFlags::values,
8079  evaluation_flag & EvaluationFlags::gradients,
8080  this->face_no,
8081  this->subface_index,
8082  this->face_orientation,
8083  this->mapping_data->descriptor[this->active_fe_index].face_orientations);
8084 }
8085 
8086 
8087 
8088 template <int dim,
8089  int fe_degree,
8090  int n_q_points_1d,
8091  int n_components_,
8092  typename Number,
8093  typename VectorizedArrayType>
8094 template <typename VectorType>
8095 inline void
8097  dim,
8098  fe_degree,
8099  n_q_points_1d,
8100  n_components_,
8101  Number,
8102  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
8103  const bool evaluate_values,
8104  const bool evaluate_gradients)
8105 {
8107  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8108  ((evaluate_gradients) ? EvaluationFlags::gradients :
8109  EvaluationFlags::nothing);
8110 
8111  gather_evaluate(input_vector, flag);
8112 }
8113 
8114 
8115 
8116 template <int dim,
8117  int fe_degree,
8118  int n_q_points_1d,
8119  int n_components_,
8120  typename Number,
8121  typename VectorizedArrayType>
8122 template <typename VectorType>
8123 inline void
8124 FEFaceEvaluation<dim,
8125  fe_degree,
8126  n_q_points_1d,
8127  n_components_,
8128  Number,
8129  VectorizedArrayType>::
8130  gather_evaluate(const VectorType & input_vector,
8131  const EvaluationFlags::EvaluationFlags evaluation_flag)
8132 {
8133  static_assert(internal::has_begin<VectorType>::value &&
8134  (std::is_same<decltype(std::declval<VectorType>().begin()),
8135  double *>::value ||
8136  std::is_same<decltype(std::declval<VectorType>().begin()),
8137  float *>::value),
8138  "This function requires a vector type with begin() function "
8139  "evaluating to a pointer to basic number (float,double). "
8140  "Use read_dof_values() followed by evaluate() instead.");
8141 
8142  Assert(
8143  (evaluation_flag &
8145  ExcMessage(
8146  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8147 
8149  fe_degree,
8150  n_q_points_1d,
8151  n_components,
8152  Number,
8153  VectorizedArrayType>::
8154  gather_evaluate(input_vector.begin(),
8155  *this->data,
8156  *this->dof_info,
8157  this->begin_values(),
8158  this->begin_gradients(),
8159  this->scratch_data,
8160  evaluation_flag & EvaluationFlags::values,
8161  evaluation_flag & EvaluationFlags::gradients,
8162  this->active_fe_index,
8163  this->first_selected_component,
8164  this->cell,
8165  this->face_no,
8166  this->subface_index,
8167  this->dof_access_index,
8168  this->face_orientation,
8169  this->mapping_data->descriptor[this->active_fe_index]
8170  .face_orientations))
8171  {
8172  this->read_dof_values(input_vector);
8173  this->evaluate(evaluation_flag);
8174  }
8175 
8176 # ifdef DEBUG
8177  if (evaluation_flag & EvaluationFlags::values)
8178  this->values_quad_initialized = true;
8179  if (evaluation_flag & EvaluationFlags::gradients)
8180  this->gradients_quad_initialized = true;
8181 # endif
8182 }
8183 
8184 
8185 
8186 template <int dim,
8187  int fe_degree,
8188  int n_q_points_1d,
8189  int n_components_,
8190  typename Number,
8191  typename VectorizedArrayType>
8192 template <typename VectorType>
8193 inline void
8195  dim,
8196  fe_degree,
8197  n_q_points_1d,
8198  n_components_,
8199  Number,
8200  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
8201  const bool integrate_gradients,
8202  VectorType &destination)
8203 {
8205  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8206  ((integrate_gradients) ? EvaluationFlags::gradients :
8207  EvaluationFlags::nothing);
8208 
8209  integrate_scatter(flag, destination);
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 template <typename VectorType>
8221 inline void
8222 FEFaceEvaluation<dim,
8223  fe_degree,
8224  n_q_points_1d,
8225  n_components_,
8226  Number,
8227  VectorizedArrayType>::
8228  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
8229  VectorType & destination)
8230 {
8231  static_assert(internal::has_begin<VectorType>::value &&
8232  (std::is_same<decltype(std::declval<VectorType>().begin()),
8233  double *>::value ||
8234  std::is_same<decltype(std::declval<VectorType>().begin()),
8235  float *>::value),
8236  "This function requires a vector type with begin() function "
8237  "evaluating to a pointer to basic number (float,double). "
8238  "Use integrate() followed by distribute_local_to_global() "
8239  "instead.");
8240 
8242  fe_degree,
8243  n_q_points_1d,
8244  n_components,
8245  Number,
8246  VectorizedArrayType>::
8247  integrate_scatter(destination.begin(),
8248  *this->data,
8249  *this->dof_info,
8250  this->begin_dof_values(),
8251  this->begin_values(),
8252  this->begin_gradients(),
8253  this->scratch_data,
8254  evaluation_flag & EvaluationFlags::values,
8255  evaluation_flag & EvaluationFlags::gradients,
8256  this->active_fe_index,
8257  this->first_selected_component,
8258  this->cell,
8259  this->face_no,
8260  this->subface_index,
8261  this->dof_access_index,
8262  this->face_orientation,
8263  this->mapping_data->descriptor[this->active_fe_index]
8264  .face_orientations))
8265  {
8266  // if we arrive here, writing into the destination vector did not succeed
8267  // because some of the assumptions in integrate_scatter were not
8268  // fulfilled (e.g. an element or degree that does not support direct
8269  // writing), so we must do it here
8270  this->distribute_local_to_global(destination);
8271  }
8272 }
8273 
8274 
8275 
8276 template <int dim,
8277  int fe_degree,
8278  int n_q_points_1d,
8279  int n_components_,
8280  typename Number,
8281  typename VectorizedArrayType>
8283 FEFaceEvaluation<dim,
8284  fe_degree,
8285  n_q_points_1d,
8286  n_components_,
8287  Number,
8288  VectorizedArrayType>::quadrature_point(const unsigned int q)
8289  const
8290 {
8291  AssertIndexRange(q, n_q_points);
8292  if (this->dof_access_index < 2)
8293  {
8294  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
8295  ExcNotImplemented());
8296  AssertIndexRange(this->cell,
8297  this->mapping_data->quadrature_point_offsets.size());
8298  return this->mapping_data->quadrature_points
8299  [this->mapping_data->quadrature_point_offsets[this->cell] + q];
8300  }
8301  else
8302  {
8303  Assert(this->matrix_info->get_mapping_info()
8304  .face_data_by_cells[this->quad_no]
8305  .quadrature_point_offsets.empty() == false,
8306  ExcNotImplemented());
8307  const unsigned int index =
8308  this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
8309  AssertIndexRange(index,
8310  this->matrix_info->get_mapping_info()
8311  .face_data_by_cells[this->quad_no]
8312  .quadrature_point_offsets.size());
8313  return this->matrix_info->get_mapping_info()
8314  .face_data_by_cells[this->quad_no]
8315  .quadrature_points[this->matrix_info->get_mapping_info()
8316  .face_data_by_cells[this->quad_no]
8317  .quadrature_point_offsets[index] +
8318  q];
8319  }
8320 }
8321 
8322 
8323 
8324 /*------------------------- end FEFaceEvaluation ------------------------- */
8325 
8326 
8327 #endif // ifndef DOXYGEN
8328 
8329 
8331 
8332 #endif
VectorizedArrayType JxW(const unsigned int q_index) const
The namespace for the EvaluationFlags enum.
Definition: fe_evaluation.h:66
const unsigned int dofs_per_component
static const unsigned int invalid_unsigned_int
Definition: types.h:191
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618