Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_evaluation.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2019 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 
23 #include <deal.II/base/array_view.h>
24 #include <deal.II/base/exceptions.h>
25 #include <deal.II/base/smartpointer.h>
27 #include <deal.II/base/template_constraints.h>
28 #include <deal.II/base/vectorization.h>
29 
30 #include <deal.II/lac/vector_operation.h>
31 
32 #include <deal.II/matrix_free/evaluation_kernels.h>
33 #include <deal.II/matrix_free/evaluation_selector.h>
34 #include <deal.II/matrix_free/mapping_data_on_the_fly.h>
35 #include <deal.II/matrix_free/matrix_free.h>
36 #include <deal.II/matrix_free/shape_info.h>
37 #include <deal.II/matrix_free/tensor_product_kernels.h>
38 #include <deal.II/matrix_free/type_traits.h>
39 #include <deal.II/matrix_free/vector_access_internal.h>
40 
41 
42 DEAL_II_NAMESPACE_OPEN
43 
44 
45 
46 namespace internal
47 {
49 }
50 
51 template <int dim,
52  int fe_degree,
53  int n_q_points_1d = fe_degree + 1,
54  int n_components_ = 1,
55  typename Number = double,
56  typename VectorizedArrayType = VectorizedArray<Number>>
58 
59 
92 template <int dim,
93  int n_components_,
94  typename Number,
95  bool is_face = false,
96  typename VectorizedArrayType = VectorizedArray<Number>>
98 {
99  static_assert(
100  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
101  "Type of Number and of VectorizedArrayType do not match.");
102 
103 public:
104  using number_type = Number;
106  using gradient_type =
108  static constexpr unsigned int dimension = dim;
109  static constexpr unsigned int n_components = n_components_;
110 
118  ~FEEvaluationBase();
119 
128  unsigned int
129  get_mapping_data_index_offset() const;
130 
138  get_cell_type() const;
139 
144  get_shape_info() const;
145 
147 
184  template <typename VectorType>
185  void
186  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
187 
216  template <typename VectorType>
217  void
218  read_dof_values_plain(const VectorType & src,
219  const unsigned int first_index = 0);
220 
251  template <typename VectorType>
252  void
253  distribute_local_to_global(
254  VectorType & dst,
255  const unsigned int first_index = 0,
256  const std::bitset<VectorizedArrayType::size()> &mask =
257  std::bitset<VectorizedArrayType::size()>().flip()) const;
258 
291  template <typename VectorType>
292  void
293  set_dof_values(VectorType & dst,
294  const unsigned int first_index = 0,
295  const std::bitset<VectorizedArrayType::size()> &mask =
296  std::bitset<VectorizedArrayType::size()>().flip()) const;
297 
299 
320  value_type
321  get_dof_value(const unsigned int dof) const;
322 
333  void
334  submit_dof_value(const value_type val_in, const unsigned int dof);
335 
347  value_type
348  get_value(const unsigned int q_point) const;
349 
361  void
362  submit_value(const value_type val_in, const unsigned int q_point);
363 
374  get_gradient(const unsigned int q_point) const;
375 
389  value_type
390  get_normal_derivative(const unsigned int q_point) const;
391 
404  void
405  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
406 
424  void
425  submit_normal_derivative(const value_type grad_in,
426  const unsigned int q_point);
427 
439  get_hessian(const unsigned int q_point) const;
440 
450  get_hessian_diagonal(const unsigned int q_point) const;
451 
462  value_type
463  get_laplacian(const unsigned int q_point) const;
464 
465 #ifdef DOXYGEN
466  // doxygen does not anyhow mention functions coming from partial template
467  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
468  // For now, hack-in those functions manually only to fix documentation:
469 
473  VectorizedArrayType
474  get_divergence(const unsigned int q_point) const;
475 
480  get_symmetric_gradient(const unsigned int q_point) const;
481 
486  get_curl(const unsigned int q_point) const;
487 
496  void
497  submit_divergence(const VectorizedArrayType div_in,
498  const unsigned int q_point);
499 
509  void
510  submit_symmetric_gradient(
512  const unsigned int q_point);
513 
523  void
524  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
525  const unsigned int q_point);
526 
527 #endif
528 
545  value_type
546  integrate_value() const;
547 
552  VectorizedArrayType
553  JxW(const unsigned int q_index) const;
554 
560  DEAL_II_DEPRECATED void
561  fill_JxW_values(AlignedVector<VectorizedArrayType> &JxW_values) const;
562 
570  inverse_jacobian(const unsigned int q_index) const;
571 
585  get_normal_vector(const unsigned int q_point) const;
586 
593  VectorizedArrayType
594  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
595 
600  template <typename T>
601  std::array<T, VectorizedArrayType::size()>
602  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
603  &array) const;
604 
609  std::array<unsigned int, VectorizedArrayType::size()>
610  get_cell_ids() const;
611 
613 
626  const VectorizedArrayType *
627  begin_dof_values() const;
628 
637  VectorizedArrayType *
638  begin_dof_values();
639 
650  const VectorizedArrayType *
651  begin_values() const;
652 
663  VectorizedArrayType *
664  begin_values();
665 
677  const VectorizedArrayType *
678  begin_gradients() const;
679 
691  VectorizedArrayType *
692  begin_gradients();
693 
706  const VectorizedArrayType *
707  begin_hessians() const;
708 
721  VectorizedArrayType *
722  begin_hessians();
723 
729  const std::vector<unsigned int> &
730  get_internal_dof_numbering() const;
731 
739  get_scratch_data() const;
740 
742 
743 protected:
754  const unsigned int dof_no,
755  const unsigned int first_selected_component,
756  const unsigned int quad_no,
757  const unsigned int fe_degree,
758  const unsigned int n_q_points,
759  const bool is_interior_face);
760 
795  template <int n_components_other>
796  FEEvaluationBase(const Mapping<dim> & mapping,
797  const FiniteElement<dim> &fe,
798  const Quadrature<1> & quadrature,
799  const UpdateFlags update_flags,
800  const unsigned int first_selected_component,
801  const FEEvaluationBase<dim,
802  n_components_other,
803  Number,
804  is_face,
805  VectorizedArrayType> *other);
806 
813  FEEvaluationBase(const FEEvaluationBase &other);
814 
821  FEEvaluationBase &
822  operator=(const FEEvaluationBase &other);
823 
830  template <typename VectorType, typename VectorOperation>
831  void
832  read_write_operation(const VectorOperation &operation,
833  VectorType * vectors[],
834  const std::bitset<VectorizedArrayType::size()> &mask,
835  const bool apply_constraints = true) const;
836 
844  template <typename VectorType, typename VectorOperation>
845  void
846  read_write_operation_contiguous(
847  const VectorOperation & operation,
848  VectorType * vectors[],
849  const std::bitset<VectorizedArrayType::size()> &mask) const;
850 
858  template <typename VectorType, typename VectorOperation>
859  void
860  read_write_operation_global(const VectorOperation &operation,
861  VectorType * vectors[]) const;
862 
867 
873  VectorizedArrayType *scratch_data;
874 
887  VectorizedArrayType *values_dofs[n_components];
888 
900  VectorizedArrayType *values_quad[n_components];
901 
915  VectorizedArrayType *gradients_quad[n_components][dim];
916 
928  VectorizedArrayType *hessians_quad[n_components][(dim * (dim + 1)) / 2];
929 
933  const unsigned int quad_no;
934 
939  const unsigned int n_fe_components;
940 
945  const unsigned int active_fe_index;
946 
951  const unsigned int active_quad_index;
952 
956  const unsigned int n_quadrature_points;
957 
962 
969 
977  (is_face ? dim - 1 : dim),
978  dim,
979  Number,
980  VectorizedArrayType> *mapping_data;
981 
989 
995 
1002  const VectorizedArrayType *J_value;
1003 
1008 
1013 
1017  const Number *quadrature_weights;
1018 
1023  unsigned int cell;
1024 
1030 
1036 
1041  unsigned int face_no;
1042 
1047  unsigned int face_orientation;
1048 
1056  unsigned int subface_index;
1057 
1065 
1072 
1079 
1086 
1093 
1100 
1107 
1112  std::shared_ptr<internal::MatrixFreeFunctions::
1113  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
1115 
1120  const unsigned int first_selected_component;
1121 
1126  mutable std::vector<types::global_dof_index> local_dof_indices;
1127 
1128 private:
1133  void
1134  set_data_pointers();
1135 
1136  // Make other FEEvaluationBase as well as FEEvaluation objects friends.
1137  template <int, int, typename, bool, typename>
1138  friend class FEEvaluationBase;
1139  template <int, int, int, int, typename, typename>
1140  friend class FEEvaluation;
1141 };
1142 
1143 
1144 
1154 template <int dim,
1155  int n_components_,
1156  typename Number,
1157  bool is_face,
1158  typename VectorizedArrayType = VectorizedArray<Number>>
1160  n_components_,
1161  Number,
1162  is_face,
1163  VectorizedArrayType>
1164 {
1165  static_assert(
1166  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1167  "Type of Number and of VectorizedArrayType do not match.");
1168 
1169 public:
1170  using number_type = Number;
1172  using gradient_type =
1174  static constexpr unsigned int dimension = dim;
1175  static constexpr unsigned int n_components = n_components_;
1176  using BaseClass =
1178 
1179 protected:
1189  const unsigned int dof_no,
1190  const unsigned int first_selected_component,
1191  const unsigned int quad_no,
1192  const unsigned int fe_degree,
1193  const unsigned int n_q_points,
1194  const bool is_interior_face = true);
1195 
1200  template <int n_components_other>
1201  FEEvaluationAccess(const Mapping<dim> & mapping,
1202  const FiniteElement<dim> &fe,
1203  const Quadrature<1> & quadrature,
1204  const UpdateFlags update_flags,
1205  const unsigned int first_selected_component,
1206  const FEEvaluationBase<dim,
1207  n_components_other,
1208  Number,
1209  is_face,
1210  VectorizedArrayType> *other);
1211 
1215  FEEvaluationAccess(const FEEvaluationAccess &other);
1216 
1220  FEEvaluationAccess &
1221  operator=(const FEEvaluationAccess &other);
1222 };
1223 
1224 
1225 
1236 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1237 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
1238  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
1239 {
1240  static_assert(
1241  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1242  "Type of Number and of VectorizedArrayType do not match.");
1243 
1244 public:
1245  using number_type = Number;
1246  using value_type = VectorizedArrayType;
1248  static constexpr unsigned int dimension = dim;
1249  using BaseClass =
1251 
1254  value_type
1255  get_dof_value(const unsigned int dof) const;
1256 
1259  void
1260  submit_dof_value(const value_type val_in, const unsigned int dof);
1261 
1264  value_type
1265  get_value(const unsigned int q_point) const;
1266 
1269  void
1270  submit_value(const value_type val_in, const unsigned int q_point);
1271 
1274  void
1275  submit_value(const Tensor<1, 1, VectorizedArrayType> val_in,
1276  const unsigned int q_point);
1277 
1281  get_gradient(const unsigned int q_point) const;
1282 
1285  value_type
1286  get_normal_derivative(const unsigned int q_point) const;
1287 
1290  void
1291  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1292 
1295  void
1296  submit_normal_derivative(const value_type grad_in,
1297  const unsigned int q_point);
1298 
1302  get_hessian(unsigned int q_point) const;
1303 
1307  get_hessian_diagonal(const unsigned int q_point) const;
1308 
1311  value_type
1312  get_laplacian(const unsigned int q_point) const;
1313 
1316  value_type
1317  integrate_value() const;
1318 
1319 protected:
1329  const unsigned int dof_no,
1330  const unsigned int first_selected_component,
1331  const unsigned int quad_no,
1332  const unsigned int fe_degree,
1333  const unsigned int n_q_points,
1334  const bool is_interior_face = true);
1335 
1340  template <int n_components_other>
1341  FEEvaluationAccess(const Mapping<dim> & mapping,
1342  const FiniteElement<dim> &fe,
1343  const Quadrature<1> & quadrature,
1344  const UpdateFlags update_flags,
1345  const unsigned int first_selected_component,
1346  const FEEvaluationBase<dim,
1347  n_components_other,
1348  Number,
1349  is_face,
1350  VectorizedArrayType> *other);
1351 
1355  FEEvaluationAccess(const FEEvaluationAccess &other);
1356 
1360  FEEvaluationAccess &
1361  operator=(const FEEvaluationAccess &other);
1362 };
1363 
1364 
1365 
1377 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1378 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1379  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1380 {
1381  static_assert(
1382  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1383  "Type of Number and of VectorizedArrayType do not match.");
1384 
1385 public:
1386  using number_type = Number;
1389  static constexpr unsigned int dimension = dim;
1390  static constexpr unsigned int n_components = dim;
1391  using BaseClass =
1393 
1397  get_gradient(const unsigned int q_point) const;
1398 
1403  VectorizedArrayType
1404  get_divergence(const unsigned int q_point) const;
1405 
1413  get_symmetric_gradient(const unsigned int q_point) const;
1414 
1420  get_curl(const unsigned int q_point) const;
1421 
1425  get_hessian(const unsigned int q_point) const;
1426 
1430  get_hessian_diagonal(const unsigned int q_point) const;
1431 
1434  void
1435  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1436 
1445  void
1446  submit_gradient(
1447  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1448  const unsigned int q_point);
1449 
1458  void
1459  submit_divergence(const VectorizedArrayType div_in,
1460  const unsigned int q_point);
1461 
1470  void
1471  submit_symmetric_gradient(
1473  const unsigned int q_point);
1474 
1479  void
1480  submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in,
1481  const unsigned int q_point);
1482 
1483 protected:
1493  const unsigned int dof_no,
1494  const unsigned int first_selected_component,
1495  const unsigned int quad_no,
1496  const unsigned int dofs_per_cell,
1497  const unsigned int n_q_points,
1498  const bool is_interior_face = true);
1499 
1504  template <int n_components_other>
1505  FEEvaluationAccess(const Mapping<dim> & mapping,
1506  const FiniteElement<dim> &fe,
1507  const Quadrature<1> & quadrature,
1508  const UpdateFlags update_flags,
1509  const unsigned int first_selected_component,
1510  const FEEvaluationBase<dim,
1511  n_components_other,
1512  Number,
1513  is_face,
1514  VectorizedArrayType> *other);
1515 
1519  FEEvaluationAccess(const FEEvaluationAccess &other);
1520 
1524  FEEvaluationAccess &
1525  operator=(const FEEvaluationAccess &other);
1526 };
1527 
1528 
1540 template <typename Number, bool is_face, typename VectorizedArrayType>
1541 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1542  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1543 {
1544  static_assert(
1545  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1546  "Type of Number and of VectorizedArrayType do not match.");
1547 
1548 public:
1549  using number_type = Number;
1550  using value_type = VectorizedArrayType;
1552  static constexpr unsigned int dimension = 1;
1553  using BaseClass =
1555 
1558  value_type
1559  get_dof_value(const unsigned int dof) const;
1560 
1563  void
1564  submit_dof_value(const value_type val_in, const unsigned int dof);
1565 
1568  value_type
1569  get_value(const unsigned int q_point) const;
1570 
1573  void
1574  submit_value(const value_type val_in, const unsigned int q_point);
1575 
1578  void
1579  submit_value(const gradient_type val_in, const unsigned int q_point);
1580 
1584  get_gradient(const unsigned int q_point) const;
1585 
1588  value_type
1589  get_normal_derivative(const unsigned int q_point) const;
1590 
1593  void
1594  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1595 
1598  void
1599  submit_gradient(const value_type grad_in, const unsigned int q_point);
1600 
1603  void
1604  submit_normal_derivative(const value_type grad_in,
1605  const unsigned int q_point);
1606 
1609  void
1610  submit_normal_derivative(const gradient_type grad_in,
1611  const unsigned int q_point);
1612 
1616  get_hessian(unsigned int q_point) const;
1617 
1621  get_hessian_diagonal(const unsigned int q_point) const;
1622 
1625  value_type
1626  get_laplacian(const unsigned int q_point) const;
1627 
1630  value_type
1631  integrate_value() const;
1632 
1633 protected:
1642  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
1643  const unsigned int dof_no,
1644  const unsigned int first_selected_component,
1645  const unsigned int quad_no,
1646  const unsigned int fe_degree,
1647  const unsigned int n_q_points,
1648  const bool is_interior_face = true);
1649 
1654  template <int n_components_other>
1655  FEEvaluationAccess(const Mapping<1> & mapping,
1656  const FiniteElement<1> &fe,
1657  const Quadrature<1> & quadrature,
1658  const UpdateFlags update_flags,
1659  const unsigned int first_selected_component,
1660  const FEEvaluationBase<1,
1661  n_components_other,
1662  Number,
1663  is_face,
1664  VectorizedArrayType> *other);
1665 
1669  FEEvaluationAccess(const FEEvaluationAccess &other);
1670 
1674  FEEvaluationAccess &
1675  operator=(const FEEvaluationAccess &other);
1676 };
1677 
1678 
1679 
2225 template <int dim,
2226  int fe_degree,
2227  int n_q_points_1d,
2228  int n_components_,
2229  typename Number,
2230  typename VectorizedArrayType>
2231 class FEEvaluation : public FEEvaluationAccess<dim,
2232  n_components_,
2233  Number,
2234  false,
2235  VectorizedArrayType>
2236 {
2237  static_assert(
2238  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2239  "Type of Number and of VectorizedArrayType do not match.");
2240 
2241 public:
2245  using BaseClass =
2247 
2251  using number_type = Number;
2252 
2259 
2266 
2270  static constexpr unsigned int dimension = dim;
2271 
2276  static constexpr unsigned int n_components = n_components_;
2277 
2284  static constexpr unsigned int static_n_q_points =
2285  Utilities::pow(n_q_points_1d, dim);
2286 
2294  static constexpr unsigned int static_dofs_per_component =
2295  Utilities::pow(fe_degree + 1, dim);
2296 
2304  static constexpr unsigned int tensor_dofs_per_cell =
2305  static_dofs_per_component * n_components;
2306 
2314  static constexpr unsigned int static_dofs_per_cell =
2315  static_dofs_per_component * n_components;
2316 
2343  const unsigned int dof_no = 0,
2344  const unsigned int quad_no = 0,
2345  const unsigned int first_selected_component = 0);
2346 
2373  FEEvaluation(const Mapping<dim> & mapping,
2374  const FiniteElement<dim> &fe,
2375  const Quadrature<1> & quadrature,
2376  const UpdateFlags update_flags,
2377  const unsigned int first_selected_component = 0);
2378 
2384  FEEvaluation(const FiniteElement<dim> &fe,
2385  const Quadrature<1> & quadrature,
2386  const UpdateFlags update_flags,
2387  const unsigned int first_selected_component = 0);
2388 
2399  template <int n_components_other>
2400  FEEvaluation(const FiniteElement<dim> & fe,
2401  const FEEvaluationBase<dim,
2402  n_components_other,
2403  Number,
2404  false,
2405  VectorizedArrayType> &other,
2406  const unsigned int first_selected_component = 0);
2407 
2414  FEEvaluation(const FEEvaluation &other);
2415 
2422  FEEvaluation &
2423  operator=(const FEEvaluation &other);
2424 
2433  void
2434  reinit(const unsigned int cell_batch_index);
2435 
2448  template <typename DoFHandlerType, bool level_dof_access>
2449  void
2451  &cell);
2452 
2463  void
2464  reinit(const typename Triangulation<dim>::cell_iterator &cell);
2465 
2475  void
2476  evaluate(const bool evaluate_values,
2477  const bool evaluate_gradients,
2478  const bool evaluate_hessians = false);
2479 
2492  void
2493  evaluate(const VectorizedArrayType *values_array,
2494  const bool evaluate_values,
2495  const bool evaluate_gradients,
2496  const bool evaluate_hessians = false);
2497 
2511  template <typename VectorType>
2512  void
2513  gather_evaluate(const VectorType &input_vector,
2514  const bool evaluate_values,
2515  const bool evaluate_gradients,
2516  const bool evaluate_hessians = false);
2517 
2528  void
2529  integrate(const bool integrate_values, const bool integrate_gradients);
2530 
2542  void
2543  integrate(const bool integrate_values,
2544  const bool integrate_gradients,
2545  VectorizedArrayType *values_array);
2546 
2560  template <typename VectorType>
2561  void
2562  integrate_scatter(const bool integrate_values,
2563  const bool integrate_gradients,
2564  VectorType &output_vector);
2565 
2571  quadrature_point(const unsigned int q_point) const;
2572 
2579  const unsigned int dofs_per_component;
2580 
2587  const unsigned int dofs_per_cell;
2588 
2596  const unsigned int n_q_points;
2597 
2598 private:
2603  void
2604  check_template_arguments(const unsigned int fe_no,
2605  const unsigned int first_selected_component);
2606 };
2607 
2608 
2609 
2647 template <int dim,
2648  int fe_degree,
2649  int n_q_points_1d = fe_degree + 1,
2650  int n_components_ = 1,
2651  typename Number = double,
2652  typename VectorizedArrayType = VectorizedArray<Number>>
2654  n_components_,
2655  Number,
2656  true,
2657  VectorizedArrayType>
2658 {
2659  static_assert(
2660  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2661  "Type of Number and of VectorizedArrayType do not match.");
2662 
2663 public:
2667  using BaseClass =
2669 
2673  using number_type = Number;
2674 
2681 
2688 
2692  static constexpr unsigned int dimension = dim;
2693 
2698  static constexpr unsigned int n_components = n_components_;
2699 
2707  static constexpr unsigned int static_n_q_points =
2708  Utilities::pow(n_q_points_1d, dim - 1);
2709 
2716  static constexpr unsigned int static_n_q_points_cell =
2717  Utilities::pow(n_q_points_1d, dim);
2718 
2725  static constexpr unsigned int static_dofs_per_component =
2726  Utilities::pow(fe_degree + 1, dim);
2727 
2734  static constexpr unsigned int tensor_dofs_per_cell =
2735  static_dofs_per_component * n_components;
2736 
2743  static constexpr unsigned int static_dofs_per_cell =
2744  static_dofs_per_component * n_components;
2745 
2777  const bool is_interior_face = true,
2778  const unsigned int dof_no = 0,
2779  const unsigned int quad_no = 0,
2780  const unsigned int first_selected_component = 0);
2781 
2792  void
2793  reinit(const unsigned int face_batch_number);
2794 
2802  void
2803  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2804 
2815  void
2816  evaluate(const bool evaluate_values, const bool evaluate_gradients);
2817 
2830  void
2831  evaluate(const VectorizedArrayType *values_array,
2832  const bool evaluate_values,
2833  const bool evaluate_gradients);
2834 
2846  template <typename VectorType>
2847  void
2848  gather_evaluate(const VectorType &input_vector,
2849  const bool evaluate_values,
2850  const bool evaluate_gradients);
2851 
2861  void
2862  integrate(const bool integrate_values, const bool integrate_gradients);
2863 
2872  void
2873  integrate(const bool integrate_values,
2874  const bool integrate_gradients,
2875  VectorizedArrayType *values_array);
2876 
2888  template <typename VectorType>
2889  void
2890  integrate_scatter(const bool integrate_values,
2891  const bool integrate_gradients,
2892  VectorType &output_vector);
2893 
2899  quadrature_point(const unsigned int q_point) const;
2900 
2907  const unsigned int dofs_per_component;
2908 
2915  const unsigned int dofs_per_cell;
2916 
2924  const unsigned int n_q_points;
2925 };
2926 
2927 
2928 
2929 namespace internal
2930 {
2931  namespace MatrixFreeFunctions
2932  {
2933  // a helper function to compute the number of DoFs of a DGP element at
2934  // compile time, depending on the degree
2935  template <int dim, int degree>
2936  struct DGP_dofs_per_component
2937  {
2938  // this division is always without remainder
2939  static constexpr unsigned int value =
2940  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2941  };
2942 
2943  // base specialization: 1d elements have 'degree+1' degrees of freedom
2944  template <int degree>
2945  struct DGP_dofs_per_component<1, degree>
2946  {
2947  static constexpr unsigned int value = degree + 1;
2948  };
2949  } // namespace MatrixFreeFunctions
2950 } // namespace internal
2951 
2952 
2953 /*----------------------- Inline functions ----------------------------------*/
2954 
2955 #ifndef DOXYGEN
2956 
2957 
2958 
2959 /*----------------------- FEEvaluationBase ----------------------------------*/
2960 
2961 template <int dim,
2962  int n_components_,
2963  typename Number,
2964  bool is_face,
2965  typename VectorizedArrayType>
2966 inline FEEvaluationBase<dim,
2967  n_components_,
2968  Number,
2969  is_face,
2970  VectorizedArrayType>::
2971  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
2972  const unsigned int dof_no,
2973  const unsigned int first_selected_component,
2974  const unsigned int quad_no_in,
2975  const unsigned int fe_degree,
2976  const unsigned int n_q_points,
2977  const bool is_interior_face)
2978  : scratch_data_array(data_in.acquire_scratch_data())
2979  , quad_no(quad_no_in)
2980  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
2981  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
2982  data_in.get_dof_info(dof_no).fe_index_from_degree(
2983  first_selected_component,
2984  fe_degree) :
2985  0)
2986  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
2987  (is_face ? data_in.get_mapping_info()
2988  .face_data[quad_no_in]
2989  .quad_index_from_n_q_points(n_q_points) :
2990  data_in.get_mapping_info()
2991  .cell_data[quad_no_in]
2992  .quad_index_from_n_q_points(n_q_points)) :
2993  0)
2994  , n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ?
2995  n_q_points :
2996  (is_face ? data_in
2997  .get_shape_info(dof_no,
2998  quad_no_in,
2999  active_fe_index,
3000  active_quad_index)
3001  .n_q_points_face :
3002  data_in
3003  .get_shape_info(dof_no,
3004  quad_no_in,
3005  active_fe_index,
3006  active_quad_index)
3007  .n_q_points))
3008  , matrix_info(&data_in)
3009  , dof_info(&data_in.get_dof_info(dof_no))
3010  , mapping_data(
3012  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
3013  data_in.get_mapping_info(),
3014  quad_no))
3015  , data(&data_in.get_shape_info(
3016  dof_no,
3017  quad_no_in,
3018  dof_info->component_to_base_index[first_selected_component],
3019  active_fe_index,
3020  active_quad_index))
3021  , jacobian(nullptr)
3022  , J_value(nullptr)
3023  , normal_vectors(nullptr)
3024  , normal_x_jacobian(nullptr)
3025  , quadrature_weights(
3026  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3028  , is_interior_face(is_interior_face)
3029  , dof_access_index(
3030  is_face ?
3031  (is_interior_face ?
3036  , dof_values_initialized(false)
3037  , values_quad_initialized(false)
3038  , gradients_quad_initialized(false)
3039  , hessians_quad_initialized(false)
3040  , values_quad_submitted(false)
3041  , gradients_quad_submitted(false)
3042  , first_selected_component(first_selected_component)
3043 {
3044  set_data_pointers();
3045  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
3047  VectorizedArrayType::size());
3048  AssertDimension((is_face ? data->n_q_points_face : data->n_q_points),
3049  n_quadrature_points);
3050  AssertDimension(n_quadrature_points,
3051  mapping_data->descriptor[active_quad_index].n_q_points);
3052  Assert(
3053  dof_info->start_components.back() == 1 ||
3054  static_cast<int>(n_components_) <=
3055  static_cast<int>(
3056  dof_info->start_components
3057  [dof_info->component_to_base_index[first_selected_component] + 1]) -
3058  first_selected_component,
3059  ExcMessage(
3060  "You tried to construct a vector-valued evaluator with " +
3061  std::to_string(n_components) +
3062  " components. However, "
3063  "the current base element has only " +
3064  std::to_string(
3065  dof_info->start_components
3066  [dof_info->component_to_base_index[first_selected_component] + 1] -
3067  first_selected_component) +
3068  " components left when starting from local element index " +
3069  std::to_string(
3070  first_selected_component -
3071  dof_info->start_components
3072  [dof_info->component_to_base_index[first_selected_component]]) +
3073  " (global index " + std::to_string(first_selected_component) + ")"));
3074 
3075  // do not check for correct dimensions of data fields here, should be done
3076  // in derived classes
3077 }
3078 
3079 
3080 
3081 template <int dim,
3082  int n_components_,
3083  typename Number,
3084  bool is_face,
3085  typename VectorizedArrayType>
3086 template <int n_components_other>
3087 inline FEEvaluationBase<dim,
3088  n_components_,
3089  Number,
3090  is_face,
3091  VectorizedArrayType>::
3092  FEEvaluationBase(const Mapping<dim> & mapping,
3093  const FiniteElement<dim> &fe,
3094  const Quadrature<1> & quadrature,
3095  const UpdateFlags update_flags,
3096  const unsigned int first_selected_component,
3097  const FEEvaluationBase<dim,
3098  n_components_other,
3099  Number,
3100  is_face,
3101  VectorizedArrayType> *other)
3102  : scratch_data_array(new AlignedVector<VectorizedArrayType>())
3104  , n_fe_components(n_components_)
3105  , active_fe_index(numbers::invalid_unsigned_int)
3106  , active_quad_index(numbers::invalid_unsigned_int)
3107  , n_quadrature_points(
3108  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3109  , matrix_info(nullptr)
3110  , dof_info(nullptr)
3111  , mapping_data(nullptr)
3112  ,
3113  // select the correct base element from the given FE component
3115  quadrature,
3116  fe,
3117  fe.component_to_base_index(first_selected_component).first))
3118  , jacobian(nullptr)
3119  , J_value(nullptr)
3120  , normal_vectors(nullptr)
3121  , normal_x_jacobian(nullptr)
3122  , quadrature_weights(nullptr)
3123  , cell(0)
3125  , is_interior_face(true)
3127  , dof_values_initialized(false)
3128  , values_quad_initialized(false)
3129  , gradients_quad_initialized(false)
3130  , hessians_quad_initialized(false)
3131  , values_quad_submitted(false)
3132  , gradients_quad_submitted(false)
3133  ,
3134  // keep the number of the selected component within the current base element
3135  // for reading dof values
3136  first_selected_component(first_selected_component)
3137 {
3138  set_data_pointers();
3139 
3140  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3141  ExcInternalError());
3142  if (other != nullptr &&
3143  other->mapped_geometry->get_quadrature() == quadrature)
3144  mapped_geometry = other->mapped_geometry;
3145  else
3146  mapped_geometry =
3147  std::make_shared<internal::MatrixFreeFunctions::
3148  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3149  mapping, quadrature, update_flags);
3150  cell = 0;
3151 
3152  mapping_data = &mapped_geometry->get_data_storage();
3153  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3154  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3155 
3156  const unsigned int base_element_number =
3157  fe.component_to_base_index(first_selected_component).first;
3158  Assert(fe.element_multiplicity(base_element_number) == 1 ||
3159  fe.element_multiplicity(base_element_number) -
3160  first_selected_component >=
3161  n_components_,
3162  ExcMessage("The underlying element must at least contain as many "
3163  "components as requested by this class"));
3164  (void)base_element_number;
3165 }
3166 
3167 
3168 
3169 template <int dim,
3170  int n_components_,
3171  typename Number,
3172  bool is_face,
3173  typename VectorizedArrayType>
3174 inline FEEvaluationBase<dim,
3175  n_components_,
3176  Number,
3177  is_face,
3178  VectorizedArrayType>::
3179  FEEvaluationBase(const FEEvaluationBase<dim,
3180  n_components_,
3181  Number,
3182  is_face,
3183  VectorizedArrayType> &other)
3184  : scratch_data_array(other.matrix_info == nullptr ?
3186  other.matrix_info->acquire_scratch_data())
3187  , quad_no(other.quad_no)
3188  , n_fe_components(other.n_fe_components)
3189  , active_fe_index(other.active_fe_index)
3190  , active_quad_index(other.active_quad_index)
3191  , n_quadrature_points(other.n_quadrature_points)
3192  , matrix_info(other.matrix_info)
3193  , dof_info(other.dof_info)
3194  , mapping_data(other.mapping_data)
3195  , data(other.matrix_info == nullptr ?
3197  *other.data) :
3198  other.data)
3199  , jacobian(nullptr)
3200  , J_value(nullptr)
3201  , normal_vectors(nullptr)
3202  , normal_x_jacobian(nullptr)
3203  , quadrature_weights(
3204  other.matrix_info == nullptr ?
3205  nullptr :
3206  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3209  , is_interior_face(other.is_interior_face)
3210  , dof_access_index(other.dof_access_index)
3211  , dof_values_initialized(false)
3212  , values_quad_initialized(false)
3213  , gradients_quad_initialized(false)
3214  , hessians_quad_initialized(false)
3215  , values_quad_submitted(false)
3216  , gradients_quad_submitted(false)
3217  , first_selected_component(other.first_selected_component)
3218 {
3219  set_data_pointers();
3220 
3221  // Create deep copy of mapped geometry for use in parallel...
3222  if (other.mapped_geometry.get() != nullptr)
3223  {
3224  mapped_geometry = std::make_shared<
3225  internal::MatrixFreeFunctions::
3226  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3227  other.mapped_geometry->get_fe_values().get_mapping(),
3228  other.mapped_geometry->get_quadrature(),
3229  other.mapped_geometry->get_fe_values().get_update_flags());
3230  mapping_data = &mapped_geometry->get_data_storage();
3231  cell = 0;
3232 
3233  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3234  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3235  }
3236 }
3237 
3238 
3239 
3240 template <int dim,
3241  int n_components_,
3242  typename Number,
3243  bool is_face,
3244  typename VectorizedArrayType>
3245 inline FEEvaluationBase<dim,
3246  n_components_,
3247  Number,
3248  is_face,
3249  VectorizedArrayType> &
3251 operator=(const FEEvaluationBase<dim,
3252  n_components_,
3253  Number,
3254  is_face,
3255  VectorizedArrayType> &other)
3256 {
3257  AssertDimension(quad_no, other.quad_no);
3258  AssertDimension(n_fe_components, other.n_fe_components);
3259  AssertDimension(active_fe_index, other.active_fe_index);
3260  AssertDimension(active_quad_index, other.active_quad_index);
3261  AssertDimension(first_selected_component, other.first_selected_component);
3262 
3263  // release old memory
3264  if (matrix_info == nullptr)
3265  {
3266  delete data;
3267  delete scratch_data_array;
3268  }
3269  else
3270  {
3271  matrix_info->release_scratch_data(scratch_data_array);
3272  }
3273 
3274  matrix_info = other.matrix_info;
3275  dof_info = other.dof_info;
3276  mapping_data = other.mapping_data;
3277  if (other.matrix_info == nullptr)
3278  {
3280  *other.data);
3281  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3282  }
3283  else
3284  {
3285  data = other.data;
3286  scratch_data_array = matrix_info->acquire_scratch_data();
3287  }
3288  set_data_pointers();
3289 
3290  quadrature_weights =
3291  (mapping_data != nullptr ?
3292  mapping_data->descriptor[active_quad_index].quadrature_weights.begin() :
3293  nullptr);
3296  is_interior_face = other.is_interior_face;
3297  dof_access_index = other.dof_access_index;
3298 
3299  // Create deep copy of mapped geometry for use in parallel...
3300  if (other.mapped_geometry.get() != nullptr)
3301  {
3302  mapped_geometry = std::make_shared<
3303  internal::MatrixFreeFunctions::
3304  MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
3305  other.mapped_geometry->get_fe_values().get_mapping(),
3306  other.mapped_geometry->get_quadrature(),
3307  other.mapped_geometry->get_fe_values().get_update_flags());
3308  cell = 0;
3309  mapping_data = &mapped_geometry->get_data_storage();
3310  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3311  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3312  }
3313 
3314  return *this;
3315 }
3316 
3317 
3318 
3319 template <int dim,
3320  int n_components_,
3321  typename Number,
3322  bool is_face,
3323  typename VectorizedArrayType>
3324 inline FEEvaluationBase<dim,
3325  n_components_,
3326  Number,
3327  is_face,
3328  VectorizedArrayType>::~FEEvaluationBase()
3329 {
3330  if (matrix_info != nullptr)
3331  {
3332  try
3333  {
3334  matrix_info->release_scratch_data(scratch_data_array);
3335  }
3336  catch (...)
3337  {}
3338  }
3339  else
3340  {
3341  delete scratch_data_array;
3342  delete data;
3343  data = nullptr;
3344  }
3345  scratch_data_array = nullptr;
3346 }
3347 
3348 
3349 
3350 template <int dim,
3351  int n_components_,
3352  typename Number,
3353  bool is_face,
3354  typename VectorizedArrayType>
3355 inline void
3358 {
3359  Assert(scratch_data_array != nullptr, ExcInternalError());
3360 
3361  const unsigned int tensor_dofs_per_component =
3362  Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
3363  const unsigned int dofs_per_component =
3364  this->data->dofs_per_component_on_cell;
3365  const unsigned int n_quadrature_points =
3366  is_face ? this->data->n_q_points_face : this->data->n_q_points;
3367 
3368  const unsigned int shift =
3369  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
3370  n_components_ * 3 +
3371  2 * n_quadrature_points;
3372  const unsigned int allocated_size =
3373  shift + n_components_ * dofs_per_component +
3374  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3375  scratch_data_array->resize_fast(allocated_size);
3376 
3377  // set the pointers to the correct position in the data array
3378  for (unsigned int c = 0; c < n_components_; ++c)
3379  {
3380  this->values_dofs[c] =
3381  scratch_data_array->begin() + c * dofs_per_component;
3382  this->values_quad[c] = scratch_data_array->begin() +
3383  n_components * dofs_per_component +
3384  c * n_quadrature_points;
3385  for (unsigned int d = 0; d < dim; ++d)
3386  this->gradients_quad[c][d] =
3387  scratch_data_array->begin() +
3388  n_components * (dofs_per_component + n_quadrature_points) +
3389  (c * dim + d) * n_quadrature_points;
3390  for (unsigned int d = 0; d < (dim * dim + dim) / 2; ++d)
3391  this->hessians_quad[c][d] =
3392  scratch_data_array->begin() +
3393  n_components *
3394  ((dim + 1) * n_quadrature_points + dofs_per_component) +
3395  (c * (dim * dim + dim) + d) * n_quadrature_points;
3396  }
3397  scratch_data =
3398  scratch_data_array->begin() + n_components_ * dofs_per_component +
3399  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3400 }
3401 
3402 
3403 
3404 template <int dim,
3405  int n_components_,
3406  typename Number,
3407  bool is_face,
3408  typename VectorizedArrayType>
3409 inline unsigned int
3412 {
3413  if (matrix_info == nullptr)
3414  return 0;
3415  else
3416  {
3417  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3418  return this->mapping_data->data_index_offsets[cell];
3419  }
3420 }
3421 
3422 
3423 
3424 template <int dim,
3425  int n_components_,
3426  typename Number,
3427  bool is_face,
3428  typename VectorizedArrayType>
3431  get_cell_type() const
3432 {
3434  return cell_type;
3435 }
3436 
3437 
3438 
3439 template <int dim,
3440  int n_components_,
3441  typename Number,
3442  bool is_face,
3443  typename VectorizedArrayType>
3446  get_shape_info() const
3447 {
3448  Assert(data != nullptr, ExcInternalError());
3449  return *data;
3450 }
3451 
3452 
3453 
3454 template <int dim,
3455  int n_components_,
3456  typename Number,
3457  bool is_face,
3458  typename VectorizedArrayType>
3459 inline void
3462 {
3463  AssertDimension(JxW_values.size(), n_quadrature_points);
3464  Assert(J_value != nullptr, ExcNotInitialized());
3465  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3466  {
3467  VectorizedArrayType J = J_value[0];
3468  for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
3469  JxW_values[q] = J * this->quadrature_weights[q];
3470  }
3471  else
3472  for (unsigned int q = 0; q < n_quadrature_points; ++q)
3473  JxW_values[q] = J_value[q];
3474 }
3475 
3476 
3477 
3478 template <int dim,
3479  int n_components_,
3480  typename Number,
3481  bool is_face,
3482  typename VectorizedArrayType>
3483 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
3485  get_normal_vector(const unsigned int q_index) const
3486 {
3487  AssertIndexRange(q_index, n_quadrature_points);
3488  Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!"));
3489  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3490  return normal_vectors[0];
3491  else
3492  return normal_vectors[q_index];
3493 }
3494 
3495 
3496 
3497 template <int dim,
3498  int n_components_,
3499  typename Number,
3500  bool is_face,
3501  typename VectorizedArrayType>
3502 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
3504  const unsigned int q_index) const
3505 {
3506  AssertIndexRange(q_index, n_quadrature_points);
3507  Assert(J_value != nullptr, ExcNotInitialized());
3508  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3509  {
3510  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3511  return J_value[0] * this->quadrature_weights[q_index];
3512  }
3513  else
3514  return J_value[q_index];
3515 }
3516 
3517 
3518 
3519 template <int dim,
3520  int n_components_,
3521  typename Number,
3522  bool is_face,
3523  typename VectorizedArrayType>
3526  inverse_jacobian(const unsigned int q_index) const
3527 {
3528  AssertIndexRange(q_index, n_quadrature_points);
3529  Assert(this->jacobian != nullptr, ExcNotImplemented());
3530  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3531  return jacobian[0];
3532  else
3533  return jacobian[q_index];
3534 }
3535 
3536 template <int dim,
3537  int n_components_,
3538  typename Number,
3539  bool is_face,
3540  typename VectorizedArrayType>
3541 std::array<unsigned int, VectorizedArrayType::size()>
3543  get_cell_ids() const
3544 {
3545  const unsigned int n_lanes = VectorizedArrayType::size();
3546  std::array<unsigned int, n_lanes> cells;
3547 
3548  // initialize array
3549  for (unsigned int i = 0; i < n_lanes; ++i)
3550  cells[i] = numbers::invalid_unsigned_int;
3551 
3552  if ((is_face == false) ||
3553  (is_face &&
3554  this->dof_access_index ==
3556  this->is_interior_face))
3557  {
3558  // cell or interior face face (element-centric loop)
3559  for (unsigned int i = 0; i < n_lanes; ++i)
3560  cells[i] = cell * n_lanes + i;
3561  }
3562  else if (is_face &&
3563  this->dof_access_index ==
3565  this->is_interior_face == false)
3566  {
3567  // exterior face (element-centric loop): for this case, we need to
3568  // look into the FaceInfo field that collects information from both
3569  // sides of a face once for the global mesh, and pick the face id that
3570  // is not the local one (cell_this).
3571  for (unsigned int i = 0; i < n_lanes; i++)
3572  {
3573  // compute actual (non vectorized) cell ID
3574  const unsigned int cell_this = this->cell * n_lanes + i;
3575  // compute face ID
3576  unsigned int face_index =
3577  this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
3578  this->face_no,
3579  i);
3580 
3581  if (face_index == numbers::invalid_unsigned_int)
3582  continue; // invalid face ID: no neighbor on boundary
3583 
3584  // get cell ID on both sides of face
3585  auto cell_m = this->matrix_info->get_face_info(face_index / n_lanes)
3586  .cells_interior[face_index % n_lanes];
3587  auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
3588  .cells_exterior[face_index % n_lanes];
3589 
3590  // compare the IDs with the given cell ID
3591  if (cell_m == cell_this)
3592  cells[i] = cell_p; // neighbor has the other ID
3593  else if (cell_p == cell_this)
3594  cells[i] = cell_m;
3595  }
3596  }
3597  else if (is_face)
3598  {
3599  // face-centric faces
3600  const unsigned int *cells_ =
3601  is_interior_face ?
3602  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3603  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3604  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3605  if (cells_[i] != numbers::invalid_unsigned_int)
3606  cells[i] = cells_[i];
3607  }
3608 
3609  return cells;
3610 }
3611 
3612 
3613 
3614 template <int dim,
3615  int n_components_,
3616  typename Number,
3617  bool is_face,
3618  typename VectorizedArrayType>
3619 inline VectorizedArrayType
3622 {
3623  Assert(matrix_info != nullptr, ExcNotImplemented());
3624  AssertDimension(array.size(),
3625  matrix_info->get_task_info().cell_partition_data.back());
3626 
3627  // 1) collect ids of cell
3628  const auto cells = this->get_cell_ids();
3629 
3630  // 2) actually gather values
3631  VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
3632  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3633  if (cells[i] != numbers::invalid_unsigned_int)
3634  out[i] = array[cells[i] / VectorizedArrayType::size()]
3635  [cells[i] % VectorizedArrayType::size()];
3636  return out;
3637 }
3638 
3639 
3640 
3641 template <int dim,
3642  int n_components_,
3643  typename Number,
3644  bool is_face,
3645  typename VectorizedArrayType>
3646 template <typename T>
3647 inline std::array<T, VectorizedArrayType::size()>
3649  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
3650  &array) const
3651 {
3652  Assert(matrix_info != nullptr, ExcNotImplemented());
3653  AssertDimension(array.size(),
3654  matrix_info->get_task_info().cell_partition_data.back());
3655 
3656  // 1) collect ids of cell
3657  const auto cells = this->get_cell_ids();
3658 
3659  // 2) actually gather values
3660  std::array<T, VectorizedArrayType::size()> out;
3661  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3662  if (cells[i] != numbers::invalid_unsigned_int)
3663  out[i] = array[cells[i] / VectorizedArrayType::size()]
3664  [cells[i] % VectorizedArrayType::size()];
3665  return out;
3666 }
3667 
3668 
3669 
3670 namespace internal
3671 {
3672  // allows to select between block vectors and non-block vectors, which
3673  // allows to use a unified interface for extracting blocks on block vectors
3674  // and doing nothing on usual vectors
3675  template <typename VectorType, bool>
3676  struct BlockVectorSelector
3677  {};
3678 
3679  template <typename VectorType>
3680  struct BlockVectorSelector<VectorType, true>
3681  {
3682  using BaseVectorType = typename VectorType::BlockType;
3683 
3684  static BaseVectorType *
3685  get_vector_component(VectorType &vec, const unsigned int component)
3686  {
3687  AssertIndexRange(component, vec.n_blocks());
3688  return &vec.block(component);
3689  }
3690  };
3691 
3692  template <typename VectorType>
3693  struct BlockVectorSelector<VectorType, false>
3694  {
3695  using BaseVectorType = VectorType;
3696 
3697  static BaseVectorType *
3698  get_vector_component(VectorType &vec, const unsigned int component)
3699  {
3700  // FEEvaluation allows to combine several vectors from a scalar
3701  // FiniteElement into a "vector-valued" FEEvaluation object with
3702  // multiple components. These components can be extracted with the other
3703  // get_vector_component functions. If we do not get a vector of vectors
3704  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3705  // must make sure that we do not duplicate the components in input
3706  // and/or duplicate the resulting integrals. In such a case, we should
3707  // only get the zeroth component in the vector contained set nullptr for
3708  // the others which allows us to catch unintended use in
3709  // read_write_operation.
3710  if (component == 0)
3711  return &vec;
3712  else
3713  return nullptr;
3714  }
3715  };
3716 
3717  template <typename VectorType>
3718  struct BlockVectorSelector<std::vector<VectorType>, false>
3719  {
3720  using BaseVectorType = VectorType;
3721 
3722  static BaseVectorType *
3723  get_vector_component(std::vector<VectorType> &vec,
3724  const unsigned int component)
3725  {
3726  AssertIndexRange(component, vec.size());
3727  return &vec[component];
3728  }
3729  };
3730 
3731  template <typename VectorType>
3732  struct BlockVectorSelector<std::vector<VectorType *>, false>
3733  {
3734  using BaseVectorType = VectorType;
3735 
3736  static BaseVectorType *
3737  get_vector_component(std::vector<VectorType *> &vec,
3738  const unsigned int component)
3739  {
3740  AssertIndexRange(component, vec.size());
3741  return vec[component];
3742  }
3743  };
3744 } // namespace internal
3745 
3746 
3747 
3748 template <int dim,
3749  int n_components_,
3750  typename Number,
3751  bool is_face,
3752  typename VectorizedArrayType>
3753 template <typename VectorType, typename VectorOperation>
3754 inline void
3756  read_write_operation(const VectorOperation &operation,
3757  VectorType * src[],
3758  const std::bitset<VectorizedArrayType::size()> &mask,
3759  const bool apply_constraints) const
3760 {
3761  // Case 1: No MatrixFree object given, simple case because we do not need to
3762  // process constraints and need not care about vectorization -> go to
3763  // separate function
3764  if (matrix_info == nullptr)
3765  {
3766  read_write_operation_global(operation, src);
3767  return;
3768  }
3769 
3770  Assert(dof_info != nullptr, ExcNotInitialized());
3771  Assert(matrix_info->indices_initialized() == true, ExcNotInitialized());
3772  if (n_fe_components == 1)
3773  for (unsigned int comp = 0; comp < n_components; ++comp)
3774  {
3775  Assert(src[comp] != nullptr,
3776  ExcMessage("The finite element underlying this FEEvaluation "
3777  "object is scalar, but you requested " +
3778  std::to_string(n_components) +
3779  " components via the template argument in "
3780  "FEEvaluation. In that case, you must pass an "
3781  "std::vector<VectorType> or a BlockVector to " +
3782  "read_dof_values and distribute_local_to_global."));
3783  internal::check_vector_compatibility(*src[comp], *dof_info);
3784  }
3785  else
3786  {
3787  internal::check_vector_compatibility(*src[0], *dof_info);
3788  }
3789 
3790  // Case 2: contiguous indices which use reduced storage of indices and can
3791  // use vectorized load/store operations -> go to separate function
3792  AssertIndexRange(cell,
3793  dof_info->index_storage_variants[dof_access_index].size());
3794  if (dof_info->index_storage_variants
3795  [is_face ? dof_access_index :
3797  [cell] >=
3799  {
3800  read_write_operation_contiguous(operation, src, mask);
3801  return;
3802  }
3803 
3804  // Case 3: standard operation with one index per degree of freedom -> go on
3805  // here
3806  constexpr unsigned int n_lanes = VectorizedArrayType::size();
3807  Assert(mask.count() == n_lanes,
3808  ExcNotImplemented("Masking currently not implemented for "
3809  "non-contiguous DoF storage"));
3810 
3811  std::integral_constant<bool,
3812  internal::is_vectorizable<VectorType, Number>::value>
3813  vector_selector;
3814 
3815  const unsigned int dofs_per_component =
3816  this->data->dofs_per_component_on_cell;
3817  if (dof_info->index_storage_variants
3818  [is_face ? dof_access_index :
3820  [cell] ==
3822  {
3823  const unsigned int *dof_indices =
3824  dof_info->dof_indices_interleaved.data() +
3825  dof_info->row_starts[cell * n_fe_components * n_lanes].first +
3826  dof_info->component_dof_indices_offset[active_fe_index]
3827  [first_selected_component] *
3828  n_lanes;
3829  if (n_components == 1 || n_fe_components == 1)
3830  for (unsigned int i = 0; i < dofs_per_component;
3831  ++i, dof_indices += n_lanes)
3832  for (unsigned int comp = 0; comp < n_components; ++comp)
3833  operation.process_dof_gather(dof_indices,
3834  *src[comp],
3835  0,
3836  values_dofs[comp][i],
3837  vector_selector);
3838  else
3839  for (unsigned int comp = 0; comp < n_components; ++comp)
3840  for (unsigned int i = 0; i < dofs_per_component;
3841  ++i, dof_indices += n_lanes)
3842  operation.process_dof_gather(
3843  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3844  return;
3845  }
3846 
3847  const unsigned int * dof_indices[n_lanes];
3848  VectorizedArrayType **values_dofs =
3849  const_cast<VectorizedArrayType **>(&this->values_dofs[0]);
3850 
3851  // Assign the appropriate cell ids for face/cell case and get the pointers
3852  // to the dof indices of the cells on all lanes
3853  unsigned int cells_copied[n_lanes];
3854  const unsigned int *cells;
3855  unsigned int n_vectorization_actual =
3856  dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
3857  bool has_constraints = false;
3858  if (is_face)
3859  {
3860  if (dof_access_index ==
3862  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3863  cells_copied[v] = cell * VectorizedArrayType::size() + v;
3864  cells = dof_access_index ==
3866  &cells_copied[0] :
3867  (is_interior_face ?
3868  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3869  &this->matrix_info->get_face_info(cell).cells_exterior[0]);
3870  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3871  {
3872  Assert(cells[v] < dof_info->row_starts.size() - 1,
3873  ExcInternalError());
3874  const std::pair<unsigned int, unsigned int> *my_index_start =
3875  &dof_info->row_starts[cells[v] * n_fe_components +
3876  first_selected_component];
3877 
3878  // check whether any of the SIMD lanes has constraints, i.e., the
3879  // constraint indicator which is the second entry of row_starts
3880  // increments on this cell
3881  if (my_index_start[n_components].second != my_index_start[0].second)
3882  has_constraints = true;
3883 
3884  dof_indices[v] =
3885  dof_info->dof_indices.data() + my_index_start[0].first;
3886  }
3887  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
3888  dof_indices[v] = nullptr;
3889  }
3890  else
3891  {
3892  AssertIndexRange((cell + 1) * n_lanes * n_fe_components,
3893  dof_info->row_starts.size());
3894  const unsigned int n_components_read =
3895  n_fe_components > 1 ? n_components : 1;
3896  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3897  {
3898  const std::pair<unsigned int, unsigned int> *my_index_start =
3899  &dof_info->row_starts[(cell * n_lanes + v) * n_fe_components +
3900  first_selected_component];
3901  if (my_index_start[n_components_read].second !=
3902  my_index_start[0].second)
3903  has_constraints = true;
3904  Assert(my_index_start[n_components_read].first ==
3905  my_index_start[0].first ||
3906  my_index_start[0].first < dof_info->dof_indices.size(),
3907  ExcIndexRange(0,
3908  my_index_start[0].first,
3909  dof_info->dof_indices.size()));
3910  dof_indices[v] =
3911  dof_info->dof_indices.data() + my_index_start[0].first;
3912  }
3913  for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
3914  dof_indices[v] = nullptr;
3915  }
3916 
3917  // Case where we have no constraints throughout the whole cell: Can go
3918  // through the list of DoFs directly
3919  if (!has_constraints)
3920  {
3921  if (n_vectorization_actual < n_lanes)
3922  for (unsigned int comp = 0; comp < n_components; ++comp)
3923  for (unsigned int i = 0; i < dofs_per_component; ++i)
3924  operation.process_empty(values_dofs[comp][i]);
3925  if (n_components == 1 || n_fe_components == 1)
3926  {
3927  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3928  for (unsigned int i = 0; i < dofs_per_component; ++i)
3929  for (unsigned int comp = 0; comp < n_components; ++comp)
3930  operation.process_dof(dof_indices[v][i],
3931  *src[comp],
3932  values_dofs[comp][i][v]);
3933  }
3934  else
3935  {
3936  for (unsigned int comp = 0; comp < n_components; ++comp)
3937  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3938  for (unsigned int i = 0; i < dofs_per_component; ++i)
3939  operation.process_dof(
3940  dof_indices[v][comp * dofs_per_component + i],
3941  *src[0],
3942  values_dofs[comp][i][v]);
3943  }
3944  return;
3945  }
3946 
3947  // In the case where there are some constraints to be resolved, loop over
3948  // all vector components that are filled and then over local dofs. ind_local
3949  // holds local number on cell, index iterates over the elements of
3950  // index_local_to_global and dof_indices points to the global indices stored
3951  // in index_local_to_global
3952  if (n_vectorization_actual < n_lanes)
3953  for (unsigned int comp = 0; comp < n_components; ++comp)
3954  for (unsigned int i = 0; i < dofs_per_component; ++i)
3955  operation.process_empty(values_dofs[comp][i]);
3956  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
3957  {
3958  unsigned int index_indicators, next_index_indicators;
3959  const unsigned int n_components_read =
3960  n_fe_components > 1 ? n_components : 1;
3961  if (is_face)
3962  {
3963  index_indicators = dof_info
3964  ->row_starts[cells[v] * n_fe_components +
3965  first_selected_component]
3966  .second;
3967  next_index_indicators = dof_info
3968  ->row_starts[cells[v] * n_fe_components +
3969  first_selected_component + 1]
3970  .second;
3971  }
3972  else
3973  {
3974  index_indicators =
3975  dof_info
3976  ->row_starts[(cell * n_lanes + v) * n_fe_components +
3977  first_selected_component]
3978  .second;
3979  next_index_indicators =
3980  dof_info
3981  ->row_starts[(cell * n_lanes + v) * n_fe_components +
3982  first_selected_component + 1]
3983  .second;
3984  }
3985 
3986  if (apply_constraints == false &&
3987  dof_info
3988  ->row_starts[(cell * n_lanes + v) * n_fe_components +
3989  first_selected_component]
3990  .second !=
3991  dof_info
3992  ->row_starts[(cell * n_lanes + v) * n_fe_components +
3993  first_selected_component + n_components_read]
3994  .second)
3995  {
3996  Assert(dof_info->row_starts_plain_indices[cell * n_lanes + v] !=
3998  ExcNotInitialized());
3999  dof_indices[v] =
4000  dof_info->plain_dof_indices.data() +
4001  dof_info->component_dof_indices_offset[active_fe_index]
4002  [first_selected_component] +
4003  (is_face ? dof_info->row_starts_plain_indices[cells[v]] :
4004  dof_info->row_starts_plain_indices[cell * n_lanes + v]);
4005  next_index_indicators = index_indicators;
4006  }
4007 
4008  if (n_components == 1 || n_fe_components == 1)
4009  {
4010  unsigned int ind_local = 0;
4011  for (; index_indicators != next_index_indicators; ++index_indicators)
4012  {
4013  const std::pair<unsigned short, unsigned short> indicator =
4014  dof_info->constraint_indicator[index_indicators];
4015  // run through values up to next constraint
4016  for (unsigned int j = 0; j < indicator.first; ++j)
4017  for (unsigned int comp = 0; comp < n_components; ++comp)
4018  operation.process_dof(dof_indices[v][j],
4019  *src[comp],
4020  values_dofs[comp][ind_local + j][v]);
4021 
4022  ind_local += indicator.first;
4023  dof_indices[v] += indicator.first;
4024 
4025  // constrained case: build the local value as a linear
4026  // combination of the global value according to constraints
4027  Number value[n_components];
4028  for (unsigned int comp = 0; comp < n_components; ++comp)
4029  operation.pre_constraints(values_dofs[comp][ind_local][v],
4030  value[comp]);
4031 
4032  const Number *data_val =
4033  matrix_info->constraint_pool_begin(indicator.second);
4034  const Number *end_pool =
4035  matrix_info->constraint_pool_end(indicator.second);
4036  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4037  for (unsigned int comp = 0; comp < n_components; ++comp)
4038  operation.process_constraint(*dof_indices[v],
4039  *data_val,
4040  *src[comp],
4041  value[comp]);
4042 
4043  for (unsigned int comp = 0; comp < n_components; ++comp)
4044  operation.post_constraints(value[comp],
4045  values_dofs[comp][ind_local][v]);
4046  ind_local++;
4047  }
4048 
4049  AssertIndexRange(ind_local, dofs_per_component + 1);
4050 
4051  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4052  for (unsigned int comp = 0; comp < n_components; ++comp)
4053  operation.process_dof(*dof_indices[v],
4054  *src[comp],
4055  values_dofs[comp][ind_local][v]);
4056  }
4057  else
4058  {
4059  // case with vector-valued finite elements where all components are
4060  // included in one single vector. Assumption: first come all entries
4061  // to the first component, then all entries to the second one, and
4062  // so on. This is ensured by the way MatrixFree reads out the
4063  // indices.
4064  for (unsigned int comp = 0; comp < n_components; ++comp)
4065  {
4066  unsigned int ind_local = 0;
4067 
4068  // check whether there is any constraint on the current cell
4069  for (; index_indicators != next_index_indicators;
4070  ++index_indicators)
4071  {
4072  const std::pair<unsigned short, unsigned short> indicator =
4073  dof_info->constraint_indicator[index_indicators];
4074 
4075  // run through values up to next constraint
4076  for (unsigned int j = 0; j < indicator.first; ++j)
4077  operation.process_dof(dof_indices[v][j],
4078  *src[0],
4079  values_dofs[comp][ind_local + j][v]);
4080  ind_local += indicator.first;
4081  dof_indices[v] += indicator.first;
4082 
4083  // constrained case: build the local value as a linear
4084  // combination of the global value according to constraints
4085  Number value;
4086  operation.pre_constraints(values_dofs[comp][ind_local][v],
4087  value);
4088 
4089  const Number *data_val =
4090  matrix_info->constraint_pool_begin(indicator.second);
4091  const Number *end_pool =
4092  matrix_info->constraint_pool_end(indicator.second);
4093 
4094  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4095  operation.process_constraint(*dof_indices[v],
4096  *data_val,
4097  *src[0],
4098  value);
4099 
4100  operation.post_constraints(value,
4101  values_dofs[comp][ind_local][v]);
4102  ind_local++;
4103  }
4104 
4105  AssertIndexRange(ind_local, dofs_per_component + 1);
4106 
4107  // get the dof values past the last constraint
4108  for (; ind_local < dofs_per_component;
4109  ++dof_indices[v], ++ind_local)
4110  {
4111  AssertIndexRange(*dof_indices[v], src[0]->size());
4112  operation.process_dof(*dof_indices[v],
4113  *src[0],
4114  values_dofs[comp][ind_local][v]);
4115  }
4116 
4117  if (apply_constraints == true && comp + 1 < n_components)
4118  {
4119  if (is_face)
4120  next_index_indicators =
4121  dof_info
4122  ->row_starts[cells[v] * n_fe_components +
4123  first_selected_component + comp + 2]
4124  .second;
4125  else
4126  next_index_indicators =
4127  dof_info
4128  ->row_starts[(cell * n_lanes + v) * n_fe_components +
4129  first_selected_component + comp + 2]
4130  .second;
4131  }
4132  }
4133  }
4134  }
4135 }
4136 
4137 
4138 
4139 template <int dim,
4140  int n_components_,
4141  typename Number,
4142  bool is_face,
4143  typename VectorizedArrayType>
4144 template <typename VectorType, typename VectorOperation>
4145 inline void
4148  VectorType * src[]) const
4149 {
4150  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4151 
4152  unsigned int index =
4153  first_selected_component * data->dofs_per_component_on_cell;
4154  for (unsigned int comp = 0; comp < n_components; ++comp)
4155  {
4156  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4157  ++i, ++index)
4158  {
4159  operation.process_empty(values_dofs[comp][i]);
4160  operation.process_dof_global(
4161  local_dof_indices[data->lexicographic_numbering[index]],
4162  *src[0],
4163  values_dofs[comp][i][0]);
4164  }
4165  }
4166 }
4167 
4168 
4169 
4170 template <int dim,
4171  int n_components_,
4172  typename Number,
4173  bool is_face,
4174  typename VectorizedArrayType>
4175 template <typename VectorType, typename VectorOperation>
4176 inline void
4179  const VectorOperation & operation,
4180  VectorType * src[],
4181  const std::bitset<VectorizedArrayType::size()> &mask) const
4182 {
4183  // This functions processes the functions read_dof_values,
4184  // distribute_local_to_global, and set_dof_values with the same code for
4185  // contiguous cell indices (DG case). The distinction between these three
4186  // cases is made by the input VectorOperation that either reads values from
4187  // a vector and puts the data into the local data field or write local data
4188  // into the vector. Certain operations are no-ops for the given use case.
4189 
4190  std::integral_constant<bool,
4191  internal::is_vectorizable<VectorType, Number>::value>
4192  vector_selector;
4194  is_face ? dof_access_index :
4196  const unsigned int n_lanes = mask.count();
4197 
4198  const std::vector<unsigned int> &dof_indices_cont =
4199  dof_info->dof_indices_contiguous[ind];
4200 
4201  // Simple case: We have contiguous storage, so we can simply copy out the
4202  // data
4203  if (dof_info->index_storage_variants[ind][cell] ==
4205  interleaved_contiguous &&
4206  n_lanes == VectorizedArrayType::size())
4207  {
4208  const unsigned int dof_index =
4209  dof_indices_cont[cell * VectorizedArrayType::size()] +
4210  dof_info->component_dof_indices_offset[active_fe_index]
4211  [first_selected_component] *
4212  VectorizedArrayType::size();
4213  if (n_components == 1 || n_fe_components == 1)
4214  for (unsigned int comp = 0; comp < n_components; ++comp)
4215  operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
4216  dof_index,
4217  *src[comp],
4218  values_dofs[comp],
4219  vector_selector);
4220  else
4221  operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
4222  n_components,
4223  dof_index,
4224  *src[0],
4225  values_dofs[0],
4226  vector_selector);
4227  return;
4228  }
4229 
4230  // More general case: Must go through the components one by one and apply
4231  // some transformations
4232  const unsigned int n_filled_lanes =
4233  dof_info->n_vectorization_lanes_filled[ind][this->cell];
4234 
4235  unsigned int dof_indices[VectorizedArrayType::size()];
4236  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4237  dof_indices[v] =
4238  dof_indices_cont[cell * VectorizedArrayType::size() + v] +
4239  dof_info->component_dof_indices_offset[active_fe_index]
4240  [first_selected_component] *
4242  [ind][cell * VectorizedArrayType::size() + v];
4243 
4244  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4245  dof_indices[v] = numbers::invalid_unsigned_int;
4246 
4247  // In the case with contiguous cell indices, we know that there are no
4248  // constraints and that the indices within each element are contiguous
4249  if (n_filled_lanes == VectorizedArrayType::size() &&
4250  n_lanes == VectorizedArrayType::size())
4251  {
4252  if (dof_info->index_storage_variants[ind][cell] ==
4254  contiguous)
4255  {
4256  if (n_components == 1 || n_fe_components == 1)
4257  for (unsigned int comp = 0; comp < n_components; ++comp)
4258  operation.process_dofs_vectorized_transpose(
4260  dof_indices,
4261  *src[comp],
4262  values_dofs[comp],
4263  vector_selector);
4264  else
4265  operation.process_dofs_vectorized_transpose(
4266  data->dofs_per_component_on_cell * n_components,
4267  dof_indices,
4268  *src[0],
4269  &values_dofs[0][0],
4270  vector_selector);
4271  }
4272  else if (dof_info->index_storage_variants[ind][cell] ==
4274  interleaved_contiguous_strided)
4275  {
4276  if (n_components == 1 || n_fe_components == 1)
4277  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4278  {
4279  for (unsigned int comp = 0; comp < n_components; ++comp)
4280  operation.process_dof_gather(dof_indices,
4281  *src[comp],
4282  i * VectorizedArrayType::size(),
4283  values_dofs[comp][i],
4284  vector_selector);
4285  }
4286  else
4287  for (unsigned int comp = 0; comp < n_components; ++comp)
4288  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4289  ++i)
4290  {
4291  operation.process_dof_gather(
4292  dof_indices,
4293  *src[0],
4294  (comp * data->dofs_per_component_on_cell + i) *
4295  VectorizedArrayType::size(),
4296  values_dofs[comp][i],
4297  vector_selector);
4298  }
4299  }
4300  else
4301  {
4302  Assert(dof_info->index_storage_variants[ind][cell] ==
4304  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4305  ExcNotImplemented());
4306  const unsigned int *offsets =
4308  [ind][VectorizedArrayType::size() * cell];
4309  if (n_components == 1 || n_fe_components == 1)
4310  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4311  {
4312  for (unsigned int comp = 0; comp < n_components; ++comp)
4313  operation.process_dof_gather(dof_indices,
4314  *src[comp],
4315  0,
4316  values_dofs[comp][i],
4317  vector_selector);
4318  DEAL_II_OPENMP_SIMD_PRAGMA
4319  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4320  dof_indices[v] += offsets[v];
4321  }
4322  else
4323  for (unsigned int comp = 0; comp < n_components; ++comp)
4324  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4325  ++i)
4326  {
4327  operation.process_dof_gather(dof_indices,
4328  *src[0],
4329  0,
4330  values_dofs[comp][i],
4331  vector_selector);
4332  DEAL_II_OPENMP_SIMD_PRAGMA
4333  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4334  dof_indices[v] += offsets[v];
4335  }
4336  }
4337  }
4338  else
4339  for (unsigned int comp = 0; comp < n_components; ++comp)
4340  {
4341  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4342  operation.process_empty(values_dofs[comp][i]);
4343  if (dof_info->index_storage_variants[ind][cell] ==
4345  contiguous)
4346  {
4347  if (n_components == 1 || n_fe_components == 1)
4348  {
4349  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4350  if (mask[v] == true)
4351  for (unsigned int i = 0;
4352  i < data->dofs_per_component_on_cell;
4353  ++i)
4354  operation.process_dof(dof_indices[v] + i,
4355  *src[comp],
4356  values_dofs[comp][i][v]);
4357  }
4358  else
4359  {
4360  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4361  if (mask[v] == true)
4362  for (unsigned int i = 0;
4363  i < data->dofs_per_component_on_cell;
4364  ++i)
4365  operation.process_dof(
4366  dof_indices[v] + i +
4367  comp * data->dofs_per_component_on_cell,
4368  *src[0],
4369  values_dofs[comp][i][v]);
4370  }
4371  }
4372  else
4373  {
4374  const unsigned int *offsets =
4376  [ind][VectorizedArrayType::size() * cell];
4377  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4378  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4379  if (n_components == 1 || n_fe_components == 1)
4380  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4381  {
4382  if (mask[v] == true)
4383  for (unsigned int i = 0;
4384  i < data->dofs_per_component_on_cell;
4385  ++i)
4386  operation.process_dof(dof_indices[v] + i * offsets[v],
4387  *src[comp],
4388  values_dofs[comp][i][v]);
4389  }
4390  else
4391  {
4392  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4393  if (mask[v] == true)
4394  for (unsigned int i = 0;
4395  i < data->dofs_per_component_on_cell;
4396  ++i)
4397  operation.process_dof(
4398  dof_indices[v] +
4399  (i + comp * data->dofs_per_component_on_cell) *
4400  offsets[v],
4401  *src[0],
4402  values_dofs[comp][i][v]);
4403  }
4404  }
4405  }
4406 }
4407 
4408 
4409 
4410 template <int dim,
4411  int n_components_,
4412  typename Number,
4413  bool is_face,
4414  typename VectorizedArrayType>
4415 template <typename VectorType>
4416 inline void
4418  read_dof_values(const VectorType &src, const unsigned int first_index)
4419 {
4420  // select between block vectors and non-block vectors. Note that the number
4421  // of components is checked in the internal data
4422  typename internal::BlockVectorSelector<
4423  VectorType,
4424  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4425  for (unsigned int d = 0; d < n_components; ++d)
4426  src_data[d] =
4427  internal::BlockVectorSelector<VectorType,
4429  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4430 
4431  internal::VectorReader<Number, VectorizedArrayType> reader;
4432  read_write_operation(reader,
4433  src_data,
4434  std::bitset<VectorizedArrayType::size()>().flip(),
4435  true);
4436 
4437 # ifdef DEBUG
4438  dof_values_initialized = true;
4439 # endif
4440 }
4441 
4442 
4443 
4444 template <int dim,
4445  int n_components_,
4446  typename Number,
4447  bool is_face,
4448  typename VectorizedArrayType>
4449 template <typename VectorType>
4450 inline void
4452  read_dof_values_plain(const VectorType &src, const unsigned int first_index)
4453 {
4454  // select between block vectors and non-block vectors. Note that the number
4455  // of components is checked in the internal data
4456  typename internal::BlockVectorSelector<
4457  VectorType,
4458  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4459  for (unsigned int d = 0; d < n_components; ++d)
4460  src_data[d] =
4461  internal::BlockVectorSelector<VectorType,
4463  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4464 
4465  internal::VectorReader<Number, VectorizedArrayType> reader;
4466  read_write_operation(reader,
4467  src_data,
4468  std::bitset<VectorizedArrayType::size()>().flip(),
4469  false);
4470 
4471 # ifdef DEBUG
4472  dof_values_initialized = true;
4473 # endif
4474 }
4475 
4476 
4477 
4478 template <int dim,
4479  int n_components_,
4480  typename Number,
4481  bool is_face,
4482  typename VectorizedArrayType>
4483 template <typename VectorType>
4484 inline void
4487  VectorType & dst,
4488  const unsigned int first_index,
4489  const std::bitset<VectorizedArrayType::size()> &mask) const
4490 {
4491 # ifdef DEBUG
4492  Assert(dof_values_initialized == true,
4494 # endif
4495 
4496  // select between block vectors and non-block vectors. Note that the number
4497  // of components is checked in the internal data
4498  typename internal::BlockVectorSelector<
4499  VectorType,
4500  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4501  for (unsigned int d = 0; d < n_components; ++d)
4502  dst_data[d] = internal::BlockVectorSelector<
4503  VectorType,
4504  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4505  d + first_index);
4506 
4507  internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
4508  distributor;
4509  read_write_operation(distributor, dst_data, mask);
4510 }
4511 
4512 
4513 
4514 template <int dim,
4515  int n_components_,
4516  typename Number,
4517  bool is_face,
4518  typename VectorizedArrayType>
4519 template <typename VectorType>
4520 inline void
4522  set_dof_values(VectorType & dst,
4523  const unsigned int first_index,
4524  const std::bitset<VectorizedArrayType::size()> &mask) const
4525 {
4526 # ifdef DEBUG
4527  Assert(dof_values_initialized == true,
4529 # endif
4530 
4531  // select between block vectors and non-block vectors. Note that the number
4532  // of components is checked in the internal data
4533  typename internal::BlockVectorSelector<
4534  VectorType,
4535  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4536  for (unsigned int d = 0; d < n_components; ++d)
4537  dst_data[d] = internal::BlockVectorSelector<
4538  VectorType,
4539  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4540  d + first_index);
4541 
4542  internal::VectorSetter<Number, VectorizedArrayType> setter;
4543  read_write_operation(setter, dst_data, mask);
4544 }
4545 
4546 
4547 
4548 /*------------------------------ access to data fields ----------------------*/
4549 
4550 template <int dim,
4551  int n_components,
4552  typename Number,
4553  bool is_face,
4554  typename VectorizedArrayType>
4555 inline const std::vector<unsigned int> &
4558 {
4559  return data->lexicographic_numbering;
4560 }
4561 
4562 
4563 
4564 template <int dim,
4565  int n_components,
4566  typename Number,
4567  bool is_face,
4568  typename VectorizedArrayType>
4571  get_scratch_data() const
4572 {
4574  const_cast<VectorizedArrayType *>(scratch_data),
4575  scratch_data_array->end() - scratch_data);
4576 }
4577 
4578 
4579 
4580 template <int dim,
4581  int n_components,
4582  typename Number,
4583  bool is_face,
4584  typename VectorizedArrayType>
4585 inline const VectorizedArrayType *
4587  begin_dof_values() const
4588 {
4589  return &values_dofs[0][0];
4590 }
4591 
4592 
4593 
4594 template <int dim,
4595  int n_components,
4596  typename Number,
4597  bool is_face,
4598  typename VectorizedArrayType>
4599 inline VectorizedArrayType *
4602 {
4603 # ifdef DEBUG
4604  dof_values_initialized = true;
4605 # endif
4606  return &values_dofs[0][0];
4607 }
4608 
4609 
4610 
4611 template <int dim,
4612  int n_components,
4613  typename Number,
4614  bool is_face,
4615  typename VectorizedArrayType>
4616 inline const VectorizedArrayType *
4618  begin_values() const
4619 {
4620 # ifdef DEBUG
4621  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
4622 # endif
4623  return &values_quad[0][0];
4624 }
4625 
4626 
4627 
4628 template <int dim,
4629  int n_components,
4630  typename Number,
4631  bool is_face,
4632  typename VectorizedArrayType>
4633 inline VectorizedArrayType *
4635  begin_values()
4636 {
4637 # ifdef DEBUG
4638  values_quad_initialized = true;
4639  values_quad_submitted = true;
4640 # endif
4641  return &values_quad[0][0];
4642 }
4643 
4644 
4645 
4646 template <int dim,
4647  int n_components,
4648  typename Number,
4649  bool is_face,
4650  typename VectorizedArrayType>
4651 inline const VectorizedArrayType *
4653  begin_gradients() const
4654 {
4655 # ifdef DEBUG
4656  Assert(gradients_quad_initialized || gradients_quad_submitted,
4657  ExcNotInitialized());
4658 # endif
4659  return &gradients_quad[0][0][0];
4660 }
4661 
4662 
4663 
4664 template <int dim,
4665  int n_components,
4666  typename Number,
4667  bool is_face,
4668  typename VectorizedArrayType>
4669 inline VectorizedArrayType *
4672 {
4673 # ifdef DEBUG
4674  gradients_quad_submitted = true;
4675  gradients_quad_initialized = true;
4676 # endif
4677  return &gradients_quad[0][0][0];
4678 }
4679 
4680 
4681 
4682 template <int dim,
4683  int n_components,
4684  typename Number,
4685  bool is_face,
4686  typename VectorizedArrayType>
4687 inline const VectorizedArrayType *
4689  begin_hessians() const
4690 {
4691 # ifdef DEBUG
4692  Assert(hessians_quad_initialized, ExcNotInitialized());
4693 # endif
4694  return &hessians_quad[0][0][0];
4695 }
4696 
4697 
4698 
4699 template <int dim,
4700  int n_components,
4701  typename Number,
4702  bool is_face,
4703  typename VectorizedArrayType>
4704 inline VectorizedArrayType *
4707 {
4708 # ifdef DEBUG
4709  hessians_quad_initialized = true;
4710 # endif
4711  return &hessians_quad[0][0][0];
4712 }
4713 
4714 
4715 
4716 template <int dim,
4717  int n_components_,
4718  typename Number,
4719  bool is_face,
4720  typename VectorizedArrayType>
4721 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType>
4723  get_dof_value(const unsigned int dof) const
4724 {
4725  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4727  for (unsigned int comp = 0; comp < n_components; comp++)
4728  return_value[comp] = this->values_dofs[comp][dof];
4729  return return_value;
4730 }
4731 
4732 
4733 
4734 template <int dim,
4735  int n_components_,
4736  typename Number,
4737  bool is_face,
4738  typename VectorizedArrayType>
4739 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType>
4741  get_value(const unsigned int q_point) const
4742 {
4743 # ifdef DEBUG
4744  Assert(this->values_quad_initialized == true,
4746 # endif
4747  AssertIndexRange(q_point, this->n_quadrature_points);
4749  for (unsigned int comp = 0; comp < n_components; comp++)
4750  return_value[comp] = this->values_quad[comp][q_point];
4751  return return_value;
4752 }
4753 
4754 
4755 
4756 template <int dim,
4757  int n_components_,
4758  typename Number,
4759  bool is_face,
4760  typename VectorizedArrayType>
4761 inline DEAL_II_ALWAYS_INLINE
4764  get_gradient(const unsigned int q_point) const
4765 {
4766 # ifdef DEBUG
4767  Assert(this->gradients_quad_initialized == true,
4769 # endif
4770  AssertIndexRange(q_point, this->n_quadrature_points);
4771 
4772  Assert(jacobian != nullptr, ExcNotInitialized());
4773 
4775 
4776  // Cartesian cell
4777  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4778  {
4779  for (unsigned int comp = 0; comp < n_components; comp++)
4780  for (unsigned int d = 0; d < dim; ++d)
4781  grad_out[comp][d] =
4782  (this->gradients_quad[comp][d][q_point] * jacobian[0][d][d]);
4783  }
4784  // cell with general/affine Jacobian
4785  else
4786  {
4788  jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4789  q_point :
4790  0];
4791  for (unsigned int comp = 0; comp < n_components; comp++)
4792  for (unsigned int d = 0; d < dim; ++d)
4793  {
4794  grad_out[comp][d] =
4795  jac[d][0] * this->gradients_quad[comp][0][q_point];
4796  for (unsigned int e = 1; e < dim; ++e)
4797  grad_out[comp][d] +=
4798  jac[d][e] * this->gradients_quad[comp][e][q_point];
4799  }
4800  }
4801  return grad_out;
4802 }
4803 
4804 
4805 
4806 template <int dim,
4807  int n_components_,
4808  typename Number,
4809  bool is_face,
4810  typename VectorizedArrayType>
4811 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType>
4813  get_normal_derivative(const unsigned int q_point) const
4814 {
4815  AssertIndexRange(q_point, this->n_quadrature_points);
4816 # ifdef DEBUG
4817  Assert(this->gradients_quad_initialized == true,
4819 # endif
4820 
4821  Assert(normal_x_jacobian != nullptr, ExcNotInitialized());
4822 
4824  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4825  for (unsigned int comp = 0; comp < n_components; comp++)
4826  grad_out[comp] = this->gradients_quad[comp][dim - 1][q_point] *
4827  (this->normal_x_jacobian[0][dim - 1]);
4828  else
4829  {
4830  const unsigned int index =
4831  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4832  for (unsigned int comp = 0; comp < n_components; comp++)
4833  {
4834  grad_out[comp] = this->gradients_quad[comp][0][q_point] *
4835  this->normal_x_jacobian[index][0];
4836  for (unsigned int d = 1; d < dim; ++d)
4837  grad_out[comp] += this->gradients_quad[comp][d][q_point] *
4838  this->normal_x_jacobian[index][d];
4839  }
4840  }
4841  return grad_out;
4842 }
4843 
4844 
4845 
4846 namespace internal
4847 {
4848  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4849  // store the lower diagonal because of symmetry
4850  template <typename VectorizedArrayType>
4851  inline void
4852  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4853  const VectorizedArrayType *const hessians_quad[1],
4854  const unsigned int q_point,
4855  VectorizedArrayType (&tmp)[1][1])
4856  {
4857  tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
4858  }
4859 
4860  template <typename VectorizedArrayType>
4861  inline void
4862  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4863  const VectorizedArrayType *const hessians_quad[3],
4864  const unsigned int q_point,
4865  VectorizedArrayType (&tmp)[2][2])
4866  {
4867  for (unsigned int d = 0; d < 2; ++d)
4868  {
4869  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
4870  jac[d][1] * hessians_quad[2][q_point]);
4871  tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
4872  jac[d][1] * hessians_quad[1][q_point]);
4873  }
4874  }
4875 
4876  template <typename VectorizedArrayType>
4877  inline void
4878  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4879  const VectorizedArrayType *const hessians_quad[6],
4880  const unsigned int q_point,
4881  VectorizedArrayType (&tmp)[3][3])
4882  {
4883  for (unsigned int d = 0; d < 3; ++d)
4884  {
4885  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
4886  jac[d][1] * hessians_quad[3][q_point] +
4887  jac[d][2] * hessians_quad[4][q_point]);
4888  tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
4889  jac[d][1] * hessians_quad[1][q_point] +
4890  jac[d][2] * hessians_quad[5][q_point]);
4891  tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
4892  jac[d][1] * hessians_quad[5][q_point] +
4893  jac[d][2] * hessians_quad[2][q_point]);
4894  }
4895  }
4896 } // namespace internal
4897 
4898 
4899 
4900 template <int dim,
4901  int n_components_,
4902  typename Number,
4903  bool is_face,
4904  typename VectorizedArrayType>
4907  get_hessian(const unsigned int q_point) const
4908 {
4909  Assert(!is_face, ExcNotImplemented());
4910 # ifdef DEBUG
4911  Assert(this->hessians_quad_initialized == true,
4913 # endif
4914  AssertIndexRange(q_point, this->n_quadrature_points);
4915 
4916  Assert(jacobian != nullptr, ExcNotImplemented());
4918  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4919  0 :
4920  q_point];
4921 
4922  Tensor<2, dim, VectorizedArrayType> hessian_out[n_components];
4923 
4924  // Cartesian cell
4925  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4926  {
4927  for (unsigned int comp = 0; comp < n_components; comp++)
4928  for (unsigned int d = 0; d < dim; ++d)
4929  {
4930  hessian_out[comp][d][d] =
4931  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
4932  switch (dim)
4933  {
4934  case 1:
4935  break;
4936  case 2:
4937  hessian_out[comp][0][1] =
4938  (this->hessians_quad[comp][2][q_point] * jac[0][0] *
4939  jac[1][1]);
4940  break;
4941  case 3:
4942  hessian_out[comp][0][1] =
4943  (this->hessians_quad[comp][3][q_point] * jac[0][0] *
4944  jac[1][1]);
4945  hessian_out[comp][0][2] =
4946  (this->hessians_quad[comp][4][q_point] * jac[0][0] *
4947  jac[2][2]);
4948  hessian_out[comp][1][2] =
4949  (this->hessians_quad[comp][5][q_point] * jac[1][1] *
4950  jac[2][2]);
4951  break;
4952  default:
4953  Assert(false, ExcNotImplemented());
4954  }
4955  for (unsigned int e = d + 1; e < dim; ++e)
4956  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4957  }
4958  }
4959  // cell with general Jacobian, but constant within the cell
4960  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4961  {
4962  for (unsigned int comp = 0; comp < n_components; comp++)
4963  {
4964  // compute laplacian before the gradient because it needs to access
4965  // unscaled gradient data
4966  VectorizedArrayType tmp[dim][dim];
4967  internal::hessian_unit_times_jac(jac,
4968  this->hessians_quad[comp],
4969  q_point,
4970  tmp);
4971 
4972  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4973  for (unsigned int d = 0; d < dim; ++d)
4974  for (unsigned int e = d; e < dim; ++e)
4975  {
4976  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4977  for (unsigned int f = 1; f < dim; ++f)
4978  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4979  }
4980 
4981  // no J' * grad(u) part here because the Jacobian is constant
4982  // throughout the cell and hence, its derivative is zero
4983 
4984  // take symmetric part
4985  for (unsigned int d = 0; d < dim; ++d)
4986  for (unsigned int e = d + 1; e < dim; ++e)
4987  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4988  }
4989  }
4990  // cell with general Jacobian
4991  else
4992  {
4993  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
4994  &jac_grad =
4995  mapping_data->jacobian_gradients
4996  [1 - this->is_interior_face]
4997  [this->mapping_data->data_index_offsets[this->cell] + q_point];
4998  for (unsigned int comp = 0; comp < n_components; comp++)
4999  {
5000  // compute laplacian before the gradient because it needs to access
5001  // unscaled gradient data
5002  VectorizedArrayType tmp[dim][dim];
5003  internal::hessian_unit_times_jac(jac,
5004  this->hessians_quad[comp],
5005  q_point,
5006  tmp);
5007 
5008  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5009  for (unsigned int d = 0; d < dim; ++d)
5010  for (unsigned int e = d; e < dim; ++e)
5011  {
5012  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5013  for (unsigned int f = 1; f < dim; ++f)
5014  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5015  }
5016 
5017  // add diagonal part of J' * grad(u)
5018  for (unsigned int d = 0; d < dim; ++d)
5019  for (unsigned int e = 0; e < dim; ++e)
5020  hessian_out[comp][d][d] +=
5021  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5022 
5023  // add off-diagonal part of J' * grad(u)
5024  for (unsigned int d = 0, count = dim; d < dim; ++d)
5025  for (unsigned int e = d + 1; e < dim; ++e, ++count)
5026  for (unsigned int f = 0; f < dim; ++f)
5027  hessian_out[comp][d][e] +=
5028  (jac_grad[count][f] * this->gradients_quad[comp][f][q_point]);
5029 
5030  // take symmetric part
5031  for (unsigned int d = 0; d < dim; ++d)
5032  for (unsigned int e = d + 1; e < dim; ++e)
5033  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5034  }
5035  }
5037  hessian_out);
5038 }
5039 
5040 
5041 
5042 template <int dim,
5043  int n_components_,
5044  typename Number,
5045  bool is_face,
5046  typename VectorizedArrayType>
5049  get_hessian_diagonal(const unsigned int q_point) const
5050 {
5051  Assert(!is_face, ExcNotImplemented());
5052 # ifdef DEBUG
5053  Assert(this->hessians_quad_initialized == true,
5055 # endif
5056  AssertIndexRange(q_point, this->n_quadrature_points);
5057 
5058  Assert(jacobian != nullptr, ExcNotImplemented());
5060  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5061  0 :
5062  q_point];
5063 
5065 
5066  // Cartesian cell
5067  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5068  {
5069  for (unsigned int comp = 0; comp < n_components; comp++)
5070  for (unsigned int d = 0; d < dim; ++d)
5071  hessian_out[comp][d] =
5072  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
5073  }
5074  // cell with general Jacobian, but constant within the cell
5075  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5076  {
5077  for (unsigned int comp = 0; comp < n_components; comp++)
5078  {
5079  // compute laplacian before the gradient because it needs to access
5080  // unscaled gradient data
5081  VectorizedArrayType tmp[dim][dim];
5082  internal::hessian_unit_times_jac(jac,
5083  this->hessians_quad[comp],
5084  q_point,
5085  tmp);
5086 
5087  // compute only the trace part of hessian, J * tmp = J *
5088  // hess_unit(u) * J^T
5089  for (unsigned int d = 0; d < dim; ++d)
5090  {
5091  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5092  for (unsigned int f = 1; f < dim; ++f)
5093  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5094  }
5095  }
5096  }
5097  // cell with general Jacobian
5098  else
5099  {
5100  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
5101  &jac_grad =
5102  mapping_data->jacobian_gradients
5103  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5104  for (unsigned int comp = 0; comp < n_components; comp++)
5105  {
5106  // compute laplacian before the gradient because it needs to access
5107  // unscaled gradient data
5108  VectorizedArrayType tmp[dim][dim];
5109  internal::hessian_unit_times_jac(jac,
5110  this->hessians_quad[comp],
5111  q_point,
5112  tmp);
5113 
5114  // compute only the trace part of hessian, J * tmp = J *
5115  // hess_unit(u) * J^T
5116  for (unsigned int d = 0; d < dim; ++d)
5117  {
5118  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5119  for (unsigned int f = 1; f < dim; ++f)
5120  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5121  }
5122 
5123  for (unsigned int d = 0; d < dim; ++d)
5124  for (unsigned int e = 0; e < dim; ++e)
5125  hessian_out[comp][d] +=
5126  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5127  }
5128  }
5129  return hessian_out;
5130 }
5131 
5132 
5133 
5134 template <int dim,
5135  int n_components_,
5136  typename Number,
5137  bool is_face,
5138  typename VectorizedArrayType>
5141  get_laplacian(const unsigned int q_point) const
5142 {
5143  Assert(is_face == false, ExcNotImplemented());
5144 # ifdef DEBUG
5145  Assert(this->hessians_quad_initialized == true,
5147 # endif
5148  AssertIndexRange(q_point, this->n_quadrature_points);
5149 
5152  hess_diag = get_hessian_diagonal(q_point);
5153  for (unsigned int comp = 0; comp < n_components; ++comp)
5154  {
5155  laplacian_out[comp] = hess_diag[comp][0];
5156  for (unsigned int d = 1; d < dim; ++d)
5157  laplacian_out[comp] += hess_diag[comp][d];
5158  }
5159  return laplacian_out;
5160 }
5161 
5162 
5163 
5164 template <int dim,
5165  int n_components_,
5166  typename Number,
5167  bool is_face,
5168  typename VectorizedArrayType>
5169 inline DEAL_II_ALWAYS_INLINE void
5172  const unsigned int dof)
5173 {
5174 # ifdef DEBUG
5175  this->dof_values_initialized = true;
5176 # endif
5177  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5178  for (unsigned int comp = 0; comp < n_components; comp++)
5179  this->values_dofs[comp][dof] = val_in[comp];
5180 }
5181 
5182 
5183 
5184 template <int dim,
5185  int n_components_,
5186  typename Number,
5187  bool is_face,
5188  typename VectorizedArrayType>
5189 inline DEAL_II_ALWAYS_INLINE void
5192  const unsigned int q_point)
5193 {
5195  AssertIndexRange(q_point, this->n_quadrature_points);
5196  Assert(this->J_value != nullptr, ExcNotInitialized());
5197 # ifdef DEBUG
5198  this->values_quad_submitted = true;
5199 # endif
5200 
5201  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5202  {
5203  const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
5204  for (unsigned int comp = 0; comp < n_components; ++comp)
5205  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5206  }
5207  else
5208  {
5209  const VectorizedArrayType JxW = J_value[q_point];
5210  for (unsigned int comp = 0; comp < n_components; ++comp)
5211  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5212  }
5213 }
5214 
5215 
5216 
5217 template <int dim,
5218  int n_components_,
5219  typename Number,
5220  bool is_face,
5221  typename VectorizedArrayType>
5222 inline DEAL_II_ALWAYS_INLINE void
5225  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5226  const unsigned int q_point)
5227 {
5229  AssertIndexRange(q_point, this->n_quadrature_points);
5230  Assert(this->J_value != nullptr, ExcNotInitialized());
5231  Assert(this->jacobian != nullptr, ExcNotInitialized());
5232 # ifdef DEBUG
5233  this->gradients_quad_submitted = true;
5234 # endif
5235 
5236  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5237  {
5238  const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
5239  for (unsigned int comp = 0; comp < n_components; comp++)
5240  for (unsigned int d = 0; d < dim; ++d)
5241  this->gradients_quad[comp][d][q_point] =
5242  (grad_in[comp][d] * jacobian[0][d][d] * JxW);
5243  }
5244  else
5245  {
5247  this->cell_type > internal::MatrixFreeFunctions::affine ?
5248  jacobian[q_point] :
5249  jacobian[0];
5250  const VectorizedArrayType JxW =
5251  this->cell_type > internal::MatrixFreeFunctions::affine ?
5252  J_value[q_point] :
5253  J_value[0] * quadrature_weights[q_point];
5254  for (unsigned int comp = 0; comp < n_components; ++comp)
5255  for (unsigned int d = 0; d < dim; ++d)
5256  {
5257  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5258  for (unsigned int e = 1; e < dim; ++e)
5259  new_val += (jac[e][d] * grad_in[comp][e]);
5260  this->gradients_quad[comp][d][q_point] = new_val * JxW;
5261  }
5262  }
5263 }
5264 
5265 
5266 
5267 template <int dim,
5268  int n_components_,
5269  typename Number,
5270  bool is_face,
5271  typename VectorizedArrayType>
5272 inline DEAL_II_ALWAYS_INLINE void
5276  const unsigned int q_point)
5277 {
5278  AssertIndexRange(q_point, this->n_quadrature_points);
5279  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
5280 # ifdef DEBUG
5281  this->gradients_quad_submitted = true;
5282 # endif
5283 
5284  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5285  for (unsigned int comp = 0; comp < n_components; comp++)
5286  {
5287  for (unsigned int d = 0; d < dim - 1; ++d)
5288  this->gradients_quad[comp][d][q_point] = VectorizedArrayType();
5289  this->gradients_quad[comp][dim - 1][q_point] =
5290  grad_in[comp] *
5291  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5292  this->quadrature_weights[q_point]);
5293  }
5294  else
5295  {
5296  const unsigned int index =
5297  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5298  for (unsigned int comp = 0; comp < n_components; comp++)
5299  {
5300  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5301  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5302  factor = factor * this->quadrature_weights[q_point];
5303  for (unsigned int d = 0; d < dim; ++d)
5304  this->gradients_quad[comp][d][q_point] =
5305  factor * this->normal_x_jacobian[index][d];
5306  }
5307  }
5308 }
5309 
5310 
5311 
5312 template <int dim,
5313  int n_components_,
5314  typename Number,
5315  bool is_face,
5316  typename VectorizedArrayType>
5319  integrate_value() const
5320 {
5322 # ifdef DEBUG
5323  Assert(this->values_quad_submitted == true,
5325 # endif
5327  for (unsigned int comp = 0; comp < n_components; ++comp)
5328  return_value[comp] = this->values_quad[comp][0];
5329  const unsigned int n_q_points = this->n_quadrature_points;
5330  for (unsigned int q = 1; q < n_q_points; ++q)
5331  for (unsigned int comp = 0; comp < n_components; ++comp)
5332  return_value[comp] += this->values_quad[comp][q];
5333  return (return_value);
5334 }
5335 
5336 
5337 
5338 /*----------------------- FEEvaluationAccess --------------------------------*/
5339 
5340 
5341 template <int dim,
5342  int n_components_,
5343  typename Number,
5344  bool is_face,
5345  typename VectorizedArrayType>
5346 inline FEEvaluationAccess<dim,
5347  n_components_,
5348  Number,
5349  is_face,
5350  VectorizedArrayType>::
5351  FEEvaluationAccess(
5353  const unsigned int dof_no,
5354  const unsigned int first_selected_component,
5355  const unsigned int quad_no_in,
5356  const unsigned int fe_degree,
5357  const unsigned int n_q_points,
5358  const bool is_interior_face)
5360  data_in,
5361  dof_no,
5362  first_selected_component,
5363  quad_no_in,
5364  fe_degree,
5365  n_q_points,
5366  is_interior_face)
5367 {}
5368 
5369 
5370 
5371 template <int dim,
5372  int n_components_,
5373  typename Number,
5374  bool is_face,
5375  typename VectorizedArrayType>
5376 template <int n_components_other>
5377 inline FEEvaluationAccess<dim,
5378  n_components_,
5379  Number,
5380  is_face,
5381  VectorizedArrayType>::
5382  FEEvaluationAccess(const Mapping<dim> & mapping,
5383  const FiniteElement<dim> &fe,
5384  const Quadrature<1> & quadrature,
5385  const UpdateFlags update_flags,
5386  const unsigned int first_selected_component,
5387  const FEEvaluationBase<dim,
5388  n_components_other,
5389  Number,
5390  is_face,
5391  VectorizedArrayType> *other)
5393  mapping,
5394  fe,
5395  quadrature,
5396  update_flags,
5397  first_selected_component,
5398  other)
5399 {}
5400 
5401 
5402 
5403 template <int dim,
5404  int n_components_,
5405  typename Number,
5406  bool is_face,
5407  typename VectorizedArrayType>
5408 inline FEEvaluationAccess<dim,
5409  n_components_,
5410  Number,
5411  is_face,
5412  VectorizedArrayType>::
5413  FEEvaluationAccess(const FEEvaluationAccess<dim,
5414  n_components_,
5415  Number,
5416  is_face,
5417  VectorizedArrayType> &other)
5419  other)
5420 {}
5421 
5422 
5423 
5424 template <int dim,
5425  int n_components_,
5426  typename Number,
5427  bool is_face,
5428  typename VectorizedArrayType>
5429 inline FEEvaluationAccess<dim,
5430  n_components_,
5431  Number,
5432  is_face,
5433  VectorizedArrayType> &
5435 operator=(const FEEvaluationAccess<dim,
5436  n_components_,
5437  Number,
5438  is_face,
5439  VectorizedArrayType> &other)
5440 {
5441  this->FEEvaluationBase<dim,
5442  n_components_,
5443  Number,
5444  is_face,
5445  VectorizedArrayType>::operator=(other);
5446  return *this;
5447 }
5448 
5449 
5450 
5451 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5452 
5453 
5454 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5458  const unsigned int dof_no,
5459  const unsigned int first_selected_component,
5460  const unsigned int quad_no_in,
5461  const unsigned int fe_degree,
5462  const unsigned int n_q_points,
5463  const bool is_interior_face)
5465  data_in,
5466  dof_no,
5467  first_selected_component,
5468  quad_no_in,
5469  fe_degree,
5470  n_q_points,
5471  is_interior_face)
5472 {}
5473 
5474 
5475 
5476 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5477 template <int n_components_other>
5479  FEEvaluationAccess(const Mapping<dim> & mapping,
5480  const FiniteElement<dim> &fe,
5481  const Quadrature<1> & quadrature,
5482  const UpdateFlags update_flags,
5483  const unsigned int first_selected_component,
5484  const FEEvaluationBase<dim,
5485  n_components_other,
5486  Number,
5487  is_face,
5488  VectorizedArrayType> *other)
5490  mapping,
5491  fe,
5492  quadrature,
5493  update_flags,
5494  first_selected_component,
5495  other)
5496 {}
5497 
5498 
5499 
5500 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5504  &other)
5506 {}
5507 
5508 
5509 
5510 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5514 {
5516  operator=(other);
5517  return *this;
5518 }
5519 
5520 
5521 
5522 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5523 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5525  const unsigned int dof) const
5526 {
5527  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5528  return this->values_dofs[0][dof];
5529 }
5530 
5531 
5532 
5533 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5534 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5536  const unsigned int q_point) const
5537 {
5538 # ifdef DEBUG
5539  Assert(this->values_quad_initialized == true,
5541 # endif
5542  AssertIndexRange(q_point, this->n_quadrature_points);
5543  return this->values_quad[0][q_point];
5544 }
5545 
5546 
5547 
5548 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5549 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5551  get_normal_derivative(const unsigned int q_point) const
5552 {
5553  return BaseClass::get_normal_derivative(q_point)[0];
5554 }
5555 
5556 
5557 
5558 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5559 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
5561  const unsigned int q_point) const
5562 {
5563  // could use the base class gradient, but that involves too many expensive
5564  // initialization operations on tensors
5565 
5566 # ifdef DEBUG
5567  Assert(this->gradients_quad_initialized == true,
5569 # endif
5570  AssertIndexRange(q_point, this->n_quadrature_points);
5571 
5572  Assert(this->jacobian != nullptr, ExcNotInitialized());
5573 
5575 
5576  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5577  {
5578  for (unsigned int d = 0; d < dim; ++d)
5579  grad_out[d] =
5580  (this->gradients_quad[0][d][q_point] * this->jacobian[0][d][d]);
5581  }
5582  // cell with general/affine Jacobian
5583  else
5584  {
5586  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5587  q_point :
5588  0];
5589  for (unsigned int d = 0; d < dim; ++d)
5590  {
5591  grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point];
5592  for (unsigned int e = 1; e < dim; ++e)
5593  grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point];
5594  }
5595  }
5596  return grad_out;
5597 }
5598 
5599 
5600 
5601 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5604  const unsigned int q_point) const
5605 {
5606  return BaseClass::get_hessian(q_point)[0];
5607 }
5608 
5609 
5610 
5611 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5614  get_hessian_diagonal(const unsigned int q_point) const
5615 {
5616  return BaseClass::get_hessian_diagonal(q_point)[0];
5617 }
5618 
5619 
5620 
5621 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5622 inline VectorizedArrayType
5624  const unsigned int q_point) const
5625 {
5626  return BaseClass::get_laplacian(q_point)[0];
5627 }
5628 
5629 
5630 
5631 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5632 inline void DEAL_II_ALWAYS_INLINE
5634  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5635 {
5636 # ifdef DEBUG
5637  this->dof_values_initialized = true;
5638  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5639 # endif
5640  this->values_dofs[0][dof] = val_in;
5641 }
5642 
5643 
5644 
5645 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5646 inline void DEAL_II_ALWAYS_INLINE
5648  const VectorizedArrayType val_in,
5649  const unsigned int q_index)
5650 {
5652  AssertIndexRange(q_index, this->n_quadrature_points);
5653  Assert(this->J_value != nullptr, ExcNotInitialized());
5654 # ifdef DEBUG
5655  this->values_quad_submitted = true;
5656 # endif
5657 
5658  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5659  {
5660  const VectorizedArrayType JxW =
5661  this->J_value[0] * this->quadrature_weights[q_index];
5662  this->values_quad[0][q_index] = val_in * JxW;
5663  }
5664  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5665  {
5666  this->values_quad[0][q_index] = val_in * this->J_value[q_index];
5667  }
5668 }
5669 
5670 
5671 
5672 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5673 inline DEAL_II_ALWAYS_INLINE void
5675  const Tensor<1, 1, VectorizedArrayType> val_in,
5676  const unsigned int q_point)
5677 {
5678  submit_value(val_in[0], q_point);
5679 }
5680 
5681 
5682 
5683 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5684 inline DEAL_II_ALWAYS_INLINE void
5686  submit_normal_derivative(const VectorizedArrayType grad_in,
5687  const unsigned int q_point)
5688 {
5690  grad[0] = grad_in;
5691  BaseClass::submit_normal_derivative(grad, q_point);
5692 }
5693 
5694 
5695 
5696 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5697 inline DEAL_II_ALWAYS_INLINE void
5700  const unsigned int q_index)
5701 {
5703  AssertIndexRange(q_index, this->n_quadrature_points);
5704  Assert(this->J_value != nullptr, ExcNotInitialized());
5705  Assert(this->jacobian != nullptr, ExcNotInitialized());
5706 # ifdef DEBUG
5707  this->gradients_quad_submitted = true;
5708 # endif
5709 
5710  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5711  {
5712  const VectorizedArrayType JxW =
5713  this->J_value[0] * this->quadrature_weights[q_index];
5714  for (unsigned int d = 0; d < dim; ++d)
5715  this->gradients_quad[0][d][q_index] =
5716  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5717  }
5718  // general/affine cell type
5719  else
5720  {
5722  this->cell_type > internal::MatrixFreeFunctions::affine ?
5723  this->jacobian[q_index] :
5724  this->jacobian[0];
5725  const VectorizedArrayType JxW =
5726  this->cell_type > internal::MatrixFreeFunctions::affine ?
5727  this->J_value[q_index] :
5728  this->J_value[0] * this->quadrature_weights[q_index];
5729  for (unsigned int d = 0; d < dim; ++d)
5730  {
5731  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5732  for (unsigned int e = 1; e < dim; ++e)
5733  new_val += jac[e][d] * grad_in[e];
5734  this->gradients_quad[0][d][q_index] = new_val * JxW;
5735  }
5736  }
5737 }
5738 
5739 
5740 
5741 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5742 inline VectorizedArrayType
5744  integrate_value() const
5745 {
5746  return BaseClass::integrate_value()[0];
5747 }
5748 
5749 
5750 
5751 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5752 
5753 
5754 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5758  const unsigned int dof_no,
5759  const unsigned int first_selected_component,
5760  const unsigned int quad_no_in,
5761  const unsigned int fe_degree,
5762  const unsigned int n_q_points,
5763  const bool is_interior_face)
5765  data_in,
5766  dof_no,
5767  first_selected_component,
5768  quad_no_in,
5769  fe_degree,
5770  n_q_points,
5771  is_interior_face)
5772 {}
5773 
5774 
5775 
5776 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5777 template <int n_components_other>
5779  FEEvaluationAccess(const Mapping<dim> & mapping,
5780  const FiniteElement<dim> &fe,
5781  const Quadrature<1> & quadrature,
5782  const UpdateFlags update_flags,
5783  const unsigned int first_selected_component,
5784  const FEEvaluationBase<dim,
5785  n_components_other,
5786  Number,
5787  is_face,
5788  VectorizedArrayType> *other)
5790  mapping,
5791  fe,
5792  quadrature,
5793  update_flags,
5794  first_selected_component,
5795  other)
5796 {}
5797 
5798 
5799 
5800 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5804  &other)
5806 {}
5807 
5808 
5809 
5810 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5814  &other)
5815 {
5817  operator=(other);
5818  return *this;
5819 }
5820 
5821 
5822 
5823 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5824 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
5826  get_gradient(const unsigned int q_point) const
5827 {
5828  return BaseClass::get_gradient(q_point);
5829 }
5830 
5831 
5832 
5833 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5834 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5836  get_divergence(const unsigned int q_point) const
5837 {
5838 # ifdef DEBUG
5839  Assert(this->gradients_quad_initialized == true,
5841 # endif
5842  AssertIndexRange(q_point, this->n_quadrature_points);
5843  Assert(this->jacobian != nullptr, ExcNotInitialized());
5844 
5845  VectorizedArrayType divergence;
5846 
5847  // Cartesian cell
5848  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5849  {
5850  divergence =
5851  (this->gradients_quad[0][0][q_point] * this->jacobian[0][0][0]);
5852  for (unsigned int d = 1; d < dim; ++d)
5853  divergence +=
5854  (this->gradients_quad[d][d][q_point] * this->jacobian[0][d][d]);
5855  }
5856  // cell with general/constant Jacobian
5857  else
5858  {
5860  this->cell_type == internal::MatrixFreeFunctions::general ?
5861  this->jacobian[q_point] :
5862  this->jacobian[0];
5863  divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
5864  for (unsigned int e = 1; e < dim; ++e)
5865  divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
5866  for (unsigned int d = 1; d < dim; ++d)
5867  for (unsigned int e = 0; e < dim; ++e)
5868  divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
5869  }
5870  return divergence;
5871 }
5872 
5873 
5874 
5875 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5876 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArrayType>
5878  get_symmetric_gradient(const unsigned int q_point) const
5879 {
5880  // copy from generic function into dim-specialization function
5881  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5882  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5883  VectorizedArrayType half = make_vectorized_array<Number>(0.5);
5884  for (unsigned int d = 0; d < dim; ++d)
5885  symmetrized[d] = grad[d][d];
5886  switch (dim)
5887  {
5888  case 1:
5889  break;
5890  case 2:
5891  symmetrized[2] = grad[0][1] + grad[1][0];
5892  symmetrized[2] *= half;
5893  break;
5894  case 3:
5895  symmetrized[3] = grad[0][1] + grad[1][0];
5896  symmetrized[3] *= half;
5897  symmetrized[4] = grad[0][2] + grad[2][0];
5898  symmetrized[4] *= half;
5899  symmetrized[5] = grad[1][2] + grad[2][1];
5900  symmetrized[5] *= half;
5901  break;
5902  default:
5903  Assert(false, ExcNotImplemented());
5904  }
5906 }
5907 
5908 
5909 
5910 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5911 inline DEAL_II_ALWAYS_INLINE
5914  const unsigned int q_point) const
5915 {
5916  // copy from generic function into dim-specialization function
5917  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5919  switch (dim)
5920  {
5921  case 1:
5922  Assert(false,
5923  ExcMessage(
5924  "Computing the curl in 1d is not a useful operation"));
5925  break;
5926  case 2:
5927  curl[0] = grad[1][0] - grad[0][1];
5928  break;
5929  case 3:
5930  curl[0] = grad[2][1] - grad[1][2];
5931  curl[1] = grad[0][2] - grad[2][0];
5932  curl[2] = grad[1][0] - grad[0][1];
5933  break;
5934  default:
5935  Assert(false, ExcNotImplemented());
5936  }
5937  return curl;
5938 }
5939 
5940 
5941 
5942 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5943 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
5945  get_hessian_diagonal(const unsigned int q_point) const
5946 {
5947  return BaseClass::get_hessian_diagonal(q_point);
5948 }
5949 
5950 
5951 
5952 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5953 inline DEAL_II_ALWAYS_INLINE Tensor<3, dim, VectorizedArrayType>
5955  const unsigned int q_point) const
5956 {
5957 # ifdef DEBUG
5958  Assert(this->hessians_quad_initialized == true,
5960 # endif
5961  AssertIndexRange(q_point, this->n_quadrature_points);
5962  return BaseClass::get_hessian(q_point);
5963 }
5964 
5965 
5966 
5967 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5968 inline DEAL_II_ALWAYS_INLINE void
5971  const unsigned int q_point)
5972 {
5973  BaseClass::submit_gradient(grad_in, q_point);
5974 }
5975 
5976 
5977 
5978 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5979 inline DEAL_II_ALWAYS_INLINE void
5982  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
5983  const unsigned int q_point)
5984 {
5985  BaseClass::submit_gradient(grad_in, q_point);
5986 }
5987 
5988 
5989 
5990 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5991 inline DEAL_II_ALWAYS_INLINE void
5993  submit_divergence(const VectorizedArrayType div_in,
5994  const unsigned int q_point)
5995 {
5997  AssertIndexRange(q_point, this->n_quadrature_points);
5998  Assert(this->J_value != nullptr, ExcNotInitialized());
5999  Assert(this->jacobian != nullptr, ExcNotInitialized());
6000 # ifdef DEBUG
6001  this->gradients_quad_submitted = true;
6002 # endif
6003 
6004  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6005  {
6006  const VectorizedArrayType fac =
6007  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6008  for (unsigned int d = 0; d < dim; ++d)
6009  {
6010  this->gradients_quad[d][d][q_point] = (fac * this->jacobian[0][d][d]);
6011  for (unsigned int e = d + 1; e < dim; ++e)
6012  {
6013  this->gradients_quad[d][e][q_point] = VectorizedArrayType();
6014  this->gradients_quad[e][d][q_point] = VectorizedArrayType();
6015  }
6016  }
6017  }
6018  else
6019  {
6021  this->cell_type == internal::MatrixFreeFunctions::general ?
6022  this->jacobian[q_point] :
6023  this->jacobian[0];
6024  const VectorizedArrayType fac =
6025  (this->cell_type == internal::MatrixFreeFunctions::general ?
6026  this->J_value[q_point] :
6027  this->J_value[0] * this->quadrature_weights[q_point]) *
6028  div_in;
6029  for (unsigned int d = 0; d < dim; ++d)
6030  {
6031  for (unsigned int e = 0; e < dim; ++e)
6032  this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
6033  }
6034  }
6035 }
6036 
6037 
6038 
6039 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6040 inline DEAL_II_ALWAYS_INLINE void
6044  const unsigned int q_point)
6045 {
6046  // could have used base class operator, but that involves some overhead
6047  // which is inefficient. it is nice to have the symmetric tensor because
6048  // that saves some operations
6050  AssertIndexRange(q_point, this->n_quadrature_points);
6051  Assert(this->J_value != nullptr, ExcNotInitialized());
6052  Assert(this->jacobian != nullptr, ExcNotInitialized());
6053 # ifdef DEBUG
6054  this->gradients_quad_submitted = true;
6055 # endif
6056 
6057  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6058  {
6059  const VectorizedArrayType JxW =
6060  this->J_value[0] * this->quadrature_weights[q_point];
6061  for (unsigned int d = 0; d < dim; ++d)
6062  this->gradients_quad[d][d][q_point] =
6063  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6064  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6065  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6066  {
6067  const VectorizedArrayType value =
6068  sym_grad.access_raw_entry(counter) * JxW;
6069  this->gradients_quad[e][d][q_point] =
6070  (value * this->jacobian[0][d][d]);
6071  this->gradients_quad[d][e][q_point] =
6072  (value * this->jacobian[0][e][e]);
6073  }
6074  }
6075  // general/affine cell type
6076  else
6077  {
6078  const VectorizedArrayType JxW =
6079  this->cell_type == internal::MatrixFreeFunctions::general ?
6080  this->J_value[q_point] :
6081  this->J_value[0] * this->quadrature_weights[q_point];
6083  this->cell_type == internal::MatrixFreeFunctions::general ?
6084  this->jacobian[q_point] :
6085  this->jacobian[0];
6086  VectorizedArrayType weighted[dim][dim];
6087  for (unsigned int i = 0; i < dim; ++i)
6088  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6089  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6090  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6091  {
6092  const VectorizedArrayType value =
6093  sym_grad.access_raw_entry(counter) * JxW;
6094  weighted[i][j] = value;
6095  weighted[j][i] = value;
6096  }
6097  for (unsigned int comp = 0; comp < dim; ++comp)
6098  for (unsigned int d = 0; d < dim; ++d)
6099  {
6100  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6101  for (unsigned int e = 1; e < dim; ++e)
6102  new_val += jac[e][d] * weighted[comp][e];
6103  this->gradients_quad[comp][d][q_point] = new_val;
6104  }
6105  }
6106 }
6107 
6108 
6109 
6110 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6111 inline DEAL_II_ALWAYS_INLINE void
6114  const unsigned int q_point)
6115 {
6117  switch (dim)
6118  {
6119  case 1:
6120  Assert(false,
6121  ExcMessage(
6122  "Testing by the curl in 1d is not a useful operation"));
6123  break;
6124  case 2:
6125  grad[1][0] = curl[0];
6126  grad[0][1] = -curl[0];
6127  break;
6128  case 3:
6129  grad[2][1] = curl[0];
6130  grad[1][2] = -curl[0];
6131  grad[0][2] = curl[1];
6132  grad[2][0] = -curl[1];
6133  grad[1][0] = curl[2];
6134  grad[0][1] = -curl[2];
6135  break;
6136  default:
6137  Assert(false, ExcNotImplemented());
6138  }
6139  submit_gradient(grad, q_point);
6140 }
6141 
6142 
6143 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6144 
6145 
6146 template <typename Number, bool is_face, typename VectorizedArrayType>
6149  const unsigned int dof_no,
6150  const unsigned int first_selected_component,
6151  const unsigned int quad_no_in,
6152  const unsigned int fe_degree,
6153  const unsigned int n_q_points,
6154  const bool is_interior_face)
6156  data_in,
6157  dof_no,
6158  first_selected_component,
6159  quad_no_in,
6160  fe_degree,
6161  n_q_points,
6162  is_interior_face)
6163 {}
6164 
6165 
6166 
6167 template <typename Number, bool is_face, typename VectorizedArrayType>
6168 template <int n_components_other>
6170  FEEvaluationAccess(const Mapping<1> & mapping,
6171  const FiniteElement<1> &fe,
6172  const Quadrature<1> & quadrature,
6173  const UpdateFlags update_flags,
6174  const unsigned int first_selected_component,
6175  const FEEvaluationBase<1,
6176  n_components_other,
6177  Number,
6178  is_face,
6179  VectorizedArrayType> *other)
6181  mapping,
6182  fe,
6183  quadrature,
6184  update_flags,
6185  first_selected_component,
6186  other)
6187 {}
6188 
6189 
6190 
6191 template <typename Number, bool is_face, typename VectorizedArrayType>
6196 {}
6197 
6198 
6199 
6200 template <typename Number, bool is_face, typename VectorizedArrayType>
6204 {
6206  other);
6207  return *this;
6208 }
6209 
6210 
6211 
6212 template <typename Number, bool is_face, typename VectorizedArrayType>
6213 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6215  const unsigned int dof) const
6216 {
6217  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6218  return this->values_dofs[0][dof];
6219 }
6220 
6221 
6222 
6223 template <typename Number, bool is_face, typename VectorizedArrayType>
6224 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6226  const unsigned int q_point) const
6227 {
6228 # ifdef DEBUG
6229  Assert(this->values_quad_initialized == true,
6231 # endif
6232  AssertIndexRange(q_point, this->n_quadrature_points);
6233  return this->values_quad[0][q_point];
6234 }
6235 
6236 
6237 
6238 template <typename Number, bool is_face, typename VectorizedArrayType>
6239 inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArrayType>
6241  const unsigned int q_point) const
6242 {
6243  // could use the base class gradient, but that involves too many inefficient
6244  // initialization operations on tensors
6245 
6246 # ifdef DEBUG
6247  Assert(this->gradients_quad_initialized == true,
6249 # endif
6250  AssertIndexRange(q_point, this->n_quadrature_points);
6251 
6253  this->cell_type == internal::MatrixFreeFunctions::general ?
6254  this->jacobian[q_point] :
6255  this->jacobian[0];
6256 
6258  grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];
6259 
6260  return grad_out;
6261 }
6262 
6263 
6264 
6265 template <typename Number, bool is_face, typename VectorizedArrayType>
6266 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6268  get_normal_derivative(const unsigned int q_point) const
6269 {
6270  return BaseClass::get_normal_derivative(q_point)[0];
6271 }
6272 
6273 
6274 
6275 template <typename Number, bool is_face, typename VectorizedArrayType>
6276 inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, VectorizedArrayType>
6278  const unsigned int q_point) const
6279 {
6280  return BaseClass::get_hessian(q_point)[0];
6281 }
6282 
6283 
6284 
6285 template <typename Number, bool is_face, typename VectorizedArrayType>
6286 inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArrayType>
6288  get_hessian_diagonal(const unsigned int q_point) const
6289 {
6290  return BaseClass::get_hessian_diagonal(q_point)[0];
6291 }
6292 
6293 
6294 
6295 template <typename Number, bool is_face, typename VectorizedArrayType>
6296 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6298  const unsigned int q_point) const
6299 {
6300  return BaseClass::get_laplacian(q_point)[0];
6301 }
6302 
6303 
6304 
6305 template <typename Number, bool is_face, typename VectorizedArrayType>
6306 inline DEAL_II_ALWAYS_INLINE void DEAL_II_ALWAYS_INLINE
6308  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6309 {
6310 # ifdef DEBUG
6311  this->dof_values_initialized = true;
6312  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6313 # endif
6314  this->values_dofs[0][dof] = val_in;
6315 }
6316 
6317 
6318 
6319 template <typename Number, bool is_face, typename VectorizedArrayType>
6320 inline DEAL_II_ALWAYS_INLINE void
6322  const VectorizedArrayType val_in,
6323  const unsigned int q_point)
6324 {
6326  AssertIndexRange(q_point, this->n_quadrature_points);
6327 # ifdef DEBUG
6328  this->values_quad_submitted = true;
6329 # endif
6330 
6331  if (this->cell_type == internal::MatrixFreeFunctions::general)
6332  {
6333  const VectorizedArrayType JxW = this->J_value[q_point];
6334  this->values_quad[0][q_point] = val_in * JxW;
6335  }
6336  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6337  {
6338  const VectorizedArrayType JxW =
6339  this->J_value[0] * this->quadrature_weights[q_point];
6340  this->values_quad[0][q_point] = val_in * JxW;
6341  }
6342 }
6343 
6344 
6345 
6346 template <typename Number, bool is_face, typename VectorizedArrayType>
6347 inline DEAL_II_ALWAYS_INLINE void
6349  const Tensor<1, 1, VectorizedArrayType> val_in,
6350  const unsigned int q_point)
6351 {
6352  submit_value(val_in[0], q_point);
6353 }
6354 
6355 
6356 
6357 template <typename Number, bool is_face, typename VectorizedArrayType>
6358 inline DEAL_II_ALWAYS_INLINE void
6360  const Tensor<1, 1, VectorizedArrayType> grad_in,
6361  const unsigned int q_point)
6362 {
6363  submit_gradient(grad_in[0], q_point);
6364 }
6365 
6366 
6367 
6368 template <typename Number, bool is_face, typename VectorizedArrayType>
6369 inline DEAL_II_ALWAYS_INLINE void
6371  const VectorizedArrayType grad_in,
6372  const unsigned int q_point)
6373 {
6375  AssertIndexRange(q_point, this->n_quadrature_points);
6376 # ifdef DEBUG
6377  this->gradients_quad_submitted = true;
6378 # endif
6379 
6381  this->cell_type == internal::MatrixFreeFunctions::general ?
6382  this->jacobian[q_point] :
6383  this->jacobian[0];
6384  const VectorizedArrayType JxW =
6385  this->cell_type == internal::MatrixFreeFunctions::general ?
6386  this->J_value[q_point] :
6387  this->J_value[0] * this->quadrature_weights[q_point];
6388 
6389  this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in * JxW;
6390 }
6391 
6392 
6393 
6394 template <typename Number, bool is_face, typename VectorizedArrayType>
6395 inline DEAL_II_ALWAYS_INLINE void
6397  submit_normal_derivative(const VectorizedArrayType grad_in,
6398  const unsigned int q_point)
6399 {
6401  grad[0] = grad_in;
6402  BaseClass::submit_normal_derivative(grad, q_point);
6403 }
6404 
6405 
6406 
6407 template <typename Number, bool is_face, typename VectorizedArrayType>
6408 inline DEAL_II_ALWAYS_INLINE void
6411  const unsigned int q_point)
6412 {
6413  BaseClass::submit_normal_derivative(grad_in, q_point);
6414 }
6415 
6416 
6417 
6418 template <typename Number, bool is_face, typename VectorizedArrayType>
6419 inline VectorizedArrayType
6421  integrate_value() const
6422 {
6423  return BaseClass::integrate_value()[0];
6424 }
6425 
6426 
6427 
6428 /*-------------------------- FEEvaluation -----------------------------------*/
6429 
6430 
6431 template <int dim,
6432  int fe_degree,
6433  int n_q_points_1d,
6434  int n_components_,
6435  typename Number,
6436  typename VectorizedArrayType>
6437 inline FEEvaluation<dim,
6438  fe_degree,
6439  n_q_points_1d,
6440  n_components_,
6441  Number,
6442  VectorizedArrayType>::
6443  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
6444  const unsigned int fe_no,
6445  const unsigned int quad_no,
6446  const unsigned int first_selected_component)
6447  : BaseClass(data_in,
6448  fe_no,
6449  first_selected_component,
6450  quad_no,
6451  fe_degree,
6452  static_n_q_points)
6453  , dofs_per_component(this->data->dofs_per_component_on_cell)
6454  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6455  , n_q_points(this->data->n_q_points)
6456 {
6457  check_template_arguments(fe_no, 0);
6458 }
6459 
6460 
6461 
6462 template <int dim,
6463  int fe_degree,
6464  int n_q_points_1d,
6465  int n_components_,
6466  typename Number,
6467  typename VectorizedArrayType>
6468 inline FEEvaluation<dim,
6469  fe_degree,
6470  n_q_points_1d,
6471  n_components_,
6472  Number,
6473  VectorizedArrayType>::
6474  FEEvaluation(const Mapping<dim> & mapping,
6475  const FiniteElement<dim> &fe,
6476  const Quadrature<1> & quadrature,
6477  const UpdateFlags update_flags,
6478  const unsigned int first_selected_component)
6479  : BaseClass(mapping,
6480  fe,
6481  quadrature,
6482  update_flags,
6483  first_selected_component,
6484  static_cast<
6486  nullptr))
6487  , dofs_per_component(this->data->dofs_per_component_on_cell)
6488  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6489  , n_q_points(this->data->n_q_points)
6490 {
6491  check_template_arguments(numbers::invalid_unsigned_int, 0);
6492 }
6493 
6494 
6495 
6496 template <int dim,
6497  int fe_degree,
6498  int n_q_points_1d,
6499  int n_components_,
6500  typename Number,
6501  typename VectorizedArrayType>
6502 inline FEEvaluation<dim,
6503  fe_degree,
6504  n_q_points_1d,
6505  n_components_,
6506  Number,
6507  VectorizedArrayType>::
6508  FEEvaluation(const FiniteElement<dim> &fe,
6509  const Quadrature<1> & quadrature,
6510  const UpdateFlags update_flags,
6511  const unsigned int first_selected_component)
6513  fe,
6514  quadrature,
6515  update_flags,
6516  first_selected_component,
6517  static_cast<
6519  nullptr))
6520  , dofs_per_component(this->data->dofs_per_component_on_cell)
6521  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6522  , n_q_points(this->data->n_q_points)
6523 {
6524  check_template_arguments(numbers::invalid_unsigned_int, 0);
6525 }
6526 
6527 
6528 
6529 template <int dim,
6530  int fe_degree,
6531  int n_q_points_1d,
6532  int n_components_,
6533  typename Number,
6534  typename VectorizedArrayType>
6535 template <int n_components_other>
6536 inline FEEvaluation<dim,
6537  fe_degree,
6538  n_q_points_1d,
6539  n_components_,
6540  Number,
6541  VectorizedArrayType>::
6542  FEEvaluation(const FiniteElement<dim> & fe,
6543  const FEEvaluationBase<dim,
6544  n_components_other,
6545  Number,
6546  false,
6547  VectorizedArrayType> &other,
6548  const unsigned int first_selected_component)
6549  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6550  fe,
6551  other.mapped_geometry->get_quadrature(),
6552  other.mapped_geometry->get_fe_values().get_update_flags(),
6553  first_selected_component,
6554  &other)
6555  , dofs_per_component(this->data->dofs_per_component_on_cell)
6556  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6557  , n_q_points(this->data->n_q_points)
6558 {
6559  check_template_arguments(numbers::invalid_unsigned_int, 0);
6560 }
6561 
6562 
6563 
6564 template <int dim,
6565  int fe_degree,
6566  int n_q_points_1d,
6567  int n_components_,
6568  typename Number,
6569  typename VectorizedArrayType>
6570 inline FEEvaluation<dim,
6571  fe_degree,
6572  n_q_points_1d,
6573  n_components_,
6574  Number,
6575  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6576  &other)
6577  : BaseClass(other)
6578  , dofs_per_component(this->data->dofs_per_component_on_cell)
6579  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6580  , n_q_points(this->data->n_q_points)
6581 {
6582  check_template_arguments(numbers::invalid_unsigned_int, 0);
6583 }
6584 
6585 
6586 
6587 template <int dim,
6588  int fe_degree,
6589  int n_q_points_1d,
6590  int n_components_,
6591  typename Number,
6592  typename VectorizedArrayType>
6593 inline FEEvaluation<dim,
6594  fe_degree,
6595  n_q_points_1d,
6596  n_components_,
6597  Number,
6598  VectorizedArrayType> &
6599 FEEvaluation<dim,
6600  fe_degree,
6601  n_q_points_1d,
6602  n_components_,
6603  Number,
6604  VectorizedArrayType>::operator=(const FEEvaluation &other)
6605 {
6606  BaseClass::operator=(other);
6607  check_template_arguments(numbers::invalid_unsigned_int, 0);
6608  return *this;
6609 }
6610 
6611 
6612 
6613 template <int dim,
6614  int fe_degree,
6615  int n_q_points_1d,
6616  int n_components_,
6617  typename Number,
6618  typename VectorizedArrayType>
6619 inline void
6620 FEEvaluation<dim,
6621  fe_degree,
6622  n_q_points_1d,
6623  n_components_,
6624  Number,
6625  VectorizedArrayType>::
6626  check_template_arguments(const unsigned int dof_no,
6627  const unsigned int first_selected_component)
6628 {
6629  (void)dof_no;
6630  (void)first_selected_component;
6631 
6632 # ifdef DEBUG
6633  // print error message when the dimensions do not match. Propose a possible
6634  // fix
6635  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6636  static_cast<unsigned int>(fe_degree) !=
6637  this->data->data.front().fe_degree) ||
6638  n_q_points != this->n_quadrature_points)
6639  {
6640  std::string message =
6641  "-------------------------------------------------------\n";
6642  message += "Illegal arguments in constructor/wrong template arguments!\n";
6643  message += " Called --> FEEvaluation<dim,";
6644  message += Utilities::int_to_string(fe_degree) + ",";
6645  message += Utilities::int_to_string(n_q_points_1d);
6646  message += "," + Utilities::int_to_string(n_components);
6647  message += ",Number>(data";
6648  if (first_selected_component != numbers::invalid_unsigned_int)
6649  {
6650  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6651  message += Utilities::int_to_string(this->quad_no) + ", ";
6652  message += Utilities::int_to_string(first_selected_component);
6653  }
6654  message += ")\n";
6655 
6656  // check whether some other vector component has the correct number of
6657  // points
6658  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6659  proposed_fe_comp = numbers::invalid_unsigned_int,
6660  proposed_quad_comp = numbers::invalid_unsigned_int;
6661  if (dof_no != numbers::invalid_unsigned_int)
6662  {
6663  if (static_cast<unsigned int>(fe_degree) ==
6664  this->data->data.front().fe_degree)
6665  {
6666  proposed_dof_comp = dof_no;
6667  proposed_fe_comp = first_selected_component;
6668  }
6669  else
6670  for (unsigned int no = 0; no < this->matrix_info->n_components();
6671  ++no)
6672  for (unsigned int nf = 0;
6673  nf < this->matrix_info->n_base_elements(no);
6674  ++nf)
6675  if (this->matrix_info
6676  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6677  .data.front()
6678  .fe_degree == static_cast<unsigned int>(fe_degree))
6679  {
6680  proposed_dof_comp = no;
6681  proposed_fe_comp = nf;
6682  break;
6683  }
6684  if (n_q_points ==
6685  this->mapping_data->descriptor[this->active_quad_index]
6686  .n_q_points)
6687  proposed_quad_comp = this->quad_no;
6688  else
6689  for (unsigned int no = 0;
6690  no < this->matrix_info->get_mapping_info().cell_data.size();
6691  ++no)
6692  if (this->matrix_info->get_mapping_info()
6693  .cell_data[no]
6694  .descriptor[this->active_quad_index]
6695  .n_q_points == n_q_points)
6696  {
6697  proposed_quad_comp = no;
6698  break;
6699  }
6700  }
6701  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6702  proposed_quad_comp != numbers::invalid_unsigned_int)
6703  {
6704  if (proposed_dof_comp != first_selected_component)
6705  message += "Wrong vector component selection:\n";
6706  else
6707  message += "Wrong quadrature formula selection:\n";
6708  message += " Did you mean FEEvaluation<dim,";
6709  message += Utilities::int_to_string(fe_degree) + ",";
6710  message += Utilities::int_to_string(n_q_points_1d);
6711  message += "," + Utilities::int_to_string(n_components);
6712  message += ",Number>(data";
6713  if (dof_no != numbers::invalid_unsigned_int)
6714  {
6715  message +=
6716  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6717  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
6718  message += Utilities::int_to_string(proposed_fe_comp);
6719  }
6720  message += ")?\n";
6721  std::string correct_pos;
6722  if (proposed_dof_comp != dof_no)
6723  correct_pos = " ^ ";
6724  else
6725  correct_pos = " ";
6726  if (proposed_quad_comp != this->quad_no)
6727  correct_pos += " ^ ";
6728  else
6729  correct_pos += " ";
6730  if (proposed_fe_comp != first_selected_component)
6731  correct_pos += " ^\n";
6732  else
6733  correct_pos += " \n";
6734  message += " " +
6735  correct_pos;
6736  }
6737  // ok, did not find the numbers specified by the template arguments in
6738  // the given list. Suggest correct template arguments
6739  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6740  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6741  message += "Wrong template arguments:\n";
6742  message += " Did you mean FEEvaluation<dim,";
6743  message +=
6744  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6745  message += Utilities::int_to_string(proposed_n_q_points_1d);
6746  message += "," + Utilities::int_to_string(n_components);
6747  message += ",Number>(data";
6748  if (dof_no != numbers::invalid_unsigned_int)
6749  {
6750  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6751  message += Utilities::int_to_string(this->quad_no);
6752  message += ", " + Utilities::int_to_string(first_selected_component);
6753  }
6754  message += ")?\n";
6755  std::string correct_pos;
6756  if (this->data->data.front().fe_degree !=
6757  static_cast<unsigned int>(fe_degree))
6758  correct_pos = " ^";
6759  else
6760  correct_pos = " ";
6761  if (proposed_n_q_points_1d != n_q_points_1d)
6762  correct_pos += " ^\n";
6763  else
6764  correct_pos += " \n";
6765  message += " " + correct_pos;
6766 
6767  Assert(static_cast<unsigned int>(fe_degree) ==
6768  this->data->data.front().fe_degree &&
6769  n_q_points == this->n_quadrature_points,
6770  ExcMessage(message));
6771  }
6772  if (dof_no != numbers::invalid_unsigned_int)
6774  n_q_points,
6775  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6776 # endif
6777 }
6778 
6779 
6780 
6781 template <int dim,
6782  int fe_degree,
6783  int n_q_points_1d,
6784  int n_components_,
6785  typename Number,
6786  typename VectorizedArrayType>
6787 inline void
6788 FEEvaluation<dim,
6789  fe_degree,
6790  n_q_points_1d,
6791  n_components_,
6792  Number,
6793  VectorizedArrayType>::reinit(const unsigned int cell_index)
6794 {
6795  Assert(this->mapped_geometry == nullptr,
6796  ExcMessage("FEEvaluation was initialized without a matrix-free object."
6797  " Integer indexing is not possible"));
6798  if (this->mapped_geometry != nullptr)
6799  return;
6800 
6801  Assert(this->dof_info != nullptr, ExcNotInitialized());
6802  Assert(this->mapping_data != nullptr, ExcNotInitialized());
6803  this->cell = cell_index;
6804  this->cell_type =
6805  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
6806 
6807  const unsigned int offsets =
6808  this->mapping_data->data_index_offsets[cell_index];
6809  this->jacobian = &this->mapping_data->jacobians[0][offsets];
6810  this->J_value = &this->mapping_data->JxW_values[offsets];
6811 
6812 # ifdef DEBUG
6813  this->dof_values_initialized = false;
6814  this->values_quad_initialized = false;
6815  this->gradients_quad_initialized = false;
6816  this->hessians_quad_initialized = false;
6817 # endif
6818 }
6819 
6820 
6821 
6822 template <int dim,
6823  int fe_degree,
6824  int n_q_points_1d,
6825  int n_components_,
6826  typename Number,
6827  typename VectorizedArrayType>
6828 template <typename DoFHandlerType, bool level_dof_access>
6829 inline void
6830 FEEvaluation<dim,
6831  fe_degree,
6832  n_q_points_1d,
6833  n_components_,
6834  Number,
6835  VectorizedArrayType>::
6836  reinit(
6838 {
6839  Assert(this->matrix_info == nullptr,
6840  ExcMessage("Cannot use initialization from cell iterator if "
6841  "initialized from MatrixFree object. Use variant for "
6842  "on the fly computation with arguments as for FEValues "
6843  "instead"));
6844  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6845  this->mapped_geometry->reinit(
6846  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6847  this->local_dof_indices.resize(cell->get_fe().dofs_per_cell);
6848  if (level_dof_access)
6849  cell->get_mg_dof_indices(this->local_dof_indices);
6850  else
6851  cell->get_dof_indices(this->local_dof_indices);
6852 }
6853 
6854 
6855 
6856 template <int dim,
6857  int fe_degree,
6858  int n_q_points_1d,
6859  int n_components_,
6860  typename Number,
6861  typename VectorizedArrayType>
6862 inline void
6863 FEEvaluation<dim,
6864  fe_degree,
6865  n_q_points_1d,
6866  n_components_,
6867  Number,
6868  VectorizedArrayType>::
6869  reinit(const typename Triangulation<dim>::cell_iterator &cell)
6870 {
6871  Assert(this->matrix_info == 0,
6872  ExcMessage("Cannot use initialization from cell iterator if "
6873  "initialized from MatrixFree object. Use variant for "
6874  "on the fly computation with arguments as for FEValues "
6875  "instead"));
6876  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6877  this->mapped_geometry->reinit(cell);
6878 }
6879 
6880 
6881 
6882 template <int dim,
6883  int fe_degree,
6884  int n_q_points_1d,
6885  int n_components_,
6886  typename Number,
6887  typename VectorizedArrayType>
6889 FEEvaluation<dim,
6890  fe_degree,
6891  n_q_points_1d,
6892  n_components_,
6893  Number,
6894  VectorizedArrayType>::quadrature_point(const unsigned int q) const
6895 {
6896  if (this->matrix_info == nullptr)
6897  {
6898  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
6900  ExcNotInitialized());
6901  }
6902  else
6903  {
6904  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
6905  ExcNotInitialized());
6906  }
6907 
6908  AssertIndexRange(q, n_q_points);
6909 
6910  const Point<dim, VectorizedArrayType> *quadrature_points =
6911  &this->mapping_data->quadrature_points
6912  [this->mapping_data->quadrature_point_offsets[this->cell]];
6913 
6914  // Cartesian/affine mesh: only first vertex of cell is stored, we must
6915  // compute it through the Jacobian (which is stored in non-inverted and
6916  // non-transposed form as index '1' in the jacobian field)
6917  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6918  {
6919  Assert(this->jacobian != nullptr, ExcNotInitialized());
6920  Point<dim, VectorizedArrayType> point = quadrature_points[0];
6921 
6922  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6923  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
6924  for (unsigned int d = 0; d < dim; ++d)
6925  point[d] += jac[d][d] *
6926  static_cast<Number>(
6927  this->mapping_data->descriptor[this->active_quad_index]
6928  .quadrature.point(q)[d]);
6929  else
6930  for (unsigned int d = 0; d < dim; ++d)
6931  for (unsigned int e = 0; e < dim; ++e)
6932  point[d] += jac[d][e] * static_cast<Number>(
6933  this->mapping_data
6934  ->descriptor[this->active_quad_index]
6935  .quadrature.point(q)[e]);
6936  return point;
6937  }
6938  else
6939  return quadrature_points[q];
6940 }
6941 
6942 
6943 
6944 template <int dim,
6945  int fe_degree,
6946  int n_q_points_1d,
6947  int n_components_,
6948  typename Number,
6949  typename VectorizedArrayType>
6950 inline void
6951 FEEvaluation<dim,
6952  fe_degree,
6953  n_q_points_1d,
6954  n_components_,
6955  Number,
6956  VectorizedArrayType>::evaluate(const bool evaluate_values,
6957  const bool evaluate_gradients,
6958  const bool evaluate_hessians)
6959 {
6960 # ifdef DEBUG
6961  Assert(this->dof_values_initialized == true,
6963 # endif
6964  evaluate(this->values_dofs[0],
6965  evaluate_values,
6966  evaluate_gradients,
6967  evaluate_hessians);
6968 }
6969 
6970 
6971 
6972 template <int dim,
6973  int fe_degree,
6974  int n_q_points_1d,
6975  int n_components_,
6976  typename Number,
6977  typename VectorizedArrayType>
6978 inline void
6979 FEEvaluation<dim,
6980  fe_degree,
6981  n_q_points_1d,
6982  n_components_,
6983  Number,
6984  VectorizedArrayType>::evaluate(const VectorizedArrayType
6985  * values_array,
6986  const bool evaluate_values,
6987  const bool evaluate_gradients,
6988  const bool evaluate_hessians)
6989 {
6991  dim,
6992  fe_degree,
6993  n_q_points_1d,
6994  n_components,
6995  VectorizedArrayType>::evaluate(*this->data,
6996  const_cast<VectorizedArrayType *>(
6997  values_array),
6998  this->values_quad[0],
6999  this->gradients_quad[0][0],
7000  this->hessians_quad[0][0],
7001  this->scratch_data,
7002  evaluate_values,
7003  evaluate_gradients,
7004  evaluate_hessians);
7005 
7006 # ifdef DEBUG
7007  if (evaluate_values == true)
7008  this->values_quad_initialized = true;
7009  if (evaluate_gradients == true)
7010  this->gradients_quad_initialized = true;
7011  if (evaluate_hessians == true)
7012  this->hessians_quad_initialized = true;
7013 # endif
7014 }
7015 
7016 
7017 
7018 template <int dim,
7019  int fe_degree,
7020  int n_q_points_1d,
7021  int n_components_,
7022  typename Number,
7023  typename VectorizedArrayType>
7024 template <typename VectorType>
7025 inline void
7026 FEEvaluation<
7027  dim,
7028  fe_degree,
7029  n_q_points_1d,
7030  n_components_,
7031  Number,
7032  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7033  const bool evaluate_values,
7034  const bool evaluate_gradients,
7035  const bool evaluate_hessians)
7036 {
7037  // If the index storage is interleaved and contiguous and the vector storage
7038  // has the correct alignment, we can directly pass the pointer into the
7039  // vector to the evaluate() call, without reading the vector entries into a
7040  // separate data field. This saves some operations.
7041  if (std::is_same<typename VectorType::value_type, Number>::value &&
7042  this->dof_info->index_storage_variants
7044  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7045  IndexStorageVariants::interleaved_contiguous &&
7046  reinterpret_cast<std::size_t>(
7047  input_vector.begin() +
7048  this->dof_info->dof_indices_contiguous
7050  [this->cell * VectorizedArrayType::size()]) %
7051  sizeof(VectorizedArrayType) ==
7052  0)
7053  {
7054  const VectorizedArrayType *vec_values =
7055  reinterpret_cast<const VectorizedArrayType *>(
7056  input_vector.begin() +
7057  this->dof_info->dof_indices_contiguous
7059  [this->cell * VectorizedArrayType::size()] +
7060  this->dof_info
7061  ->component_dof_indices_offset[this->active_fe_index]
7062  [this->first_selected_component] *
7063  VectorizedArrayType::size());
7064 
7065  evaluate(vec_values,
7066  evaluate_values,
7067  evaluate_gradients,
7068  evaluate_hessians);
7069  }
7070  else
7071  {
7072  this->read_dof_values(input_vector);
7073  evaluate(this->begin_dof_values(),
7074  evaluate_values,
7075  evaluate_gradients,
7076  evaluate_hessians);
7077  }
7078 }
7079 
7080 
7081 
7082 template <int dim,
7083  int fe_degree,
7084  int n_q_points_1d,
7085  int n_components_,
7086  typename Number,
7087  typename VectorizedArrayType>
7088 inline void
7089 FEEvaluation<dim,
7090  fe_degree,
7091  n_q_points_1d,
7092  n_components_,
7093  Number,
7094  VectorizedArrayType>::integrate(const bool integrate_values,
7095  const bool integrate_gradients)
7096 {
7097  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7098 
7099 # ifdef DEBUG
7100  this->dof_values_initialized = true;
7101 # endif
7102 }
7103 
7104 
7105 
7106 template <int dim,
7107  int fe_degree,
7108  int n_q_points_1d,
7109  int n_components_,
7110  typename Number,
7111  typename VectorizedArrayType>
7112 inline void
7113 FEEvaluation<dim,
7114  fe_degree,
7115  n_q_points_1d,
7116  n_components_,
7117  Number,
7118  VectorizedArrayType>::integrate(const bool integrate_values,
7119  const bool integrate_gradients,
7120  VectorizedArrayType *values_array)
7121 {
7122 # ifdef DEBUG
7123  if (integrate_values == true)
7124  Assert(this->values_quad_submitted == true,
7126  if (integrate_gradients == true)
7127  Assert(this->gradients_quad_submitted == true,
7129 # endif
7130  Assert(this->matrix_info != nullptr ||
7131  this->mapped_geometry->is_initialized(),
7132  ExcNotInitialized());
7133 
7134  SelectEvaluator<dim,
7135  fe_degree,
7136  n_q_points_1d,
7137  n_components,
7138  VectorizedArrayType>::integrate(*this->data,
7139  values_array,
7140  this->values_quad[0],
7141  this->gradients_quad[0][0],
7142  this->scratch_data,
7143  integrate_values,
7144  integrate_gradients,
7145  false);
7146 
7147 # ifdef DEBUG
7148  this->dof_values_initialized = true;
7149 # endif
7150 }
7151 
7152 
7153 
7154 template <int dim,
7155  int fe_degree,
7156  int n_q_points_1d,
7157  int n_components_,
7158  typename Number,
7159  typename VectorizedArrayType>
7160 template <typename VectorType>
7161 inline void
7162 FEEvaluation<
7163  dim,
7164  fe_degree,
7165  n_q_points_1d,
7166  n_components_,
7167  Number,
7168  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7169  const bool integrate_gradients,
7170  VectorType &destination)
7171 {
7172  // If the index storage is interleaved and contiguous and the vector storage
7173  // has the correct alignment, we can directly pass the pointer into the
7174  // vector to the integrate() call, without writing temporary results into a
7175  // separate data field that will later be added into the vector. This saves
7176  // some operations.
7177  if (std::is_same<typename VectorType::value_type, Number>::value &&
7178  this->dof_info->index_storage_variants
7180  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7181  IndexStorageVariants::interleaved_contiguous &&
7182  reinterpret_cast<std::size_t>(
7183  destination.begin() +
7184  this->dof_info->dof_indices_contiguous
7186  [this->cell * VectorizedArrayType::size()]) %
7187  sizeof(VectorizedArrayType) ==
7188  0)
7189  {
7190  VectorizedArrayType *vec_values = reinterpret_cast<VectorizedArrayType *>(
7191  destination.begin() +
7192  this->dof_info->dof_indices_contiguous
7194  [this->cell * VectorizedArrayType::size()] +
7195  this->dof_info
7196  ->component_dof_indices_offset[this->active_fe_index]
7197  [this->first_selected_component] *
7198  VectorizedArrayType::size());
7199  SelectEvaluator<dim,
7200  fe_degree,
7201  n_q_points_1d,
7202  n_components,
7203  VectorizedArrayType>::integrate(*this->data,
7204  vec_values,
7205  this->values_quad[0],
7206  this
7207  ->gradients_quad[0][0],
7208  this->scratch_data,
7209  integrate_values,
7210  integrate_gradients,
7211  true);
7212  }
7213  else
7214  {
7215  integrate(integrate_values,
7216  integrate_gradients,
7217  this->begin_dof_values());
7218  this->distribute_local_to_global(destination);
7219  }
7220 }
7221 
7222 
7223 
7224 /*-------------------------- FEFaceEvaluation ---------------------------*/
7225 
7226 
7227 
7228 template <int dim,
7229  int fe_degree,
7230  int n_q_points_1d,
7231  int n_components_,
7232  typename Number,
7233  typename VectorizedArrayType>
7234 inline FEFaceEvaluation<dim,
7235  fe_degree,
7236  n_q_points_1d,
7237  n_components_,
7238  Number,
7239  VectorizedArrayType>::
7240  FEFaceEvaluation(
7242  const bool is_interior_face,
7243  const unsigned int dof_no,
7244  const unsigned int quad_no,
7245  const unsigned int first_selected_component)
7246  : BaseClass(matrix_free,
7247  dof_no,
7248  first_selected_component,
7249  quad_no,
7250  fe_degree,
7251  static_n_q_points,
7252  is_interior_face)
7253  , dofs_per_component(this->data->dofs_per_component_on_cell)
7254  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7255  , n_q_points(this->data->n_q_points_face)
7256 {}
7257 
7258 
7259 
7260 template <int dim,
7261  int fe_degree,
7262  int n_q_points_1d,
7263  int n_components_,
7264  typename Number,
7265  typename VectorizedArrayType>
7266 inline void
7267 FEFaceEvaluation<dim,
7268  fe_degree,
7269  n_q_points_1d,
7270  n_components_,
7271  Number,
7272  VectorizedArrayType>::reinit(const unsigned int face_index)
7273 {
7274  Assert(this->mapped_geometry == nullptr,
7275  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7276  " Integer indexing is not possible"));
7277  if (this->mapped_geometry != nullptr)
7278  return;
7279 
7280  this->cell = face_index;
7281  this->dof_access_index =
7282  this->is_interior_face ?
7285  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7287  VectorizedArrayType::size()> &faces =
7288  this->matrix_info->get_face_info(face_index);
7289  if (face_index >=
7290  this->matrix_info->get_task_info().face_partition_data.back() &&
7291  face_index <
7292  this->matrix_info->get_task_info().boundary_partition_data.back())
7293  Assert(this->is_interior_face,
7294  ExcMessage("Boundary faces do not have a neighbor"));
7295 
7296  this->face_no =
7297  (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
7298  this->subface_index = faces.subface_index;
7299  if (this->is_interior_face == true)
7300  {
7301  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7302  if (faces.face_orientation > 8)
7303  this->face_orientation = faces.face_orientation - 8;
7304  else
7305  this->face_orientation = 0;
7306  }
7307  else
7308  {
7309  if (faces.face_orientation < 8)
7310  this->face_orientation = faces.face_orientation;
7311  else
7312  this->face_orientation = 0;
7313  }
7314 
7315  this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
7316  const unsigned int offsets =
7317  this->mapping_data->data_index_offsets[face_index];
7318  this->J_value = &this->mapping_data->JxW_values[offsets];
7319  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7320  this->jacobian =
7321  &this->mapping_data->jacobians[!this->is_interior_face][offsets];
7322  this->normal_x_jacobian =
7323  &this->mapping_data
7324  ->normals_times_jacobians[!this->is_interior_face][offsets];
7325 
7326 # ifdef DEBUG
7327  this->dof_values_initialized = false;
7328  this->values_quad_initialized = false;
7329  this->gradients_quad_initialized = false;
7330  this->hessians_quad_initialized = false;
7331 # endif
7332 }
7333 
7334 
7335 
7336 template <int dim,
7337  int fe_degree,
7338  int n_q_points_1d,
7339  int n_components_,
7340  typename Number,
7341  typename VectorizedArrayType>
7342 inline void
7343 FEFaceEvaluation<dim,
7344  fe_degree,
7345  n_q_points_1d,
7346  n_components_,
7347  Number,
7348  VectorizedArrayType>::reinit(const unsigned int cell_index,
7349  const unsigned int face_number)
7350 {
7351  Assert(
7352  this->quad_no <
7353  this->matrix_info->get_mapping_info().face_data_by_cells.size(),
7354  ExcMessage(
7355  "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7357  AssertIndexRange(cell_index,
7358  this->matrix_info->get_mapping_info().cell_type.size());
7359  Assert(this->mapped_geometry == nullptr,
7360  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7361  " Integer indexing is not possible"));
7362  if (this->mapped_geometry != nullptr)
7363  return;
7364  Assert(this->matrix_info != nullptr, ExcNotInitialized());
7365 
7366  this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
7367  this->cell = cell_index;
7368  this->face_orientation = 0;
7369  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7370  this->face_no = face_number;
7371  this->dof_access_index =
7373 
7374  const unsigned int offsets =
7375  this->matrix_info->get_mapping_info()
7376  .face_data_by_cells[this->quad_no]
7377  .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
7378  face_number];
7379  AssertIndexRange(offsets,
7380  this->matrix_info->get_mapping_info()
7381  .face_data_by_cells[this->quad_no]
7382  .JxW_values.size());
7383  this->J_value = &this->matrix_info->get_mapping_info()
7384  .face_data_by_cells[this->quad_no]
7385  .JxW_values[offsets];
7386  this->normal_vectors = &this->matrix_info->get_mapping_info()
7387  .face_data_by_cells[this->quad_no]
7388  .normal_vectors[offsets];
7389  this->jacobian = &this->matrix_info->get_mapping_info()
7390  .face_data_by_cells[this->quad_no]
7391  .jacobians[!this->is_interior_face][offsets];
7392  this->normal_x_jacobian =
7393  &this->matrix_info->get_mapping_info()
7394  .face_data_by_cells[this->quad_no]
7395  .normals_times_jacobians[!this->is_interior_face][offsets];
7396 
7397 # ifdef DEBUG
7398  this->dof_values_initialized = false;
7399  this->values_quad_initialized = false;
7400  this->gradients_quad_initialized = false;
7401  this->hessians_quad_initialized = false;
7402 # endif
7403 }
7404 
7405 
7406 
7407 template <int dim,
7408  int fe_degree,
7409  int n_q_points_1d,
7410  int n_components,
7411  typename Number,
7412  typename VectorizedArrayType>
7413 inline void
7414 FEFaceEvaluation<dim,
7415  fe_degree,
7416  n_q_points_1d,
7417  n_components,
7418  Number,
7419  VectorizedArrayType>::evaluate(const bool evaluate_values,
7420  const bool evaluate_gradients)
7421 {
7422 # ifdef DEBUG
7423  Assert(this->dof_values_initialized, ExcNotInitialized());
7424 # endif
7425 
7426  evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients);
7427 }
7428 
7429 
7430 
7431 template <int dim,
7432  int fe_degree,
7433  int n_q_points_1d,
7434  int n_components,
7435  typename Number,
7436  typename VectorizedArrayType>
7437 inline void
7438 FEFaceEvaluation<dim,
7439  fe_degree,
7440  n_q_points_1d,
7441  n_components,
7442  Number,
7443  VectorizedArrayType>::evaluate(const VectorizedArrayType
7444  * values_array,
7445  const bool evaluate_values,
7446  const bool evaluate_gradients)
7447 {
7448  if (!(evaluate_values + evaluate_gradients))
7449  return;
7450 
7451  internal::FEFaceEvaluationSelector<
7452  dim,
7453  fe_degree,
7454  n_q_points_1d,
7455  n_components,
7456  Number,
7457  VectorizedArrayType>::evaluate(*this->data,
7458  values_array,
7459  this->begin_values(),
7460  this->begin_gradients(),
7461  this->scratch_data,
7462  evaluate_values,
7463  evaluate_gradients,
7464  this->face_no,
7465  this->subface_index,
7466  this->face_orientation,
7467  this->mapping_data
7468  ->descriptor[this->active_fe_index]
7469  .face_orientations);
7470 
7471 # ifdef DEBUG
7472  if (evaluate_values == true)
7473  this->values_quad_initialized = true;
7474  if (evaluate_gradients == true)
7475  this->gradients_quad_initialized = true;
7476 # endif
7477 }
7478 
7479 
7480 
7481 template <int dim,
7482  int fe_degree,
7483  int n_q_points_1d,
7484  int n_components,
7485  typename Number,
7486  typename VectorizedArrayType>
7487 inline void
7488 FEFaceEvaluation<dim,
7489  fe_degree,
7490  n_q_points_1d,
7491  n_components,
7492  Number,
7493  VectorizedArrayType>::integrate(const bool integrate_values,
7494  const bool integrate_gradients)
7495 {
7496  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7497 
7498 # ifdef DEBUG
7499  this->dof_values_initialized = true;
7500 # endif
7501 }
7502 
7503 
7504 
7505 template <int dim,
7506  int fe_degree,
7507  int n_q_points_1d,
7508  int n_components,
7509  typename Number,
7510  typename VectorizedArrayType>
7511 inline void
7512 FEFaceEvaluation<dim,
7513  fe_degree,
7514  n_q_points_1d,
7515  n_components,
7516  Number,
7517  VectorizedArrayType>::integrate(const bool integrate_values,
7518  const bool integrate_gradients,
7519  VectorizedArrayType
7520  *values_array)
7521 {
7522  if (!(integrate_values + integrate_gradients))
7523  return;
7524 
7525  internal::FEFaceEvaluationSelector<
7526  dim,
7527  fe_degree,
7528  n_q_points_1d,
7529  n_components,
7530  Number,
7531  VectorizedArrayType>::integrate(*this->data,
7532  values_array,
7533  this->begin_values(),
7534  this->begin_gradients(),
7535  this->scratch_data,
7536  integrate_values,
7537  integrate_gradients,
7538  this->face_no,
7539  this->subface_index,
7540  this->face_orientation,
7541  this->mapping_data
7542  ->descriptor[this->active_fe_index]
7543  .face_orientations);
7544 }
7545 
7546 
7547 
7548 template <int dim,
7549  int fe_degree,
7550  int n_q_points_1d,
7551  int n_components_,
7552  typename Number,
7553  typename VectorizedArrayType>
7554 template <typename VectorType>
7555 inline void
7557  dim,
7558  fe_degree,
7559  n_q_points_1d,
7560  n_components_,
7561  Number,
7562  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7563  const bool evaluate_values,
7564  const bool evaluate_gradients)
7565 {
7566  static_assert(internal::has_begin<VectorType>::value &&
7567  (std::is_same<decltype(std::declval<VectorType>().begin()),
7568  double *>::value ||
7569  std::is_same<decltype(std::declval<VectorType>().begin()),
7570  float *>::value),
7571  "This function requires a vector type with begin() function "
7572  "evaluating to a pointer to basic number (float,double). "
7573  "Use read_dof_values() followed by evaluate() instead.");
7574 
7575  if (!internal::FEFaceEvaluationSelector<dim,
7576  fe_degree,
7577  n_q_points_1d,
7578  n_components,
7579  Number,
7580  VectorizedArrayType>::
7581  gather_evaluate(input_vector.begin(),
7582  *this->data,
7583  *this->dof_info,
7584  this->begin_values(),
7585  this->begin_gradients(),
7586  this->scratch_data,
7587  evaluate_values,
7588  evaluate_gradients,
7589  this->active_fe_index,
7590  this->first_selected_component,
7591  this->cell,
7592  this->face_no,
7593  this->subface_index,
7594  this->dof_access_index,
7595  this->face_orientation,
7596  this->mapping_data->descriptor[this->active_fe_index]
7597  .face_orientations))
7598  {
7599  this->read_dof_values(input_vector);
7600  this->evaluate(evaluate_values, evaluate_gradients);
7601  }
7602 
7603 # ifdef DEBUG
7604  if (evaluate_values == true)
7605  this->values_quad_initialized = true;
7606  if (evaluate_gradients == true)
7607  this->gradients_quad_initialized = true;
7608 # endif
7609 }
7610 
7611 
7612 
7613 template <int dim,
7614  int fe_degree,
7615  int n_q_points_1d,
7616  int n_components_,
7617  typename Number,
7618  typename VectorizedArrayType>
7619 template <typename VectorType>
7620 inline void
7622  dim,
7623  fe_degree,
7624  n_q_points_1d,
7625  n_components_,
7626  Number,
7627  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7628  const bool integrate_gradients,
7629  VectorType &destination)
7630 {
7631  static_assert(internal::has_begin<VectorType>::value &&
7632  (std::is_same<decltype(std::declval<VectorType>().begin()),
7633  double *>::value ||
7634  std::is_same<decltype(std::declval<VectorType>().begin()),
7635  float *>::value),
7636  "This function requires a vector type with begin() function "
7637  "evaluating to a pointer to basic number (float,double). "
7638  "Use integrate() followed by distribute_local_to_global() "
7639  "instead.");
7640 
7641  if (!internal::FEFaceEvaluationSelector<dim,
7642  fe_degree,
7643  n_q_points_1d,
7644  n_components,
7645  Number,
7646  VectorizedArrayType>::
7647  integrate_scatter(destination.begin(),
7648  *this->data,
7649  *this->dof_info,
7650  this->begin_dof_values(),
7651  this->begin_values(),
7652  this->begin_gradients(),
7653  this->scratch_data,
7654  integrate_values,
7655  integrate_gradients,
7656  this->active_fe_index,
7657  this->first_selected_component,
7658  this->cell,
7659  this->face_no,
7660  this->subface_index,
7661  this->dof_access_index,
7662  this->face_orientation,
7663  this->mapping_data->descriptor[this->active_fe_index]
7664  .face_orientations))
7665  {
7666  // if we arrive here, writing into the destination vector did not succeed
7667  // because some of the assumptions in integrate_scatter were not
7668  // fulfilled (e.g. an element or degree that does not support direct
7669  // writing), so we must do it here
7670  this->distribute_local_to_global(destination);
7671  }
7672 }
7673 
7674 
7675 
7676 template <int dim,
7677  int fe_degree,
7678  int n_q_points_1d,
7679  int n_components_,
7680  typename Number,
7681  typename VectorizedArrayType>
7683 FEFaceEvaluation<dim,
7684  fe_degree,
7685  n_q_points_1d,
7686  n_components_,
7687  Number,
7688  VectorizedArrayType>::quadrature_point(const unsigned int q)
7689  const
7690 {
7691  AssertIndexRange(q, n_q_points);
7692  if (this->dof_access_index < 2)
7693  {
7694  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
7695  ExcNotImplemented());
7696  AssertIndexRange(this->cell,
7697  this->mapping_data->quadrature_point_offsets.size());
7698  return this->mapping_data->quadrature_points
7699  [this->mapping_data->quadrature_point_offsets[this->cell] + q];
7700  }
7701  else
7702  {
7703  Assert(this->matrix_info->get_mapping_info()
7704  .face_data_by_cells[this->quad_no]
7705  .quadrature_point_offsets.empty() == false,
7706  ExcNotImplemented());
7707  const unsigned int index =
7708  this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
7709  AssertIndexRange(index,
7710  this->matrix_info->get_mapping_info()
7711  .face_data_by_cells[this->quad_no]
7712  .quadrature_point_offsets.size());
7713  return this->matrix_info->get_mapping_info()
7714  .face_data_by_cells[this->quad_no]
7715  .quadrature_points[this->matrix_info->get_mapping_info()
7716  .face_data_by_cells[this->quad_no]
7717  .quadrature_point_offsets[index] +
7718  q];
7719  }
7720 }
7721 
7722 
7723 
7724 /*------------------------- end FEFaceEvaluation ------------------------- */
7725 
7726 
7727 #endif // ifndef DOXYGEN
7728 
7729 
7730 DEAL_II_NAMESPACE_CLOSE
7731 
7732 #endif
VectorizedArrayType JxW(const unsigned int q_index) const
std::vector< unsigned int > plain_dof_indices
Definition: dof_info.h:519
const unsigned int dofs_per_component
static const unsigned int invalid_unsigned_int
Definition: types.h:190
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
const unsigned int n_q_points
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
unsigned int face_no
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
unsigned int face_orientation
const unsigned int first_selected_component
static ::ExceptionBase & ExcAccessToUninitializedField()
std::vector< unsigned int > component_to_base_index
Definition: dof_info.h:549
value_type get_dof_value(const unsigned int dof) const
std::vector< IndexStorageVariants > index_storage_variants[3]
Definition: dof_info.h:395
unsigned int n_components() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
std::vector< unsigned int > dof_indices_interleave_strides[3]
Definition: dof_info.h:457
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
std::vector< types::global_dof_index > local_dof_indices
const Tensor< 1, dim, VectorizedArrayType > * normal_vectors
const Point< dim > & point(const unsigned int i) const
const unsigned int n_quadrature_points
bool mapping_initialized() const
value_type get_laplacian(const unsigned int q_point) const
void read_write_operation_contiguous(const VectorOperation &operation, VectorType *vectors[], const std::bitset< VectorizedArrayType::size()> &mask) const
VectorizedArrayType read_cell_data(const AlignedVector< VectorizedArrayType > &array) const
Transformed quadrature points.
static ::ExceptionBase & ExcNotInitialized()
unsigned int cell
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:468
VectorizedArrayType * scratch_data
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_info
void set_data_pointers()
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > * data
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Tensor< 1, dim, VectorizedArrayType > get_normal_vector(const unsigned int q_point) const
std::vector< unsigned int > dof_indices
Definition: dof_info.h:420
Definition: point.h:111
const Number * constraint_pool_end(const unsigned int pool_index) const
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
std::vector< unsigned char > n_vectorization_lanes_filled[3]
Definition: dof_info.h:468
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
unsigned int subface_index
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > get_hessian(const unsigned int q_point) const
const Tensor< 2, dim, VectorizedArrayType > * jacobian
const Number * constraint_pool_begin(const unsigned int pool_index) const
const VectorizedArrayType * J_value
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
static ::ExceptionBase & ExcMessage(std::string arg1)
bool indices_initialized() const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
Definition: fe.h:3231
std::vector< unsigned int > dof_indices_interleaved
Definition: dof_info.h:437
#define Assert(cond, exc)
Definition: exceptions.h:1419
const unsigned int n_q_points
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3133
UpdateFlags
value_type integrate_value() const
std::vector< unsigned int > row_starts_plain_indices
Definition: dof_info.h:509
const unsigned int active_quad_index
const unsigned int active_fe_index
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:486
unsigned int get_mapping_data_index_offset() const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition: dof_info.h:403
const VectorizedArrayType * begin_values() const
const Number * quadrature_weights
const unsigned int dofs_per_component
VectorizedArrayType get_divergence(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition: dof_info.h:432
internal::MatrixFreeFunctions::GeometryType get_cell_type() const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:472
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true)
const std::vector< unsigned int > & get_internal_dof_numbering() const
unsigned int size() const
std::array< unsigned int, VectorizedArrayType::size()> get_cell_ids() const
ArrayView< VectorizedArrayType > get_scratch_data() const
Tensor< 2, dim, VectorizedArrayType > inverse_jacobian(const unsigned int q_index) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
bool hessians_quad_initialized
AlignedVector< VectorizedArrayType > * scratch_data_array
std::vector< unsigned int > face_partition_data
Definition: task_info.h:477
gradient_type get_gradient(const unsigned int q_point) const
value_type get_normal_derivative(const unsigned int q_point) const
void fill_JxW_values(AlignedVector< VectorizedArrayType > &JxW_values) const
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
const unsigned int dofs_per_cell
const unsigned int n_fe_components
internal::MatrixFreeFunctions::GeometryType cell_type
const unsigned int quad_no
Definition: tensor.h:417
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
std::vector< unsigned int > dof_indices_contiguous[3]
Definition: dof_info.h:447
const VectorizedArrayType * begin_dof_values() const
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition: dof_info.h:562
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info() const
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0)
const VectorizedArrayType * begin_hessians() const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void submit_value(const value_type val_in, const unsigned int q_point)
size_type size() const
const VectorizedArrayType * begin_gradients() const
FEEvaluationBase & operator=(const FEEvaluationBase &other)
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number, VectorizedArrayType > > mapped_geometry
static ::ExceptionBase & ExcNotImplemented()
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
std::vector< UnivariateShapeData< Number > > data
Definition: shape_info.h:355
const Tensor< 1, dim, VectorizedArrayType > * normal_x_jacobian
bool gradients_quad_initialized
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
std::vector< unsigned int > start_components
Definition: dof_info.h:543
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:349
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_number) const
const internal::MatrixFreeFunctions::DoFInfo * dof_info
void read_write_operation_global(const VectorOperation &operation, VectorType *vectors[]) const
void read_write_operation(const VectorOperation &operation, VectorType *vectors[], const std::bitset< VectorizedArrayType::size()> &mask, const bool apply_constraints=true) const
const unsigned int dofs_per_cell
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
value_type get_value(const unsigned int q_point) const
static ::ExceptionBase & ExcInternalError()
gradient_type get_hessian_diagonal(const unsigned int q_point) const