Reference documentation for deal.II version Git 040c6ad7d4 2020-09-26 18:01:03 +0200
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_fe_evaluation_h
18 #define dealii_matrix_free_fe_evaluation_h
19 
20 
21 #include <deal.II/base/config.h>
22 
29 
31 
42 
43 
45 
46 
47 
48 namespace internal
49 {
51 }
52 
53 template <int dim,
54  int fe_degree,
55  int n_q_points_1d = fe_degree + 1,
56  int n_components_ = 1,
57  typename Number = double,
58  typename VectorizedArrayType = VectorizedArray<Number>>
60 
61 
62 
89 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
91 {
92  static_assert(
93  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
94  "Type of Number and of VectorizedArrayType do not match.");
95 
96 public:
97  static constexpr unsigned int dimension = dim;
98 
103 
112  unsigned int
113  get_mapping_data_index_offset() const;
114 
122  get_cell_type() const;
123 
128  get_shape_info() const;
129 
134  VectorizedArrayType
135  JxW(const unsigned int q_point) const;
136 
144  inverse_jacobian(const unsigned int q_point) const;
145 
159  get_normal_vector(const unsigned int q_point) const;
160 
167  VectorizedArrayType
168  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
169 
176  void
177  set_cell_data(AlignedVector<VectorizedArrayType> &array,
178  const VectorizedArrayType & value) const;
179 
184  template <typename T>
185  std::array<T, VectorizedArrayType::size()>
186  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
187  &array) const;
188 
193  template <typename T>
194  void
195  set_cell_data(
196  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
197  const std::array<T, VectorizedArrayType::size()> & value) const;
198 
203  std::array<unsigned int, VectorizedArrayType::size()>
204  get_cell_ids() const;
205 
210  std::array<unsigned int, VectorizedArrayType::size()>
211  get_cell_or_face_ids() const;
212 
213 
219  const std::vector<unsigned int> &
220  get_internal_dof_numbering() const;
221 
229  get_scratch_data() const;
230 
234  unsigned int
235  get_quadrature_index() const;
236 
240  unsigned int
241  get_current_cell_index() const;
242 
243 protected:
252  const unsigned int dof_no,
253  const unsigned int first_selected_component,
254  const unsigned int quad_no,
255  const unsigned int fe_degree,
256  const unsigned int n_q_points,
257  const bool is_interior_face);
258 
264  const Mapping<dim> & mapping,
265  const FiniteElement<dim> &fe,
266  const Quadrature<1> & quadrature,
267  const UpdateFlags update_flags,
268  const unsigned int first_selected_component,
270  *other);
271 
279 
287  operator=(const FEEvaluationBaseData &other);
288 
293 
299  VectorizedArrayType *scratch_data;
300 
304  const unsigned int quad_no;
305 
309  const unsigned int active_fe_index;
310 
315  const unsigned int active_quad_index;
316 
320  const unsigned int n_quadrature_points;
321 
326 
333 
341  (is_face ? dim - 1 : dim),
342  dim,
343  Number,
344  VectorizedArrayType> *mapping_data;
345 
353 
359 
366  const VectorizedArrayType *J_value;
367 
372 
377 
381  const Number *quadrature_weights;
382 
387  unsigned int cell;
388 
394 
400 
405  unsigned int face_no;
406 
411  unsigned int face_orientation;
412 
420  unsigned int subface_index;
421 
429 
434  std::shared_ptr<internal::MatrixFreeFunctions::
435  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
437 
438  // Make FEEvaluation objects friends for access to protected member
439  // mapped_geometry.
440  template <int, int, int, int, typename, typename>
441  friend class FEEvaluation;
442 };
443 
444 
445 
483 template <int dim,
484  int n_components_,
485  typename Number,
486  bool is_face = false,
487  typename VectorizedArrayType = VectorizedArray<Number>>
489  : public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
490 {
491 public:
492  using number_type = Number;
494  using gradient_type =
496  static constexpr unsigned int dimension = dim;
497  static constexpr unsigned int n_components = n_components_;
498 
535  template <typename VectorType>
536  void
537  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
538 
567  template <typename VectorType>
568  void
569  read_dof_values_plain(const VectorType & src,
570  const unsigned int first_index = 0);
571 
603  template <typename VectorType>
604  void
605  distribute_local_to_global(
606  VectorType & dst,
607  const unsigned int first_index = 0,
608  const std::bitset<VectorizedArrayType::size()> &mask =
609  std::bitset<VectorizedArrayType::size()>().flip()) const;
610 
649  template <typename VectorType>
650  void
651  set_dof_values(VectorType & dst,
652  const unsigned int first_index = 0,
653  const std::bitset<VectorizedArrayType::size()> &mask =
654  std::bitset<VectorizedArrayType::size()>().flip()) const;
655 
659  template <typename VectorType>
660  void
661  set_dof_values_plain(
662  VectorType & dst,
663  const unsigned int first_index = 0,
664  const std::bitset<VectorizedArrayType::size()> &mask =
665  std::bitset<VectorizedArrayType::size()>().flip()) const;
666 
668 
689  value_type
690  get_dof_value(const unsigned int dof) const;
691 
702  void
703  submit_dof_value(const value_type val_in, const unsigned int dof);
704 
717  value_type
718  get_value(const unsigned int q_point) const;
719 
732  void
733  submit_value(const value_type val_in, const unsigned int q_point);
734 
746  get_gradient(const unsigned int q_point) const;
747 
762  value_type
763  get_normal_derivative(const unsigned int q_point) const;
764 
777  void
778  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
779 
798  void
799  submit_normal_derivative(const value_type grad_in,
800  const unsigned int q_point);
801 
814  get_hessian(const unsigned int q_point) const;
815 
826  get_hessian_diagonal(const unsigned int q_point) const;
827 
839  value_type
840  get_laplacian(const unsigned int q_point) const;
841 
842 #ifdef DOXYGEN
843  // doxygen does not anyhow mention functions coming from partial template
844  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
845  // For now, hack in those functions manually only to fix documentation:
846 
853  VectorizedArrayType
854  get_divergence(const unsigned int q_point) const;
855 
865  get_symmetric_gradient(const unsigned int q_point) const;
866 
874  get_curl(const unsigned int q_point) const;
875 
891  void
892  submit_divergence(const VectorizedArrayType div_in,
893  const unsigned int q_point);
894 
911  void
912  submit_symmetric_gradient(
914  const unsigned int q_point);
915 
928  void
929  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
930  const unsigned int q_point);
931 
932 #endif
933 
950  value_type
951  integrate_value() const;
952 
954 
967  const VectorizedArrayType *
968  begin_dof_values() const;
969 
978  VectorizedArrayType *
979  begin_dof_values();
980 
991  const VectorizedArrayType *
992  begin_values() const;
993 
1004  VectorizedArrayType *
1005  begin_values();
1006 
1018  const VectorizedArrayType *
1019  begin_gradients() const;
1020 
1032  VectorizedArrayType *
1033  begin_gradients();
1034 
1047  const VectorizedArrayType *
1048  begin_hessians() const;
1049 
1062  VectorizedArrayType *
1063  begin_hessians();
1064 
1066 
1067 protected:
1078  const unsigned int dof_no,
1079  const unsigned int first_selected_component,
1080  const unsigned int quad_no,
1081  const unsigned int fe_degree,
1082  const unsigned int n_q_points,
1083  const bool is_interior_face);
1084 
1121  const Mapping<dim> & mapping,
1122  const FiniteElement<dim> &fe,
1123  const Quadrature<1> & quadrature,
1124  const UpdateFlags update_flags,
1125  const unsigned int first_selected_component,
1127  *other);
1128 
1135  FEEvaluationBase(const FEEvaluationBase &other);
1136 
1144  operator=(const FEEvaluationBase &other);
1145 
1152  template <typename VectorType, typename VectorOperation>
1153  void
1154  read_write_operation(const VectorOperation &operation,
1155  VectorType * vectors[],
1156  const std::bitset<VectorizedArrayType::size()> &mask,
1157  const bool apply_constraints = true) const;
1158 
1166  template <typename VectorType, typename VectorOperation>
1167  void
1168  read_write_operation_contiguous(
1169  const VectorOperation & operation,
1170  VectorType * vectors[],
1171  const std::bitset<VectorizedArrayType::size()> &mask) const;
1172 
1180  template <typename VectorType, typename VectorOperation>
1181  void
1182  read_write_operation_global(const VectorOperation &operation,
1183  VectorType * vectors[]) const;
1184 
1197  VectorizedArrayType *values_dofs[n_components];
1198 
1210  VectorizedArrayType *values_quad;
1211 
1225  VectorizedArrayType *gradients_quad;
1226 
1238  VectorizedArrayType *hessians_quad;
1239 
1244  const unsigned int n_fe_components;
1245 
1252 
1259 
1266 
1273 
1280 
1287 
1292  const unsigned int first_selected_component;
1293 
1298  mutable std::vector<types::global_dof_index> local_dof_indices;
1299 
1300 private:
1305  void
1306  set_data_pointers();
1307 };
1308 
1309 
1310 
1318 template <int dim,
1319  int n_components_,
1320  typename Number,
1321  bool is_face,
1322  typename VectorizedArrayType = VectorizedArray<Number>>
1324  n_components_,
1325  Number,
1326  is_face,
1327  VectorizedArrayType>
1328 {
1329  static_assert(
1330  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1331  "Type of Number and of VectorizedArrayType do not match.");
1332 
1333 public:
1334  using number_type = Number;
1336  using gradient_type =
1338  static constexpr unsigned int dimension = dim;
1339  static constexpr unsigned int n_components = n_components_;
1340  using BaseClass =
1342 
1343 protected:
1353  const unsigned int dof_no,
1354  const unsigned int first_selected_component,
1355  const unsigned int quad_no,
1356  const unsigned int fe_degree,
1357  const unsigned int n_q_points,
1358  const bool is_interior_face = true);
1359 
1365  const Mapping<dim> & mapping,
1366  const FiniteElement<dim> &fe,
1367  const Quadrature<1> & quadrature,
1368  const UpdateFlags update_flags,
1369  const unsigned int first_selected_component,
1371  *other);
1372 
1376  FEEvaluationAccess(const FEEvaluationAccess &other);
1377 
1382  operator=(const FEEvaluationAccess &other);
1383 };
1384 
1385 
1386 
1395 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1396 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1397  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1398 {
1399  static_assert(
1400  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1401  "Type of Number and of VectorizedArrayType do not match.");
1402 
1403 public:
1404  using number_type = Number;
1405  using value_type = VectorizedArrayType;
1407  static constexpr unsigned int dimension = dim;
1408  using BaseClass =
1410 
1413  value_type
1414  get_dof_value(const unsigned int dof) const;
1415 
1418  void
1419  submit_dof_value(const value_type val_in, const unsigned int dof);
1420 
1423  value_type
1424  get_value(const unsigned int q_point) const;
1425 
1428  void
1429  submit_value(const value_type val_in, const unsigned int q_point);
1430 
1433  void
1434  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
1435  const unsigned int q_point);
1436 
1440  get_gradient(const unsigned int q_point) const;
1441 
1444  value_type
1445  get_normal_derivative(const unsigned int q_point) const;
1446 
1449  void
1450  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1451 
1454  void
1455  submit_normal_derivative(const value_type grad_in,
1456  const unsigned int q_point);
1457 
1461  get_hessian(unsigned int q_point) const;
1462 
1466  get_hessian_diagonal(const unsigned int q_point) const;
1467 
1470  value_type
1471  get_laplacian(const unsigned int q_point) const;
1472 
1475  value_type
1476  integrate_value() const;
1477 
1478 protected:
1488  const unsigned int dof_no,
1489  const unsigned int first_selected_component,
1490  const unsigned int quad_no,
1491  const unsigned int fe_degree,
1492  const unsigned int n_q_points,
1493  const bool is_interior_face = true);
1494 
1500  const Mapping<dim> & mapping,
1501  const FiniteElement<dim> &fe,
1502  const Quadrature<1> & quadrature,
1503  const UpdateFlags update_flags,
1504  const unsigned int first_selected_component,
1506  *other);
1507 
1511  FEEvaluationAccess(const FEEvaluationAccess &other);
1512 
1517  operator=(const FEEvaluationAccess &other);
1518 };
1519 
1520 
1521 
1531 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1532 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1533  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1534 {
1535  static_assert(
1536  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1537  "Type of Number and of VectorizedArrayType do not match.");
1538 
1539 public:
1540  using number_type = Number;
1543  static constexpr unsigned int dimension = dim;
1544  static constexpr unsigned int n_components = dim;
1545  using BaseClass =
1547 
1551  get_gradient(const unsigned int q_point) const;
1552 
1557  VectorizedArrayType
1558  get_divergence(const unsigned int q_point) const;
1559 
1567  get_symmetric_gradient(const unsigned int q_point) const;
1568 
1574  get_curl(const unsigned int q_point) const;
1575 
1579  get_hessian(const unsigned int q_point) const;
1580 
1584  get_hessian_diagonal(const unsigned int q_point) const;
1585 
1588  void
1589  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1590 
1599  void
1600  submit_gradient(
1601  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1602  const unsigned int q_point);
1603 
1612  void
1613  submit_divergence(const VectorizedArrayType div_in,
1614  const unsigned int q_point);
1615 
1624  void
1625  submit_symmetric_gradient(
1627  const unsigned int q_point);
1628 
1633  void
1634  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1635  const unsigned int q_point);
1636 
1637 protected:
1647  const unsigned int dof_no,
1648  const unsigned int first_selected_component,
1649  const unsigned int quad_no,
1650  const unsigned int dofs_per_cell,
1651  const unsigned int n_q_points,
1652  const bool is_interior_face = true);
1653 
1659  const Mapping<dim> & mapping,
1660  const FiniteElement<dim> &fe,
1661  const Quadrature<1> & quadrature,
1662  const UpdateFlags update_flags,
1663  const unsigned int first_selected_component,
1665  *other);
1666 
1670  FEEvaluationAccess(const FEEvaluationAccess &other);
1671 
1676  operator=(const FEEvaluationAccess &other);
1677 };
1678 
1679 
1690 template <typename Number, bool is_face, typename VectorizedArrayType>
1691 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1692  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1693 {
1694  static_assert(
1695  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1696  "Type of Number and of VectorizedArrayType do not match.");
1697 
1698 public:
1699  using number_type = Number;
1700  using value_type = VectorizedArrayType;
1702  static constexpr unsigned int dimension = 1;
1703  using BaseClass =
1705 
1708  value_type
1709  get_dof_value(const unsigned int dof) const;
1710 
1713  void
1714  submit_dof_value(const value_type val_in, const unsigned int dof);
1715 
1718  value_type
1719  get_value(const unsigned int q_point) const;
1720 
1723  void
1724  submit_value(const value_type val_in, const unsigned int q_point);
1725 
1728  void
1729  submit_value(const gradient_type val_in, const unsigned int q_point);
1730 
1734  get_gradient(const unsigned int q_point) const;
1735 
1738  value_type
1739  get_normal_derivative(const unsigned int q_point) const;
1740 
1743  void
1744  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1745 
1748  void
1749  submit_gradient(const value_type grad_in, const unsigned int q_point);
1750 
1753  void
1754  submit_normal_derivative(const value_type grad_in,
1755  const unsigned int q_point);
1756 
1759  void
1760  submit_normal_derivative(const gradient_type grad_in,
1761  const unsigned int q_point);
1762 
1766  get_hessian(unsigned int q_point) const;
1767 
1771  get_hessian_diagonal(const unsigned int q_point) const;
1772 
1775  value_type
1776  get_laplacian(const unsigned int q_point) const;
1777 
1780  value_type
1781  integrate_value() const;
1782 
1783 protected:
1792  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
1793  const unsigned int dof_no,
1794  const unsigned int first_selected_component,
1795  const unsigned int quad_no,
1796  const unsigned int fe_degree,
1797  const unsigned int n_q_points,
1798  const bool is_interior_face = true);
1799 
1805  const Mapping<1> & mapping,
1806  const FiniteElement<1> &fe,
1807  const Quadrature<1> & quadrature,
1808  const UpdateFlags update_flags,
1809  const unsigned int first_selected_component,
1811 
1815  FEEvaluationAccess(const FEEvaluationAccess &other);
1816 
1821  operator=(const FEEvaluationAccess &other);
1822 };
1823 
1824 
1825 
2377 template <int dim,
2378  int fe_degree,
2379  int n_q_points_1d,
2380  int n_components_,
2381  typename Number,
2382  typename VectorizedArrayType>
2383 class FEEvaluation : public FEEvaluationAccess<dim,
2384  n_components_,
2385  Number,
2386  false,
2387  VectorizedArrayType>
2388 {
2389  static_assert(
2390  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2391  "Type of Number and of VectorizedArrayType do not match.");
2392 
2393 public:
2397  using BaseClass =
2399 
2403  using number_type = Number;
2404 
2411 
2418 
2422  static constexpr unsigned int dimension = dim;
2423 
2428  static constexpr unsigned int n_components = n_components_;
2429 
2436  static constexpr unsigned int static_n_q_points =
2437  Utilities::pow(n_q_points_1d, dim);
2438 
2446  static constexpr unsigned int static_dofs_per_component =
2447  Utilities::pow(fe_degree + 1, dim);
2448 
2456  static constexpr unsigned int tensor_dofs_per_cell =
2457  static_dofs_per_component * n_components;
2458 
2466  static constexpr unsigned int static_dofs_per_cell =
2467  static_dofs_per_component * n_components;
2468 
2495  const unsigned int dof_no = 0,
2496  const unsigned int quad_no = 0,
2497  const unsigned int first_selected_component = 0);
2498 
2525  FEEvaluation(const Mapping<dim> & mapping,
2526  const FiniteElement<dim> &fe,
2527  const Quadrature<1> & quadrature,
2528  const UpdateFlags update_flags,
2529  const unsigned int first_selected_component = 0);
2530 
2536  FEEvaluation(const FiniteElement<dim> &fe,
2537  const Quadrature<1> & quadrature,
2538  const UpdateFlags update_flags,
2539  const unsigned int first_selected_component = 0);
2540 
2551  FEEvaluation(
2552  const FiniteElement<dim> & fe,
2554  const unsigned int first_selected_component = 0);
2555 
2562  FEEvaluation(const FEEvaluation &other);
2563 
2570  FEEvaluation &
2571  operator=(const FEEvaluation &other);
2572 
2581  void
2582  reinit(const unsigned int cell_batch_index);
2583 
2596  template <bool level_dof_access>
2597  void
2599 
2610  void
2611  reinit(const typename Triangulation<dim>::cell_iterator &cell);
2612 
2622  void
2623  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2624 
2629  void
2630  evaluate(const bool evaluate_values,
2631  const bool evaluate_gradients,
2632  const bool evaluate_hessians = false);
2633 
2646  void
2647  evaluate(const VectorizedArrayType * values_array,
2648  const EvaluationFlags::EvaluationFlags evaluation_flag);
2649 
2654  void
2655  evaluate(const VectorizedArrayType *values_array,
2656  const bool evaluate_values,
2657  const bool evaluate_gradients,
2658  const bool evaluate_hessians = false);
2659 
2673  template <typename VectorType>
2674  void
2675  gather_evaluate(const VectorType & input_vector,
2676  const EvaluationFlags::EvaluationFlags evaluation_flag);
2677 
2681  template <typename VectorType>
2682  void
2683  gather_evaluate(const VectorType &input_vector,
2684  const bool evaluate_values,
2685  const bool evaluate_gradients,
2686  const bool evaluate_hessians = false);
2687 
2698  void
2699  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2700 
2701 
2705  void
2706  integrate(const bool integrate_values, const bool integrate_gradients);
2707 
2719  void
2720  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2721  VectorizedArrayType * values_array);
2722 
2726  void
2727  integrate(const bool integrate_values,
2728  const bool integrate_gradients,
2729  VectorizedArrayType *values_array);
2730 
2744  template <typename VectorType>
2745  void
2746  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
2747  VectorType & output_vector);
2748 
2752  template <typename VectorType>
2753  void
2754  integrate_scatter(const bool integrate_values,
2755  const bool integrate_gradients,
2756  VectorType &output_vector);
2757 
2763  quadrature_point(const unsigned int q_point) const;
2764 
2771  const unsigned int dofs_per_component;
2772 
2779  const unsigned int dofs_per_cell;
2780 
2788  const unsigned int n_q_points;
2789 
2790 private:
2795  void
2796  check_template_arguments(const unsigned int fe_no,
2797  const unsigned int first_selected_component);
2798 };
2799 
2800 
2801 
2837 template <int dim,
2838  int fe_degree,
2839  int n_q_points_1d = fe_degree + 1,
2840  int n_components_ = 1,
2841  typename Number = double,
2842  typename VectorizedArrayType = VectorizedArray<Number>>
2844  n_components_,
2845  Number,
2846  true,
2847  VectorizedArrayType>
2848 {
2849  static_assert(
2850  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2851  "Type of Number and of VectorizedArrayType do not match.");
2852 
2853 public:
2857  using BaseClass =
2859 
2863  using number_type = Number;
2864 
2871 
2878 
2882  static constexpr unsigned int dimension = dim;
2883 
2888  static constexpr unsigned int n_components = n_components_;
2889 
2897  static constexpr unsigned int static_n_q_points =
2898  Utilities::pow(n_q_points_1d, dim - 1);
2899 
2906  static constexpr unsigned int static_n_q_points_cell =
2907  Utilities::pow(n_q_points_1d, dim);
2908 
2915  static constexpr unsigned int static_dofs_per_component =
2916  Utilities::pow(fe_degree + 1, dim);
2917 
2924  static constexpr unsigned int tensor_dofs_per_cell =
2925  static_dofs_per_component * n_components;
2926 
2933  static constexpr unsigned int static_dofs_per_cell =
2934  static_dofs_per_component * n_components;
2935 
2967  const bool is_interior_face = true,
2968  const unsigned int dof_no = 0,
2969  const unsigned int quad_no = 0,
2970  const unsigned int first_selected_component = 0);
2971 
2982  void
2983  reinit(const unsigned int face_batch_number);
2984 
2992  void
2993  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2994 
3005  void
3006  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3007 
3011  void
3012  evaluate(const bool evaluate_values, const bool evaluate_gradients);
3013 
3026  void
3027  evaluate(const VectorizedArrayType * values_array,
3028  const EvaluationFlags::EvaluationFlags evaluation_flag);
3029 
3033  void
3034  evaluate(const VectorizedArrayType *values_array,
3035  const bool evaluate_values,
3036  const bool evaluate_gradients);
3037 
3049  template <typename VectorType>
3050  void
3051  gather_evaluate(const VectorType & input_vector,
3052  const EvaluationFlags::EvaluationFlags evaluation_flag);
3053 
3057  template <typename VectorType>
3058  void
3059  gather_evaluate(const VectorType &input_vector,
3060  const bool evaluate_values,
3061  const bool evaluate_gradients);
3062 
3072  void
3073  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag);
3074 
3078  void
3079  integrate(const bool integrate_values, const bool integrate_gradients);
3080 
3089  void
3090  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag,
3091  VectorizedArrayType * values_array);
3092 
3096  void
3097  integrate(const bool integrate_values,
3098  const bool integrate_gradients,
3099  VectorizedArrayType *values_array);
3100 
3112  template <typename VectorType>
3113  void
3114  integrate_scatter(const EvaluationFlags::EvaluationFlags evaluation_flag,
3115  VectorType & output_vector);
3116 
3120  template <typename VectorType>
3121  void
3122  integrate_scatter(const bool integrate_values,
3123  const bool integrate_gradients,
3124  VectorType &output_vector);
3125 
3131  quadrature_point(const unsigned int q_point) const;
3132 
3139  const unsigned int dofs_per_component;
3140 
3147  const unsigned int dofs_per_cell;
3148 
3156  const unsigned int n_q_points;
3157 
3158 
3159 private:
3163  std::array<unsigned int, VectorizedArrayType::size()>
3164  compute_face_no_data();
3165 
3169  std::array<unsigned int, VectorizedArrayType::size()>
3170  compute_face_orientations();
3171 };
3172 
3173 
3174 
3175 namespace internal
3176 {
3177  namespace MatrixFreeFunctions
3178  {
3179  // a helper function to compute the number of DoFs of a DGP element at
3180  // compile time, depending on the degree
3181  template <int dim, int degree>
3183  {
3184  // this division is always without remainder
3185  static constexpr unsigned int value =
3186  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
3187  };
3188 
3189  // base specialization: 1d elements have 'degree+1' degrees of freedom
3190  template <int degree>
3191  struct DGP_dofs_per_component<1, degree>
3192  {
3193  static constexpr unsigned int value = degree + 1;
3194  };
3195  } // namespace MatrixFreeFunctions
3196 } // namespace internal
3197 
3198 
3199 /*----------------------- Inline functions ----------------------------------*/
3200 
3201 #ifndef DOXYGEN
3202 
3203 
3204 /*----------------------- FEEvaluationBaseData ------------------------*/
3205 
3206 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3210  const unsigned int dof_no,
3211  const unsigned int first_selected_component,
3212  const unsigned int quad_no_in,
3213  const unsigned int fe_degree,
3214  const unsigned int n_q_points,
3215  const bool is_interior_face)
3216  : scratch_data_array(data_in.acquire_scratch_data())
3217  , quad_no(quad_no_in)
3218  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
3219  data_in.get_dof_info(dof_no).fe_index_from_degree(
3220  first_selected_component,
3221  fe_degree) :
3222  0)
3223  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
3224  (is_face ? data_in.get_mapping_info()
3225  .face_data[quad_no_in]
3226  .quad_index_from_n_q_points(n_q_points) :
3227  data_in.get_mapping_info()
3228  .cell_data[quad_no_in]
3229  .quad_index_from_n_q_points(n_q_points)) :
3230  0)
3231  , n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ?
3232  n_q_points :
3233  (is_face ? data_in
3234  .get_shape_info(dof_no,
3235  quad_no_in,
3236  active_fe_index,
3237  active_quad_index)
3238  .n_q_points_face :
3239  data_in
3240  .get_shape_info(dof_no,
3241  quad_no_in,
3242  active_fe_index,
3243  active_quad_index)
3244  .n_q_points))
3245  , matrix_info(&data_in)
3246  , dof_info(&data_in.get_dof_info(dof_no))
3247  , mapping_data(
3248  internal::MatrixFreeFunctions::
3249  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3250  data_in.get_mapping_info(),
3251  quad_no))
3252  , data(&data_in.get_shape_info(
3253  dof_no,
3254  quad_no_in,
3255  dof_info->component_to_base_index[first_selected_component],
3256  active_fe_index,
3257  active_quad_index))
3258  , jacobian(nullptr)
3259  , J_value(nullptr)
3260  , normal_vectors(nullptr)
3261  , normal_x_jacobian(nullptr)
3262  , quadrature_weights(
3263  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3264  , cell(numbers::invalid_unsigned_int)
3265  , is_interior_face(is_interior_face)
3266  , dof_access_index(
3267  is_face ?
3268  (is_interior_face ?
3269  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
3270  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
3271  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3272  , cell_type(internal::MatrixFreeFunctions::general)
3273 {
3274  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3275  AssertDimension(matrix_info->get_task_info().vectorization_length,
3276  VectorizedArrayType::size());
3277  AssertDimension((is_face ? data->n_q_points_face : data->n_q_points),
3278  n_quadrature_points);
3279  AssertDimension(n_quadrature_points,
3280  mapping_data->descriptor[active_quad_index].n_q_points);
3281 }
3282 
3283 
3284 
3285 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3288  const Mapping<dim> & mapping,
3289  const FiniteElement<dim> &fe,
3290  const Quadrature<1> & quadrature,
3291  const UpdateFlags update_flags,
3292  const unsigned int first_selected_component,
3294  *other)
3295  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3297  , active_fe_index(numbers::invalid_unsigned_int)
3298  , active_quad_index(numbers::invalid_unsigned_int)
3299  , n_quadrature_points(
3300  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3301  , matrix_info(nullptr)
3302  , dof_info(nullptr)
3303  , mapping_data(nullptr)
3304  ,
3305  // select the correct base element from the given FE component
3307  quadrature,
3308  fe,
3309  fe.component_to_base_index(first_selected_component).first))
3310  , jacobian(nullptr)
3311  , J_value(nullptr)
3312  , normal_vectors(nullptr)
3313  , normal_x_jacobian(nullptr)
3314  , quadrature_weights(nullptr)
3315  , cell(0)
3317  , is_interior_face(true)
3319 {
3320  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3321  ExcInternalError());
3322  if (other != nullptr &&
3323  other->mapped_geometry->get_quadrature() == quadrature)
3324  mapped_geometry = other->mapped_geometry;
3325  else
3326  mapped_geometry =
3327  std::make_shared<internal::MatrixFreeFunctions::
3329  mapping, quadrature, update_flags);
3330  cell = 0;
3331 
3332  mapping_data = &mapped_geometry->get_data_storage();
3333  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3334  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3335 }
3336 
3337 
3338 
3339 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3343  &other)
3344  : scratch_data_array(other.matrix_info == nullptr ?
3346  other.matrix_info->acquire_scratch_data())
3347  , quad_no(other.quad_no)
3348  , active_fe_index(other.active_fe_index)
3349  , active_quad_index(other.active_quad_index)
3350  , n_quadrature_points(other.n_quadrature_points)
3351  , matrix_info(other.matrix_info)
3352  , dof_info(other.dof_info)
3353  , mapping_data(other.mapping_data)
3354  , data(other.matrix_info == nullptr ?
3356  *other.data) :
3357  other.data)
3358  , jacobian(nullptr)
3359  , J_value(nullptr)
3360  , normal_vectors(nullptr)
3361  , normal_x_jacobian(nullptr)
3362  , quadrature_weights(
3363  other.matrix_info == nullptr ?
3364  nullptr :
3365  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3368  , is_interior_face(other.is_interior_face)
3369  , dof_access_index(other.dof_access_index)
3370 {
3371  // Create deep copy of mapped geometry for use in parallel...
3372  if (other.mapped_geometry.get() != nullptr)
3373  {
3374  mapped_geometry = std::make_shared<
3375  internal::MatrixFreeFunctions::
3376  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3377  other.mapped_geometry->get_fe_values().get_mapping(),
3378  other.mapped_geometry->get_quadrature(),
3379  other.mapped_geometry->get_fe_values().get_update_flags());
3380  mapping_data = &mapped_geometry->get_data_storage();
3381  cell = 0;
3382 
3383  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3384  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3385  }
3386 }
3387 
3388 
3389 
3390 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3394 {
3395  AssertDimension(quad_no, other.quad_no);
3396  AssertDimension(active_fe_index, other.active_fe_index);
3397  AssertDimension(active_quad_index, other.active_quad_index);
3398 
3399  // release old memory
3400  if (matrix_info == nullptr)
3401  {
3402  delete data;
3403  delete scratch_data_array;
3404  }
3405  else
3406  {
3407  matrix_info->release_scratch_data(scratch_data_array);
3408  }
3409 
3410  matrix_info = other.matrix_info;
3411  dof_info = other.dof_info;
3412  mapping_data = other.mapping_data;
3413  if (other.matrix_info == nullptr)
3414  {
3416  *other.data);
3417  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3418  }
3419  else
3420  {
3421  data = other.data;
3422  scratch_data_array = matrix_info->acquire_scratch_data();
3423  }
3424 
3425  quadrature_weights =
3426  (mapping_data != nullptr ?
3427  mapping_data->descriptor[active_quad_index].quadrature_weights.begin() :
3428  nullptr);
3431  is_interior_face = other.is_interior_face;
3432  dof_access_index = other.dof_access_index;
3433 
3434  // Create deep copy of mapped geometry for use in parallel...
3435  if (other.mapped_geometry.get() != nullptr)
3436  {
3437  mapped_geometry = std::make_shared<
3438  internal::MatrixFreeFunctions::
3439  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3440  other.mapped_geometry->get_fe_values().get_mapping(),
3441  other.mapped_geometry->get_quadrature(),
3442  other.mapped_geometry->get_fe_values().get_update_flags());
3443  cell = 0;
3444  mapping_data = &mapped_geometry->get_data_storage();
3445  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3446  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3447  }
3448 
3449  return *this;
3450 }
3451 
3452 
3453 
3454 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3457 {
3458  if (matrix_info != nullptr)
3459  {
3460  try
3461  {
3462  matrix_info->release_scratch_data(scratch_data_array);
3463  }
3464  catch (...)
3465  {}
3466  }
3467  else
3468  {
3469  delete scratch_data_array;
3470  delete data;
3471  data = nullptr;
3472  }
3473  scratch_data_array = nullptr;
3474 }
3475 
3476 
3477 
3478 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3479 inline unsigned int
3482 {
3483  if (matrix_info == nullptr)
3484  return 0;
3485  else
3486  {
3487  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3488  return this->mapping_data->data_index_offsets[cell];
3489  }
3490 }
3491 
3492 
3493 
3494 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3497  const
3498 {
3500  return cell_type;
3501 }
3502 
3503 
3504 
3505 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3508  get_shape_info() const
3509 {
3510  Assert(data != nullptr, ExcInternalError());
3511  return *data;
3512 }
3513 
3514 
3515 
3516 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3519  get_normal_vector(const unsigned int q_point) const
3520 {
3521  AssertIndexRange(q_point, n_quadrature_points);
3522  Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!"));
3523  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3524  return normal_vectors[0];
3525  else
3526  return normal_vectors[q_point];
3527 }
3528 
3529 
3530 
3531 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3532 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3534  const unsigned int q_point) const
3535 {
3536  AssertIndexRange(q_point, n_quadrature_points);
3537  Assert(J_value != nullptr, ExcNotInitialized());
3538  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3539  {
3540  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3541  return J_value[0] * this->quadrature_weights[q_point];
3542  }
3543  else
3544  return J_value[q_point];
3545 }
3546 
3547 
3548 
3549 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3552  inverse_jacobian(const unsigned int q_point) const
3553 {
3554  AssertIndexRange(q_point, n_quadrature_points);
3555  Assert(this->jacobian != nullptr, ExcNotImplemented());
3556  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3557  return jacobian[0];
3558  else
3559  return jacobian[q_point];
3560 }
3561 
3562 
3563 
3564 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3565 inline std::array<unsigned int, VectorizedArrayType::size()>
3567  const
3568 {
3569  Assert(this->matrix_info != nullptr, ExcNotInitialized());
3570 
3571  const unsigned int n_lanes = VectorizedArrayType::size();
3572  std::array<unsigned int, n_lanes> cells;
3573 
3574  // initialize array
3575  for (unsigned int i = 0; i < n_lanes; ++i)
3576  cells[i] = numbers::invalid_unsigned_int;
3577 
3578  if ((is_face == false) ||
3579  (is_face &&
3580  this->dof_access_index ==
3582  this->is_interior_face))
3583  {
3584  // cell or interior face face (element-centric loop)
3585  for (unsigned int i = 0; i < n_lanes; ++i)
3586  cells[i] = cell * n_lanes + i;
3587  }
3588  else if (is_face &&
3589  this->dof_access_index ==
3591  this->is_interior_face == false)
3592  {
3593  // exterior face (element-centric loop): for this case, we need to
3594  // look into the FaceInfo field that collects information from both
3595  // sides of a face once for the global mesh, and pick the face id that
3596  // is not the local one (cell_this).
3597  for (unsigned int i = 0; i < n_lanes; i++)
3598  {
3599  // compute actual (non vectorized) cell ID
3600  const unsigned int cell_this = this->cell * n_lanes + i;
3601  // compute face ID
3602  unsigned int face_index =
3603  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3604  this->face_no,
3605  i);
3606 
3607  if (face_index == numbers::invalid_unsigned_int)
3608  continue; // invalid face ID: no neighbor on boundary
3609 
3610  // get cell ID on both sides of face
3611  auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3612  .cells_interior[face_index % n_lanes];
3613  auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3614  .cells_exterior[face_index % n_lanes];
3615 
3616  // compare the IDs with the given cell ID
3617  if (cell_m == cell_this)
3618  cells[i] = cell_p; // neighbor has the other ID
3619  else if (cell_p == cell_this)
3620  cells[i] = cell_m;
3621  }
3622  }
3623  else if (is_face)
3624  {
3625  // face-centric faces
3626  const unsigned int *cells_ =
3627  is_interior_face ?
3628  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3629  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3630  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3631  if (cells_[i] != numbers::invalid_unsigned_int)
3632  cells[i] = cells_[i];
3633  }
3634 
3635  return cells;
3636 }
3637 
3638 
3639 namespace internal
3640 {
3641  template <int dim,
3642  typename Number,
3643  bool is_face,
3644  typename VectorizedArrayType,
3645  typename VectorizedArrayType2,
3646  typename GlobalVectorType,
3647  typename FU>
3648  inline void
3649  process_cell_data(
3652  GlobalVectorType & array,
3653  VectorizedArrayType2 & out,
3654  const FU & fu)
3655  {
3656  (void)matrix_info;
3657  Assert(matrix_info != nullptr, ExcNotImplemented());
3658  AssertDimension(array.size(),
3659  matrix_info->get_task_info().cell_partition_data.back());
3660 
3661  // 1) collect ids of cell
3662  const auto cells = phi.get_cell_ids();
3663 
3664  // 2) actually gather values
3665  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3666  if (cells[i] != numbers::invalid_unsigned_int)
3667  fu(out[i],
3668  array[cells[i] / VectorizedArrayType::size()]
3669  [cells[i] % VectorizedArrayType::size()]);
3670  }
3671 } // namespace internal
3672 
3673 
3674 
3675 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3676 std::array<unsigned int, VectorizedArrayType::size()>
3678  get_cell_or_face_ids() const
3679 {
3680  const unsigned int v_len = VectorizedArrayType::size();
3681  std::array<unsigned int, VectorizedArrayType::size()> cells;
3682 
3683  // initialize array
3684  for (unsigned int i = 0; i < v_len; ++i)
3685  cells[i] = numbers::invalid_unsigned_int;
3686 
3687  if (is_face &&
3688  this->dof_access_index ==
3690  this->is_interior_face == false)
3691  {
3692  // cell-based face-loop: plus face
3693  for (unsigned int i = 0; i < v_len; i++)
3694  {
3695  // compute actual (non vectorized) cell ID
3696  const unsigned int cell_this = this->cell * v_len + i;
3697  // compute face ID
3698  unsigned int fn =
3699  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3700  this->face_no,
3701  i);
3702 
3704  continue; // invalid face ID: no neighbor on boundary
3705 
3706  // get cell ID on both sides of face
3707  auto cell_m = this->matrix_info->get_face_info(fn / v_len)
3708  .cells_interior[fn % v_len];
3709  auto cell_p = this->matrix_info->get_face_info(fn / v_len)
3710  .cells_exterior[fn % v_len];
3711 
3712  // compare the IDs with the given cell ID
3713  if (cell_m == cell_this)
3714  cells[i] = cell_p; // neighbor has the other ID
3715  else if (cell_p == cell_this)
3716  cells[i] = cell_m;
3717  }
3718  }
3719  else
3720  {
3721  for (unsigned int i = 0; i < v_len; ++i)
3722  cells[i] = cell * v_len + i;
3723  }
3724 
3725  return cells;
3726 }
3727 
3728 
3729 
3730 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3731 inline VectorizedArrayType
3733  const AlignedVector<VectorizedArrayType> &array) const
3734 {
3735  VectorizedArrayType out = Number(1.);
3736  internal::process_cell_data(
3737  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3738  local = global;
3739  });
3740  return out;
3741 }
3742 
3743 
3744 
3745 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3746 inline void
3749  const VectorizedArrayType & in) const
3750 {
3751  internal::process_cell_data(
3752  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3753  global = local;
3754  });
3755 }
3756 
3757 
3758 
3759 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3760 template <typename T>
3761 inline std::array<T, VectorizedArrayType::size()>
3763  const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
3764 {
3765  std::array<T, VectorizedArrayType::size()> out;
3766  internal::process_cell_data(
3767  *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
3768  local = global;
3769  });
3770  return out;
3771 }
3772 
3773 
3774 
3775 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3776 template <typename T>
3777 inline void
3779  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
3780  const std::array<T, VectorizedArrayType::size()> & in) const
3781 {
3782  internal::process_cell_data(
3783  *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
3784  global = local;
3785  });
3786 }
3787 
3788 
3789 
3790 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3791 inline const std::vector<unsigned int> &
3794 {
3795  return data->lexicographic_numbering;
3796 }
3797 
3798 
3799 
3800 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3803  get_scratch_data() const
3804 {
3806  const_cast<VectorizedArrayType *>(scratch_data),
3807  scratch_data_array->end() - scratch_data);
3808 }
3809 
3810 
3811 
3812 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3813 inline unsigned int
3815  get_quadrature_index() const
3816 {
3817  return this->quad_no;
3818 }
3819 
3820 
3821 
3822 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
3823 inline unsigned int
3826 {
3827  if (is_face && this->dof_access_index ==
3829  return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
3830  else
3831  return this->cell;
3832 }
3833 
3834 
3835 /*----------------------- FEEvaluationBase ----------------------------------*/
3836 
3837 template <int dim,
3838  int n_components_,
3839  typename Number,
3840  bool is_face,
3841  typename VectorizedArrayType>
3842 inline FEEvaluationBase<dim,
3843  n_components_,
3844  Number,
3845  is_face,
3846  VectorizedArrayType>::
3847  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
3848  const unsigned int dof_no,
3849  const unsigned int first_selected_component,
3850  const unsigned int quad_no_in,
3851  const unsigned int fe_degree,
3852  const unsigned int n_q_points,
3853  const bool is_interior_face)
3855  data_in,
3856  dof_no,
3857  first_selected_component,
3858  quad_no_in,
3859  fe_degree,
3860  n_q_points,
3861  is_interior_face)
3862  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
3863  , dof_values_initialized(false)
3864  , values_quad_initialized(false)
3865  , gradients_quad_initialized(false)
3866  , hessians_quad_initialized(false)
3867  , values_quad_submitted(false)
3868  , gradients_quad_submitted(false)
3869  , first_selected_component(first_selected_component)
3870 {
3871  set_data_pointers();
3872  Assert(
3873  this->dof_info->start_components.back() == 1 ||
3874  static_cast<int>(n_components_) <=
3875  static_cast<int>(
3876  this->dof_info->start_components
3877  [this->dof_info->component_to_base_index[first_selected_component] +
3878  1]) -
3879  first_selected_component,
3880  ExcMessage(
3881  "You tried to construct a vector-valued evaluator with " +
3882  std::to_string(n_components) +
3883  " components. However, "
3884  "the current base element has only " +
3886  this->dof_info->start_components
3887  [this->dof_info->component_to_base_index[first_selected_component] +
3888  1] -
3889  first_selected_component) +
3890  " components left when starting from local element index " +
3892  first_selected_component -
3893  this->dof_info->start_components
3894  [this->dof_info->component_to_base_index[first_selected_component]]) +
3895  " (global index " + std::to_string(first_selected_component) + ")"));
3896 
3897  // do not check for correct dimensions of data fields here, should be done
3898  // in derived classes
3899 }
3900 
3901 
3902 
3903 template <int dim,
3904  int n_components_,
3905  typename Number,
3906  bool is_face,
3907  typename VectorizedArrayType>
3908 inline FEEvaluationBase<dim,
3909  n_components_,
3910  Number,
3911  is_face,
3912  VectorizedArrayType>::
3913  FEEvaluationBase(
3914  const Mapping<dim> & mapping,
3915  const FiniteElement<dim> &fe,
3916  const Quadrature<1> & quadrature,
3917  const UpdateFlags update_flags,
3918  const unsigned int first_selected_component,
3920  *other)
3922  mapping,
3923  fe,
3924  quadrature,
3925  update_flags,
3926  first_selected_component,
3927  other)
3928  , n_fe_components(n_components_)
3929  , dof_values_initialized(false)
3930  , values_quad_initialized(false)
3931  , gradients_quad_initialized(false)
3932  , hessians_quad_initialized(false)
3933  , values_quad_submitted(false)
3934  , gradients_quad_submitted(false)
3935  // keep the number of the selected component within the current base element
3936  // for reading dof values
3937  , first_selected_component(first_selected_component)
3938 {
3939  set_data_pointers();
3940 
3941  const unsigned int base_element_number =
3942  fe.component_to_base_index(first_selected_component).first;
3943  Assert(fe.element_multiplicity(base_element_number) == 1 ||
3944  fe.element_multiplicity(base_element_number) -
3945  first_selected_component >=
3946  n_components_,
3947  ExcMessage("The underlying element must at least contain as many "
3948  "components as requested by this class"));
3949  (void)base_element_number;
3950 }
3951 
3952 
3953 
3954 template <int dim,
3955  int n_components_,
3956  typename Number,
3957  bool is_face,
3958  typename VectorizedArrayType>
3959 inline FEEvaluationBase<dim,
3960  n_components_,
3961  Number,
3962  is_face,
3963  VectorizedArrayType>::
3964  FEEvaluationBase(const FEEvaluationBase<dim,
3965  n_components_,
3966  Number,
3967  is_face,
3968  VectorizedArrayType> &other)
3970  , n_fe_components(other.n_fe_components)
3971  , dof_values_initialized(false)
3972  , values_quad_initialized(false)
3973  , gradients_quad_initialized(false)
3974  , hessians_quad_initialized(false)
3975  , values_quad_submitted(false)
3976  , gradients_quad_submitted(false)
3977  , first_selected_component(other.first_selected_component)
3978 {
3979  set_data_pointers();
3980 }
3981 
3982 
3983 
3984 template <int dim,
3985  int n_components_,
3986  typename Number,
3987  bool is_face,
3988  typename VectorizedArrayType>
3989 inline FEEvaluationBase<dim,
3990  n_components_,
3991  Number,
3992  is_face,
3993  VectorizedArrayType> &
3995 operator=(const FEEvaluationBase<dim,
3996  n_components_,
3997  Number,
3998  is_face,
3999  VectorizedArrayType> &other)
4000 {
4002  operator=(other);
4003  AssertDimension(n_fe_components, other.n_fe_components);
4004  AssertDimension(first_selected_component, other.first_selected_component);
4005 
4006  return *this;
4007 }
4008 
4009 
4010 
4011 template <int dim,
4012  int n_components_,
4013  typename Number,
4014  bool is_face,
4015  typename VectorizedArrayType>
4016 inline void
4019 {
4020  Assert(this->scratch_data_array != nullptr, ExcInternalError());
4021 
4022  const unsigned int tensor_dofs_per_component =
4023  Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
4024  const unsigned int dofs_per_component =
4025  this->data->dofs_per_component_on_cell;
4026  const unsigned int n_quadrature_points =
4027  is_face ? this->data->n_q_points_face : this->data->n_q_points;
4028 
4029  const unsigned int shift =
4030  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
4031  n_components_ * 3 +
4032  2 * n_quadrature_points;
4033  const unsigned int allocated_size =
4034  shift + n_components_ * dofs_per_component +
4035  (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4036  this->scratch_data_array->resize_fast(allocated_size);
4037 
4038  // set the pointers to the correct position in the data array
4039  for (unsigned int c = 0; c < n_components_; ++c)
4040  {
4041  values_dofs[c] =
4042  this->scratch_data_array->begin() + c * dofs_per_component;
4043  }
4044  values_quad =
4045  this->scratch_data_array->begin() + n_components * dofs_per_component;
4046  gradients_quad = this->scratch_data_array->begin() +
4047  n_components * (dofs_per_component + n_quadrature_points);
4048  hessians_quad =
4049  this->scratch_data_array->begin() +
4050  n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
4051  this->scratch_data =
4052  this->scratch_data_array->begin() + n_components_ * dofs_per_component +
4053  (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
4054 }
4055 
4056 
4057 
4058 namespace internal
4059 {
4060  // allows to select between block vectors and non-block vectors, which
4061  // allows to use a unified interface for extracting blocks on block vectors
4062  // and doing nothing on usual vectors
4063  template <typename VectorType, bool>
4064  struct BlockVectorSelector
4065  {};
4066 
4067  template <typename VectorType>
4068  struct BlockVectorSelector<VectorType, true>
4069  {
4070  using BaseVectorType = typename VectorType::BlockType;
4071 
4072  static BaseVectorType *
4073  get_vector_component(VectorType &vec, const unsigned int component)
4074  {
4075  AssertIndexRange(component, vec.n_blocks());
4076  return &vec.block(component);
4077  }
4078  };
4079 
4080  template <typename VectorType>
4081  struct BlockVectorSelector<VectorType, false>
4082  {
4083  using BaseVectorType = VectorType;
4084 
4085  static BaseVectorType *
4086  get_vector_component(VectorType &vec, const unsigned int component)
4087  {
4088  // FEEvaluation allows to combine several vectors from a scalar
4089  // FiniteElement into a "vector-valued" FEEvaluation object with
4090  // multiple components. These components can be extracted with the other
4091  // get_vector_component functions. If we do not get a vector of vectors
4092  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
4093  // must make sure that we do not duplicate the components in input
4094  // and/or duplicate the resulting integrals. In such a case, we should
4095  // only get the zeroth component in the vector contained set nullptr for
4096  // the others which allows us to catch unintended use in
4097  // read_write_operation.
4098  if (component == 0)
4099  return &vec;
4100  else
4101  return nullptr;
4102  }
4103  };
4104 
4105  template <typename VectorType>
4106  struct BlockVectorSelector<std::vector<VectorType>, false>
4107  {
4108  using BaseVectorType = VectorType;
4109 
4110  static BaseVectorType *
4111  get_vector_component(std::vector<VectorType> &vec,
4112  const unsigned int component)
4113  {
4114  AssertIndexRange(component, vec.size());
4115  return &vec[component];
4116  }
4117  };
4118 
4119  template <typename VectorType>
4120  struct BlockVectorSelector<std::vector<VectorType *>, false>
4121  {
4122  using BaseVectorType = VectorType;
4123 
4124  static BaseVectorType *
4125  get_vector_component(std::vector<VectorType *> &vec,
4126  const unsigned int component)
4127  {
4128  AssertIndexRange(component, vec.size());
4129  return vec[component];
4130  }
4131  };
4132 } // namespace internal
4133 
4134 
4135 
4136 template <int dim,
4137  int n_components_,
4138  typename Number,
4139  bool is_face,
4140  typename VectorizedArrayType>
4141 template <typename VectorType, typename VectorOperation>
4142 inline void
4144  read_write_operation(const VectorOperation &operation,
4145  VectorType * src[],
4146  const std::bitset<VectorizedArrayType::size()> &mask,
4147  const bool apply_constraints) const
4148 {
4149  // Case 1: No MatrixFree object given, simple case because we do not need to
4150  // process constraints and need not care about vectorization -> go to
4151  // separate function
4152  if (this->matrix_info == nullptr)
4153  {
4154  read_write_operation_global(operation, src);
4155  return;
4156  }
4157 
4158  Assert(this->dof_info != nullptr, ExcNotInitialized());
4159  Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
4160  if (n_fe_components == 1)
4161  for (unsigned int comp = 0; comp < n_components; ++comp)
4162  {
4163  Assert(src[comp] != nullptr,
4164  ExcMessage("The finite element underlying this FEEvaluation "
4165  "object is scalar, but you requested " +
4166  std::to_string(n_components) +
4167  " components via the template argument in "
4168  "FEEvaluation. In that case, you must pass an "
4169  "std::vector<VectorType> or a BlockVector to " +
4170  "read_dof_values and distribute_local_to_global."));
4171  internal::check_vector_compatibility(*src[comp], *this->dof_info);
4172  }
4173  else
4174  {
4175  internal::check_vector_compatibility(*src[0], *this->dof_info);
4176  }
4177 
4178  // Case 2: contiguous indices which use reduced storage of indices and can
4179  // use vectorized load/store operations -> go to separate function
4181  this->cell,
4182  this->dof_info->index_storage_variants[this->dof_access_index].size());
4183  if (this->dof_info->index_storage_variants
4184  [is_face ? this->dof_access_index :
4186  [this->cell] >=
4188  {
4189  read_write_operation_contiguous(operation, src, mask);
4190  return;
4191  }
4192 
4193  // Case 3: standard operation with one index per degree of freedom -> go on
4194  // here
4195  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4196  Assert(mask.count() == n_lanes,
4197  ExcNotImplemented("Masking currently not implemented for "
4198  "non-contiguous DoF storage"));
4199 
4200  std::integral_constant<bool,
4202  vector_selector;
4203 
4204  const unsigned int dofs_per_component =
4205  this->data->dofs_per_component_on_cell;
4206  if (this->dof_info->index_storage_variants
4207  [is_face ? this->dof_access_index :
4209  [this->cell] ==
4211  {
4212  const unsigned int *dof_indices =
4213  this->dof_info->dof_indices_interleaved.data() +
4214  this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
4215  .first +
4216  this->dof_info
4217  ->component_dof_indices_offset[this->active_fe_index]
4218  [this->first_selected_component] *
4219  n_lanes;
4220  if (n_components == 1 || n_fe_components == 1)
4221  for (unsigned int i = 0; i < dofs_per_component;
4222  ++i, dof_indices += n_lanes)
4223  for (unsigned int comp = 0; comp < n_components; ++comp)
4224  operation.process_dof_gather(dof_indices,
4225  *src[comp],
4226  0,
4227  values_dofs[comp][i],
4228  vector_selector);
4229  else
4230  for (unsigned int comp = 0; comp < n_components; ++comp)
4231  for (unsigned int i = 0; i < dofs_per_component;
4232  ++i, dof_indices += n_lanes)
4233  operation.process_dof_gather(
4234  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4235  return;
4236  }
4237 
4238  const unsigned int * dof_indices[n_lanes];
4239  VectorizedArrayType **values_dofs =
4240  const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
4241 
4242  // Assign the appropriate cell ids for face/cell case and get the pointers
4243  // to the dof indices of the cells on all lanes
4244  unsigned int cells_copied[n_lanes];
4245  const unsigned int *cells;
4246  unsigned int n_vectorization_actual =
4247  this->dof_info
4248  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
4249  bool has_constraints = false;
4250  const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
4251  if (is_face)
4252  {
4253  if (this->dof_access_index ==
4255  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4256  cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
4257  cells =
4258  this->dof_access_index ==
4260  &cells_copied[0] :
4261  (this->is_interior_face ?
4262  &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
4263  &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
4264  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4265  {
4266  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
4267  ExcInternalError());
4268  const std::pair<unsigned int, unsigned int> *my_index_start =
4269  &this->dof_info->row_starts[cells[v] * n_fe_components +
4270  first_selected_component];
4271 
4272  // check whether any of the SIMD lanes has constraints, i.e., the
4273  // constraint indicator which is the second entry of row_starts
4274  // increments on this cell
4275  if (my_index_start[n_components_read].second !=
4276  my_index_start[0].second)
4277  has_constraints = true;
4278 
4279  dof_indices[v] =
4280  this->dof_info->dof_indices.data() + my_index_start[0].first;
4281  }
4282  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4283  dof_indices[v] = nullptr;
4284  }
4285  else
4286  {
4287  AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
4288  this->dof_info->row_starts.size());
4289  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4290  {
4291  const std::pair<unsigned int, unsigned int> *my_index_start =
4292  &this->dof_info
4293  ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
4294  first_selected_component];
4295  if (my_index_start[n_components_read].second !=
4296  my_index_start[0].second)
4297  has_constraints = true;
4298  Assert(my_index_start[n_components_read].first ==
4299  my_index_start[0].first ||
4300  my_index_start[0].first < this->dof_info->dof_indices.size(),
4301  ExcIndexRange(0,
4302  my_index_start[0].first,
4303  this->dof_info->dof_indices.size()));
4304  dof_indices[v] =
4305  this->dof_info->dof_indices.data() + my_index_start[0].first;
4306  }
4307  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
4308  dof_indices[v] = nullptr;
4309  }
4310 
4311  // Case where we have no constraints throughout the whole cell: Can go
4312  // through the list of DoFs directly
4313  if (!has_constraints)
4314  {
4315  if (n_vectorization_actual < n_lanes)
4316  for (unsigned int comp = 0; comp < n_components; ++comp)
4317  for (unsigned int i = 0; i < dofs_per_component; ++i)
4318  operation.process_empty(values_dofs[comp][i]);
4319  if (n_components == 1 || n_fe_components == 1)
4320  {
4321  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4322  for (unsigned int i = 0; i < dofs_per_component; ++i)
4323  for (unsigned int comp = 0; comp < n_components; ++comp)
4324  operation.process_dof(dof_indices[v][i],
4325  *src[comp],
4326  values_dofs[comp][i][v]);
4327  }
4328  else
4329  {
4330  for (unsigned int comp = 0; comp < n_components; ++comp)
4331  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4332  for (unsigned int i = 0; i < dofs_per_component; ++i)
4333  operation.process_dof(
4334  dof_indices[v][comp * dofs_per_component + i],
4335  *src[0],
4336  values_dofs[comp][i][v]);
4337  }
4338  return;
4339  }
4340 
4341  // In the case where there are some constraints to be resolved, loop over
4342  // all vector components that are filled and then over local dofs. ind_local
4343  // holds local number on cell, index iterates over the elements of
4344  // index_local_to_global and dof_indices points to the global indices stored
4345  // in index_local_to_global
4346  if (n_vectorization_actual < n_lanes)
4347  for (unsigned int comp = 0; comp < n_components; ++comp)
4348  for (unsigned int i = 0; i < dofs_per_component; ++i)
4349  operation.process_empty(values_dofs[comp][i]);
4350  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4351  {
4352  const unsigned int cell_index =
4353  is_face ? cells[v] : this->cell * n_lanes + v;
4354  const unsigned int cell_dof_index =
4355  cell_index * n_fe_components + first_selected_component;
4356  const unsigned int n_components_read =
4357  n_fe_components > 1 ? n_components : 1;
4358  unsigned int index_indicators =
4359  this->dof_info->row_starts[cell_dof_index].second;
4360  unsigned int next_index_indicators =
4361  this->dof_info->row_starts[cell_dof_index + 1].second;
4362 
4363  // For read_dof_values_plain, redirect the dof_indices field to the
4364  // unconstrained indices
4365  if (apply_constraints == false &&
4366  this->dof_info->row_starts[cell_dof_index].second !=
4367  this->dof_info->row_starts[cell_dof_index + n_components_read]
4368  .second)
4369  {
4370  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
4372  ExcNotInitialized());
4373  dof_indices[v] =
4374  this->dof_info->plain_dof_indices.data() +
4375  this->dof_info
4376  ->component_dof_indices_offset[this->active_fe_index]
4377  [this->first_selected_component] +
4378  this->dof_info->row_starts_plain_indices[cell_index];
4379  next_index_indicators = index_indicators;
4380  }
4381 
4382  if (n_components == 1 || n_fe_components == 1)
4383  {
4384  unsigned int ind_local = 0;
4385  for (; index_indicators != next_index_indicators; ++index_indicators)
4386  {
4387  const std::pair<unsigned short, unsigned short> indicator =
4388  this->dof_info->constraint_indicator[index_indicators];
4389  // run through values up to next constraint
4390  for (unsigned int j = 0; j < indicator.first; ++j)
4391  for (unsigned int comp = 0; comp < n_components; ++comp)
4392  operation.process_dof(dof_indices[v][j],
4393  *src[comp],
4394  values_dofs[comp][ind_local + j][v]);
4395 
4396  ind_local += indicator.first;
4397  dof_indices[v] += indicator.first;
4398 
4399  // constrained case: build the local value as a linear
4400  // combination of the global value according to constraints
4401  Number value[n_components];
4402  for (unsigned int comp = 0; comp < n_components; ++comp)
4403  operation.pre_constraints(values_dofs[comp][ind_local][v],
4404  value[comp]);
4405 
4406  const Number *data_val =
4407  this->matrix_info->constraint_pool_begin(indicator.second);
4408  const Number *end_pool =
4409  this->matrix_info->constraint_pool_end(indicator.second);
4410  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4411  for (unsigned int comp = 0; comp < n_components; ++comp)
4412  operation.process_constraint(*dof_indices[v],
4413  *data_val,
4414  *src[comp],
4415  value[comp]);
4416 
4417  for (unsigned int comp = 0; comp < n_components; ++comp)
4418  operation.post_constraints(value[comp],
4419  values_dofs[comp][ind_local][v]);
4420  ind_local++;
4421  }
4422 
4423  AssertIndexRange(ind_local, dofs_per_component + 1);
4424 
4425  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4426  for (unsigned int comp = 0; comp < n_components; ++comp)
4427  operation.process_dof(*dof_indices[v],
4428  *src[comp],
4429  values_dofs[comp][ind_local][v]);
4430  }
4431  else
4432  {
4433  // case with vector-valued finite elements where all components are
4434  // included in one single vector. Assumption: first come all entries
4435  // to the first component, then all entries to the second one, and
4436  // so on. This is ensured by the way MatrixFree reads out the
4437  // indices.
4438  for (unsigned int comp = 0; comp < n_components; ++comp)
4439  {
4440  unsigned int ind_local = 0;
4441 
4442  // check whether there is any constraint on the current cell
4443  for (; index_indicators != next_index_indicators;
4444  ++index_indicators)
4445  {
4446  const std::pair<unsigned short, unsigned short> indicator =
4447  this->dof_info->constraint_indicator[index_indicators];
4448 
4449  // run through values up to next constraint
4450  for (unsigned int j = 0; j < indicator.first; ++j)
4451  operation.process_dof(dof_indices[v][j],
4452  *src[0],
4453  values_dofs[comp][ind_local + j][v]);
4454  ind_local += indicator.first;
4455  dof_indices[v] += indicator.first;
4456 
4457  // constrained case: build the local value as a linear
4458  // combination of the global value according to constraints
4459  Number value;
4460  operation.pre_constraints(values_dofs[comp][ind_local][v],
4461  value);
4462 
4463  const Number *data_val =
4464  this->matrix_info->constraint_pool_begin(indicator.second);
4465  const Number *end_pool =
4466  this->matrix_info->constraint_pool_end(indicator.second);
4467 
4468  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4469  operation.process_constraint(*dof_indices[v],
4470  *data_val,
4471  *src[0],
4472  value);
4473 
4474  operation.post_constraints(value,
4475  values_dofs[comp][ind_local][v]);
4476  ind_local++;
4477  }
4478 
4479  AssertIndexRange(ind_local, dofs_per_component + 1);
4480 
4481  // get the dof values past the last constraint
4482  for (; ind_local < dofs_per_component;
4483  ++dof_indices[v], ++ind_local)
4484  {
4485  AssertIndexRange(*dof_indices[v], src[0]->size());
4486  operation.process_dof(*dof_indices[v],
4487  *src[0],
4488  values_dofs[comp][ind_local][v]);
4489  }
4490 
4491  if (apply_constraints == true && comp + 1 < n_components)
4492  next_index_indicators =
4493  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
4494  }
4495  }
4496  }
4497 }
4498 
4499 
4500 
4501 template <int dim,
4502  int n_components_,
4503  typename Number,
4504  bool is_face,
4505  typename VectorizedArrayType>
4506 template <typename VectorType, typename VectorOperation>
4507 inline void
4510  VectorType * src[]) const
4511 {
4512  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4513 
4514  unsigned int index =
4515  first_selected_component * this->data->dofs_per_component_on_cell;
4516  for (unsigned int comp = 0; comp < n_components; ++comp)
4517  {
4518  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4519  ++i, ++index)
4520  {
4521  operation.process_empty(values_dofs[comp][i]);
4522  operation.process_dof_global(
4523  local_dof_indices[this->data->lexicographic_numbering[index]],
4524  *src[0],
4525  values_dofs[comp][i][0]);
4526  }
4527  }
4528 }
4529 
4530 
4531 
4532 template <int dim,
4533  int n_components_,
4534  typename Number,
4535  bool is_face,
4536  typename VectorizedArrayType>
4537 template <typename VectorType, typename VectorOperation>
4538 inline void
4541  const VectorOperation & operation,
4542  VectorType * src[],
4543  const std::bitset<VectorizedArrayType::size()> &mask) const
4544 {
4545  // This functions processes the functions read_dof_values,
4546  // distribute_local_to_global, and set_dof_values with the same code for
4547  // contiguous cell indices (DG case). The distinction between these three
4548  // cases is made by the input VectorOperation that either reads values from
4549  // a vector and puts the data into the local data field or write local data
4550  // into the vector. Certain operations are no-ops for the given use case.
4551 
4552  std::integral_constant<bool,
4554  vector_selector;
4556  is_face ? this->dof_access_index :
4558  const unsigned int n_lanes = mask.count();
4559 
4560  const std::vector<unsigned int> &dof_indices_cont =
4561  this->dof_info->dof_indices_contiguous[ind];
4562 
4563  // Simple case: We have contiguous storage, so we can simply copy out the
4564  // data
4565  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4567  interleaved_contiguous &&
4568  n_lanes == VectorizedArrayType::size())
4569  {
4570  const unsigned int dof_index =
4571  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
4572  this->dof_info->component_dof_indices_offset[this->active_fe_index]
4573  [first_selected_component] *
4574  VectorizedArrayType::size();
4575  if (n_components == 1 || n_fe_components == 1)
4576  for (unsigned int comp = 0; comp < n_components; ++comp)
4577  operation.process_dofs_vectorized(
4578  this->data->dofs_per_component_on_cell,
4579  dof_index,
4580  *src[comp],
4581  values_dofs[comp],
4582  vector_selector);
4583  else
4584  operation.process_dofs_vectorized(
4585  this->data->dofs_per_component_on_cell * n_components,
4586  dof_index,
4587  *src[0],
4588  values_dofs[0],
4589  vector_selector);
4590  return;
4591  }
4592 
4593  // More general case: Must go through the components one by one and apply
4594  // some transformations
4595  const unsigned int n_filled_lanes =
4596  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
4597 
4598  unsigned int dof_indices[VectorizedArrayType::size()];
4599  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4600  dof_indices[v] =
4601  dof_indices_cont[this->cell * VectorizedArrayType::size() + v] +
4602  this->dof_info
4603  ->component_dof_indices_offset[this->active_fe_index]
4604  [this->first_selected_component] *
4605  this->dof_info->dof_indices_interleave_strides
4606  [ind][this->cell * VectorizedArrayType::size() + v];
4607 
4608  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4609  dof_indices[v] = numbers::invalid_unsigned_int;
4610 
4611  // In the case with contiguous cell indices, we know that there are no
4612  // constraints and that the indices within each element are contiguous
4613  if (n_filled_lanes == VectorizedArrayType::size() &&
4614  n_lanes == VectorizedArrayType::size())
4615  {
4616  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4618  contiguous)
4619  {
4620  if (n_components == 1 || n_fe_components == 1)
4621  for (unsigned int comp = 0; comp < n_components; ++comp)
4622  operation.process_dofs_vectorized_transpose(
4623  this->data->dofs_per_component_on_cell,
4624  dof_indices,
4625  *src[comp],
4626  values_dofs[comp],
4627  vector_selector);
4628  else
4629  operation.process_dofs_vectorized_transpose(
4630  this->data->dofs_per_component_on_cell * n_components,
4631  dof_indices,
4632  *src[0],
4633  &values_dofs[0][0],
4634  vector_selector);
4635  }
4636  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
4638  interleaved_contiguous_strided)
4639  {
4640  if (n_components == 1 || n_fe_components == 1)
4641  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4642  ++i)
4643  {
4644  for (unsigned int comp = 0; comp < n_components; ++comp)
4645  operation.process_dof_gather(dof_indices,
4646  *src[comp],
4647  i * VectorizedArrayType::size(),
4648  values_dofs[comp][i],
4649  vector_selector);
4650  }
4651  else
4652  for (unsigned int comp = 0; comp < n_components; ++comp)
4653  for (unsigned int i = 0;
4654  i < this->data->dofs_per_component_on_cell;
4655  ++i)
4656  {
4657  operation.process_dof_gather(
4658  dof_indices,
4659  *src[0],
4660  (comp * this->data->dofs_per_component_on_cell + i) *
4661  VectorizedArrayType::size(),
4662  values_dofs[comp][i],
4663  vector_selector);
4664  }
4665  }
4666  else
4667  {
4668  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
4670  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4671  ExcNotImplemented());
4672  const unsigned int *offsets =
4673  &this->dof_info->dof_indices_interleave_strides
4674  [ind][VectorizedArrayType::size() * this->cell];
4675  if (n_components == 1 || n_fe_components == 1)
4676  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4677  ++i)
4678  {
4679  for (unsigned int comp = 0; comp < n_components; ++comp)
4680  operation.process_dof_gather(dof_indices,
4681  *src[comp],
4682  0,
4683  values_dofs[comp][i],
4684  vector_selector);
4686  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4687  dof_indices[v] += offsets[v];
4688  }
4689  else
4690  for (unsigned int comp = 0; comp < n_components; ++comp)
4691  for (unsigned int i = 0;
4692  i < this->data->dofs_per_component_on_cell;
4693  ++i)
4694  {
4695  operation.process_dof_gather(dof_indices,
4696  *src[0],
4697  0,
4698  values_dofs[comp][i],
4699  vector_selector);
4701  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4702  dof_indices[v] += offsets[v];
4703  }
4704  }
4705  }
4706  else
4707  for (unsigned int comp = 0; comp < n_components; ++comp)
4708  {
4709  for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
4710  ++i)
4711  operation.process_empty(values_dofs[comp][i]);
4712  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4714  contiguous)
4715  {
4716  if (n_components == 1 || n_fe_components == 1)
4717  {
4718  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4719  if (mask[v] == true)
4720  for (unsigned int i = 0;
4721  i < this->data->dofs_per_component_on_cell;
4722  ++i)
4723  operation.process_dof(dof_indices[v] + i,
4724  *src[comp],
4725  values_dofs[comp][i][v]);
4726  }
4727  else
4728  {
4729  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4730  if (mask[v] == true)
4731  for (unsigned int i = 0;
4732  i < this->data->dofs_per_component_on_cell;
4733  ++i)
4734  operation.process_dof(
4735  dof_indices[v] + i +
4736  comp * this->data->dofs_per_component_on_cell,
4737  *src[0],
4738  values_dofs[comp][i][v]);
4739  }
4740  }
4741  else
4742  {
4743  const unsigned int *offsets =
4744  &this->dof_info->dof_indices_interleave_strides
4745  [ind][VectorizedArrayType::size() * this->cell];
4746  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4747  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4748  if (n_components == 1 || n_fe_components == 1)
4749  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4750  {
4751  if (mask[v] == true)
4752  for (unsigned int i = 0;
4753  i < this->data->dofs_per_component_on_cell;
4754  ++i)
4755  operation.process_dof(dof_indices[v] + i * offsets[v],
4756  *src[comp],
4757  values_dofs[comp][i][v]);
4758  }
4759  else
4760  {
4761  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4762  if (mask[v] == true)
4763  for (unsigned int i = 0;
4764  i < this->data->dofs_per_component_on_cell;
4765  ++i)
4766  operation.process_dof(
4767  dof_indices[v] +
4768  (i + comp * this->data->dofs_per_component_on_cell) *
4769  offsets[v],
4770  *src[0],
4771  values_dofs[comp][i][v]);
4772  }
4773  }
4774  }
4775 }
4776 
4777 
4778 
4779 template <int dim,
4780  int n_components_,
4781  typename Number,
4782  bool is_face,
4783  typename VectorizedArrayType>
4784 template <typename VectorType>
4785 inline void
4787  read_dof_values(const VectorType &src, const unsigned int first_index)
4788 {
4789  // select between block vectors and non-block vectors. Note that the number
4790  // of components is checked in the internal data
4791  typename internal::BlockVectorSelector<
4792  VectorType,
4793  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4794  for (unsigned int d = 0; d < n_components; ++d)
4795  src_data[d] =
4796  internal::BlockVectorSelector<VectorType,
4798  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4799 
4801  read_write_operation(reader,
4802  src_data,
4803  std::bitset<VectorizedArrayType::size()>().flip(),
4804  true);
4805 
4806 # ifdef DEBUG
4807  dof_values_initialized = true;
4808 # endif
4809 }
4810 
4811 
4812 
4813 template <int dim,
4814  int n_components_,
4815  typename Number,
4816  bool is_face,
4817  typename VectorizedArrayType>
4818 template <typename VectorType>
4819 inline void
4821  read_dof_values_plain(const VectorType &src, const unsigned int first_index)
4822 {
4823  // select between block vectors and non-block vectors. Note that the number
4824  // of components is checked in the internal data
4825  typename internal::BlockVectorSelector<
4826  VectorType,
4827  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4828  for (unsigned int d = 0; d < n_components; ++d)
4829  src_data[d] =
4830  internal::BlockVectorSelector<VectorType,
4832  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4833 
4835  read_write_operation(reader,
4836  src_data,
4837  std::bitset<VectorizedArrayType::size()>().flip(),
4838  false);
4839 
4840 # ifdef DEBUG
4841  dof_values_initialized = true;
4842 # endif
4843 }
4844 
4845 
4846 
4847 template <int dim,
4848  int n_components_,
4849  typename Number,
4850  bool is_face,
4851  typename VectorizedArrayType>
4852 template <typename VectorType>
4853 inline void
4856  VectorType & dst,
4857  const unsigned int first_index,
4858  const std::bitset<VectorizedArrayType::size()> &mask) const
4859 {
4860 # ifdef DEBUG
4861  Assert(dof_values_initialized == true,
4863 # endif
4864 
4865  // select between block vectors and non-block vectors. Note that the number
4866  // of components is checked in the internal data
4867  typename internal::BlockVectorSelector<
4868  VectorType,
4869  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4870  for (unsigned int d = 0; d < n_components; ++d)
4871  dst_data[d] = internal::BlockVectorSelector<
4872  VectorType,
4873  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4874  d + first_index);
4875 
4877  distributor;
4878  read_write_operation(distributor, dst_data, mask);
4879 }
4880 
4881 
4882 
4883 template <int dim,
4884  int n_components_,
4885  typename Number,
4886  bool is_face,
4887  typename VectorizedArrayType>
4888 template <typename VectorType>
4889 inline void
4892  const unsigned int first_index,
4893  const std::bitset<VectorizedArrayType::size()> &mask) const
4894 {
4895 # ifdef DEBUG
4896  Assert(dof_values_initialized == true,
4898 # endif
4899 
4900  // select between block vectors and non-block vectors. Note that the number
4901  // of components is checked in the internal data
4902  typename internal::BlockVectorSelector<
4903  VectorType,
4904  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4905  for (unsigned int d = 0; d < n_components; ++d)
4906  dst_data[d] = internal::BlockVectorSelector<
4907  VectorType,
4908  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4909  d + first_index);
4910 
4912  read_write_operation(setter, dst_data, mask);
4913 }
4914 
4915 
4916 
4917 template <int dim,
4918  int n_components_,
4919  typename Number,
4920  bool is_face,
4921  typename VectorizedArrayType>
4922 template <typename VectorType>
4923 inline void
4926  VectorType & dst,
4927  const unsigned int first_index,
4928  const std::bitset<VectorizedArrayType::size()> &mask) const
4929 {
4930 # ifdef DEBUG
4931  Assert(dof_values_initialized == true,
4933 # endif
4934 
4935  // select between block vectors and non-block vectors. Note that the number
4936  // of components is checked in the internal data
4937  typename internal::BlockVectorSelector<
4938  VectorType,
4939  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4940  for (unsigned int d = 0; d < n_components; ++d)
4941  dst_data[d] = internal::BlockVectorSelector<
4942  VectorType,
4943  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4944  d + first_index);
4945 
4947  read_write_operation(setter, dst_data, mask, false);
4948 }
4949 
4950 
4951 
4952 /*------------------------------ access to data fields ----------------------*/
4953 
4954 
4955 
4956 template <int dim,
4957  int n_components,
4958  typename Number,
4959  bool is_face,
4960  typename VectorizedArrayType>
4961 inline const VectorizedArrayType *
4963  begin_dof_values() const
4964 {
4965  return &values_dofs[0][0];
4966 }
4967 
4968 
4969 
4970 template <int dim,
4971  int n_components,
4972  typename Number,
4973  bool is_face,
4974  typename VectorizedArrayType>
4975 inline VectorizedArrayType *
4978 {
4979 # ifdef DEBUG
4980  dof_values_initialized = true;
4981 # endif
4982  return &values_dofs[0][0];
4983 }
4984 
4985 
4986 
4987 template <int dim,
4988  int n_components,
4989  typename Number,
4990  bool is_face,
4991  typename VectorizedArrayType>
4992 inline const VectorizedArrayType *
4994  begin_values() const
4995 {
4996 # ifdef DEBUG
4997  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
4998 # endif
4999  return values_quad;
5000 }
5001 
5002 
5003 
5004 template <int dim,
5005  int n_components,
5006  typename Number,
5007  bool is_face,
5008  typename VectorizedArrayType>
5009 inline VectorizedArrayType *
5011  begin_values()
5012 {
5013 # ifdef DEBUG
5014  values_quad_initialized = true;
5015  values_quad_submitted = true;
5016 # endif
5017  return values_quad;
5018 }
5019 
5020 
5021 
5022 template <int dim,
5023  int n_components,
5024  typename Number,
5025  bool is_face,
5026  typename VectorizedArrayType>
5027 inline const VectorizedArrayType *
5029  begin_gradients() const
5030 {
5031 # ifdef DEBUG
5032  Assert(gradients_quad_initialized || gradients_quad_submitted,
5033  ExcNotInitialized());
5034 # endif
5035  return gradients_quad;
5036 }
5037 
5038 
5039 
5040 template <int dim,
5041  int n_components,
5042  typename Number,
5043  bool is_face,
5044  typename VectorizedArrayType>
5045 inline VectorizedArrayType *
5048 {
5049 # ifdef DEBUG
5050  gradients_quad_submitted = true;
5051  gradients_quad_initialized = true;
5052 # endif
5053  return gradients_quad;
5054 }
5055 
5056 
5057 
5058 template <int dim,
5059  int n_components,
5060  typename Number,
5061  bool is_face,
5062  typename VectorizedArrayType>
5063 inline const VectorizedArrayType *
5065  begin_hessians() const
5066 {
5067 # ifdef DEBUG
5068  Assert(hessians_quad_initialized, ExcNotInitialized());
5069 # endif
5070  return hessians_quad;
5071 }
5072 
5073 
5074 
5075 template <int dim,
5076  int n_components,
5077  typename Number,
5078  bool is_face,
5079  typename VectorizedArrayType>
5080 inline VectorizedArrayType *
5083 {
5084 # ifdef DEBUG
5085  hessians_quad_initialized = true;
5086 # endif
5087  return hessians_quad;
5088 }
5089 
5090 
5091 
5092 template <int dim,
5093  int n_components_,
5094  typename Number,
5095  bool is_face,
5096  typename VectorizedArrayType>
5099  get_dof_value(const unsigned int dof) const
5100 {
5101  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5103  for (unsigned int comp = 0; comp < n_components; comp++)
5104  return_value[comp] = this->values_dofs[comp][dof];
5105  return return_value;
5106 }
5107 
5108 
5109 
5110 template <int dim,
5111  int n_components_,
5112  typename Number,
5113  bool is_face,
5114  typename VectorizedArrayType>
5117  get_value(const unsigned int q_point) const
5118 {
5119 # ifdef DEBUG
5120  Assert(this->values_quad_initialized == true,
5122 # endif
5123 
5124  AssertIndexRange(q_point, this->n_quadrature_points);
5125  const std::size_t nqp = this->n_quadrature_points;
5127  for (unsigned int comp = 0; comp < n_components; comp++)
5128  return_value[comp] = values_quad[comp * nqp + q_point];
5129  return return_value;
5130 }
5131 
5132 
5133 
5134 template <int dim,
5135  int n_components_,
5136  typename Number,
5137  bool is_face,
5138  typename VectorizedArrayType>
5139 inline DEAL_II_ALWAYS_INLINE
5142  get_gradient(const unsigned int q_point) const
5143 {
5144 # ifdef DEBUG
5145  Assert(this->gradients_quad_initialized == true,
5147 # endif
5148 
5149  AssertIndexRange(q_point, this->n_quadrature_points);
5150  Assert(this->jacobian != nullptr, ExcNotInitialized());
5151  const std::size_t nqp = this->n_quadrature_points;
5153 
5154  // Cartesian cell
5155  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5156  {
5157  for (unsigned int d = 0; d < dim; ++d)
5158  for (unsigned int comp = 0; comp < n_components; comp++)
5159  grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] *
5160  this->jacobian[0][d][d];
5161  }
5162  // cell with general/affine Jacobian
5163  else
5164  {
5166  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5167  q_point :
5168  0];
5169  for (unsigned int comp = 0; comp < n_components; comp++)
5170  for (unsigned int d = 0; d < dim; ++d)
5171  {
5172  grad_out[comp][d] =
5173  jac[d][0] * gradients_quad[(comp * dim) * nqp + q_point];
5174  for (unsigned int e = 1; e < dim; ++e)
5175  grad_out[comp][d] +=
5176  jac[d][e] * gradients_quad[(comp * dim + e) * nqp + q_point];
5177  }
5178  }
5179  return grad_out;
5180 }
5181 
5182 
5183 
5184 template <int dim,
5185  int n_components_,
5186  typename Number,
5187  bool is_face,
5188  typename VectorizedArrayType>
5191  get_normal_derivative(const unsigned int q_point) const
5192 {
5193  AssertIndexRange(q_point, this->n_quadrature_points);
5194 # ifdef DEBUG
5195  Assert(this->gradients_quad_initialized == true,
5197 # endif
5198 
5199  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
5200 
5201  const std::size_t nqp = this->n_quadrature_points;
5203 
5204  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5205  for (unsigned int comp = 0; comp < n_components; comp++)
5206  grad_out[comp] = gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
5207  (this->normal_x_jacobian[0][dim - 1]);
5208  else
5209  {
5210  const std::size_t index =
5211  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5212  for (unsigned int comp = 0; comp < n_components; comp++)
5213  {
5214  grad_out[comp] = gradients_quad[comp * dim * nqp + q_point] *
5215  this->normal_x_jacobian[index][0];
5216  for (unsigned int d = 1; d < dim; ++d)
5217  grad_out[comp] += gradients_quad[(comp * dim + d) * nqp + q_point] *
5218  this->normal_x_jacobian[index][d];
5219  }
5220  }
5221  return grad_out;
5222 }
5223 
5224 
5225 
5226 namespace internal
5227 {
5228  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5229  // store the lower diagonal because of symmetry
5230  template <typename VectorizedArrayType>
5231  inline void
5232  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
5233  const VectorizedArrayType *const hessians,
5234  const unsigned int,
5235  VectorizedArrayType (&tmp)[1][1])
5236  {
5237  tmp[0][0] = jac[0][0] * hessians[0];
5238  }
5239 
5240  template <typename VectorizedArrayType>
5241  inline void
5242  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
5243  const VectorizedArrayType *const hessians,
5244  const unsigned int nqp,
5245  VectorizedArrayType (&tmp)[2][2])
5246  {
5247  for (unsigned int d = 0; d < 2; ++d)
5248  {
5249  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
5250  tmp[1][d] =
5251  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
5252  }
5253  }
5254 
5255  template <typename VectorizedArrayType>
5256  inline void
5257  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
5258  const VectorizedArrayType *const hessians,
5259  const unsigned int nqp,
5260  VectorizedArrayType (&tmp)[3][3])
5261  {
5262  for (unsigned int d = 0; d < 3; ++d)
5263  {
5264  tmp[0][d] =
5265  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
5266  jac[d][2] * hessians[4 * nqp]);
5267  tmp[1][d] =
5268  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
5269  jac[d][2] * hessians[5 * nqp]);
5270  tmp[2][d] =
5271  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
5272  jac[d][2] * hessians[2 * nqp]);
5273  }
5274  }
5275 } // namespace internal
5276 
5277 
5278 
5279 template <int dim,
5280  int n_components_,
5281  typename Number,
5282  bool is_face,
5283  typename VectorizedArrayType>
5286  get_hessian(const unsigned int q_point) const
5287 {
5288  Assert(!is_face, ExcNotImplemented());
5289 # ifdef DEBUG
5290  Assert(this->hessians_quad_initialized == true,
5292 # endif
5293  AssertIndexRange(q_point, this->n_quadrature_points);
5294 
5295  Assert(this->jacobian != nullptr, ExcNotImplemented());
5297  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5298  0 :
5299  q_point];
5300 
5302 
5303  const std::size_t nqp = this->n_quadrature_points;
5304  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5305 
5306  // Cartesian cell
5307  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5308  {
5309  for (unsigned int comp = 0; comp < n_components; comp++)
5310  {
5311  for (unsigned int d = 0; d < dim; ++d)
5312  hessian_out[comp][d][d] =
5313  hessians_quad[(comp * hdim + d) * nqp + q_point] *
5314  (jac[d][d] * jac[d][d]);
5315  switch (dim)
5316  {
5317  case 1:
5318  break;
5319  case 2:
5320  hessian_out[comp][0][1] =
5321  hessians_quad[(comp * hdim + 2) * nqp + q_point] *
5322  (jac[0][0] * jac[1][1]);
5323  break;
5324  case 3:
5325  hessian_out[comp][0][1] =
5326  hessians_quad[(comp * hdim + 3) * nqp + q_point] *
5327  (jac[0][0] * jac[1][1]);
5328  hessian_out[comp][0][2] =
5329  hessians_quad[(comp * hdim + 4) * nqp + q_point] *
5330  (jac[0][0] * jac[2][2]);
5331  hessian_out[comp][1][2] =
5332  hessians_quad[(comp * hdim + 5) * nqp + q_point] *
5333  (jac[1][1] * jac[2][2]);
5334  break;
5335  default:
5336  Assert(false, ExcNotImplemented());
5337  }
5338  for (unsigned int d = 0; d < dim; ++d)
5339  for (unsigned int e = d + 1; e < dim; ++e)
5340  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5341  }
5342  }
5343  // cell with general Jacobian, but constant within the cell
5344  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5345  {
5346  for (unsigned int comp = 0; comp < n_components; comp++)
5347  {
5348  VectorizedArrayType tmp[dim][dim];
5349  internal::hessian_unit_times_jac(
5350  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5351 
5352  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5353  for (unsigned int d = 0; d < dim; ++d)
5354  for (unsigned int e = d; e < dim; ++e)
5355  {
5356  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5357  for (unsigned int f = 1; f < dim; ++f)
5358  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5359  }
5360 
5361  // no J' * grad(u) part here because the Jacobian is constant
5362  // throughout the cell and hence, its derivative is zero
5363 
5364  // take symmetric part
5365  for (unsigned int d = 0; d < dim; ++d)
5366  for (unsigned int e = d + 1; e < dim; ++e)
5367  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5368  }
5369  }
5370  // cell with general Jacobian
5371  else
5372  {
5373  const auto &jac_grad =
5374  this->mapping_data->jacobian_gradients
5375  [1 - this->is_interior_face]
5376  [this->mapping_data->data_index_offsets[this->cell] + q_point];
5377  for (unsigned int comp = 0; comp < n_components; comp++)
5378  {
5379  // compute laplacian before the gradient because it needs to access
5380  // unscaled gradient data
5381  VectorizedArrayType tmp[dim][dim];
5382  internal::hessian_unit_times_jac(
5383  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5384 
5385  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5386  for (unsigned int d = 0; d < dim; ++d)
5387  for (unsigned int e = d; e < dim; ++e)
5388  {
5389  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5390  for (unsigned int f = 1; f < dim; ++f)
5391  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5392  }
5393 
5394  // add diagonal part of J' * grad(u)
5395  for (unsigned int d = 0; d < dim; ++d)
5396  for (unsigned int e = 0; e < dim; ++e)
5397  hessian_out[comp][d][d] +=
5398  jac_grad[d][e] *
5399  gradients_quad[(comp * dim + e) * nqp + q_point];
5400 
5401  // add off-diagonal part of J' * grad(u)
5402  for (unsigned int d = 0, count = dim; d < dim; ++d)
5403  for (unsigned int e = d + 1; e < dim; ++e, ++count)
5404  for (unsigned int f = 0; f < dim; ++f)
5405  hessian_out[comp][d][e] +=
5406  jac_grad[count][f] *
5407  gradients_quad[(comp * dim + f) * nqp + q_point];
5408 
5409  // take symmetric part
5410  for (unsigned int d = 0; d < dim; ++d)
5411  for (unsigned int e = d + 1; e < dim; ++e)
5412  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5413  }
5414  }
5415  return hessian_out;
5416 }
5417 
5418 
5419 
5420 template <int dim,
5421  int n_components_,
5422  typename Number,
5423  bool is_face,
5424  typename VectorizedArrayType>
5427  get_hessian_diagonal(const unsigned int q_point) const
5428 {
5429  Assert(!is_face, ExcNotImplemented());
5430 # ifdef DEBUG
5431  Assert(this->hessians_quad_initialized == true,
5433 # endif
5434  AssertIndexRange(q_point, this->n_quadrature_points);
5435 
5436  Assert(this->jacobian != nullptr, ExcNotImplemented());
5438  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5439  0 :
5440  q_point];
5441 
5442  const std::size_t nqp = this->n_quadrature_points;
5443  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5445 
5446  // Cartesian cell
5447  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5448  {
5449  for (unsigned int comp = 0; comp < n_components; comp++)
5450  for (unsigned int d = 0; d < dim; ++d)
5451  hessian_out[comp][d] =
5452  hessians_quad[(comp * hdim + d) * nqp + q_point] *
5453  (jac[d][d] * jac[d][d]);
5454  }
5455  // cell with general Jacobian, but constant within the cell
5456  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5457  {
5458  for (unsigned int comp = 0; comp < n_components; comp++)
5459  {
5460  // compute laplacian before the gradient because it needs to access
5461  // unscaled gradient data
5462  VectorizedArrayType tmp[dim][dim];
5463  internal::hessian_unit_times_jac(
5464  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5465 
5466  // compute only the trace part of hessian, J * tmp = J *
5467  // hess_unit(u) * J^T
5468  for (unsigned int d = 0; d < dim; ++d)
5469  {
5470  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5471  for (unsigned int f = 1; f < dim; ++f)
5472  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5473  }
5474  }
5475  }
5476  // cell with general Jacobian
5477  else
5478  {
5479  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5480  &jac_grad =
5481  this->mapping_data->jacobian_gradients
5482  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5483  for (unsigned int comp = 0; comp < n_components; comp++)
5484  {
5485  // compute laplacian before the gradient because it needs to access
5486  // unscaled gradient data
5487  VectorizedArrayType tmp[dim][dim];
5488  internal::hessian_unit_times_jac(
5489  jac, hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
5490 
5491  // compute only the trace part of hessian, J * tmp = J *
5492  // hess_unit(u) * J^T
5493  for (unsigned int d = 0; d < dim; ++d)
5494  {
5495  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5496  for (unsigned int f = 1; f < dim; ++f)
5497  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5498  }
5499 
5500  for (unsigned int d = 0; d < dim; ++d)
5501  for (unsigned int e = 0; e < dim; ++e)
5502  hessian_out[comp][d] +=
5503  jac_grad[d][e] *
5504  gradients_quad[(comp * dim + e) * nqp + q_point];
5505  }
5506  }
5507  return hessian_out;
5508 }
5509 
5510 
5511 
5512 template <int dim,
5513  int n_components_,
5514  typename Number,
5515  bool is_face,
5516  typename VectorizedArrayType>
5519  get_laplacian(const unsigned int q_point) const
5520 {
5521  Assert(is_face == false, ExcNotImplemented());
5522 # ifdef DEBUG
5523  Assert(this->hessians_quad_initialized == true,
5525 # endif
5526  AssertIndexRange(q_point, this->n_quadrature_points);
5527 
5529  const auto hess_diag = get_hessian_diagonal(q_point);
5530  for (unsigned int comp = 0; comp < n_components; ++comp)
5531  {
5532  laplacian_out[comp] = hess_diag[comp][0];
5533  for (unsigned int d = 1; d < dim; ++d)
5534  laplacian_out[comp] += hess_diag[comp][d];
5535  }
5536  return laplacian_out;
5537 }
5538 
5539 
5540 
5541 template <int dim,
5542  int n_components_,
5543  typename Number,
5544  bool is_face,
5545  typename VectorizedArrayType>
5546 inline DEAL_II_ALWAYS_INLINE void
5549  const unsigned int dof)
5550 {
5551 # ifdef DEBUG
5552  this->dof_values_initialized = true;
5553 # endif
5554  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5555  for (unsigned int comp = 0; comp < n_components; comp++)
5556  this->values_dofs[comp][dof] = val_in[comp];
5557 }
5558 
5559 
5560 
5561 template <int dim,
5562  int n_components_,
5563  typename Number,
5564  bool is_face,
5565  typename VectorizedArrayType>
5566 inline DEAL_II_ALWAYS_INLINE void
5569  const unsigned int q_point)
5570 {
5572  AssertIndexRange(q_point, this->n_quadrature_points);
5573  Assert(this->J_value != nullptr, ExcNotInitialized());
5574 # ifdef DEBUG
5575  this->values_quad_submitted = true;
5576 # endif
5577 
5578  const std::size_t nqp = this->n_quadrature_points;
5579  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5580  {
5581  const VectorizedArrayType JxW =
5582  this->J_value[0] * this->quadrature_weights[q_point];
5583  for (unsigned int comp = 0; comp < n_components; ++comp)
5584  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5585  }
5586  else
5587  {
5588  const VectorizedArrayType JxW = this->J_value[q_point];
5589  for (unsigned int comp = 0; comp < n_components; ++comp)
5590  values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5591  }
5592 }
5593 
5594 
5595 
5596 template <int dim,
5597  int n_components_,
5598  typename Number,
5599  bool is_face,
5600  typename VectorizedArrayType>
5601 inline DEAL_II_ALWAYS_INLINE void
5604  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5605  const unsigned int q_point)
5606 {
5608  AssertIndexRange(q_point, this->n_quadrature_points);
5609  Assert(this->J_value != nullptr, ExcNotInitialized());
5610  Assert(this->jacobian != nullptr, ExcNotInitialized());
5611 # ifdef DEBUG
5612  this->gradients_quad_submitted = true;
5613 # endif
5614 
5615  const std::size_t nqp = this->n_quadrature_points;
5616  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5617  {
5618  const VectorizedArrayType JxW =
5619  this->J_value[0] * this->quadrature_weights[q_point];
5620  for (unsigned int d = 0; d < dim; ++d)
5621  {
5622  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
5623  for (unsigned int comp = 0; comp < n_components; comp++)
5624  gradients_quad[(comp * dim + d) * nqp + q_point] =
5625  grad_in[comp][d] * factor;
5626  }
5627  }
5628  else
5629  {
5631  this->cell_type > internal::MatrixFreeFunctions::affine ?
5632  this->jacobian[q_point] :
5633  this->jacobian[0];
5634  const VectorizedArrayType JxW =
5635  this->cell_type > internal::MatrixFreeFunctions::affine ?
5636  this->J_value[q_point] :
5637  this->J_value[0] * this->quadrature_weights[q_point];
5638  for (unsigned int comp = 0; comp < n_components; ++comp)
5639  for (unsigned int d = 0; d < dim; ++d)
5640  {
5641  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5642  for (unsigned int e = 1; e < dim; ++e)
5643  new_val += (jac[e][d] * grad_in[comp][e]);
5644  gradients_quad[(comp * dim + d) * nqp + q_point] = new_val * JxW;
5645  }
5646  }
5647 }
5648 
5649 
5650 
5651 template <int dim,
5652  int n_components_,
5653  typename Number,
5654  bool is_face,
5655  typename VectorizedArrayType>
5656 inline DEAL_II_ALWAYS_INLINE void
5660  const unsigned int q_point)
5661 {
5662  AssertIndexRange(q_point, this->n_quadrature_points);
5663  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
5664 # ifdef DEBUG
5665  this->gradients_quad_submitted = true;
5666 # endif
5667 
5668  const std::size_t nqp = this->n_quadrature_points;
5669  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5670  for (unsigned int comp = 0; comp < n_components; comp++)
5671  {
5672  for (unsigned int d = 0; d < dim - 1; ++d)
5673  gradients_quad[(comp * dim + d) * nqp + q_point] =
5674  VectorizedArrayType();
5675  gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5676  grad_in[comp] *
5677  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5678  this->quadrature_weights[q_point]);
5679  }
5680  else
5681  {
5682  const unsigned int index =
5683  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5685  this->normal_x_jacobian[index];
5686  for (unsigned int comp = 0; comp < n_components; comp++)
5687  {
5688  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5689  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5690  factor = factor * this->quadrature_weights[q_point];
5691  for (unsigned int d = 0; d < dim; ++d)
5692  gradients_quad[(comp * dim + d) * nqp + q_point] = factor * jac[d];
5693  }
5694  }
5695 }
5696 
5697 
5698 
5699 template <int dim,
5700  int n_components_,
5701  typename Number,
5702  bool is_face,
5703  typename VectorizedArrayType>
5706  integrate_value() const
5707 {
5709 # ifdef DEBUG
5710  Assert(this->values_quad_submitted == true,
5712 # endif
5713 
5715  const std::size_t nqp = this->n_quadrature_points;
5716  for (unsigned int q = 0; q < nqp; ++q)
5717  for (unsigned int comp = 0; comp < n_components; ++comp)
5718  return_value[comp] += this->values_quad[comp * nqp + q];
5719  return (return_value);
5720 }
5721 
5722 
5723 
5724 /*----------------------- FEEvaluationAccess --------------------------------*/
5725 
5726 
5727 template <int dim,
5728  int n_components_,
5729  typename Number,
5730  bool is_face,
5731  typename VectorizedArrayType>
5732 inline FEEvaluationAccess<dim,
5733  n_components_,
5734  Number,
5735  is_face,
5736  VectorizedArrayType>::
5737  FEEvaluationAccess(
5739  const unsigned int dof_no,
5740  const unsigned int first_selected_component,
5741  const unsigned int quad_no_in,
5742  const unsigned int fe_degree,
5743  const unsigned int n_q_points,
5744  const bool is_interior_face)
5746  data_in,
5747  dof_no,
5748  first_selected_component,
5749  quad_no_in,
5750  fe_degree,
5751  n_q_points,
5752  is_interior_face)
5753 {}
5754 
5755 
5756 
5757 template <int dim,
5758  int n_components_,
5759  typename Number,
5760  bool is_face,
5761  typename VectorizedArrayType>
5762 inline FEEvaluationAccess<dim,
5763  n_components_,
5764  Number,
5765  is_face,
5766  VectorizedArrayType>::
5767  FEEvaluationAccess(
5768  const Mapping<dim> & mapping,
5769  const FiniteElement<dim> &fe,
5770  const Quadrature<1> & quadrature,
5771  const UpdateFlags update_flags,
5772  const unsigned int first_selected_component,
5774  *other)
5776  mapping,
5777  fe,
5778  quadrature,
5779  update_flags,
5780  first_selected_component,
5781  other)
5782 {}
5783 
5784 
5785 
5786 template <int dim,
5787  int n_components_,
5788  typename Number,
5789  bool is_face,
5790  typename VectorizedArrayType>
5791 inline FEEvaluationAccess<dim,
5792  n_components_,
5793  Number,
5794  is_face,
5795  VectorizedArrayType>::
5796  FEEvaluationAccess(const FEEvaluationAccess<dim,
5797  n_components_,
5798  Number,
5799  is_face,
5800  VectorizedArrayType> &other)
5802  other)
5803 {}
5804 
5805 
5806 
5807 template <int dim,
5808  int n_components_,
5809  typename Number,
5810  bool is_face,
5811  typename VectorizedArrayType>
5812 inline FEEvaluationAccess<dim,
5813  n_components_,
5814  Number,
5815  is_face,
5816  VectorizedArrayType> &
5818 operator=(const FEEvaluationAccess<dim,
5819  n_components_,
5820  Number,
5821  is_face,
5822  VectorizedArrayType> &other)
5823 {
5824  this->FEEvaluationBase<dim,
5825  n_components_,
5826  Number,
5827  is_face,
5828  VectorizedArrayType>::operator=(other);
5829  return *this;
5830 }
5831 
5832 
5833 
5834 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5835 
5836 
5837 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5841  const unsigned int dof_no,
5842  const unsigned int first_selected_component,
5843  const unsigned int quad_no_in,
5844  const unsigned int fe_degree,
5845  const unsigned int n_q_points,
5846  const bool is_interior_face)
5848  data_in,
5849  dof_no,
5850  first_selected_component,
5851  quad_no_in,
5852  fe_degree,
5853  n_q_points,
5854  is_interior_face)
5855 {}
5856 
5857 
5858 
5859 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5862  const Mapping<dim> & mapping,
5863  const FiniteElement<dim> &fe,
5864  const Quadrature<1> & quadrature,
5865  const UpdateFlags update_flags,
5866  const unsigned int first_selected_component,
5868  *other)
5870  mapping,
5871  fe,
5872  quadrature,
5873  update_flags,
5874  first_selected_component,
5875  other)
5876 {}
5877 
5878 
5879 
5880 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5884  &other)
5886 {}
5887 
5888 
5889 
5890 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5894 {
5896  operator=(other);
5897  return *this;
5898 }
5899 
5900 
5901 
5902 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5903 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5905  const unsigned int dof) const
5906 {
5907  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5908  return this->values_dofs[0][dof];
5909 }
5910 
5911 
5912 
5913 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5914 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5916  const unsigned int q_point) const
5917 {
5918 # ifdef DEBUG
5919  Assert(this->values_quad_initialized == true,
5921 # endif
5922  AssertIndexRange(q_point, this->n_quadrature_points);
5923  return this->values_quad[q_point];
5924 }
5925 
5926 
5927 
5928 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5929 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5931  get_normal_derivative(const unsigned int q_point) const
5932 {
5933  return BaseClass::get_normal_derivative(q_point)[0];
5934 }
5935 
5936 
5937 
5938 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5941  const unsigned int q_point) const
5942 {
5943  // could use the base class gradient, but that involves too many expensive
5944  // initialization operations on tensors
5945 
5946 # ifdef DEBUG
5947  Assert(this->gradients_quad_initialized == true,
5949 # endif
5950  AssertIndexRange(q_point, this->n_quadrature_points);
5951 
5952  Assert(this->jacobian != nullptr, ExcNotInitialized());
5953 
5955 
5956  const std::size_t nqp = this->n_quadrature_points;
5957  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5958  {
5959  for (unsigned int d = 0; d < dim; ++d)
5960  grad_out[d] =
5961  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5962  }
5963  // cell with general/affine Jacobian
5964  else
5965  {
5967  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5968  q_point :
5969  0];
5970  for (unsigned int d = 0; d < dim; ++d)
5971  {
5972  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5973  for (unsigned int e = 1; e < dim; ++e)
5974  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5975  }
5976  }
5977  return grad_out;
5978 }
5979 
5980 
5981 
5982 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5985  const unsigned int q_point) const
5986 {
5987  return BaseClass::get_hessian(q_point)[0];
5988 }
5989 
5990 
5991 
5992 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5995  get_hessian_diagonal(const unsigned int q_point) const
5996 {
5997  return BaseClass::get_hessian_diagonal(q_point)[0];
5998 }
5999 
6000 
6001 
6002 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6003 inline VectorizedArrayType
6005  const unsigned int q_point) const
6006 {
6007  return BaseClass::get_laplacian(q_point)[0];
6008 }
6009 
6010 
6011 
6012 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6013 inline void DEAL_II_ALWAYS_INLINE
6015  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6016 {
6017 # ifdef DEBUG
6018  this->dof_values_initialized = true;
6019  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6020 # endif
6021  this->values_dofs[0][dof] = val_in;
6022 }
6023 
6024 
6025 
6026 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6027 inline void DEAL_II_ALWAYS_INLINE
6029  const VectorizedArrayType val_in,
6030  const unsigned int q_point)
6031 {
6033  AssertIndexRange(q_point, this->n_quadrature_points);
6034  Assert(this->J_value != nullptr, ExcNotInitialized());
6035 # ifdef DEBUG
6036  this->values_quad_submitted = true;
6037 # endif
6038 
6039  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6040  {
6041  const VectorizedArrayType JxW =
6042  this->J_value[0] * this->quadrature_weights[q_point];
6043  this->values_quad[q_point] = val_in * JxW;
6044  }
6045  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
6046  {
6047  this->values_quad[q_point] = val_in * this->J_value[q_point];
6048  }
6049 }
6050 
6051 
6052 
6053 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6054 inline DEAL_II_ALWAYS_INLINE void
6056  const Tensor<1, 1, VectorizedArrayType> val_in,
6057  const unsigned int q_point)
6058 {
6059  submit_value(val_in[0], q_point);
6060 }
6061 
6062 
6063 
6064 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6065 inline DEAL_II_ALWAYS_INLINE void
6067  submit_normal_derivative(const VectorizedArrayType grad_in,
6068  const unsigned int q_point)
6069 {
6071  grad[0] = grad_in;
6072  BaseClass::submit_normal_derivative(grad, q_point);
6073 }
6074 
6075 
6076 
6077 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6078 inline DEAL_II_ALWAYS_INLINE void
6081  const unsigned int q_point)
6082 {
6084  AssertIndexRange(q_point, this->n_quadrature_points);
6085  Assert(this->J_value != nullptr, ExcNotInitialized());
6086  Assert(this->jacobian != nullptr, ExcNotInitialized());
6087 # ifdef DEBUG
6088  this->gradients_quad_submitted = true;
6089 # endif
6090 
6091  const std::size_t nqp = this->n_quadrature_points;
6092  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6093  {
6094  const VectorizedArrayType JxW =
6095  this->J_value[0] * this->quadrature_weights[q_point];
6096  for (unsigned int d = 0; d < dim; ++d)
6097  this->gradients_quad[d * nqp + q_point] =
6098  (grad_in[d] * this->jacobian[0][d][d] * JxW);
6099  }
6100  // general/affine cell type
6101  else
6102  {
6104  this->cell_type > internal::MatrixFreeFunctions::affine ?
6105  this->jacobian[q_point] :
6106  this->jacobian[0];
6107  const VectorizedArrayType JxW =
6108  this->cell_type > internal::MatrixFreeFunctions::affine ?
6109  this->J_value[q_point] :
6110  this->J_value[0] * this->quadrature_weights[q_point];
6111  for (unsigned int d = 0; d < dim; ++d)
6112  {
6113  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
6114  for (unsigned int e = 1; e < dim; ++e)
6115  new_val += jac[e][d] * grad_in[e];
6116  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
6117  }
6118  }
6119 }
6120 
6121 
6122 
6123 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6124 inline VectorizedArrayType
6126  integrate_value() const
6127 {
6128  return BaseClass::integrate_value()[0];
6129 }
6130 
6131 
6132 
6133 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
6134 
6135 
6136 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6140  const unsigned int dof_no,
6141  const unsigned int first_selected_component,
6142  const unsigned int quad_no_in,
6143  const unsigned int fe_degree,
6144  const unsigned int n_q_points,
6145  const bool is_interior_face)
6147  data_in,
6148  dof_no,
6149  first_selected_component,
6150  quad_no_in,
6151  fe_degree,
6152  n_q_points,
6153  is_interior_face)
6154 {}
6155 
6156 
6157 
6158 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6161  const Mapping<dim> & mapping,
6162  const FiniteElement<dim> &fe,
6163  const Quadrature<1> & quadrature,
6164  const UpdateFlags update_flags,
6165  const unsigned int first_selected_component,
6167  *other)
6169  mapping,
6170  fe,
6171  quadrature,
6172  update_flags,
6173  first_selected_component,
6174  other)
6175 {}
6176 
6177 
6178 
6179 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6183  &other)
6185 {}
6186 
6187 
6188 
6189 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6193  &other)
6194 {
6196  operator=(other);
6197  return *this;
6198 }
6199 
6200 
6201 
6202 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6205  get_gradient(const unsigned int q_point) const
6206 {
6207  return BaseClass::get_gradient(q_point);
6208 }
6209 
6210 
6211 
6212 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6213 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6215  get_divergence(const unsigned int q_point) const
6216 {
6217 # ifdef DEBUG
6218  Assert(this->gradients_quad_initialized == true,
6220 # endif
6221  AssertIndexRange(q_point, this->n_quadrature_points);
6222  Assert(this->jacobian != nullptr, ExcNotInitialized());
6223 
6224  VectorizedArrayType divergence;
6225  const std::size_t nqp = this->n_quadrature_points;
6226 
6227  // Cartesian cell
6228  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6229  {
6230  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6231  for (unsigned int d = 1; d < dim; ++d)
6232  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
6233  this->jacobian[0][d][d];
6234  }
6235  // cell with general/constant Jacobian
6236  else
6237  {
6239  this->cell_type == internal::MatrixFreeFunctions::general ?
6240  this->jacobian[q_point] :
6241  this->jacobian[0];
6242  divergence = jac[0][0] * this->gradients_quad[q_point];
6243  for (unsigned int e = 1; e < dim; ++e)
6244  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
6245  for (unsigned int d = 1; d < dim; ++d)
6246  for (unsigned int e = 0; e < dim; ++e)
6247  divergence +=
6248  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
6249  }
6250  return divergence;
6251 }
6252 
6253 
6254 
6255 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6258  get_symmetric_gradient(const unsigned int q_point) const
6259 {
6260  // copy from generic function into dim-specialization function
6261  const auto grad = get_gradient(q_point);
6262  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6263  VectorizedArrayType half = Number(0.5);
6264  for (unsigned int d = 0; d < dim; ++d)
6265  symmetrized[d] = grad[d][d];
6266  switch (dim)
6267  {
6268  case 1:
6269  break;
6270  case 2:
6271  symmetrized[2] = grad[0][1] + grad[1][0];
6272  symmetrized[2] *= half;
6273  break;
6274  case 3:
6275  symmetrized[3] = grad[0][1] + grad[1][0];
6276  symmetrized[3] *= half;
6277  symmetrized[4] = grad[0][2] + grad[2][0];
6278  symmetrized[4] *= half;
6279  symmetrized[5] = grad[1][2] + grad[2][1];
6280  symmetrized[5] *= half;
6281  break;
6282  default:
6283  Assert(false, ExcNotImplemented());
6284  }
6286 }
6287 
6288 
6289 
6290 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6291 inline DEAL_II_ALWAYS_INLINE
6294  const unsigned int q_point) const
6295 {
6296  // copy from generic function into dim-specialization function
6297  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6299  switch (dim)
6300  {
6301  case 1:
6302  Assert(false,
6303  ExcMessage(
6304  "Computing the curl in 1d is not a useful operation"));
6305  break;
6306  case 2:
6307  curl[0] = grad[1][0] - grad[0][1];
6308  break;
6309  case 3:
6310  curl[0] = grad[2][1] - grad[1][2];
6311  curl[1] = grad[0][2] - grad[2][0];
6312  curl[2] = grad[1][0] - grad[0][1];
6313  break;
6314  default:
6315  Assert(false, ExcNotImplemented());
6316  }
6317  return curl;
6318 }
6319 
6320 
6321 
6322 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6325  get_hessian_diagonal(const unsigned int q_point) const
6326 {
6327  return BaseClass::get_hessian_diagonal(q_point);
6328 }
6329 
6330 
6331 
6332 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6335  const unsigned int q_point) const
6336 {
6337 # ifdef DEBUG
6338  Assert(this->hessians_quad_initialized == true,
6340 # endif
6341  AssertIndexRange(q_point, this->n_quadrature_points);
6342  return BaseClass::get_hessian(q_point);
6343 }
6344 
6345 
6346 
6347 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6348 inline DEAL_II_ALWAYS_INLINE void
6351  const unsigned int q_point)
6352 {
6353  BaseClass::submit_gradient(grad_in, q_point);
6354 }
6355 
6356 
6357 
6358 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6359 inline DEAL_II_ALWAYS_INLINE void
6362  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6363  const unsigned int q_point)
6364 {
6365  BaseClass::submit_gradient(grad_in, q_point);
6366 }
6367 
6368 
6369 
6370 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6371 inline DEAL_II_ALWAYS_INLINE void
6373  submit_divergence(const VectorizedArrayType div_in,
6374  const unsigned int q_point)
6375 {
6377  AssertIndexRange(q_point, this->n_quadrature_points);
6378  Assert(this->J_value != nullptr, ExcNotInitialized());
6379  Assert(this->jacobian != nullptr, ExcNotInitialized());
6380 # ifdef DEBUG
6381  this->gradients_quad_submitted = true;
6382 # endif
6383 
6384  const std::size_t nqp = this->n_quadrature_points;
6385  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6386  {
6387  const VectorizedArrayType fac =
6388  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6389  for (unsigned int d = 0; d < dim; ++d)
6390  {
6391  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6392  (fac * this->jacobian[0][d][d]);
6393  for (unsigned int e = d + 1; e < dim; ++e)
6394  {
6395  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6396  VectorizedArrayType();
6397  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6398  VectorizedArrayType();
6399  }
6400  }
6401  }
6402  else
6403  {
6405  this->cell_type == internal::MatrixFreeFunctions::general ?
6406  this->jacobian[q_point] :
6407  this->jacobian[0];
6408  const VectorizedArrayType fac =
6409  (this->cell_type == internal::MatrixFreeFunctions::general ?
6410  this->J_value[q_point] :
6411  this->J_value[0] * this->quadrature_weights[q_point]) *
6412  div_in;
6413  for (unsigned int d = 0; d < dim; ++d)
6414  {
6415  for (unsigned int e = 0; e < dim; ++e)
6416  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6417  jac[d][e] * fac;
6418  }
6419  }
6420 }
6421 
6422 
6423 
6424 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6425 inline DEAL_II_ALWAYS_INLINE void
6429  const unsigned int q_point)
6430 {
6431  // could have used base class operator, but that involves some overhead
6432  // which is inefficient. it is nice to have the symmetric tensor because
6433  // that saves some operations
6435  AssertIndexRange(q_point, this->n_quadrature_points);
6436  Assert(this->J_value != nullptr, ExcNotInitialized());
6437  Assert(this->jacobian != nullptr, ExcNotInitialized());
6438 # ifdef DEBUG
6439  this->gradients_quad_submitted = true;
6440 # endif
6441 
6442  const std::size_t nqp = this->n_quadrature_points;
6443  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6444  {
6445  const VectorizedArrayType JxW =
6446  this->J_value[0] * this->quadrature_weights[q_point];
6447  for (unsigned int d = 0; d < dim; ++d)
6448  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6449  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6450  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6451  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6452  {
6453  const VectorizedArrayType value =
6454  sym_grad.access_raw_entry(counter) * JxW;
6455  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6456  value * this->jacobian[0][d][d];
6457  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6458  value * this->jacobian[0][e][e];
6459  }
6460  }
6461  // general/affine cell type
6462  else
6463  {
6464  const VectorizedArrayType JxW =
6465  this->cell_type == internal::MatrixFreeFunctions::general ?
6466  this->J_value[q_point] :
6467  this->J_value[0] * this->quadrature_weights[q_point];
6469  this->cell_type == internal::MatrixFreeFunctions::general ?
6470  this->jacobian[q_point] :
6471  this->jacobian[0];
6472  VectorizedArrayType weighted[dim][dim];
6473  for (unsigned int i = 0; i < dim; ++i)
6474  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6475  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6476  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6477  {
6478  const VectorizedArrayType value =
6479  sym_grad.access_raw_entry(counter) * JxW;
6480  weighted[i][j] = value;
6481  weighted[j][i] = value;
6482  }
6483  for (unsigned int comp = 0; comp < dim; ++comp)
6484  for (unsigned int d = 0; d < dim; ++d)
6485  {
6486  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6487  for (unsigned int e = 1; e < dim; ++e)
6488  new_val += jac[e][d] * weighted[comp][e];
6489  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6490  }
6491  }
6492 }
6493 
6494 
6495 
6496 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6497 inline DEAL_II_ALWAYS_INLINE void
6500  const unsigned int q_point)
6501 {
6503  switch (dim)
6504  {
6505  case 1:
6506  Assert(false,
6507  ExcMessage(
6508  "Testing by the curl in 1d is not a useful operation"));
6509  break;
6510  case 2:
6511  grad[1][0] = curl[0];
6512  grad[0][1] = -curl[0];
6513  break;
6514  case 3:
6515  grad[2][1] = curl[0];
6516  grad[1][2] = -curl[0];
6517  grad[0][2] = curl[1];
6518  grad[2][0] = -curl[1];
6519  grad[1][0] = curl[2];
6520  grad[0][1] = -curl[2];
6521  break;
6522  default:
6523  Assert(false, ExcNotImplemented());
6524  }
6525  submit_gradient(grad, q_point);
6526 }
6527 
6528 
6529 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6530 
6531 
6532 template <typename Number, bool is_face, typename VectorizedArrayType>
6535  const unsigned int dof_no,
6536  const unsigned int first_selected_component,
6537  const unsigned int quad_no_in,
6538  const unsigned int fe_degree,
6539  const unsigned int n_q_points,
6540  const bool is_interior_face)
6542  data_in,
6543  dof_no,
6544  first_selected_component,
6545  quad_no_in,
6546  fe_degree,
6547  n_q_points,
6548  is_interior_face)
6549 {}
6550 
6551 
6552 
6553 template <typename Number, bool is_face, typename VectorizedArrayType>
6556  const Mapping<1> & mapping,
6557  const FiniteElement<1> &fe,
6558  const Quadrature<1> & quadrature,
6559  const UpdateFlags update_flags,
6560  const unsigned int first_selected_component,
6563  mapping,
6564  fe,
6565  quadrature,
6566  update_flags,
6567  first_selected_component,
6568  other)
6569 {}
6570 
6571 
6572 
6573 template <typename Number, bool is_face, typename VectorizedArrayType>
6578 {}
6579 
6580 
6581 
6582 template <typename Number, bool is_face, typename VectorizedArrayType>
6586 {
6588  other);
6589  return *this;
6590 }
6591 
6592 
6593 
6594 template <typename Number, bool is_face, typename VectorizedArrayType>
6595 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6597  const unsigned int dof) const
6598 {
6599  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6600  return this->values_dofs[0][dof];
6601 }
6602 
6603 
6604 
6605 template <typename Number, bool is_face, typename VectorizedArrayType>
6606 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6608  const unsigned int q_point) const
6609 {
6610 # ifdef DEBUG
6611  Assert(this->values_quad_initialized == true,
6613 # endif
6614  AssertIndexRange(q_point, this->n_quadrature_points);
6615  return this->values_quad[q_point];
6616 }
6617 
6618 
6619 
6620 template <typename Number, bool is_face, typename VectorizedArrayType>
6623  const unsigned int q_point) const
6624 {
6625  // could use the base class gradient, but that involves too many inefficient
6626  // initialization operations on tensors
6627 
6628 # ifdef DEBUG
6629  Assert(this->gradients_quad_initialized == true,
6631 # endif
6632  AssertIndexRange(q_point, this->n_quadrature_points);
6633 
6635  this->cell_type == internal::MatrixFreeFunctions::general ?
6636  this->jacobian[q_point] :
6637  this->jacobian[0];
6638 
6640  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6641 
6642  return grad_out;
6643 }
6644 
6645 
6646 
6647 template <typename Number, bool is_face, typename VectorizedArrayType>
6648 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6650  get_normal_derivative(const unsigned int q_point) const
6651 {
6652  return BaseClass::get_normal_derivative(q_point)[0];
6653 }
6654 
6655 
6656 
6657 template <typename Number, bool is_face, typename VectorizedArrayType>
6660  const unsigned int q_point) const
6661 {
6662  return BaseClass::get_hessian(q_point)[0];
6663 }
6664 
6665 
6666 
6667 template <typename Number, bool is_face, typename VectorizedArrayType>
6670  get_hessian_diagonal(const unsigned int q_point) const
6671 {
6672  return BaseClass::get_hessian_diagonal(q_point)[0];
6673 }
6674 
6675 
6676 
6677 template <typename Number, bool is_face, typename VectorizedArrayType>
6678 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6680  const unsigned int q_point) const
6681 {
6682  return BaseClass::get_laplacian(q_point)[0];
6683 }
6684 
6685 
6686 
6687 template <typename Number, bool is_face, typename VectorizedArrayType>
6690  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6691 {
6692 # ifdef DEBUG
6693  this->dof_values_initialized = true;
6694  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6695 # endif
6696  this->values_dofs[0][dof] = val_in;
6697 }
6698 
6699 
6700 
6701 template <typename Number, bool is_face, typename VectorizedArrayType>
6702 inline DEAL_II_ALWAYS_INLINE void
6704  const VectorizedArrayType val_in,
6705  const unsigned int q_point)
6706 {
6708  AssertIndexRange(q_point, this->n_quadrature_points);
6709 # ifdef DEBUG
6710  this->values_quad_submitted = true;
6711 # endif
6712 
6713  if (this->cell_type == internal::MatrixFreeFunctions::general)
6714  {
6715  const VectorizedArrayType JxW = this->J_value[q_point];
6716  this->values_quad[q_point] = val_in * JxW;
6717  }
6718  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6719  {
6720  const VectorizedArrayType JxW =
6721  this->J_value[0] * this->quadrature_weights[q_point];
6722  this->values_quad[q_point] = val_in * JxW;
6723  }
6724 }
6725 
6726 
6727 
6728 template <typename Number, bool is_face, typename VectorizedArrayType>
6729 inline DEAL_II_ALWAYS_INLINE void
6731  const Tensor<1, 1, VectorizedArrayType> val_in,
6732  const unsigned int q_point)
6733 {
6734  submit_value(val_in[0], q_point);
6735 }
6736 
6737 
6738 
6739 template <typename Number, bool is_face, typename VectorizedArrayType>
6740 inline DEAL_II_ALWAYS_INLINE void
6742  const Tensor<1, 1, VectorizedArrayType> grad_in,
6743  const unsigned int q_point)
6744 {
6745  submit_gradient(grad_in[0], q_point);
6746 }
6747 
6748 
6749 
6750 template <typename Number, bool is_face, typename VectorizedArrayType>
6751 inline DEAL_II_ALWAYS_INLINE void
6753  const VectorizedArrayType grad_in,
6754  const unsigned int q_point)
6755 {
6757  AssertIndexRange(q_point, this->n_quadrature_points);
6758 # ifdef DEBUG
6759  this->gradients_quad_submitted = true;
6760 # endif
6761 
6763  this->cell_type == internal::MatrixFreeFunctions::general ?
6764  this->jacobian[q_point] :
6765  this->jacobian[0];
6766  const VectorizedArrayType JxW =
6767  this->cell_type == internal::MatrixFreeFunctions::general ?
6768  this->J_value[q_point] :
6769  this->J_value[0] * this->quadrature_weights[q_point];
6770 
6771  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6772 }
6773 
6774 
6775 
6776 template <typename Number, bool is_face, typename VectorizedArrayType>
6777 inline DEAL_II_ALWAYS_INLINE void
6779  submit_normal_derivative(const VectorizedArrayType grad_in,
6780  const unsigned int q_point)
6781 {
6783  grad[0] = grad_in;
6784  BaseClass::submit_normal_derivative(grad, q_point);
6785 }
6786 
6787 
6788 
6789 template <typename Number, bool is_face, typename VectorizedArrayType>
6790 inline DEAL_II_ALWAYS_INLINE void
6793  const unsigned int q_point)
6794 {
6795  BaseClass::submit_normal_derivative(grad_in, q_point);
6796 }
6797 
6798 
6799 
6800 template <typename Number, bool is_face, typename VectorizedArrayType>
6801 inline VectorizedArrayType
6803  integrate_value() const
6804 {
6805  return BaseClass::integrate_value()[0];
6806 }
6807 
6808 
6809 
6810 /*-------------------------- FEEvaluation -----------------------------------*/
6811 
6812 
6813 template <int dim,
6814  int fe_degree,
6815  int n_q_points_1d,
6816  int n_components_,
6817  typename Number,
6818  typename VectorizedArrayType>
6819 inline FEEvaluation<dim,
6820  fe_degree,
6821  n_q_points_1d,
6822  n_components_,
6823  Number,
6824  VectorizedArrayType>::
6825  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
6826  const unsigned int fe_no,
6827  const unsigned int quad_no,
6828  const unsigned int first_selected_component)
6829  : BaseClass(data_in,
6830  fe_no,
6831  first_selected_component,
6832  quad_no,
6833  fe_degree,
6834  static_n_q_points)
6835  , dofs_per_component(this->data->dofs_per_component_on_cell)
6836  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6837  , n_q_points(this->data->n_q_points)
6838 {
6839  check_template_arguments(fe_no, 0);
6840 }
6841 
6842 
6843 
6844 template <int dim,
6845  int fe_degree,
6846  int n_q_points_1d,
6847  int n_components_,
6848  typename Number,
6849  typename VectorizedArrayType>
6850 inline FEEvaluation<dim,
6851  fe_degree,
6852  n_q_points_1d,
6853  n_components_,
6854  Number,
6855  VectorizedArrayType>::
6856  FEEvaluation(const Mapping<dim> & mapping,
6857  const FiniteElement<dim> &fe,
6858  const Quadrature<1> & quadrature,
6859  const UpdateFlags update_flags,
6860  const unsigned int first_selected_component)
6861  : BaseClass(mapping,
6862  fe,
6863  quadrature,
6864  update_flags,
6865  first_selected_component,
6866  nullptr)
6867  , dofs_per_component(this->data->dofs_per_component_on_cell)
6868  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6869  , n_q_points(this->data->n_q_points)
6870 {
6871  check_template_arguments(numbers::invalid_unsigned_int, 0);
6872 }
6873 
6874 
6875 
6876 template <int dim,
6877  int fe_degree,
6878  int n_q_points_1d,
6879  int n_components_,
6880  typename Number,
6881  typename VectorizedArrayType>
6882 inline FEEvaluation<dim,
6883  fe_degree,
6884  n_q_points_1d,
6885  n_components_,
6886  Number,
6887  VectorizedArrayType>::
6888  FEEvaluation(const FiniteElement<dim> &fe,
6889  const Quadrature<1> & quadrature,
6890  const UpdateFlags update_flags,
6891  const unsigned int first_selected_component)
6892  : BaseClass(StaticMappingQ1<dim>::mapping,
6893  fe,
6894  quadrature,
6895  update_flags,
6896  first_selected_component,
6897  nullptr)
6898  , dofs_per_component(this->data->dofs_per_component_on_cell)
6899  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6900  , n_q_points(this->data->n_q_points)
6901 {
6902  check_template_arguments(numbers::invalid_unsigned_int, 0);
6903 }
6904 
6905 
6906 
6907 template <int dim,
6908  int fe_degree,
6909  int n_q_points_1d,
6910  int n_components_,
6911  typename Number,
6912  typename VectorizedArrayType>
6913 inline FEEvaluation<dim,
6914  fe_degree,
6915  n_q_points_1d,
6916  n_components_,
6917  Number,
6918  VectorizedArrayType>::
6919  FEEvaluation(
6920  const FiniteElement<dim> & fe,
6922  const unsigned int first_selected_component)
6923  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6924  fe,
6925  other.mapped_geometry->get_quadrature(),
6926  other.mapped_geometry->get_fe_values().get_update_flags(),
6927  first_selected_component,
6928  &other)
6929  , dofs_per_component(this->data->dofs_per_component_on_cell)
6930  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6931  , n_q_points(this->data->n_q_points)
6932 {
6933  check_template_arguments(numbers::invalid_unsigned_int, 0);
6934 }
6935 
6936 
6937 
6938 template <int dim,
6939  int fe_degree,
6940  int n_q_points_1d,
6941  int n_components_,
6942  typename Number,
6943  typename VectorizedArrayType>
6944 inline FEEvaluation<dim,
6945  fe_degree,
6946  n_q_points_1d,
6947  n_components_,
6948  Number,
6949  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6950  &other)
6951  : BaseClass(other)
6952  , dofs_per_component(this->data->dofs_per_component_on_cell)
6953  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6954  , n_q_points(this->data->n_q_points)
6955 {
6956  check_template_arguments(numbers::invalid_unsigned_int, 0);
6957 }
6958 
6959 
6960 
6961 template <int dim,
6962  int fe_degree,
6963  int n_q_points_1d,
6964  int n_components_,
6965  typename Number,
6966  typename VectorizedArrayType>
6967 inline FEEvaluation<dim,
6968  fe_degree,
6969  n_q_points_1d,
6970  n_components_,
6971  Number,
6972  VectorizedArrayType> &
6973 FEEvaluation<dim,
6974  fe_degree,
6975  n_q_points_1d,
6976  n_components_,
6977  Number,
6978  VectorizedArrayType>::operator=(const FEEvaluation &other)
6979 {
6980  BaseClass::operator=(other);
6981  check_template_arguments(numbers::invalid_unsigned_int, 0);
6982  return *this;
6983 }
6984 
6985 
6986 
6987 template <int dim,
6988  int fe_degree,
6989  int n_q_points_1d,
6990  int n_components_,
6991  typename Number,
6992  typename VectorizedArrayType>
6993 inline void
6994 FEEvaluation<dim,
6995  fe_degree,
6996  n_q_points_1d,
6997  n_components_,
6998  Number,
6999  VectorizedArrayType>::
7000  check_template_arguments(const unsigned int dof_no,
7001  const unsigned int first_selected_component)
7002 {
7003  (void)dof_no;
7004  (void)first_selected_component;
7005 
7006 # ifdef DEBUG
7007  // print error message when the dimensions do not match. Propose a possible
7008  // fix
7009  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7010  static_cast<unsigned int>(fe_degree) !=
7011  this->data->data.front().fe_degree) ||
7012  n_q_points != this->n_quadrature_points)
7013  {
7014  std::string message =
7015  "-------------------------------------------------------\n";
7016  message += "Illegal arguments in constructor/wrong template arguments!\n";
7017  message += " Called --> FEEvaluation<dim,";
7018  message += Utilities::int_to_string(fe_degree) + ",";
7019  message += Utilities::int_to_string(n_q_points_1d);
7020  message += "," + Utilities::int_to_string(n_components);
7021  message += ",Number>(data";
7022  if (first_selected_component != numbers::invalid_unsigned_int)
7023  {
7024  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7025  message += Utilities::int_to_string(this->quad_no) + ", ";
7026  message += Utilities::int_to_string(first_selected_component);
7027  }
7028  message += ")\n";
7029 
7030  // check whether some other vector component has the correct number of
7031  // points
7032  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7033  proposed_fe_comp = numbers::invalid_unsigned_int,
7034  proposed_quad_comp = numbers::invalid_unsigned_int;
7035  if (dof_no != numbers::invalid_unsigned_int)
7036  {
7037  if (static_cast<unsigned int>(fe_degree) ==
7038  this->data->data.front().fe_degree)
7039  {
7040  proposed_dof_comp = dof_no;
7041  proposed_fe_comp = first_selected_component;
7042  }
7043  else
7044  for (unsigned int no = 0; no < this->matrix_info->n_components();
7045  ++no)
7046  for (unsigned int nf = 0;
7047  nf < this->matrix_info->n_base_elements(no);
7048  ++nf)
7049  if (this->matrix_info
7050  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7051  .data.front()
7052  .fe_degree == static_cast<unsigned int>(fe_degree))
7053  {
7054  proposed_dof_comp = no;
7055  proposed_fe_comp = nf;
7056  break;
7057  }
7058  if (n_q_points ==
7059  this->mapping_data->descriptor[this->active_quad_index]
7060  .n_q_points)
7061  proposed_quad_comp = this->quad_no;
7062  else
7063  for (unsigned int no = 0;
7064  no < this->matrix_info->get_mapping_info().cell_data.size();
7065  ++no)
7066  if (this->matrix_info->get_mapping_info()
7067  .cell_data[no]
7068  .descriptor[this->active_quad_index]
7069  .n_q_points == n_q_points)
7070  {
7071  proposed_quad_comp = no;
7072  break;
7073  }
7074  }
7075  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7076  proposed_quad_comp != numbers::invalid_unsigned_int)
7077  {
7078  if (proposed_dof_comp != first_selected_component)
7079  message += "Wrong vector component selection:\n";
7080  else
7081  message += "Wrong quadrature formula selection:\n";
7082  message += " Did you mean FEEvaluation<dim,";
7083  message += Utilities::int_to_string(fe_degree) + ",";
7084  message += Utilities::int_to_string(n_q_points_1d);
7085  message += "," + Utilities::int_to_string(n_components);
7086  message += ",Number>(data";
7087  if (dof_no != numbers::invalid_unsigned_int)
7088  {
7089  message +=
7090  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7091  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7092  message += Utilities::int_to_string(proposed_fe_comp);
7093  }
7094  message += ")?\n";
7095  std::string correct_pos;
7096  if (proposed_dof_comp != dof_no)
7097  correct_pos = " ^ ";
7098  else
7099  correct_pos = " ";
7100  if (proposed_quad_comp != this->quad_no)
7101  correct_pos += " ^ ";
7102  else
7103  correct_pos += " ";
7104  if (proposed_fe_comp != first_selected_component)
7105  correct_pos += " ^\n";
7106  else
7107  correct_pos += " \n";
7108  message += " " +
7109  correct_pos;
7110  }
7111  // ok, did not find the numbers specified by the template arguments in
7112  // the given list. Suggest correct template arguments
7113  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7114  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7115  message += "Wrong template arguments:\n";
7116  message += " Did you mean FEEvaluation<dim,";
7117  message +=
7118  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7119  message += Utilities::int_to_string(proposed_n_q_points_1d);
7120  message += "," + Utilities::int_to_string(n_components);
7121  message += ",Number>(data";
7122  if (dof_no != numbers::invalid_unsigned_int)
7123  {
7124  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7125  message += Utilities::int_to_string(this->quad_no);
7126  message += ", " + Utilities::int_to_string(first_selected_component);
7127  }
7128  message += ")?\n";
7129  std::string correct_pos;
7130  if (this->data->data.front().fe_degree !=
7131  static_cast<unsigned int>(fe_degree))
7132  correct_pos = " ^";
7133  else
7134  correct_pos = " ";
7135  if (proposed_n_q_points_1d != n_q_points_1d)
7136  correct_pos += " ^\n";
7137  else
7138  correct_pos += " \n";
7139  message += " " + correct_pos;
7140 
7141  Assert(static_cast<unsigned int>(fe_degree) ==
7142  this->data->data.front().fe_degree &&
7143  n_q_points == this->n_quadrature_points,
7144  ExcMessage(message));
7145  }
7146  if (dof_no != numbers::invalid_unsigned_int)
7148  n_q_points,
7149  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7150 # endif
7151 }
7152 
7153 
7154 
7155 template <int dim,
7156  int fe_degree,
7157  int n_q_points_1d,
7158  int n_components_,
7159  typename Number,
7160  typename VectorizedArrayType>
7161 inline void
7162 FEEvaluation<dim,
7163  fe_degree,
7164  n_q_points_1d,
7165  n_components_,
7166  Number,
7167  VectorizedArrayType>::reinit(const unsigned int cell_index)
7168 {
7169  Assert(this->mapped_geometry == nullptr,
7170  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7171  " Integer indexing is not possible"));
7172  if (this->mapped_geometry != nullptr)
7173  return;
7174 
7175  Assert(this->dof_info != nullptr, ExcNotInitialized());
7176  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7177  this->cell = cell_index;
7178  this->cell_type =
7179  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
7180 
7181  const unsigned int offsets =
7182  this->mapping_data->data_index_offsets[cell_index];
7183  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7184  this->J_value = &this->mapping_data->JxW_values[offsets];
7185 
7186 # ifdef DEBUG
7187  this->dof_values_initialized = false;
7188  this->values_quad_initialized = false;
7189  this->gradients_quad_initialized = false;
7190  this->hessians_quad_initialized = false;
7191 # endif
7192 }
7193 
7194 
7195 
7196 template <int dim,
7197  int fe_degree,
7198  int n_q_points_1d,
7199  int n_components_,
7200  typename Number,
7201  typename VectorizedArrayType>
7202 template <bool level_dof_access>
7203 inline void
7204 FEEvaluation<dim,
7205  fe_degree,
7206  n_q_points_1d,
7207  n_components_,
7208  Number,
7209  VectorizedArrayType>::
7211 {
7212  Assert(this->matrix_info == nullptr,
7213  ExcMessage("Cannot use initialization from cell iterator if "
7214  "initialized from MatrixFree object. Use variant for "
7215  "on the fly computation with arguments as for FEValues "
7216  "instead"));
7217  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7218  this->mapped_geometry->reinit(
7219  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7220  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7221  if (level_dof_access)
7222  cell->get_mg_dof_indices(this->local_dof_indices);
7223  else
7224  cell->get_dof_indices(this->local_dof_indices);
7225 }
7226 
7227 
7228 
7229 template <int dim,
7230  int fe_degree,
7231  int n_q_points_1d,
7232  int n_components_,
7233  typename Number,
7234  typename VectorizedArrayType>
7235 inline void
7236 FEEvaluation<dim,
7237  fe_degree,
7238  n_q_points_1d,
7239  n_components_,
7240  Number,
7241  VectorizedArrayType>::
7242  reinit(const typename Triangulation<dim>::cell_iterator &cell)
7243 {
7244  Assert(this->matrix_info == 0,
7245  ExcMessage("Cannot use initialization from cell iterator if "
7246  "initialized from MatrixFree object. Use variant for "
7247  "on the fly computation with arguments as for FEValues "
7248  "instead"));
7249  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7250  this->mapped_geometry->reinit(cell);
7251 }
7252 
7253 
7254 
7255 template <int dim,
7256  int fe_degree,
7257  int n_q_points_1d,
7258  int n_components_,
7259  typename Number,
7260  typename VectorizedArrayType>
7262 FEEvaluation<dim,
7263  fe_degree,
7264  n_q_points_1d,
7265  n_components_,
7266  Number,
7267  VectorizedArrayType>::quadrature_point(const unsigned int q) const
7268 {
7269  if (this->matrix_info == nullptr)
7270  {
7271  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7273  ExcNotInitialized());
7274  }
7275  else
7276  {
7277  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
7278  ExcNotInitialized());
7279  }
7280 
7281  AssertIndexRange(q, n_q_points);
7282 
7284  &this->mapping_data->quadrature_points
7285  [this->mapping_data->quadrature_point_offsets[this->cell]];
7286 
7287  // Cartesian/affine mesh: only first vertex of cell is stored, we must
7288  // compute it through the Jacobian (which is stored in non-inverted and
7289  // non-transposed form as index '1' in the jacobian field)
7290  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
7291  {
7292  Assert(this->jacobian != nullptr, ExcNotInitialized());
7293  Point<dim, VectorizedArrayType> point = quadrature_points[0];
7294 
7295  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
7296  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
7297  for (unsigned int d = 0; d < dim; ++d)
7298  point[d] += jac[d][d] *
7299  static_cast<Number>(
7300  this->mapping_data->descriptor[this->active_quad_index]
7301  .quadrature.point(q)[d]);
7302  else
7303  for (unsigned int d = 0; d < dim; ++d)
7304  for (unsigned int e = 0; e < dim; ++e)
7305  point[d] += jac[d][e] * static_cast<Number>(
7306  this->mapping_data
7307  ->descriptor[this->active_quad_index]
7308  .quadrature.point(q)[e]);
7309  return point;
7310  }
7311  else
7312  return quadrature_points[q];
7313 }
7314 
7315 
7316 
7317 template <int dim,
7318  int fe_degree,
7319  int n_q_points_1d,
7320  int n_components_,
7321  typename Number,
7322  typename VectorizedArrayType>
7323 inline void
7324 FEEvaluation<dim,
7325  fe_degree,
7326  n_q_points_1d,
7327  n_components_,
7328  Number,
7329  VectorizedArrayType>::evaluate(const bool evaluate_values,
7330  const bool evaluate_gradients,
7331  const bool evaluate_hessians)
7332 {
7333 # ifdef DEBUG
7334  Assert(this->dof_values_initialized == true,
7336 # endif
7337  evaluate(this->values_dofs[0],
7338  evaluate_values,
7339  evaluate_gradients,
7340  evaluate_hessians);
7341 }
7342 
7343 
7344 template <int dim,
7345  int fe_degree,
7346  int n_q_points_1d,
7347  int n_components_,
7348  typename Number,
7349  typename VectorizedArrayType>
7350 inline void
7351 FEEvaluation<dim,
7352  fe_degree,
7353  n_q_points_1d,
7354  n_components_,
7355  Number,
7356  VectorizedArrayType>::
7357  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7358 {
7359 # ifdef DEBUG
7360  Assert(this->dof_values_initialized == true,
7362 # endif
7363  evaluate(this->values_dofs[0], evaluation_flags);
7364 }
7365 
7366 
7367 
7368 template <int dim,
7369  int fe_degree,
7370  int n_q_points_1d,
7371  int n_components_,
7372  typename Number,
7373  typename VectorizedArrayType>
7374 inline void
7375 FEEvaluation<dim,
7376  fe_degree,
7377  n_q_points_1d,
7378  n_components_,
7379  Number,
7380  VectorizedArrayType>::evaluate(const VectorizedArrayType
7381  * values_array,
7382  const bool evaluate_values,
7383  const bool evaluate_gradients,
7384  const bool evaluate_hessians)
7385 {
7387  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7388  ((evaluate_gradients) ? EvaluationFlags::gradients :
7389  EvaluationFlags::nothing) |
7390  ((evaluate_hessians) ? EvaluationFlags::hessians :
7391  EvaluationFlags::nothing);
7392 
7393  evaluate(values_array, flag);
7394 }
7395 
7396 
7397 
7398 template <int dim,
7399  int fe_degree,
7400  int n_q_points_1d,
7401  int n_components_,
7402  typename Number,
7403  typename VectorizedArrayType>
7404 inline void
7405 FEEvaluation<dim,
7406  fe_degree,
7407  n_q_points_1d,
7408  n_components_,
7409  Number,
7410  VectorizedArrayType>::
7411  evaluate(const VectorizedArrayType * values_array,
7412  const EvaluationFlags::EvaluationFlags evaluation_flags)
7413 {
7414  if (fe_degree > -1)
7416  evaluate(n_components,
7417  evaluation_flags,
7418  *this->data,
7419  const_cast<VectorizedArrayType *>(values_array),
7420  this->values_quad,
7421  this->gradients_quad,
7422  this->hessians_quad,
7423  this->scratch_data);
7424  else
7426  n_components,
7427  evaluation_flags,
7428  *this->data,
7429  const_cast<VectorizedArrayType *>(values_array),
7430  this->values_quad,
7431  this->gradients_quad,
7432  this->hessians_quad,
7433  this->scratch_data);
7434 
7435 # ifdef DEBUG
7436  if (evaluation_flags & EvaluationFlags::values)
7437  this->values_quad_initialized = true;
7438  if (evaluation_flags & EvaluationFlags::gradients)
7439  this->gradients_quad_initialized = true;
7440  if (evaluation_flags & EvaluationFlags::hessians)
7441  this->hessians_quad_initialized = true;
7442 # endif
7443 }
7444 
7445 
7446 
7447 template <int dim,
7448  int fe_degree,
7449  int n_q_points_1d,
7450  int n_components_,
7451  typename Number,
7452  typename VectorizedArrayType>
7453 template <typename VectorType>
7454 inline void
7455 FEEvaluation<
7456  dim,
7457  fe_degree,
7458  n_q_points_1d,
7459  n_components_,
7460  Number,
7461  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7462  const bool evaluate_values,
7463  const bool evaluate_gradients,
7464  const bool evaluate_hessians)
7465 {
7467  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7468  ((evaluate_gradients) ? EvaluationFlags::gradients :
7469  EvaluationFlags::nothing) |
7470  ((evaluate_hessians) ? EvaluationFlags::hessians :
7471  EvaluationFlags::nothing);
7472 
7473  gather_evaluate(input_vector, flag);
7474 }
7475 
7476 
7477 namespace internal
7478 {
7482  template <typename Number,
7483  typename VectorizedArrayType,
7484  typename VectorType,
7485  typename T,
7486  typename std::enable_if<
7488  std::is_same<decltype(std::declval<VectorType>().begin()),
7489  Number *>::value,
7490  VectorType>::type * = nullptr>
7491  bool
7492  try_gather_evaluate_inplace(
7493  T phi,
7494  const VectorType & input_vector,
7495  const unsigned int cell,
7496  const unsigned int active_fe_index,
7497  const unsigned int first_selected_component,
7499  const EvaluationFlags::EvaluationFlags evaluation_flag)
7500  {
7501  // If the index storage is interleaved and contiguous and the vector storage
7502  // has the correct alignment, we can directly pass the pointer into the
7503  // vector to the evaluate() call, without reading the vector entries into a
7504  // separate data field. This saves some operations.
7505  if (std::is_same<typename VectorType::value_type, Number>::value &&
7506  dof_info->index_storage_variants
7509  interleaved_contiguous &&
7510  reinterpret_cast<std::size_t>(
7511  input_vector.begin() +
7512  dof_info->dof_indices_contiguous
7514  [cell * VectorizedArrayType::size()]) %
7515  sizeof(VectorizedArrayType) ==
7516  0)
7517  {
7518  const VectorizedArrayType *vec_values =
7519  reinterpret_cast<const VectorizedArrayType *>(
7520  input_vector.begin() +
7521  dof_info->dof_indices_contiguous
7523  [cell * VectorizedArrayType::size()] +
7524  dof_info->component_dof_indices_offset[active_fe_index]
7525  [first_selected_component] *
7526  VectorizedArrayType::size());
7527 
7528  phi->evaluate(vec_values, evaluation_flag);
7529 
7530  return true;
7531  }
7532 
7533  return false;
7534  }
7535 
7539  template <typename Number,
7540  typename VectorizedArrayType,
7541  typename VectorType,
7542  typename T,
7543  typename std::enable_if<
7545  !std::is_same<decltype(std::declval<VectorType>().begin()),
7546  Number *>::value,
7547  VectorType>::type * = nullptr>
7548  bool
7549  try_gather_evaluate_inplace(T,
7550  const VectorType &,
7551  const unsigned int,
7552  const unsigned int,
7553  const unsigned int,
7556  {
7557  return false;
7558  }
7559 
7563  template <int dim,
7564  int fe_degree,
7565  int n_q_points_1d,
7566  typename Number,
7567  typename VectorizedArrayType,
7568  typename VectorType,
7569  typename std::enable_if<
7571  std::is_same<decltype(std::declval<VectorType>().begin()),
7572  Number *>::value,
7573  VectorType>::type * = nullptr>
7574  bool
7575  try_integrate_scatter_inplace(
7576  VectorType & destination,
7577  const unsigned int cell,
7578  const unsigned int n_components,
7579  const unsigned int active_fe_index,
7580  const unsigned int first_selected_component,
7582  VectorizedArrayType * values_quad,
7583  VectorizedArrayType * gradients_quad,
7584  VectorizedArrayType * scratch_data,
7586  const EvaluationFlags::EvaluationFlags integration_flag)
7587  {
7588  // If the index storage is interleaved and contiguous and the vector storage
7589  // has the correct alignment, we can directly pass the pointer into the
7590  // vector to the integrate() call, without writing temporary results into a
7591  // separate data field that will later be added into the vector. This saves
7592  // some operations.
7593  if (std::is_same<typename VectorType::value_type, Number>::value &&
7594  dof_info->index_storage_variants
7597  interleaved_contiguous &&
7598  reinterpret_cast<std::size_t>(
7599  destination.begin() +
7600  dof_info->dof_indices_contiguous
7602  [cell * VectorizedArrayType::size()]) %
7603  sizeof(VectorizedArrayType) ==
7604  0)
7605  {
7606  VectorizedArrayType *vec_values =
7607  reinterpret_cast<VectorizedArrayType *>(
7608  destination.begin() +
7609  dof_info->dof_indices_contiguous
7611  [cell * VectorizedArrayType::size()] +
7612  dof_info->component_dof_indices_offset[active_fe_index]
7613  [first_selected_component] *
7614  VectorizedArrayType::size());
7615  if (fe_degree > -1)
7617  integrate(n_components,
7618  integration_flag,
7619  *data,
7620  vec_values,
7621  values_quad,
7622  gradients_quad,
7623  scratch_data,
7624  true);
7625  else
7627  n_components,
7628  integration_flag,
7629  *data,
7630  vec_values,
7631  values_quad,
7632  gradients_quad,
7633  scratch_data,
7634  true);
7635 
7636  return true;
7637  }
7638 
7639  return false;
7640  }
7641 
7645  template <int dim,
7646  int fe_degree,
7647  int n_q_points_1d,
7648  typename Number,
7649  typename VectorizedArrayType,
7650  typename VectorType,
7651  typename std::enable_if<
7653  !std::is_same<decltype(std::declval<VectorType>().begin()),
7654  Number *>::value,
7655  VectorType>::type * = nullptr>
7656  bool
7657  try_integrate_scatter_inplace(
7658  VectorType &,
7659  const unsigned int,
7660  const unsigned int,
7661  const unsigned int,
7662  const unsigned int,
7664  const VectorizedArrayType *,
7665  const VectorizedArrayType *,
7666  const VectorizedArrayType *,
7669  {
7670  return false;
7671  }
7672 
7673  template <typename Number,
7674  typename VectorType,
7675  typename std::enable_if<
7677  (std::is_same<decltype(std::declval<VectorType>().begin()),
7678  const double *>::value ||
7679  std::is_same<decltype(std::declval<VectorType>().begin()),
7680  double *>::value ||
7681  std::is_same<decltype(std::declval<VectorType>().begin()),
7682  const float *>::value ||
7683  std::is_same<decltype(std::declval<VectorType>().begin()),
7684  float *>::value),
7685  VectorType>::type * = nullptr>
7686  decltype(std::declval<VectorType>().begin())
7687  get_beginning(VectorType &vec)
7688  {
7689  return vec.begin();
7690  }
7691 
7692  template <typename Number,
7693  typename VectorType,
7694  typename std::enable_if<
7696  !(std::is_same<decltype(std::declval<VectorType>().begin()),
7697  const double *>::value ||
7698  std::is_same<decltype(std::declval<VectorType>().begin()),
7699  double *>::value ||
7700  std::is_same<decltype(std::declval<VectorType>().begin()),
7701  const float *>::value ||
7702  std::is_same<decltype(std::declval<VectorType>().begin()),
7703  float *>::value),
7704  VectorType>::type * = nullptr>
7705  typename VectorType::value_type *
7706  get_beginning(VectorType &)
7707  {
7708  return nullptr;
7709  }
7710 } // namespace internal
7711 
7712 
7713 
7714 template <int dim,
7715  int fe_degree,
7716  int n_q_points_1d,
7717  int n_components_,
7718  typename Number,
7719  typename VectorizedArrayType>
7720 template <typename VectorType>
7721 inline void
7722 FEEvaluation<dim,
7723  fe_degree,
7724  n_q_points_1d,
7725  n_components_,
7726  Number,
7727  VectorizedArrayType>::
7728  gather_evaluate(const VectorType & input_vector,
7729  const EvaluationFlags::EvaluationFlags evaluation_flag)
7730 {
7731  if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
7732  this,
7733  input_vector,
7734  this->cell,
7735  this->active_fe_index,
7736  this->first_selected_component,
7737  this->dof_info,
7738  evaluation_flag) == false)
7739  {
7740  this->read_dof_values(input_vector);
7741  evaluate(this->begin_dof_values(), evaluation_flag);
7742  }
7743 }
7744 
7745 
7746 
7747 template <int dim,
7748  int fe_degree,
7749  int n_q_points_1d,
7750  int n_components_,
7751  typename Number,
7752  typename VectorizedArrayType>
7753 inline void
7754 FEEvaluation<dim,
7755  fe_degree,
7756  n_q_points_1d,
7757  n_components_,
7758  Number,
7759  VectorizedArrayType>::integrate(const bool integrate_values,
7760  const bool integrate_gradients)
7761 {
7762  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7763 
7764 # ifdef DEBUG
7765  this->dof_values_initialized = true;
7766 # endif
7767 }
7768 
7769 
7770 
7771 template <int dim,
7772  int fe_degree,
7773  int n_q_points_1d,
7774  int n_components_,
7775  typename Number,
7776  typename VectorizedArrayType>
7777 inline void
7778 FEEvaluation<dim,
7779  fe_degree,
7780  n_q_points_1d,
7781  n_components_,
7782  Number,
7783  VectorizedArrayType>::
7784  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7785 {
7786  integrate(integration_flag, this->values_dofs[0]);
7787 
7788 # ifdef DEBUG
7789  this->dof_values_initialized = true;
7790 # endif
7791 }
7792 
7793 
7794 
7795 template <int dim,
7796  int fe_degree,
7797  int n_q_points_1d,
7798  int n_components_,
7799  typename Number,
7800  typename VectorizedArrayType>
7801 inline void
7802 FEEvaluation<dim,
7803  fe_degree,
7804  n_q_points_1d,
7805  n_components_,
7806  Number,
7807  VectorizedArrayType>::integrate(const bool integrate_values,
7808  const bool integrate_gradients,
7809  VectorizedArrayType *values_array)
7810 {
7812  (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
7813  (integrate_gradients ? EvaluationFlags::gradients :
7815  integrate(flag, values_array);
7816 }
7817 
7818 
7819 
7820 template <int dim,
7821  int fe_degree,
7822  int n_q_points_1d,
7823  int n_components_,
7824  typename Number,
7825  typename VectorizedArrayType>
7826 inline void
7827 FEEvaluation<dim,
7828  fe_degree,
7829  n_q_points_1d,
7830  n_components_,
7831  Number,
7832  VectorizedArrayType>::
7833  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7834  VectorizedArrayType * values_array)
7835 {
7836 # ifdef DEBUG
7837  if (integration_flag & EvaluationFlags::values)
7838  Assert(this->values_quad_submitted == true,
7840  if (integration_flag & EvaluationFlags::gradients)
7841  Assert(this->gradients_quad_submitted == true,
7843 # endif
7844  Assert(this->matrix_info != nullptr ||
7845  this->mapped_geometry->is_initialized(),
7846  ExcNotInitialized());
7847 
7848  Assert(
7849  (integration_flag &
7850  ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0,
7851  ExcMessage(
7852  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
7853 
7854  if (fe_degree > -1)
7856  integrate(n_components,
7857  integration_flag,
7858  *this->data,
7859  values_array,
7860  this->values_quad,
7861  this->gradients_quad,
7862  this->scratch_data,
7863  false);
7864  else
7866  n_components,
7867  integration_flag,
7868  *this->data,
7869  values_array,
7870  this->values_quad,
7871  this->gradients_quad,
7872  this->scratch_data,
7873  false);
7874 
7875 # ifdef DEBUG
7876  this->dof_values_initialized = true;
7877 # endif
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 template <typename VectorType>
7889 inline void
7890 FEEvaluation<
7891  dim,
7892  fe_degree,
7893  n_q_points_1d,
7894  n_components_,
7895  Number,
7896  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7897  const bool integrate_gradients,
7898  VectorType &destination)
7899 {
7901  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7902  ((integrate_gradients) ? EvaluationFlags::gradients :
7903  EvaluationFlags::nothing);
7904 
7905  integrate_scatter(flag, destination);
7906 }
7907 
7908 
7909 
7910 template <int dim,
7911  int fe_degree,
7912  int n_q_points_1d,
7913  int n_components_,
7914  typename Number,
7915  typename VectorizedArrayType>
7916 template <typename VectorType>
7917 inline void
7918 FEEvaluation<dim,
7919  fe_degree,
7920  n_q_points_1d,
7921  n_components_,
7922  Number,
7923  VectorizedArrayType>::
7924  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
7925  VectorType & destination)
7926 {
7927  if (internal::try_integrate_scatter_inplace<dim,
7928  fe_degree,
7929  n_q_points_1d,
7930  Number,
7931  VectorizedArrayType>(
7932  destination,
7933  this->cell,
7934  n_components,
7935  this->active_fe_index,
7936  this->first_selected_component,
7937  this->dof_info,
7938  this->values_quad,
7939  this->gradients_quad,
7940  this->scratch_data,
7941  this->data,
7942  integration_flag) == false)
7943  {
7944  integrate(integration_flag, this->begin_dof_values());
7945  this->distribute_local_to_global(destination);
7946  }
7947 }
7948 
7949 
7950 
7951 /*-------------------------- FEFaceEvaluation ---------------------------*/
7952 
7953 
7954 
7955 template <int dim,
7956  int fe_degree,
7957  int n_q_points_1d,
7958  int n_components_,
7959  typename Number,
7960  typename VectorizedArrayType>
7961 inline FEFaceEvaluation<dim,
7962  fe_degree,
7963  n_q_points_1d,
7964  n_components_,
7965  Number,
7966  VectorizedArrayType>::
7967  FEFaceEvaluation(
7969  const bool is_interior_face,
7970  const unsigned int dof_no,
7971  const unsigned int quad_no,
7972  const unsigned int first_selected_component)
7973  : BaseClass(matrix_free,
7974  dof_no,
7975  first_selected_component,
7976  quad_no,
7977  fe_degree,
7978  static_n_q_points,
7979  is_interior_face)
7980  , dofs_per_component(this->data->dofs_per_component_on_cell)
7981  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7982  , n_q_points(this->data->n_q_points_face)
7983 {}
7984 
7985 
7986 
7987 template <int dim,
7988  int fe_degree,
7989  int n_q_points_1d,
7990  int n_components_,
7991  typename Number,
7992  typename VectorizedArrayType>
7993 inline void
7994 FEFaceEvaluation<dim,
7995  fe_degree,
7996  n_q_points_1d,
7997  n_components_,
7998  Number,
7999  VectorizedArrayType>::reinit(const unsigned int face_index)
8000 {
8001  Assert(this->mapped_geometry == nullptr,
8002  ExcMessage("FEEvaluation was initialized without a matrix-free object."
8003  " Integer indexing is not possible"));
8004  if (this->mapped_geometry != nullptr)
8005  return;
8006 
8007  this->cell = face_index;
8008  this->dof_access_index =
8009  this->is_interior_face ?
8012  Assert(this->mapping_data != nullptr, ExcNotInitialized());
8014  VectorizedArrayType::size()> &faces =
8015  this->matrix_info->get_face_info(face_index);
8016  if (face_index >=
8017  this->matrix_info->get_task_info().face_partition_data.back() &&
8018  face_index <
8019  this->matrix_info->get_task_info().boundary_partition_data.back())
8020  Assert(this->is_interior_face,
8021  ExcMessage("Boundary faces do not have a neighbor"));
8022 
8023  this->face_no =
8024  (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
8025  this->subface_index = this->is_interior_face == true ?
8027  faces.subface_index;
8028 
8029  // First check if interior or exterior cell has non-standard orientation
8030  // (i.e. the third bit is one or not). Then set zero if this cell has
8031  // standard-orientation else copy the first three bits
8032  // (which is equivalent to modulo 8). See also the documentation of
8033  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
8034  this->face_orientation =
8035  (this->is_interior_face == (faces.face_orientation >= 8)) ?
8036  (faces.face_orientation % 8) :
8037  0;
8038 
8039  this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
8040  const unsigned int offsets =
8041  this->mapping_data->data_index_offsets[face_index];
8042  this->J_value = &this->mapping_data->JxW_values[offsets];
8043  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8044  this->jacobian =
8045  &this->mapping_data->jacobians[!this->is_interior_face][offsets];
8046  this->normal_x_jacobian =
8047  &this->mapping_data
8048  ->normals_times_jacobians[!this->is_interior_face][offsets];
8049 
8050 # ifdef DEBUG
8051  this->dof_values_initialized = false;
8052  this->values_quad_initialized = false;
8053  this->gradients_quad_initialized = false;
8054  this->hessians_quad_initialized = false;
8055 # endif
8056 }
8057 
8058 
8059 
8060 template <int dim,
8061  int fe_degree,
8062  int n_q_points_1d,
8063  int n_components_,
8064  typename Number,
8065  typename VectorizedArrayType>
8066 inline void
8067 FEFaceEvaluation<dim,
8068  fe_degree,
8069  n_q_points_1d,
8070  n_components_,
8071  Number,
8072  VectorizedArrayType>::reinit(const unsigned int cell_index,
8073  const unsigned int face_number)
8074 {
8075  Assert(
8076  this->quad_no <
8077  this->matrix_info->get_mapping_info().face_data_by_cells.size(),
8078  ExcMessage(
8079  "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8081  AssertIndexRange(cell_index,
8082  this->matrix_info->get_mapping_info().cell_type.size());
8083  Assert(this->mapped_geometry == nullptr,
8084  ExcMessage("FEEvaluation was initialized without a matrix-free object."
8085  " Integer indexing is not possible"));
8086  if (this->mapped_geometry != nullptr)
8087  return;
8088  Assert(this->matrix_info != nullptr, ExcNotInitialized());
8089 
8090  this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
8091  this->cell = cell_index;
8092  this->face_orientation = 0;
8093  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
8094  this->face_no = face_number;
8095  this->dof_access_index =
8097 
8098  const unsigned int offsets =
8099  this->matrix_info->get_mapping_info()
8100  .face_data_by_cells[this->quad_no]
8101  .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
8102  face_number];
8103  AssertIndexRange(offsets,
8104  this->matrix_info->get_mapping_info()
8105  .face_data_by_cells[this->quad_no]
8106  .JxW_values.size());
8107  this->J_value = &this->matrix_info->get_mapping_info()
8108  .face_data_by_cells[this->quad_no]
8109  .JxW_values[offsets];
8110  this->normal_vectors = &this->matrix_info->get_mapping_info()
8111  .face_data_by_cells[this->quad_no]
8112  .normal_vectors[offsets];
8113  this->jacobian = &this->matrix_info->get_mapping_info()
8114  .face_data_by_cells[this->quad_no]
8115  .jacobians[!this->is_interior_face][offsets];
8116  this->normal_x_jacobian =
8117  &this->matrix_info->get_mapping_info()
8118  .face_data_by_cells[this->quad_no]
8119  .normals_times_jacobians[!this->is_interior_face][offsets];
8120 
8121 # ifdef DEBUG
8122  this->dof_values_initialized = false;
8123  this->values_quad_initialized = false;
8124  this->gradients_quad_initialized = false;
8125  this->hessians_quad_initialized = false;
8126 # endif
8127 }
8128 
8129 
8130 
8131 template <int dim,
8132  int fe_degree,
8133  int n_q_points_1d,
8134  int n_components,
8135  typename Number,
8136  typename VectorizedArrayType>
8137 inline void
8138 FEFaceEvaluation<dim,
8139  fe_degree,
8140  n_q_points_1d,
8141  n_components,
8142  Number,
8143  VectorizedArrayType>::evaluate(const bool evaluate_values,
8144  const bool evaluate_gradients)
8145 {
8146 # ifdef DEBUG
8147  Assert(this->dof_values_initialized, ExcNotInitialized());
8148 # endif
8149 
8150  evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients);
8151 }
8152 
8153 
8154 
8155 template <int dim,
8156  int fe_degree,
8157  int n_q_points_1d,
8158  int n_components,
8159  typename Number,
8160  typename VectorizedArrayType>
8161 inline void
8162 FEFaceEvaluation<dim,
8163  fe_degree,
8164  n_q_points_1d,
8165  n_components,
8166  Number,
8167  VectorizedArrayType>::
8168  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
8169 {
8170 # ifdef DEBUG
8171  Assert(this->dof_values_initialized, ExcNotInitialized());
8172 # endif
8173 
8174  evaluate(this->values_dofs[0], evaluation_flag);
8175 }
8176 
8177 
8178 
8179 template <int dim,
8180  int fe_degree,
8181  int n_q_points_1d,
8182  int n_components,
8183  typename Number,
8184  typename VectorizedArrayType>
8185 inline void
8186 FEFaceEvaluation<dim,
8187  fe_degree,
8188  n_q_points_1d,
8189  n_components,
8190  Number,
8191  VectorizedArrayType>::evaluate(const VectorizedArrayType
8192  * values_array,
8193  const bool evaluate_values,
8194  const bool evaluate_gradients)
8195 {
8197  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8198  ((evaluate_gradients) ? EvaluationFlags::gradients :
8199  EvaluationFlags::nothing);
8200 
8201  evaluate(values_array, flag);
8202 }
8203 
8204 
8205 
8206 template <int dim,
8207  int fe_degree,
8208  int n_q_points_1d,
8209  int n_components,
8210  typename Number,
8211  typename VectorizedArrayType>
8212 inline void
8213 FEFaceEvaluation<dim,
8214  fe_degree,
8215  n_q_points_1d,
8216  n_components,
8217  Number,
8218  VectorizedArrayType>::
8219  evaluate(const VectorizedArrayType * values_array,
8220  const EvaluationFlags::EvaluationFlags evaluation_flag)
8221 {
8222  Assert(
8223  (evaluation_flag &
8225  ExcMessage(
8226  "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
8227 
8228  if (!(evaluation_flag & EvaluationFlags::values) &&
8229  !(evaluation_flag & EvaluationFlags::gradients))
8230  return;
8231 
8232  if (fe_degree > -1)
8234  template run<fe_degree, n_q_points_1d>(
8235  n_components,
8236  *this->data,
8237  values_array,
8238  this->begin_values(),
8239  this->begin_gradients(),
8240  this->scratch_data,
8241  evaluation_flag & EvaluationFlags::values,
8242  evaluation_flag & EvaluationFlags::gradients,
8243  this->face_no,
8244  this->subface_index,
8245  this->face_orientation,
8246  this->mapping_data->descriptor[this->active_fe_index]
8247  .face_orientations);
8248  else
8250  evaluate(n_components,
8251  *this->data,
8252  values_array,
8253  this->begin_values(),
8254  this->begin_gradients(),
8255  this->scratch_data,
8256  evaluation_flag & EvaluationFlags::values,
8257  evaluation_flag & EvaluationFlags::gradients,
8258  this->face_no,
8259  this->subface_index,
8260  this->face_orientation,
8261  this->mapping_data->descriptor[this->active_fe_index]
8262  .face_orientations);
8263 
8264 # ifdef DEBUG
8265  if (evaluation_flag & EvaluationFlags::values)
8266  this->values_quad_initialized = true;
8267  if (evaluation_flag & EvaluationFlags::gradients)
8268  this->gradients_quad_initialized = true;
8269 # endif
8270 }
8271 
8272 
8273 
8274 template <int dim,
8275  int fe_degree,
8276  int n_q_points_1d,
8277  int n_components,
8278  typename Number,
8279  typename VectorizedArrayType>
8280 inline void
8281 FEFaceEvaluation<dim,
8282  fe_degree,
8283  n_q_points_1d,
8284  n_components,
8285  Number,
8286  VectorizedArrayType>::
8287  integrate(const EvaluationFlags::EvaluationFlags evaluation_flag)
8288 {
8289  integrate(evaluation_flag, this->values_dofs[0]);
8290 
8291 # ifdef DEBUG
8292  this->dof_values_initialized = true;
8293 # endif
8294 }
8295 
8296 
8297 
8298 template <int dim,
8299  int fe_degree,
8300  int n_q_points_1d,
8301  int n_components,
8302  typename Number,
8303  typename VectorizedArrayType>
8304 inline void
8305 FEFaceEvaluation<dim,
8306  fe_degree,
8307  n_q_points_1d,
8308  n_components,
8309  Number,
8310  VectorizedArrayType>::integrate(const bool integrate_values,
8311  const bool integrate_gradients)
8312 {
8313  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
8314 
8315 # ifdef DEBUG
8316  this->dof_values_initialized = true;
8317 # endif
8318 }
8319 
8320 
8321 
8322 template <int dim,
8323  int fe_degree,
8324  int n_q_points_1d,
8325  int n_components,
8326  typename Number,
8327  typename VectorizedArrayType>
8328 inline void
8329 FEFaceEvaluation<dim,
8330  fe_degree,
8331  n_q_points_1d,
8332  n_components,
8333  Number,
8334  VectorizedArrayType>::integrate(const bool integrate_values,
8335  const