Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
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 
40 
41 DEAL_II_NAMESPACE_OPEN
42 
43 
44 
45 namespace internal
46 {
48 }
49 
50 template <int dim,
51  int fe_degree,
52  int n_q_points_1d = fe_degree + 1,
53  int n_components_ = 1,
54  typename Number = double>
56 
57 
83 template <int dim, int n_components_, typename Number, bool is_face = false>
85 {
86 public:
87  using number_type = Number;
89  using gradient_type =
91  static constexpr unsigned int dimension = dim;
92  static constexpr unsigned int n_components = n_components_;
93 
102 
106  DEAL_II_DEPRECATED
107  unsigned int
108  get_cell_data_number() const;
109 
118  unsigned int
120 
128  get_cell_type() const;
129 
134  get_shape_info() const;
135 
137 
174  template <typename VectorType>
175  void
176  read_dof_values(const VectorType &src, const unsigned int first_index = 0);
177 
206  template <typename VectorType>
207  void
208  read_dof_values_plain(const VectorType & src,
209  const unsigned int first_index = 0);
210 
241  template <typename VectorType>
242  void
244  VectorType & dst,
245  const unsigned int first_index = 0,
246  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask =
247  std::bitset<VectorizedArray<Number>::n_array_elements>().flip()) const;
248 
281  template <typename VectorType>
282  void
284  VectorType & dst,
285  const unsigned int first_index = 0,
286  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask =
287  std::bitset<VectorizedArray<Number>::n_array_elements>().flip()) const;
288 
290 
311  value_type
312  get_dof_value(const unsigned int dof) const;
313 
324  void
325  submit_dof_value(const value_type val_in, const unsigned int dof);
326 
338  value_type
339  get_value(const unsigned int q_point) const;
340 
352  void
353  submit_value(const value_type val_in, const unsigned int q_point);
354 
365  get_gradient(const unsigned int q_point) const;
366 
380  value_type
381  get_normal_derivative(const unsigned int q_point) const;
382 
395  void
396  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
397 
415  void
416  submit_normal_derivative(const value_type grad_in,
417  const unsigned int q_point);
418 
430  get_hessian(const unsigned int q_point) const;
431 
441  get_hessian_diagonal(const unsigned int q_point) const;
442 
453  value_type
454  get_laplacian(const unsigned int q_point) const;
455 
456 #ifdef DOXYGEN
457  // doxygen does not anyhow mention functions coming from partial template
458  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
459  // For now, hack-in those functions manually only to fix documentation:
460 
465  get_divergence(const unsigned int q_point) const;
466 
471  get_symmetric_gradient(const unsigned int q_point) const;
472 
477  get_curl(const unsigned int q_point) const;
478 
487  void
489  const unsigned int q_point);
490 
500  void
502  const SymmetricTensor<2, dim, VectorizedArray<Number>> grad_in,
503  const unsigned int q_point);
504 
514  void
515  submit_curl(
516  const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl_in,
517  const unsigned int q_point);
518 
519 #endif
520 
537  value_type
538  integrate_value() const;
539 
545  JxW(const unsigned int q_index) const;
546 
552  DEAL_II_DEPRECATED void
554 
562  inverse_jacobian(const unsigned int q_index) const;
563 
577  get_normal_vector(const unsigned int q_point) const;
578 
587 
589 
603  begin_dof_values() const;
604 
615 
627  begin_values() const;
628 
640  begin_values();
641 
654  begin_gradients() const;
655 
668  begin_gradients();
669 
683  begin_hessians() const;
684 
698  begin_hessians();
699 
705  const std::vector<unsigned int> &
707 
715  get_scratch_data() const;
716 
718 
719 protected:
728  FEEvaluationBase(const MatrixFree<dim, Number> &matrix_free,
729  const unsigned int dof_no,
730  const unsigned int first_selected_component,
731  const unsigned int quad_no,
732  const unsigned int fe_degree,
733  const unsigned int n_q_points,
734  const bool is_interior_face);
735 
770  template <int n_components_other>
772  const Mapping<dim> & mapping,
773  const FiniteElement<dim> &fe,
774  const Quadrature<1> & quadrature,
775  const UpdateFlags update_flags,
776  const unsigned int first_selected_component,
778 
785  FEEvaluationBase(const FEEvaluationBase &other);
786 
794  operator=(const FEEvaluationBase &other);
795 
802  template <typename VectorType, typename VectorOperation>
803  void
805  const VectorOperation & operation,
806  VectorType * vectors[],
807  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask,
808  const bool apply_constraints = true) const;
809 
817  template <typename VectorType, typename VectorOperation>
818  void
820  const VectorOperation & operation,
821  VectorType * vectors[],
822  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask) const;
823 
831  template <typename VectorType, typename VectorOperation>
832  void
834  VectorType * vectors[]) const;
835 
840 
847 
861 
874 
889 
901  VectorizedArray<Number> *hessians_quad[n_components][(dim * (dim + 1)) / 2];
902 
906  const unsigned int quad_no;
907 
912  const unsigned int n_fe_components;
913 
918  const unsigned int active_fe_index;
919 
924  const unsigned int active_quad_index;
925 
929  const unsigned int n_quadrature_points;
930 
935 
942 
949  const internal::MatrixFreeFunctions::
950  MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number> *mapping_data;
951 
959 
965 
973 
978 
983 
987  const Number *quadrature_weights;
988 
993  unsigned int cell;
994 
1000 
1006 
1011  unsigned int face_no;
1012 
1017  unsigned int face_orientation;
1018 
1026  unsigned int subface_index;
1027 
1035 
1042 
1049 
1056 
1063 
1070 
1077 
1082  std::shared_ptr<
1085 
1090  const unsigned int first_selected_component;
1091 
1096  mutable std::vector<types::global_dof_index> local_dof_indices;
1097 
1098 private:
1103  void
1105 
1109  template <int, int, typename, bool>
1110  friend class FEEvaluationBase;
1111  template <int, int, int, int, typename>
1112  friend class FEEvaluation;
1113 };
1114 
1115 
1116 
1126 template <int dim, int n_components_, typename Number, bool is_face>
1128  : public FEEvaluationBase<dim, n_components_, Number, is_face>
1129 {
1130 public:
1131  using number_type = Number;
1133  using gradient_type =
1135  static constexpr unsigned int dimension = dim;
1136  static constexpr unsigned int n_components = n_components_;
1138 
1139 protected:
1147  FEEvaluationAccess(const MatrixFree<dim, Number> &matrix_free,
1148  const unsigned int dof_no,
1149  const unsigned int first_selected_component,
1150  const unsigned int quad_no,
1151  const unsigned int fe_degree,
1152  const unsigned int n_q_points,
1153  const bool is_interior_face = true);
1154 
1159  template <int n_components_other>
1161  const Mapping<dim> & mapping,
1162  const FiniteElement<dim> &fe,
1163  const Quadrature<1> & quadrature,
1164  const UpdateFlags update_flags,
1165  const unsigned int first_selected_component,
1167 
1171  FEEvaluationAccess(const FEEvaluationAccess &other);
1172 
1177  operator=(const FEEvaluationAccess &other);
1178 };
1179 
1180 
1181 
1192 template <int dim, typename Number, bool is_face>
1193 class FEEvaluationAccess<dim, 1, Number, is_face>
1194  : public FEEvaluationBase<dim, 1, Number, is_face>
1195 {
1196 public:
1197  using number_type = Number;
1200  static constexpr unsigned int dimension = dim;
1202 
1205  value_type
1206  get_dof_value(const unsigned int dof) const;
1207 
1210  void
1211  submit_dof_value(const value_type val_in, const unsigned int dof);
1212 
1215  value_type
1216  get_value(const unsigned int q_point) const;
1217 
1220  void
1221  submit_value(const value_type val_in, const unsigned int q_point);
1222 
1225  void
1226  submit_value(const Tensor<1, 1, VectorizedArray<Number>> val_in,
1227  const unsigned int q_point);
1228 
1232  get_gradient(const unsigned int q_point) const;
1233 
1236  value_type
1237  get_normal_derivative(const unsigned int q_point) const;
1238 
1241  void
1242  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1243 
1246  void
1247  submit_normal_derivative(const value_type grad_in,
1248  const unsigned int q_point);
1249 
1253  get_hessian(unsigned int q_point) const;
1254 
1258  get_hessian_diagonal(const unsigned int q_point) const;
1259 
1262  value_type
1263  get_laplacian(const unsigned int q_point) const;
1264 
1267  value_type
1268  integrate_value() const;
1269 
1270 protected:
1278  FEEvaluationAccess(const MatrixFree<dim, Number> &matrix_free,
1279  const unsigned int dof_no,
1280  const unsigned int first_selected_component,
1281  const unsigned int quad_no,
1282  const unsigned int fe_degree,
1283  const unsigned int n_q_points,
1284  const bool is_interior_face = true);
1285 
1290  template <int n_components_other>
1292  const Mapping<dim> & mapping,
1293  const FiniteElement<dim> &fe,
1294  const Quadrature<1> & quadrature,
1295  const UpdateFlags update_flags,
1296  const unsigned int first_selected_component,
1298 
1302  FEEvaluationAccess(const FEEvaluationAccess &other);
1303 
1308  operator=(const FEEvaluationAccess &other);
1309 };
1310 
1311 
1312 
1324 template <int dim, typename Number, bool is_face>
1325 class FEEvaluationAccess<dim, dim, Number, is_face>
1326  : public FEEvaluationBase<dim, dim, Number, is_face>
1327 {
1328 public:
1329  using number_type = Number;
1332  static constexpr unsigned int dimension = dim;
1333  static constexpr unsigned int n_components = dim;
1335 
1339  get_gradient(const unsigned int q_point) const;
1340 
1346  get_divergence(const unsigned int q_point) const;
1347 
1355  get_symmetric_gradient(const unsigned int q_point) const;
1356 
1362  get_curl(const unsigned int q_point) const;
1363 
1367  get_hessian(const unsigned int q_point) const;
1368 
1372  get_hessian_diagonal(const unsigned int q_point) const;
1373 
1376  void
1377  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1378 
1387  void
1389  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> grad_in,
1390  const unsigned int q_point);
1391 
1400  void
1402  const unsigned int q_point);
1403 
1412  void
1414  const SymmetricTensor<2, dim, VectorizedArray<Number>> grad_in,
1415  const unsigned int q_point);
1416 
1421  void
1422  submit_curl(
1423  const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl_in,
1424  const unsigned int q_point);
1425 
1426 protected:
1434  FEEvaluationAccess(const MatrixFree<dim, Number> &matrix_free,
1435  const unsigned int dof_no,
1436  const unsigned int first_selected_component,
1437  const unsigned int quad_no,
1438  const unsigned int dofs_per_cell,
1439  const unsigned int n_q_points,
1440  const bool is_interior_face = true);
1441 
1446  template <int n_components_other>
1448  const Mapping<dim> & mapping,
1449  const FiniteElement<dim> &fe,
1450  const Quadrature<1> & quadrature,
1451  const UpdateFlags update_flags,
1452  const unsigned int first_selected_component,
1454 
1458  FEEvaluationAccess(const FEEvaluationAccess &other);
1459 
1464  operator=(const FEEvaluationAccess &other);
1465 };
1466 
1467 
1479 template <typename Number, bool is_face>
1480 class FEEvaluationAccess<1, 1, Number, is_face>
1481  : public FEEvaluationBase<1, 1, Number, is_face>
1482 {
1483 public:
1484  using number_type = Number;
1487  static constexpr unsigned int dimension = 1;
1489 
1492  value_type
1493  get_dof_value(const unsigned int dof) const;
1494 
1497  void
1498  submit_dof_value(const value_type val_in, const unsigned int dof);
1499 
1502  value_type
1503  get_value(const unsigned int q_point) const;
1504 
1507  void
1508  submit_value(const value_type val_in, const unsigned int q_point);
1509 
1512  void
1513  submit_value(const gradient_type val_in, const unsigned int q_point);
1514 
1518  get_gradient(const unsigned int q_point) const;
1519 
1522  value_type
1523  get_normal_derivative(const unsigned int q_point) const;
1524 
1527  void
1528  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1529 
1532  void
1533  submit_gradient(const value_type grad_in, const unsigned int q_point);
1534 
1537  void
1538  submit_normal_derivative(const value_type grad_in,
1539  const unsigned int q_point);
1540 
1543  void
1545  const unsigned int q_point);
1546 
1550  get_hessian(unsigned int q_point) const;
1551 
1555  get_hessian_diagonal(const unsigned int q_point) const;
1556 
1559  value_type
1560  get_laplacian(const unsigned int q_point) const;
1561 
1564  value_type
1565  integrate_value() const;
1566 
1567 protected:
1575  FEEvaluationAccess(const MatrixFree<1, Number> &matrix_free,
1576  const unsigned int dof_no,
1577  const unsigned int first_selected_component,
1578  const unsigned int quad_no,
1579  const unsigned int fe_degree,
1580  const unsigned int n_q_points,
1581  const bool is_interior_face = true);
1582 
1587  template <int n_components_other>
1589  const Mapping<1> & mapping,
1590  const FiniteElement<1> &fe,
1591  const Quadrature<1> & quadrature,
1592  const UpdateFlags update_flags,
1593  const unsigned int first_selected_component,
1595 
1599  FEEvaluationAccess(const FEEvaluationAccess &other);
1600 
1605  operator=(const FEEvaluationAccess &other);
1606 };
1607 
1608 
1609 
2159 template <int dim,
2160  int fe_degree,
2161  int n_q_points_1d,
2162  int n_components_,
2163  typename Number>
2164 class FEEvaluation
2165  : public FEEvaluationAccess<dim, n_components_, Number, false>
2166 {
2167 public:
2172 
2176  using number_type = Number;
2177 
2184 
2191 
2195  static constexpr unsigned int dimension = dim;
2196 
2201  static constexpr unsigned int n_components = n_components_;
2202 
2209  static constexpr unsigned int static_n_q_points =
2210  Utilities::pow(n_q_points_1d, dim);
2211 
2219  static constexpr unsigned int static_dofs_per_component =
2220  Utilities::pow(fe_degree + 1, dim);
2221 
2229  static constexpr unsigned int tensor_dofs_per_cell =
2231 
2239  static constexpr unsigned int static_dofs_per_cell =
2241 
2267  FEEvaluation(const MatrixFree<dim, Number> &matrix_free,
2268  const unsigned int dof_no = 0,
2269  const unsigned int quad_no = 0,
2270  const unsigned int first_selected_component = 0);
2271 
2298  FEEvaluation(const Mapping<dim> & mapping,
2299  const FiniteElement<dim> &fe,
2300  const Quadrature<1> & quadrature,
2301  const UpdateFlags update_flags,
2302  const unsigned int first_selected_component = 0);
2303 
2309  FEEvaluation(const FiniteElement<dim> &fe,
2310  const Quadrature<1> & quadrature,
2311  const UpdateFlags update_flags,
2312  const unsigned int first_selected_component = 0);
2313 
2324  template <int n_components_other>
2325  FEEvaluation(const FiniteElement<dim> & fe,
2327  const unsigned int first_selected_component = 0);
2328 
2335  FEEvaluation(const FEEvaluation &other);
2336 
2343  FEEvaluation &
2344  operator=(const FEEvaluation &other);
2345 
2354  void
2355  reinit(const unsigned int cell_batch_index);
2356 
2369  template <typename DoFHandlerType, bool level_dof_access>
2370  void
2372  &cell);
2373 
2384  void
2386 
2396  void
2397  evaluate(const bool evaluate_values,
2398  const bool evaluate_gradients,
2399  const bool evaluate_hessians = false);
2400 
2413  void
2414  evaluate(const VectorizedArray<Number> *values_array,
2415  const bool evaluate_values,
2416  const bool evaluate_gradients,
2417  const bool evaluate_hessians = false);
2418 
2432  template <typename VectorType>
2433  void
2434  gather_evaluate(const VectorType &input_vector,
2435  const bool evaluate_values,
2436  const bool evaluate_gradients,
2437  const bool evaluate_hessians = false);
2438 
2449  void
2450  integrate(const bool integrate_values, const bool integrate_gradients);
2451 
2463  void
2464  integrate(const bool integrate_values,
2465  const bool integrate_gradients,
2466  VectorizedArray<Number> *values_array);
2467 
2481  template <typename VectorType>
2482  void
2483  integrate_scatter(const bool integrate_values,
2484  const bool integrate_gradients,
2485  VectorType &output_vector);
2486 
2492  quadrature_point(const unsigned int q_point) const;
2493 
2500  const unsigned int dofs_per_component;
2501 
2508  const unsigned int dofs_per_cell;
2509 
2517  const unsigned int n_q_points;
2518 
2519 private:
2524  void
2525  check_template_arguments(const unsigned int fe_no,
2526  const unsigned int first_selected_component);
2527 };
2528 
2529 
2530 
2562 template <int dim,
2563  int fe_degree,
2564  int n_q_points_1d = fe_degree + 1,
2565  int n_components_ = 1,
2566  typename Number = double>
2568  : public FEEvaluationAccess<dim, n_components_, Number, true>
2569 {
2570 public:
2575 
2579  using number_type = Number;
2580 
2587 
2594 
2598  static constexpr unsigned int dimension = dim;
2599 
2604  static constexpr unsigned int n_components = n_components_;
2605 
2613  static constexpr unsigned int static_n_q_points =
2614  Utilities::pow(n_q_points_1d, dim - 1);
2615 
2622  static constexpr unsigned int static_n_q_points_cell =
2623  Utilities::pow(n_q_points_1d, dim);
2624 
2631  static constexpr unsigned int static_dofs_per_component =
2632  Utilities::pow(fe_degree + 1, dim);
2633 
2640  static constexpr unsigned int tensor_dofs_per_cell =
2642 
2649  static constexpr unsigned int static_dofs_per_cell =
2651 
2681  FEFaceEvaluation(const MatrixFree<dim, Number> &matrix_free,
2682  const bool is_interior_face = true,
2683  const unsigned int dof_no = 0,
2684  const unsigned int quad_no = 0,
2685  const unsigned int first_selected_component = 0);
2686 
2697  void
2698  reinit(const unsigned int face_batch_number);
2699 
2707  void
2708  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2709 
2720  void
2721  evaluate(const bool evaluate_values, const bool evaluate_gradients);
2722 
2735  void
2736  evaluate(const VectorizedArray<Number> *values_array,
2737  const bool evaluate_values,
2738  const bool evaluate_gradients);
2739 
2751  template <typename VectorType>
2752  void
2753  gather_evaluate(const VectorType &input_vector,
2754  const bool evaluate_values,
2755  const bool evaluate_gradients);
2756 
2766  void
2767  integrate(const bool integrate_values, const bool integrate_gradients);
2768 
2777  void
2778  integrate(const bool integrate_values,
2779  const bool integrate_gradients,
2780  VectorizedArray<Number> *values_array);
2781 
2793  template <typename VectorType>
2794  void
2795  integrate_scatter(const bool integrate_values,
2796  const bool integrate_gradients,
2797  VectorType &output_vector);
2798 
2804  quadrature_point(const unsigned int q_point) const;
2805 
2812  const unsigned int dofs_per_component;
2813 
2820  const unsigned int dofs_per_cell;
2821 
2829  const unsigned int n_q_points;
2830 
2831 protected:
2837  void
2839  const bool values,
2840  const bool gradients);
2841 };
2842 
2843 
2844 
2845 namespace internal
2846 {
2847  namespace MatrixFreeFunctions
2848  {
2849  // a helper function to compute the number of DoFs of a DGP element at
2850  // compile time, depending on the degree
2851  template <int dim, int degree>
2852  struct DGP_dofs_per_component
2853  {
2854  // this division is always without remainder
2855  static constexpr unsigned int value =
2856  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2857  };
2858 
2859  // base specialization: 1d elements have 'degree+1' degrees of freedom
2860  template <int degree>
2861  struct DGP_dofs_per_component<1, degree>
2862  {
2863  static constexpr unsigned int value = degree + 1;
2864  };
2865  } // namespace MatrixFreeFunctions
2866 } // namespace internal
2867 
2868 
2869 /*----------------------- Inline functions ----------------------------------*/
2870 
2871 #ifndef DOXYGEN
2872 
2873 
2874 
2875 /*----------------------- FEEvaluationBase ----------------------------------*/
2876 
2877 template <int dim, int n_components_, typename Number, bool is_face>
2879  const MatrixFree<dim, Number> &data_in,
2880  const unsigned int dof_no,
2881  const unsigned int first_selected_component,
2882  const unsigned int quad_no_in,
2883  const unsigned int fe_degree,
2884  const unsigned int n_q_points,
2885  const bool is_interior_face)
2886  : scratch_data_array(data_in.acquire_scratch_data())
2887  , quad_no(quad_no_in)
2888  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
2889  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
2890  data_in.get_dof_info(dof_no).fe_index_from_degree(
2891  first_selected_component,
2892  fe_degree) :
2893  0)
2894  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
2895  (is_face ? data_in.get_mapping_info()
2896  .face_data[quad_no_in]
2897  .quad_index_from_n_q_points(n_q_points) :
2898  data_in.get_mapping_info()
2899  .cell_data[quad_no_in]
2900  .quad_index_from_n_q_points(n_q_points)) :
2901  0)
2902  , n_quadrature_points(fe_degree != numbers::invalid_unsigned_int ?
2903  n_q_points :
2904  (is_face ? data_in
2905  .get_shape_info(dof_no,
2906  quad_no_in,
2907  active_fe_index,
2908  active_quad_index)
2909  .n_q_points_face :
2910  data_in
2911  .get_shape_info(dof_no,
2912  quad_no_in,
2913  active_fe_index,
2914  active_quad_index)
2915  .n_q_points))
2916  , matrix_info(&data_in)
2917  , dof_info(&data_in.get_dof_info(dof_no))
2918  , mapping_data(internal::MatrixFreeFunctions::
2919  MappingInfoCellsOrFaces<dim, Number, is_face>::get(
2920  data_in.get_mapping_info(),
2921  quad_no))
2922  , data(&data_in.get_shape_info(
2923  dof_no,
2924  quad_no_in,
2925  dof_info->component_to_base_index[first_selected_component],
2926  active_fe_index,
2927  active_quad_index))
2928  , jacobian(nullptr)
2929  , J_value(nullptr)
2930  , normal_vectors(nullptr)
2931  , normal_x_jacobian(nullptr)
2932  , quadrature_weights(
2933  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
2934  , cell(numbers::invalid_unsigned_int)
2935  , is_interior_face(is_interior_face)
2936  , dof_access_index(
2937  is_face ?
2938  (is_interior_face ?
2939  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
2940  internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
2941  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
2942  , cell_type(internal::MatrixFreeFunctions::general)
2943  , dof_values_initialized(false)
2944  , values_quad_initialized(false)
2945  , gradients_quad_initialized(false)
2946  , hessians_quad_initialized(false)
2947  , values_quad_submitted(false)
2948  , gradients_quad_submitted(false)
2949  , first_selected_component(first_selected_component)
2950 {
2952  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
2953  AssertDimension(matrix_info->get_size_info().vectorization_length,
2959  Assert(
2960  dof_info->start_components.back() == 1 ||
2961  static_cast<int>(n_components_) <=
2962  static_cast<int>(
2966  ExcMessage(
2967  "You tried to construct a vector-valued evaluator with " +
2969  " components. However, "
2970  "the current base element has only " +
2975  " components left when starting from local element index " +
2980  " (global index " + Utilities::to_string(first_selected_component) +
2981  ")"));
2982 
2983  // do not check for correct dimensions of data fields here, should be done
2984  // in derived classes
2985 }
2986 
2987 
2988 
2989 template <int dim, int n_components_, typename Number, bool is_face>
2990 template <int n_components_other>
2992  const Mapping<dim> & mapping,
2993  const FiniteElement<dim> &fe,
2994  const Quadrature<1> & quadrature,
2995  const UpdateFlags update_flags,
2996  const unsigned int first_selected_component,
2998  : scratch_data_array(new AlignedVector<VectorizedArray<Number>>())
2999  , quad_no(numbers::invalid_unsigned_int)
3000  , n_fe_components(n_components_)
3001  , active_fe_index(numbers::invalid_unsigned_int)
3002  , active_quad_index(numbers::invalid_unsigned_int)
3003  , n_quadrature_points(
3004  Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
3005  , matrix_info(nullptr)
3006  , dof_info(nullptr)
3007  , mapping_data(nullptr)
3008  ,
3009  // select the correct base element from the given FE component
3010  data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(
3011  quadrature,
3012  fe,
3013  fe.component_to_base_index(first_selected_component).first))
3014  , jacobian(nullptr)
3015  , J_value(nullptr)
3016  , normal_vectors(nullptr)
3017  , normal_x_jacobian(nullptr)
3018  , quadrature_weights(nullptr)
3019  , cell(0)
3020  , cell_type(internal::MatrixFreeFunctions::general)
3021  , is_interior_face(true)
3022  , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
3023  , dof_values_initialized(false)
3024  , values_quad_initialized(false)
3025  , gradients_quad_initialized(false)
3026  , hessians_quad_initialized(false)
3027  , values_quad_submitted(false)
3028  , gradients_quad_submitted(false)
3029  ,
3030  // keep the number of the selected component within the current base element
3031  // for reading dof values
3032  first_selected_component(first_selected_component)
3033 {
3035 
3036  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
3037  ExcInternalError());
3038  if (other != nullptr &&
3039  other->mapped_geometry->get_quadrature() == quadrature)
3041  else
3042  mapped_geometry = std::make_shared<
3044  mapping, quadrature, update_flags);
3045  cell = 0;
3046 
3047  mapping_data = &mapped_geometry->get_data_storage();
3048  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3049  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3050 
3051  const unsigned int base_element_number =
3053  Assert(fe.element_multiplicity(base_element_number) == 1 ||
3054  fe.element_multiplicity(base_element_number) -
3056  n_components_,
3057  ExcMessage("The underlying element must at least contain as many "
3058  "components as requested by this class"));
3059  (void)base_element_number;
3060 }
3061 
3062 
3063 
3064 template <int dim, int n_components_, typename Number, bool is_face>
3067  : scratch_data_array(other.matrix_info == nullptr ?
3068  new AlignedVector<VectorizedArray<Number>>() :
3069  other.matrix_info->acquire_scratch_data())
3070  , quad_no(other.quad_no)
3071  , n_fe_components(other.n_fe_components)
3072  , active_fe_index(other.active_fe_index)
3073  , active_quad_index(other.active_quad_index)
3074  , n_quadrature_points(other.n_quadrature_points)
3075  , matrix_info(other.matrix_info)
3076  , dof_info(other.dof_info)
3077  , mapping_data(other.mapping_data)
3078  , data(
3079  other.matrix_info == nullptr ?
3080  new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(
3081  *other.data) :
3082  other.data)
3083  , jacobian(nullptr)
3084  , J_value(nullptr)
3085  , normal_vectors(nullptr)
3086  , normal_x_jacobian(nullptr)
3087  , quadrature_weights(
3088  other.matrix_info == nullptr ?
3089  nullptr :
3090  mapping_data->descriptor[active_quad_index].quadrature_weights.begin())
3091  , cell(numbers::invalid_unsigned_int)
3092  , cell_type(internal::MatrixFreeFunctions::general)
3093  , is_interior_face(other.is_interior_face)
3094  , dof_access_index(other.dof_access_index)
3095  , dof_values_initialized(false)
3096  , values_quad_initialized(false)
3097  , gradients_quad_initialized(false)
3098  , hessians_quad_initialized(false)
3099  , values_quad_submitted(false)
3100  , gradients_quad_submitted(false)
3101  , first_selected_component(other.first_selected_component)
3102 {
3104 
3105  // Create deep copy of mapped geometry for use in parallel...
3106  if (other.mapped_geometry.get() != nullptr)
3107  {
3108  mapped_geometry.reset(
3110  other.mapped_geometry->get_fe_values().get_mapping(),
3111  other.mapped_geometry->get_quadrature(),
3112  other.mapped_geometry->get_fe_values().get_update_flags()));
3113  mapping_data = &mapped_geometry->get_data_storage();
3114  cell = 0;
3115 
3116  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3117  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3118  }
3119 }
3120 
3121 
3122 
3123 template <int dim, int n_components_, typename Number, bool is_face>
3127 {
3128  AssertDimension(quad_no, other.quad_no);
3129  AssertDimension(n_fe_components, other.n_fe_components);
3130  AssertDimension(active_fe_index, other.active_fe_index);
3131  AssertDimension(active_quad_index, other.active_quad_index);
3132  AssertDimension(first_selected_component, other.first_selected_component);
3133 
3134  // release old memory
3135  if (matrix_info == nullptr)
3136  {
3137  delete data;
3138  delete scratch_data_array;
3139  }
3140  else
3141  {
3142  matrix_info->release_scratch_data(scratch_data_array);
3143  }
3144 
3145  matrix_info = other.matrix_info;
3146  dof_info = other.dof_info;
3147  mapping_data = other.mapping_data;
3148  if (other.matrix_info == nullptr)
3149  {
3150  data =
3152  *other.data);
3153  scratch_data_array = new AlignedVector<VectorizedArray<Number>>();
3154  }
3155  else
3156  {
3157  data = other.data;
3158  scratch_data_array = matrix_info->acquire_scratch_data();
3159  }
3160  set_data_pointers();
3161 
3162  quadrature_weights =
3163  (mapping_data != nullptr ?
3164  mapping_data->descriptor[active_quad_index].quadrature_weights.begin() :
3165  nullptr);
3168  is_interior_face = other.is_interior_face;
3169  dof_access_index = other.dof_access_index;
3170 
3171  // Create deep copy of mapped geometry for use in parallel...
3172  if (other.mapped_geometry.get() != nullptr)
3173  {
3174  mapped_geometry.reset(
3176  other.mapped_geometry->get_fe_values().get_mapping(),
3177  other.mapped_geometry->get_quadrature(),
3178  other.mapped_geometry->get_fe_values().get_update_flags()));
3179  cell = 0;
3180  mapping_data = &mapped_geometry->get_data_storage();
3181  jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
3182  J_value = mapped_geometry->get_data_storage().JxW_values.begin();
3183  }
3184 
3185  return *this;
3186 }
3187 
3188 
3189 
3190 template <int dim, int n_components_, typename Number, bool is_face>
3193 {
3194  if (matrix_info != nullptr)
3195  {
3196  try
3197  {
3198  matrix_info->release_scratch_data(scratch_data_array);
3199  }
3200  catch (...)
3201  {}
3202  }
3203  else
3204  {
3205  delete scratch_data_array;
3206  delete data;
3207  data = nullptr;
3208  }
3209  scratch_data_array = nullptr;
3210 }
3211 
3212 
3213 
3214 template <int dim, int n_components_, typename Number, bool is_face>
3215 inline void
3217 {
3218  Assert(scratch_data_array != nullptr, ExcInternalError());
3219 
3220  const unsigned int tensor_dofs_per_component =
3221  Utilities::fixed_power<dim>(this->data->fe_degree + 1);
3222  const unsigned int dofs_per_component =
3223  this->data->dofs_per_component_on_cell;
3224  const unsigned int n_quadrature_points =
3225  is_face ? this->data->n_q_points_face : this->data->n_q_points;
3226 
3227  const unsigned int shift =
3228  std::max(tensor_dofs_per_component + 1, dofs_per_component) *
3229  n_components_ * 3 +
3230  2 * n_quadrature_points;
3231  const unsigned int allocated_size =
3232  shift + n_components_ * dofs_per_component +
3233  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3234  scratch_data_array->resize_fast(allocated_size);
3235 
3236  // set the pointers to the correct position in the data array
3237  for (unsigned int c = 0; c < n_components_; ++c)
3238  {
3239  this->values_dofs[c] =
3240  scratch_data_array->begin() + c * dofs_per_component;
3241  this->values_quad[c] = scratch_data_array->begin() +
3242  n_components * dofs_per_component +
3243  c * n_quadrature_points;
3244  for (unsigned int d = 0; d < dim; ++d)
3245  this->gradients_quad[c][d] =
3246  scratch_data_array->begin() +
3247  n_components * (dofs_per_component + n_quadrature_points) +
3248  (c * dim + d) * n_quadrature_points;
3249  for (unsigned int d = 0; d < (dim * dim + dim) / 2; ++d)
3250  this->hessians_quad[c][d] =
3251  scratch_data_array->begin() +
3252  n_components *
3253  ((dim + 1) * n_quadrature_points + dofs_per_component) +
3254  (c * (dim * dim + dim) + d) * n_quadrature_points;
3255  }
3256  scratch_data =
3257  scratch_data_array->begin() + n_components_ * dofs_per_component +
3258  (n_components_ * (dim * dim + 2 * dim + 1) * n_quadrature_points);
3259 }
3260 
3261 
3262 
3263 template <int dim, int n_components_, typename Number, bool is_face>
3264 inline unsigned int
3266  const
3267 {
3268  return get_mapping_data_index_offset();
3269 }
3270 
3271 
3272 
3273 template <int dim, int n_components_, typename Number, bool is_face>
3274 inline unsigned int
3277 {
3278  if (matrix_info == nullptr)
3279  return 0;
3280  else
3281  {
3282  AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
3283  return this->mapping_data->data_index_offsets[cell];
3284  }
3285 }
3286 
3287 
3288 
3289 template <int dim, int n_components_, typename Number, bool is_face>
3292 {
3294  return cell_type;
3295 }
3296 
3297 
3298 
3299 template <int dim, int n_components_, typename Number, bool is_face>
3302 {
3303  Assert(data != nullptr, ExcInternalError());
3304  return *data;
3305 }
3306 
3307 
3308 
3309 template <int dim, int n_components_, typename Number, bool is_face>
3310 inline void
3312  AlignedVector<VectorizedArray<Number>> &JxW_values) const
3313 {
3314  AssertDimension(JxW_values.size(), n_quadrature_points);
3315  Assert(J_value != nullptr, ExcNotInitialized());
3316  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3317  {
3318  VectorizedArray<Number> J = J_value[0];
3319  for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
3320  JxW_values[q] = J * this->quadrature_weights[q];
3321  }
3322  else
3323  for (unsigned int q = 0; q < n_quadrature_points; ++q)
3324  JxW_values[q] = J_value[q];
3325 }
3326 
3327 
3328 
3329 template <int dim, int n_components_, typename Number, bool is_face>
3330 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArray<Number>>
3332  const unsigned int q_index) const
3333 {
3334  AssertIndexRange(q_index, n_quadrature_points);
3335  Assert(normal_vectors != nullptr, ExcMessage("Did not call reinit()!"));
3336  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
3337  return normal_vectors[0];
3338  else
3339  return normal_vectors[q_index];
3340 }
3341 
3342 
3343 
3344 template <int dim, int n_components_, typename Number, bool is_face>
3345 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
3347  const unsigned int q_index) const
3348 {
3349  AssertIndexRange(q_index, n_quadrature_points);
3350  Assert(J_value != nullptr, ExcNotInitialized());
3351  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3352  {
3353  Assert(this->quadrature_weights != nullptr, ExcInternalError());
3354  return J_value[0] * this->quadrature_weights[q_index];
3355  }
3356  else
3357  return J_value[q_index];
3358 }
3359 
3360 
3361 
3362 template <int dim, int n_components_, typename Number, bool is_face>
3365  const unsigned int q_index) const
3366 {
3367  AssertIndexRange(q_index, n_quadrature_points);
3368  Assert(this->jacobian != nullptr, ExcNotImplemented());
3369  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
3370  return jacobian[0];
3371  else
3372  return jacobian[q_index];
3373 }
3374 
3375 
3376 
3377 template <int dim, int n_components_, typename Number, bool is_face>
3380  const AlignedVector<VectorizedArray<Number>> &array) const
3381 {
3382  Assert(matrix_info != nullptr, ExcNotImplemented());
3383  AssertDimension(array.size(),
3384  matrix_info->get_task_info().cell_partition_data.back());
3385  if (is_face)
3386  {
3387  VectorizedArray<Number> out = make_vectorized_array<Number>(Number(1.));
3388  const unsigned int * cells =
3389  is_interior_face ?
3390  &this->matrix_info->get_face_info(cell).cells_interior[0] :
3391  &this->matrix_info->get_face_info(cell).cells_exterior[0];
3392  for (unsigned int i = 0; i < VectorizedArray<Number>::n_array_elements;
3393  ++i)
3394  if (cells[i] != numbers::invalid_unsigned_int)
3395  out[i] = array[cells[i] / VectorizedArray<Number>::n_array_elements]
3397  return out;
3398  }
3399  else
3400  return array[cell];
3401 }
3402 
3403 
3404 
3405 namespace internal
3406 {
3407  // below we use type-traits from matrix-free/type_traits.h
3408 
3409  // access to generic const vectors that have operator ().
3410  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
3411  // where we should first do Partitioner::local_to_global()
3412  template <typename VectorType,
3413  typename std::enable_if<!has_local_element<VectorType>::value,
3414  VectorType>::type * = nullptr>
3415  inline typename VectorType::value_type
3416  vector_access(const VectorType &vec, const unsigned int entry)
3417  {
3418  return vec(entry);
3419  }
3420 
3421 
3422 
3423  // access to generic non-const vectors that have operator ().
3424  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
3425  // where we should first do Partitioner::local_to_global()
3426  template <typename VectorType,
3427  typename std::enable_if<!has_local_element<VectorType>::value,
3428  VectorType>::type * = nullptr>
3429  inline typename VectorType::value_type &
3430  vector_access(VectorType &vec, const unsigned int entry)
3431  {
3432  return vec(entry);
3433  }
3434 
3435 
3436 
3437  // access to distributed MPI vectors that have a local_element(uint)
3438  // method to access data in local index space, which is what we use in
3439  // DoFInfo and hence in read_dof_values etc.
3440  template <typename VectorType,
3441  typename std::enable_if<has_local_element<VectorType>::value,
3442  VectorType>::type * = nullptr>
3443  inline typename VectorType::value_type &
3444  vector_access(VectorType &vec, const unsigned int entry)
3445  {
3446  return vec.local_element(entry);
3447  }
3448 
3449 
3450 
3451  // same for const access
3452  template <typename VectorType,
3453  typename std::enable_if<has_local_element<VectorType>::value,
3454  VectorType>::type * = nullptr>
3455  inline typename VectorType::value_type
3456  vector_access(const VectorType &vec, const unsigned int entry)
3457  {
3458  return vec.local_element(entry);
3459  }
3460 
3461 
3462 
3463  template <typename VectorType,
3464  typename std::enable_if<has_add_local_element<VectorType>::value,
3465  VectorType>::type * = nullptr>
3466  inline void
3467  vector_access_add(VectorType & vec,
3468  const unsigned int entry,
3469  const typename VectorType::value_type &val)
3470  {
3471  vec.add_local_element(entry, val);
3472  }
3473 
3474 
3475 
3476  template <typename VectorType,
3477  typename std::enable_if<!has_add_local_element<VectorType>::value,
3478  VectorType>::type * = nullptr>
3479  inline void
3480  vector_access_add(VectorType & vec,
3481  const unsigned int entry,
3482  const typename VectorType::value_type &val)
3483  {
3484  vector_access(vec, entry) += val;
3485  }
3486 
3487 
3488 
3489  template <typename VectorType,
3490  typename std::enable_if<has_add_local_element<VectorType>::value,
3491  VectorType>::type * = nullptr>
3492  inline void
3493  vector_access_add_global(VectorType & vec,
3494  const types::global_dof_index entry,
3495  const typename VectorType::value_type &val)
3496  {
3497  vec.add(entry, val);
3498  }
3499 
3500 
3501 
3502  template <typename VectorType,
3503  typename std::enable_if<!has_add_local_element<VectorType>::value,
3504  VectorType>::type * = nullptr>
3505  inline void
3506  vector_access_add_global(VectorType & vec,
3507  const types::global_dof_index entry,
3508  const typename VectorType::value_type &val)
3509  {
3510  vec(entry) += val;
3511  }
3512 
3513 
3514 
3515  template <typename VectorType,
3516  typename std::enable_if<has_set_local_element<VectorType>::value,
3517  VectorType>::type * = nullptr>
3518  inline void
3519  vector_access_set(VectorType & vec,
3520  const unsigned int entry,
3521  const typename VectorType::value_type &val)
3522  {
3523  vec.set_local_element(entry, val);
3524  }
3525 
3526 
3527 
3528  template <typename VectorType,
3529  typename std::enable_if<!has_set_local_element<VectorType>::value,
3530  VectorType>::type * = nullptr>
3531  inline void
3532  vector_access_set(VectorType & vec,
3533  const unsigned int entry,
3534  const typename VectorType::value_type &val)
3535  {
3536  vector_access(vec, entry) = val;
3537  }
3538 
3539 
3540 
3541  // this is to make sure that the parallel partitioning in VectorType
3542  // is really the same as stored in MatrixFree.
3543  // version below is when has_partitioners_are_compatible == false
3544  // FIXME: this is incorrect for PETSc/Trilinos MPI vectors
3545  template <
3546  typename VectorType,
3547  typename std::enable_if<!has_partitioners_are_compatible<VectorType>::value,
3548  VectorType>::type * = nullptr>
3549  inline void
3550  check_vector_compatibility(
3551  const VectorType & vec,
3553  {
3554  (void)vec;
3555  (void)dof_info;
3556 
3557  AssertDimension(vec.size(), dof_info.vector_partitioner->size());
3558  }
3559 
3560 
3561 
3562  // same as above for has_partitioners_are_compatible == true
3563  template <
3564  typename VectorType,
3565  typename std::enable_if<has_partitioners_are_compatible<VectorType>::value,
3566  VectorType>::type * = nullptr>
3567  inline void
3568  check_vector_compatibility(
3569  const VectorType & vec,
3571  {
3572  (void)vec;
3573  (void)dof_info;
3574  Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
3575  ExcMessage(
3576  "The parallel layout of the given vector is not "
3577  "compatible with the parallel partitioning in MatrixFree. "
3578  "Use MatrixFree::initialize_dof_vector to get a "
3579  "compatible vector."));
3580  }
3581 
3582 
3583 
3584  // Below, three classes (VectorReader, VectorSetter,
3585  // VectorDistributorLocalToGlobal) implement the same interface and can be
3586  // used to to read from vector, set elements of a vector and add to elements
3587  // of the vector.
3588 
3589  // 1. A class to read data from vector
3590  template <typename Number>
3591  struct VectorReader
3592  {
3593  template <typename VectorType>
3594  void
3595  process_dof(const unsigned int index,
3596  const VectorType & vec,
3597  Number & res) const
3598  {
3599  res = vector_access(vec, index);
3600  }
3601 
3602 
3603 
3604  template <typename VectorType>
3605  void
3606  process_dofs_vectorized(const unsigned int dofs_per_cell,
3607  const unsigned int dof_index,
3608  VectorType & vec,
3609  VectorizedArray<Number> *dof_values,
3610  std::integral_constant<bool, true>) const
3611  {
3612  const Number *vec_ptr = vec.begin() + dof_index;
3613  for (unsigned int i = 0; i < dofs_per_cell;
3615  dof_values[i].load(vec_ptr);
3616  }
3617 
3618 
3619 
3620  template <typename VectorType>
3621  void
3622  process_dofs_vectorized(const unsigned int dofs_per_cell,
3623  const unsigned int dof_index,
3624  const VectorType & vec,
3625  VectorizedArray<Number> *dof_values,
3626  std::integral_constant<bool, false>) const
3627  {
3628  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3629  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3630  ++v)
3631  dof_values[i][v] = vector_access(
3632  vec, dof_index + v + i * VectorizedArray<Number>::n_array_elements);
3633  }
3634 
3635 
3636 
3637  template <typename VectorType>
3638  void
3639  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3640  const unsigned int * dof_indices,
3641  VectorType & vec,
3642  VectorizedArray<Number> *dof_values,
3643  std::integral_constant<bool, true>) const
3644  {
3645  ::vectorized_load_and_transpose(dofs_per_cell,
3646  vec.begin(),
3647  dof_indices,
3648  dof_values);
3649  }
3650 
3651 
3652 
3653  template <typename VectorType>
3654  void
3655  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3656  const unsigned int * dof_indices,
3657  const VectorType & vec,
3658  VectorizedArray<Number> *dof_values,
3659  std::integral_constant<bool, false>) const
3660  {
3661  for (unsigned int d = 0; d < dofs_per_cell; ++d)
3662  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3663  ++v)
3664  dof_values[d][v] = vector_access(vec, dof_indices[v] + d);
3665  }
3666 
3667 
3668 
3669  // variant where VectorType::value_type is the same as Number -> can call
3670  // gather
3671  template <typename VectorType>
3672  void
3673  process_dof_gather(const unsigned int * indices,
3674  VectorType & vec,
3675  const unsigned int constant_offset,
3677  std::integral_constant<bool, true>) const
3678  {
3679  res.gather(vec.begin() + constant_offset, indices);
3680  }
3681 
3682 
3683 
3684  // variant where VectorType::value_type is not the same as Number -> must
3685  // manually load the data
3686  template <typename VectorType>
3687  void
3688  process_dof_gather(const unsigned int * indices,
3689  const VectorType & vec,
3690  const unsigned int constant_offset,
3692  std::integral_constant<bool, false>) const
3693  {
3694  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3695  ++v)
3696  res[v] = vector_access(vec, indices[v] + constant_offset);
3697  }
3698 
3699 
3700 
3701  template <typename VectorType>
3702  void
3703  process_dof_global(const types::global_dof_index index,
3704  const VectorType & vec,
3705  Number & res) const
3706  {
3707  res = vec(index);
3708  }
3709 
3710 
3711 
3712  void
3713  pre_constraints(const Number &, Number &res) const
3714  {
3715  res = Number();
3716  }
3717 
3718 
3719 
3720  template <typename VectorType>
3721  void
3722  process_constraint(const unsigned int index,
3723  const Number weight,
3724  const VectorType & vec,
3725  Number & res) const
3726  {
3727  res += weight * vector_access(vec, index);
3728  }
3729 
3730 
3731 
3732  void
3733  post_constraints(const Number &sum, Number &write_pos) const
3734  {
3735  write_pos = sum;
3736  }
3737 
3738 
3739 
3740  void
3741  process_empty(VectorizedArray<Number> &res) const
3742  {
3743  res = VectorizedArray<Number>();
3744  }
3745  };
3746 
3747 
3748 
3749  // 2. A class to add values to the vector during
3750  // FEEvaluation::distribute_local_to_global() call
3751  template <typename Number>
3752  struct VectorDistributorLocalToGlobal
3753  {
3754  template <typename VectorType>
3755  void
3756  process_dof(const unsigned int index, VectorType &vec, Number &res) const
3757  {
3758  vector_access_add(vec, index, res);
3759  }
3760 
3761 
3762 
3763  template <typename VectorType>
3764  void
3765  process_dofs_vectorized(const unsigned int dofs_per_cell,
3766  const unsigned int dof_index,
3767  VectorType & vec,
3768  VectorizedArray<Number> *dof_values,
3769  std::integral_constant<bool, true>) const
3770  {
3771  Number *vec_ptr = vec.begin() + dof_index;
3772  for (unsigned int i = 0; i < dofs_per_cell;
3774  {
3776  tmp.load(vec_ptr);
3777  tmp += dof_values[i];
3778  tmp.store(vec_ptr);
3779  }
3780  }
3781 
3782 
3783 
3784  template <typename VectorType>
3785  void
3786  process_dofs_vectorized(const unsigned int dofs_per_cell,
3787  const unsigned int dof_index,
3788  VectorType & vec,
3789  VectorizedArray<Number> *dof_values,
3790  std::integral_constant<bool, false>) const
3791  {
3792  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3793  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3794  ++v)
3795  vector_access_add(vec,
3796  dof_index + v +
3798  dof_values[i][v]);
3799  }
3800 
3801 
3802 
3803  template <typename VectorType>
3804  void
3805  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3806  const unsigned int * dof_indices,
3807  VectorType & vec,
3808  VectorizedArray<Number> *dof_values,
3809  std::integral_constant<bool, true>) const
3810  {
3811  vectorized_transpose_and_store(
3812  true, dofs_per_cell, dof_values, dof_indices, vec.begin());
3813  }
3814 
3815 
3816 
3817  template <typename VectorType>
3818  void
3819  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3820  const unsigned int * dof_indices,
3821  VectorType & vec,
3822  VectorizedArray<Number> *dof_values,
3823  std::integral_constant<bool, false>) const
3824  {
3825  for (unsigned int d = 0; d < dofs_per_cell; ++d)
3826  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3827  ++v)
3828  vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]);
3829  }
3830 
3831 
3832 
3833  // variant where VectorType::value_type is the same as Number -> can call
3834  // scatter
3835  template <typename VectorType>
3836  void
3837  process_dof_gather(const unsigned int * indices,
3838  VectorType & vec,
3839  const unsigned int constant_offset,
3841  std::integral_constant<bool, true>) const
3842  {
3843 # if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
3844  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3845  ++v)
3846  vector_access(vec, indices[v] + constant_offset) += res[v];
3847 # else
3848  // only use gather in case there is also scatter.
3850  tmp.gather(vec.begin() + constant_offset, indices);
3851  tmp += res;
3852  tmp.scatter(indices, vec.begin() + constant_offset);
3853 # endif
3854  }
3855 
3856 
3857 
3858  // variant where VectorType::value_type is not the same as Number -> must
3859  // manually append all data
3860  template <typename VectorType>
3861  void
3862  process_dof_gather(const unsigned int * indices,
3863  VectorType & vec,
3864  const unsigned int constant_offset,
3866  std::integral_constant<bool, false>) const
3867  {
3868  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3869  ++v)
3870  vector_access_add(vec, indices[v] + constant_offset, res[v]);
3871  }
3872 
3873 
3874 
3875  template <typename VectorType>
3876  void
3877  process_dof_global(const types::global_dof_index index,
3878  VectorType & vec,
3879  Number & res) const
3880  {
3881  vector_access_add_global(vec, index, res);
3882  }
3883 
3884 
3885 
3886  void
3887  pre_constraints(const Number &input, Number &res) const
3888  {
3889  res = input;
3890  }
3891 
3892 
3893 
3894  template <typename VectorType>
3895  void
3896  process_constraint(const unsigned int index,
3897  const Number weight,
3898  VectorType & vec,
3899  Number & res) const
3900  {
3901  vector_access_add(vec, index, weight * res);
3902  }
3903 
3904 
3905 
3906  void
3907  post_constraints(const Number &, Number &) const
3908  {}
3909 
3910 
3911 
3912  void
3913  process_empty(VectorizedArray<Number> &) const
3914  {}
3915  };
3916 
3917 
3918 
3919  // 3. A class to set elements of the vector
3920  template <typename Number>
3921  struct VectorSetter
3922  {
3923  template <typename VectorType>
3924  void
3925  process_dof(const unsigned int index, VectorType &vec, Number &res) const
3926  {
3927  vector_access(vec, index) = res;
3928  }
3929 
3930 
3931 
3932  template <typename VectorType>
3933  void
3934  process_dofs_vectorized(const unsigned int dofs_per_cell,
3935  const unsigned int dof_index,
3936  VectorType & vec,
3937  VectorizedArray<Number> *dof_values,
3938  std::integral_constant<bool, true>) const
3939  {
3940  Number *vec_ptr = vec.begin() + dof_index;
3941  for (unsigned int i = 0; i < dofs_per_cell;
3943  dof_values[i].store(vec_ptr);
3944  }
3945 
3946 
3947 
3948  template <typename VectorType>
3949  void
3950  process_dofs_vectorized(const unsigned int dofs_per_cell,
3951  const unsigned int dof_index,
3952  VectorType & vec,
3953  VectorizedArray<Number> *dof_values,
3954  std::integral_constant<bool, false>) const
3955  {
3956  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3957  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3958  ++v)
3959  vector_access(vec,
3960  dof_index + v +
3962  dof_values[i][v];
3963  }
3964 
3965 
3966 
3967  template <typename VectorType>
3968  void
3969  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3970  const unsigned int * dof_indices,
3971  VectorType & vec,
3972  VectorizedArray<Number> *dof_values,
3973  std::integral_constant<bool, true>) const
3974  {
3975  vectorized_transpose_and_store(
3976  false, dofs_per_cell, dof_values, dof_indices, vec.begin());
3977  }
3978 
3979 
3980 
3981  template <typename VectorType, bool booltype>
3982  void
3983  process_dofs_vectorized_transpose(const unsigned int dofs_per_cell,
3984  const unsigned int * dof_indices,
3985  VectorType & vec,
3986  VectorizedArray<Number> *dof_values,
3987  std::integral_constant<bool, false>) const
3988  {
3989  for (unsigned int i = 0; i < dofs_per_cell; ++i)
3990  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
3991  ++v)
3992  vector_access(vec, dof_indices[v] + i) = dof_values[i][v];
3993  }
3994 
3995 
3996 
3997  template <typename VectorType>
3998  void
3999  process_dof_gather(const unsigned int * indices,
4000  VectorType & vec,
4001  const unsigned int constant_offset,
4003  std::integral_constant<bool, true>) const
4004  {
4005  res.scatter(indices, vec.begin() + constant_offset);
4006  }
4007 
4008 
4009 
4010  template <typename VectorType>
4011  void
4012  process_dof_gather(const unsigned int * indices,
4013  VectorType & vec,
4014  const unsigned int constant_offset,
4016  std::integral_constant<bool, false>) const
4017  {
4018  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
4019  ++v)
4020  vector_access(vec, indices[v] + constant_offset) = res[v];
4021  }
4022 
4023 
4024 
4025  template <typename VectorType>
4026  void
4027  process_dof_global(const types::global_dof_index index,
4028  VectorType & vec,
4029  Number & res) const
4030  {
4031  vec(index) = res;
4032  }
4033 
4034 
4035 
4036  void
4037  pre_constraints(const Number &, Number &) const
4038  {}
4039 
4040 
4041 
4042  template <typename VectorType>
4043  void
4044  process_constraint(const unsigned int,
4045  const Number,
4046  VectorType &,
4047  Number &) const
4048  {}
4049 
4050 
4051 
4052  void
4053  post_constraints(const Number &, Number &) const
4054  {}
4055 
4056 
4057 
4058  void
4059  process_empty(VectorizedArray<Number> &) const
4060  {}
4061  };
4062 
4063 
4064 
4065  // allows to select between block vectors and non-block vectors, which
4066  // allows to use a unified interface for extracting blocks on block vectors
4067  // and doing nothing on usual vectors
4068  template <typename VectorType, bool>
4069  struct BlockVectorSelector
4070  {};
4071 
4072  template <typename VectorType>
4073  struct BlockVectorSelector<VectorType, true>
4074  {
4075  using BaseVectorType = typename VectorType::BlockType;
4076 
4077  static BaseVectorType *
4078  get_vector_component(VectorType &vec, const unsigned int component)
4079  {
4080  AssertIndexRange(component, vec.n_blocks());
4081  return &vec.block(component);
4082  }
4083  };
4084 
4085  template <typename VectorType>
4086  struct BlockVectorSelector<VectorType, false>
4087  {
4088  using BaseVectorType = VectorType;
4089 
4090  static BaseVectorType *
4091  get_vector_component(VectorType &vec, const unsigned int component)
4092  {
4093  // FEEvaluation allows to combine several vectors from a scalar
4094  // FiniteElement into a "vector-valued" FEEvaluation object with
4095  // multiple components. These components can be extracted with the other
4096  // get_vector_component functions. If we do not get a vector of vectors
4097  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
4098  // must make sure that we do not duplicate the components in input
4099  // and/or duplicate the resulting integrals. In such a case, we should
4100  // only get the zeroth component in the vector contained set nullptr for
4101  // the others which allows us to catch unintended use in
4102  // read_write_operation.
4103  if (component == 0)
4104  return &vec;
4105  else
4106  return nullptr;
4107  }
4108  };
4109 
4110  template <typename VectorType>
4111  struct BlockVectorSelector<std::vector<VectorType>, false>
4112  {
4113  using BaseVectorType = VectorType;
4114 
4115  static BaseVectorType *
4116  get_vector_component(std::vector<VectorType> &vec,
4117  const unsigned int component)
4118  {
4119  AssertIndexRange(component, vec.size());
4120  return &vec[component];
4121  }
4122  };
4123 
4124  template <typename VectorType>
4125  struct BlockVectorSelector<std::vector<VectorType *>, false>
4126  {
4127  using BaseVectorType = VectorType;
4128 
4129  static BaseVectorType *
4130  get_vector_component(std::vector<VectorType *> &vec,
4131  const unsigned int component)
4132  {
4133  AssertIndexRange(component, vec.size());
4134  return vec[component];
4135  }
4136  };
4137 } // namespace internal
4138 
4139 
4140 
4141 template <int dim, int n_components_, typename Number, bool is_face>
4142 template <typename VectorType, typename VectorOperation>
4143 inline void
4145  const VectorOperation & operation,
4146  VectorType * src[],
4147  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask,
4148  const bool apply_constraints) const
4149 {
4150  // Case 1: No MatrixFree object given, simple case because we do not need to
4151  // process constraints and need not care about vectorization -> go to
4152  // separate function
4153  if (matrix_info == nullptr)
4154  {
4155  read_write_operation_global(operation, src);
4156  return;
4157  }
4158 
4159  Assert(dof_info != nullptr, ExcNotInitialized());
4160  Assert(matrix_info->indices_initialized() == true, ExcNotInitialized());
4161  if (n_fe_components == 1)
4162  for (unsigned int comp = 0; comp < n_components; ++comp)
4163  internal::check_vector_compatibility(*src[comp], *dof_info);
4164  else
4165  {
4166  internal::check_vector_compatibility(*src[0], *dof_info);
4167  }
4168 
4169  // Case 2: contiguous indices which use reduced storage of indices and can
4170  // use vectorized load/store operations -> go to separate function
4171  AssertIndexRange(cell,
4172  dof_info->index_storage_variants[dof_access_index].size());
4173  if (dof_info->index_storage_variants
4174  [is_face ? dof_access_index :
4176  [cell] >=
4178  {
4179  read_write_operation_contiguous(operation, src, mask);
4180  return;
4181  }
4182 
4183  // Case 3: standard operation with one index per degree of freedom -> go on
4184  // here
4185  constexpr unsigned int n_vectorization =
4187  Assert(mask.count() == n_vectorization,
4188  ExcNotImplemented("Masking currently not implemented for "
4189  "non-contiguous DoF storage"));
4190 
4191  std::integral_constant<bool,
4192  internal::is_vectorizable<VectorType, Number>::value>
4193  vector_selector;
4194 
4195  const unsigned int dofs_per_component =
4196  this->data->dofs_per_component_on_cell;
4197  if (dof_info->index_storage_variants
4198  [is_face ? dof_access_index :
4200  [cell] ==
4202  {
4203  const unsigned int *dof_indices =
4204  dof_info->dof_indices_interleaved.data() +
4205  dof_info->row_starts[cell * n_fe_components * n_vectorization].first +
4206  dof_info->component_dof_indices_offset[active_fe_index]
4207  [first_selected_component] *
4208  n_vectorization;
4209  if (n_components == 1 || n_fe_components == 1)
4210  for (unsigned int i = 0; i < dofs_per_component;
4211  ++i, dof_indices += n_vectorization)
4212  for (unsigned int comp = 0; comp < n_components; ++comp)
4213  operation.process_dof_gather(dof_indices,
4214  *src[comp],
4215  0,
4216  values_dofs[comp][i],
4217  vector_selector);
4218  else
4219  for (unsigned int comp = 0; comp < n_components; ++comp)
4220  for (unsigned int i = 0; i < dofs_per_component;
4221  ++i, dof_indices += n_vectorization)
4222  operation.process_dof_gather(
4223  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
4224  return;
4225  }
4226 
4227  const unsigned int * dof_indices[n_vectorization];
4228  VectorizedArray<Number> **values_dofs =
4229  const_cast<VectorizedArray<Number> **>(&this->values_dofs[0]);
4230 
4231  unsigned int cells_copied[n_vectorization];
4232  const unsigned int *cells;
4233  unsigned int n_vectorization_actual =
4234  dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
4235  bool has_constraints = false;
4236  if (is_face)
4237  {
4238  if (dof_access_index ==
4240  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4241  cells_copied[v] =
4243  cells = dof_access_index ==
4245  &cells_copied[0] :
4246  (is_interior_face ?
4247  &this->matrix_info->get_face_info(cell).cells_interior[0] :
4248  &this->matrix_info->get_face_info(cell).cells_exterior[0]);
4249  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4250  {
4251  Assert(cells[v] < dof_info->row_starts.size() - 1,
4252  ExcInternalError());
4253  has_constraints =
4254  has_constraints &&
4255  dof_info
4256  ->row_starts[cells[v] * n_fe_components +
4257  first_selected_component + n_components]
4258  .second != dof_info
4259  ->row_starts[cells[v] * n_fe_components +
4260  first_selected_component]
4261  .second;
4262  dof_indices[v] = dof_info->dof_indices.data() +
4263  dof_info
4264  ->row_starts[cells[v] * n_fe_components +
4265  first_selected_component]
4266  .first;
4267  }
4268  for (unsigned int v = n_vectorization_actual; v < n_vectorization; ++v)
4269  dof_indices[v] = nullptr;
4270  }
4271  else
4272  {
4273  AssertIndexRange((cell + 1) * n_vectorization * n_fe_components,
4274  dof_info->row_starts.size());
4275  const unsigned int n_components_read =
4276  n_fe_components > 1 ? n_components : 1;
4277  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4278  {
4279  if (dof_info
4280  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4281  first_selected_component + n_components_read]
4282  .second !=
4283  dof_info
4284  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4285  first_selected_component]
4286  .second)
4287  has_constraints = true;
4288  Assert(
4289  dof_info
4290  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4291  first_selected_component + n_components_read]
4292  .first ==
4293  dof_info
4294  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4295  first_selected_component]
4296  .first ||
4297  dof_info
4298  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4299  first_selected_component]
4300  .first < dof_info->dof_indices.size(),
4301  ExcIndexRange(
4302  0,
4303  dof_info
4304  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4305  first_selected_component]
4306  .first,
4307  dof_info->dof_indices.size()));
4308  dof_indices[v] =
4309  dof_info->dof_indices.data() +
4310  dof_info
4311  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4312  first_selected_component]
4313  .first;
4314  }
4315  for (unsigned int v = n_vectorization_actual; v < n_vectorization; ++v)
4316  dof_indices[v] = nullptr;
4317  }
4318 
4319  // Case where we have no constraints throughout the whole cell: Can go
4320  // through the list of DoFs directly
4321  if (!has_constraints)
4322  {
4323  if (n_vectorization_actual < n_vectorization)
4324  for (unsigned int comp = 0; comp < n_components; ++comp)
4325  for (unsigned int i = 0; i < dofs_per_component; ++i)
4326  operation.process_empty(values_dofs[comp][i]);
4327  if (n_components == 1 || n_fe_components == 1)
4328  {
4329  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4330  for (unsigned int i = 0; i < dofs_per_component; ++i)
4331  for (unsigned int comp = 0; comp < n_components; ++comp)
4332  operation.process_dof(dof_indices[v][i],
4333  *src[comp],
4334  values_dofs[comp][i][v]);
4335  }
4336  else
4337  {
4338  for (unsigned int comp = 0; comp < n_components; ++comp)
4339  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4340  for (unsigned int i = 0; i < dofs_per_component; ++i)
4341  operation.process_dof(
4342  dof_indices[v][comp * dofs_per_component + i],
4343  *src[0],
4344  values_dofs[comp][i][v]);
4345  }
4346  return;
4347  }
4348 
4349  // In the case where there are some constraints to be resolved, loop over
4350  // all vector components that are filled and then over local dofs. ind_local
4351  // holds local number on cell, index iterates over the elements of
4352  // index_local_to_global and dof_indices points to the global indices stored
4353  // in index_local_to_global
4354  if (n_vectorization_actual < n_vectorization)
4355  for (unsigned int comp = 0; comp < n_components; ++comp)
4356  for (unsigned int i = 0; i < dofs_per_component; ++i)
4357  operation.process_empty(values_dofs[comp][i]);
4358  for (unsigned int v = 0; v < n_vectorization_actual; ++v)
4359  {
4360  unsigned int index_indicators, next_index_indicators;
4361  const unsigned int n_components_read =
4362  n_fe_components > 1 ? n_components : 1;
4363  if (is_face)
4364  {
4365  index_indicators = dof_info
4366  ->row_starts[cells[v] * n_fe_components +
4367  first_selected_component]
4368  .second;
4369  next_index_indicators = dof_info
4370  ->row_starts[cells[v] * n_fe_components +
4371  first_selected_component + 1]
4372  .second;
4373  }
4374  else
4375  {
4376  index_indicators =
4377  dof_info
4378  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4379  first_selected_component]
4380  .second;
4381  next_index_indicators =
4382  dof_info
4383  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4384  first_selected_component + 1]
4385  .second;
4386  }
4387 
4388  if (apply_constraints == false &&
4389  dof_info
4390  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4391  first_selected_component]
4392  .second !=
4393  dof_info
4394  ->row_starts[(cell * n_vectorization + v) * n_fe_components +
4395  first_selected_component + n_components_read]
4396  .second)
4397  {
4398  Assert(
4399  dof_info->row_starts_plain_indices[cell * n_vectorization + v] !=
4401  ExcNotInitialized());
4402  dof_indices[v] =
4403  dof_info->plain_dof_indices.data() +
4404  dof_info->component_dof_indices_offset[active_fe_index]
4405  [first_selected_component] +
4406  (is_face ?
4407  dof_info->row_starts_plain_indices[cells[v]] :
4408  dof_info->row_starts_plain_indices[cell * n_vectorization + v]);
4409  next_index_indicators = index_indicators;
4410  }
4411 
4412  if (n_components == 1 || n_fe_components == 1)
4413  {
4414  for (unsigned int c = 0; c < n_components; ++c)
4415  Assert(src[c] != nullptr,
4416  ExcMessage(
4417  "The finite element underlying this FEEvaluation "
4418  "object is scalar, but you requested " +
4419  std::to_string(n_components) +
4420  " components via the template argument in "
4421  "FEEvaluation. In that case, you must pass an "
4422  "std::vector<VectorType> or a BlockVector to " +
4423  "read_dof_values and distribute_local_to_global."));
4424 
4425  unsigned int ind_local = 0;
4426  for (; index_indicators != next_index_indicators; ++index_indicators)
4427  {
4428  const std::pair<unsigned short, unsigned short> indicator =
4429  dof_info->constraint_indicator[index_indicators];
4430  // run through values up to next constraint
4431  for (unsigned int j = 0; j < indicator.first; ++j)
4432  for (unsigned int comp = 0; comp < n_components; ++comp)
4433  operation.process_dof(dof_indices[v][j],
4434  *src[comp],
4435  values_dofs[comp][ind_local + j][v]);
4436 
4437  ind_local += indicator.first;
4438  dof_indices[v] += indicator.first;
4439 
4440  // constrained case: build the local value as a linear
4441  // combination of the global value according to constraints
4442  Number value[n_components];
4443  for (unsigned int comp = 0; comp < n_components; ++comp)
4444  operation.pre_constraints(values_dofs[comp][ind_local][v],
4445  value[comp]);
4446 
4447  const Number *data_val =
4448  matrix_info->constraint_pool_begin(indicator.second);
4449  const Number *end_pool =
4450  matrix_info->constraint_pool_end(indicator.second);
4451  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4452  for (unsigned int comp = 0; comp < n_components; ++comp)
4453  operation.process_constraint(*dof_indices[v],
4454  *data_val,
4455  *src[comp],
4456  value[comp]);
4457 
4458  for (unsigned int comp = 0; comp < n_components; ++comp)
4459  operation.post_constraints(value[comp],
4460  values_dofs[comp][ind_local][v]);
4461  ind_local++;
4462  }
4463 
4464  AssertIndexRange(ind_local, dofs_per_component + 1);
4465 
4466  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
4467  for (unsigned int comp = 0; comp < n_components; ++comp)
4468  operation.process_dof(*dof_indices[v],
4469  *src[comp],
4470  values_dofs[comp][ind_local][v]);
4471  }
4472  else
4473  {
4474  // case with vector-valued finite elements where all components are
4475  // included in one single vector. Assumption: first come all entries
4476  // to the first component, then all entries to the second one, and
4477  // so on. This is ensured by the way MatrixFree reads out the
4478  // indices.
4479  for (unsigned int comp = 0; comp < n_components; ++comp)
4480  {
4481  unsigned int ind_local = 0;
4482 
4483  // check whether there is any constraint on the current cell
4484  for (; index_indicators != next_index_indicators;
4485  ++index_indicators)
4486  {
4487  const std::pair<unsigned short, unsigned short> indicator =
4488  dof_info->constraint_indicator[index_indicators];
4489 
4490  // run through values up to next constraint
4491  for (unsigned int j = 0; j < indicator.first; ++j)
4492  operation.process_dof(dof_indices[v][j],
4493  *src[0],
4494  values_dofs[comp][ind_local + j][v]);
4495  ind_local += indicator.first;
4496  dof_indices[v] += indicator.first;
4497 
4498  // constrained case: build the local value as a linear
4499  // combination of the global value according to constraints
4500  Number value;
4501  operation.pre_constraints(values_dofs[comp][ind_local][v],
4502  value);
4503 
4504  const Number *data_val =
4505  matrix_info->constraint_pool_begin(indicator.second);
4506  const Number *end_pool =
4507  matrix_info->constraint_pool_end(indicator.second);
4508 
4509  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
4510  operation.process_constraint(*dof_indices[v],
4511  *data_val,
4512  *src[0],
4513  value);
4514 
4515  operation.post_constraints(value,
4516  values_dofs[comp][ind_local][v]);
4517  ind_local++;
4518  }
4519 
4520  AssertIndexRange(ind_local, dofs_per_component + 1);
4521 
4522  // get the dof values past the last constraint
4523  for (; ind_local < dofs_per_component;
4524  ++dof_indices[v], ++ind_local)
4525  {
4526  AssertIndexRange(*dof_indices[v], src[0]->size());
4527  operation.process_dof(*dof_indices[v],
4528  *src[0],
4529  values_dofs[comp][ind_local][v]);
4530  }
4531 
4532  if (apply_constraints == true && comp + 1 < n_components)
4533  {
4534  if (is_face)
4535  next_index_indicators =
4536  dof_info
4537  ->row_starts[cells[v] * n_fe_components +
4538  first_selected_component + comp + 2]
4539  .second;
4540  else
4541  next_index_indicators =
4542  dof_info
4543  ->row_starts[(cell * n_vectorization + v) *
4544  n_fe_components +
4545  first_selected_component + comp + 2]
4546  .second;
4547  }
4548  }
4549  }
4550  }
4551 }
4552 
4553 
4554 
4555 template <int dim, int n_components_, typename Number, bool is_face>
4556 template <typename VectorType, typename VectorOperation>
4557 inline void
4560  VectorType * src[]) const
4561 {
4562  Assert(!local_dof_indices.empty(), ExcNotInitialized());
4563 
4564  unsigned int index =
4565  first_selected_component * data->dofs_per_component_on_cell;
4566  for (unsigned int comp = 0; comp < n_components; ++comp)
4567  {
4568  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4569  ++i, ++index)
4570  {
4571  operation.process_empty(values_dofs[comp][i]);
4572  operation.process_dof_global(
4573  local_dof_indices[data->lexicographic_numbering[index]],
4574  *src[0],
4575  values_dofs[comp][i][0]);
4576  }
4577  }
4578 }
4579 
4580 
4581 
4582 template <int dim, int n_components_, typename Number, bool is_face>
4583 template <typename VectorType, typename VectorOperation>
4584 inline void
4587  const VectorOperation & operation,
4588  VectorType * src[],
4589  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask) const
4590 {
4591  // This functions processes the functions read_dof_values,
4592  // distribute_local_to_global, and set_dof_values with the same code for
4593  // contiguous cell indices (DG case). The distinction between these three
4594  // cases is made by the input VectorOperation that either reads values from
4595  // a vector and puts the data into the local data field or write local data
4596  // into the vector. Certain operations are no-ops for the given use case.
4597 
4598  std::integral_constant<bool,
4599  internal::is_vectorizable<VectorType, Number>::value>
4600  vector_selector;
4602  is_face ? dof_access_index :
4604  const unsigned int n_lanes = mask.count();
4605 
4606  const std::vector<unsigned int> &dof_indices_cont =
4607  dof_info->dof_indices_contiguous[ind];
4608 
4609  // Simple case: We have contiguous storage, so we can simply copy out the
4610  // data
4611  if (dof_info->index_storage_variants[ind][cell] ==
4613  interleaved_contiguous &&
4615  {
4616  const unsigned int dof_index =
4617  dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements] +
4618  dof_info->component_dof_indices_offset[active_fe_index]
4619  [first_selected_component] *
4621  if (n_components == 1 || n_fe_components == 1)
4622  for (unsigned int comp = 0; comp < n_components; ++comp)
4623  operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
4624  dof_index,
4625  *src[comp],
4626  values_dofs[comp],
4627  vector_selector);
4628  else
4629  operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
4630  n_components,
4631  dof_index,
4632  *src[0],
4633  values_dofs[0],
4634  vector_selector);
4635  return;
4636  }
4637 
4638  // More general case: Must go through the components one by one and apply
4639  // some transformations
4640  const unsigned int vectorization_populated =
4641  dof_info->n_vectorization_lanes_filled[ind][this->cell];
4642 
4643  unsigned int dof_indices[VectorizedArray<Number>::n_array_elements];
4644  for (unsigned int v = 0; v < vectorization_populated; ++v)
4645  dof_indices[v] =
4646  dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements + v] +
4647  dof_info->component_dof_indices_offset[active_fe_index]
4648  [first_selected_component] *
4651 
4652  for (unsigned int v = vectorization_populated;
4653  v < VectorizedArray<Number>::n_array_elements;
4654  ++v)
4655  dof_indices[v] = numbers::invalid_unsigned_int;
4656 
4657  // In the case with contiguous cell indices, we know that there are no
4658  // constraints and that the indices within each element are contiguous
4659  if (vectorization_populated == VectorizedArray<Number>::n_array_elements &&
4661  {
4662  if (dof_info->index_storage_variants[ind][cell] ==
4664  contiguous)
4665  {
4666  if (n_components == 1 || n_fe_components == 1)
4667  for (unsigned int comp = 0; comp < n_components; ++comp)
4668  operation.process_dofs_vectorized_transpose(
4670  dof_indices,
4671  *src[comp],
4672  values_dofs[comp],
4673  vector_selector);
4674  else
4675  operation.process_dofs_vectorized_transpose(
4677  dof_indices,
4678  *src[0],
4679  &values_dofs[0][0],
4680  vector_selector);
4681  }
4682  else if (dof_info->index_storage_variants[ind][cell] ==
4684  interleaved_contiguous_strided)
4685  {
4686  if (n_components == 1 || n_fe_components == 1)
4687  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4688  {
4689  for (unsigned int comp = 0; comp < n_components; ++comp)
4690  operation.process_dof_gather(
4691  dof_indices,
4692  *src[comp],
4694  values_dofs[comp][i],
4695  vector_selector);
4696  }
4697  else
4698  for (unsigned int comp = 0; comp < n_components; ++comp)
4699  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4700  ++i)
4701  {
4702  operation.process_dof_gather(
4703  dof_indices,
4704  *src[0],
4705  (comp * data->dofs_per_component_on_cell + i) *
4707  values_dofs[comp][i],
4708  vector_selector);
4709  }
4710  }
4711  else
4712  {
4713  Assert(dof_info->index_storage_variants[ind][cell] ==
4715  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4716  ExcNotImplemented());
4717  const unsigned int *offsets =
4720  if (n_components == 1 || n_fe_components == 1)
4721  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4722  {
4723  for (unsigned int comp = 0; comp < n_components; ++comp)
4724  operation.process_dof_gather(dof_indices,
4725  *src[comp],
4726  0,
4727  values_dofs[comp][i],
4728  vector_selector);
4729  DEAL_II_OPENMP_SIMD_PRAGMA
4730  for (unsigned int v = 0;
4731  v < VectorizedArray<Number>::n_array_elements;
4732  ++v)
4733  dof_indices[v] += offsets[v];
4734  }
4735  else
4736  for (unsigned int comp = 0; comp < n_components; ++comp)
4737  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
4738  ++i)
4739  {
4740  operation.process_dof_gather(dof_indices,
4741  *src[0],
4742  0,
4743  values_dofs[comp][i],
4744  vector_selector);
4745  DEAL_II_OPENMP_SIMD_PRAGMA
4746  for (unsigned int v = 0;
4747  v < VectorizedArray<Number>::n_array_elements;
4748  ++v)
4749  dof_indices[v] += offsets[v];
4750  }
4751  }
4752  }
4753  else
4754  for (unsigned int comp = 0; comp < n_components; ++comp)
4755  {
4756  for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
4757  operation.process_empty(values_dofs[comp][i]);
4758  if (dof_info->index_storage_variants[ind][cell] ==
4760  contiguous)
4761  {
4762  if (n_components == 1 || n_fe_components == 1)
4763  {
4764  for (unsigned int v = 0; v < vectorization_populated; ++v)
4765  if (mask[v] == true)
4766  for (unsigned int i = 0;
4767  i < data->dofs_per_component_on_cell;
4768  ++i)
4769  operation.process_dof(dof_indices[v] + i,
4770  *src[comp],
4771  values_dofs[comp][i][v]);
4772  }
4773  else
4774  {
4775  for (unsigned int v = 0; v < vectorization_populated; ++v)
4776  if (mask[v] == true)
4777  for (unsigned int i = 0;
4778  i < data->dofs_per_component_on_cell;
4779  ++i)
4780  operation.process_dof(
4781  dof_indices[v] + i +
4782  comp * data->dofs_per_component_on_cell,
4783  *src[0],
4784  values_dofs[comp][i][v]);
4785  }
4786  }
4787  else
4788  {
4789  const unsigned int *offsets =
4792  for (unsigned int v = 0; v < vectorization_populated; ++v)
4793  AssertIndexRange(offsets[v],
4795  if (n_components == 1 || n_fe_components == 1)
4796  for (unsigned int v = 0; v < vectorization_populated; ++v)
4797  {
4798  if (mask[v] == true)
4799  for (unsigned int i = 0;
4800  i < data->dofs_per_component_on_cell;
4801  ++i)
4802  operation.process_dof(dof_indices[v] + i * offsets[v],
4803  *src[comp],
4804  values_dofs[comp][i][v]);
4805  }
4806  else
4807  {
4808  for (unsigned int v = 0; v < vectorization_populated; ++v)
4809  if (mask[v] == true)
4810  for (unsigned int i = 0;
4811  i < data->dofs_per_component_on_cell;
4812  ++i)
4813  operation.process_dof(
4814  dof_indices[v] +
4815  (i + comp * data->dofs_per_component_on_cell) *
4816  offsets[v],
4817  *src[0],
4818  values_dofs[comp][i][v]);
4819  }
4820  }
4821  }
4822 }
4823 
4824 
4825 
4826 template <int dim, int n_components_, typename Number, bool is_face>
4827 template <typename VectorType>
4828 inline void
4830  const VectorType & src,
4831  const unsigned int first_index)
4832 {
4833  // select between block vectors and non-block vectors. Note that the number
4834  // of components is checked in the internal data
4835  typename internal::BlockVectorSelector<
4836  VectorType,
4837  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4838  for (unsigned int d = 0; d < n_components; ++d)
4839  src_data[d] =
4840  internal::BlockVectorSelector<VectorType,
4842  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4843 
4844  internal::VectorReader<Number> reader;
4845  read_write_operation(
4846  reader,
4847  src_data,
4848  std::bitset<VectorizedArray<Number>::n_array_elements>().flip(),
4849  true);
4850 
4851 # ifdef DEBUG
4852  dof_values_initialized = true;
4853 # endif
4854 }
4855 
4856 
4857 
4858 template <int dim, int n_components_, typename Number, bool is_face>
4859 template <typename VectorType>
4860 inline void
4862  const VectorType & src,
4863  const unsigned int first_index)
4864 {
4865  // select between block vectors and non-block vectors. Note that the number
4866  // of components is checked in the internal data
4867  typename internal::BlockVectorSelector<
4868  VectorType,
4869  IsBlockVector<VectorType>::value>::BaseVectorType *src_data[n_components];
4870  for (unsigned int d = 0; d < n_components; ++d)
4871  src_data[d] =
4872  internal::BlockVectorSelector<VectorType,
4874  get_vector_component(const_cast<VectorType &>(src), d + first_index);
4875 
4876  internal::VectorReader<Number> reader;
4877  read_write_operation(
4878  reader,
4879  src_data,
4880  std::bitset<VectorizedArray<Number>::n_array_elements>().flip(),
4881  false);
4882 
4883 # ifdef DEBUG
4884  dof_values_initialized = true;
4885 # endif
4886 }
4887 
4888 
4889 
4890 template <int dim, int n_components_, typename Number, bool is_face>
4891 template <typename VectorType>
4892 inline void
4895  VectorType & dst,
4896  const unsigned int first_index,
4897  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask) const
4898 {
4899  Assert(dof_values_initialized == true,
4901 
4902  // select between block vectors and non-block vectors. Note that the number
4903  // of components is checked in the internal data
4904  typename internal::BlockVectorSelector<
4905  VectorType,
4906  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4907  for (unsigned int d = 0; d < n_components; ++d)
4908  dst_data[d] = internal::BlockVectorSelector<
4909  VectorType,
4910  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4911  d + first_index);
4912 
4913  internal::VectorDistributorLocalToGlobal<Number> distributor;
4914  read_write_operation(distributor, dst_data, mask);
4915 }
4916 
4917 
4918 
4919 template <int dim, int n_components_, typename Number, bool is_face>
4920 template <typename VectorType>
4921 inline void
4923  VectorType & dst,
4924  const unsigned int first_index,
4925  const std::bitset<VectorizedArray<Number>::n_array_elements> &mask) const
4926 {
4927  Assert(dof_values_initialized == true,
4929 
4930  // select between block vectors and non-block vectors. Note that the number
4931  // of components is checked in the internal data
4932  typename internal::BlockVectorSelector<
4933  VectorType,
4934  IsBlockVector<VectorType>::value>::BaseVectorType *dst_data[n_components];
4935  for (unsigned int d = 0; d < n_components; ++d)
4936  dst_data[d] = internal::BlockVectorSelector<
4937  VectorType,
4938  IsBlockVector<VectorType>::value>::get_vector_component(dst,
4939  d + first_index);
4940 
4941  internal::VectorSetter<Number> setter;
4942  read_write_operation(setter, dst_data, mask);
4943 }
4944 
4945 
4946 
4947 /*------------------------------ access to data fields ----------------------*/
4948 
4949 template <int dim, int n_components, typename Number, bool is_face>
4950 inline const std::vector<unsigned int> &
4953 {
4954  return data->lexicographic_numbering;
4955 }
4956 
4957 
4958 
4959 template <int dim, int n_components, typename Number, bool is_face>
4962 {
4964  const_cast<VectorizedArray<Number> *>(scratch_data),
4965  scratch_data_array->end() - scratch_data);
4966 }
4967 
4968 
4969 
4970 template <int dim, int n_components, typename Number, bool is_face>
4971 inline const VectorizedArray<Number> *
4973 {
4974  return &values_dofs[0][0];
4975 }
4976 
4977 
4978 
4979 template <int dim, int n_components, typename Number, bool is_face>
4980 inline VectorizedArray<Number> *
4982 {
4983 # ifdef DEBUG
4984  dof_values_initialized = true;
4985 # endif
4986  return &values_dofs[0][0];
4987 }
4988 
4989 
4990 
4991 template <int dim, int n_components, typename Number, bool is_face>
4992 inline const VectorizedArray<Number> *
4994 {
4995  Assert(values_quad_initialized || values_quad_submitted, ExcNotInitialized());
4996  return &values_quad[0][0];
4997 }
4998 
4999 
5000 
5001 template <int dim, int n_components, typename Number, bool is_face>
5002 inline VectorizedArray<Number> *
5004 {
5005 # ifdef DEBUG
5006  values_quad_initialized = true;
5007  values_quad_submitted = true;
5008 # endif
5009  return &values_quad[0][0];
5010 }
5011 
5012 
5013 
5014 template <int dim, int n_components, typename Number, bool is_face>
5015 inline const VectorizedArray<Number> *
5017 {
5018  Assert(gradients_quad_initialized || gradients_quad_submitted,
5019  ExcNotInitialized());
5020  return &gradients_quad[0][0][0];
5021 }
5022 
5023 
5024 
5025 template <int dim, int n_components, typename Number, bool is_face>
5026 inline VectorizedArray<Number> *
5028 {
5029 # ifdef DEBUG
5030  gradients_quad_submitted = true;
5031  gradients_quad_initialized = true;
5032 # endif
5033  return &gradients_quad[0][0][0];
5034 }
5035 
5036 
5037 
5038 template <int dim, int n_components, typename Number, bool is_face>
5039 inline const VectorizedArray<Number> *
5041 {
5042  Assert(hessians_quad_initialized, ExcNotInitialized());
5043  return &hessians_quad[0][0][0];
5044 }
5045 
5046 
5047 
5048 template <int dim, int n_components, typename Number, bool is_face>
5049 inline VectorizedArray<Number> *
5051 {
5052 # ifdef DEBUG
5053  hessians_quad_initialized = true;
5054 # endif
5055  return &hessians_quad[0][0][0];
5056 }
5057 
5058 
5059 
5060 template <int dim, int n_components_, typename Number, bool is_face>
5061 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
5063  const unsigned int dof) const
5064 {
5065  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5067  for (unsigned int comp = 0; comp < n_components; comp++)
5068  return_value[comp] = this->values_dofs[comp][dof];
5069  return return_value;
5070 }
5071 
5072 
5073 
5074 template <int dim, int n_components_, typename Number, bool is_face>
5075 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
5077  const unsigned int q_point) const
5078 {
5079  Assert(this->values_quad_initialized == true,
5081  AssertIndexRange(q_point, this->n_quadrature_points);
5083  for (unsigned int comp = 0; comp < n_components; comp++)
5084  return_value[comp] = this->values_quad[comp][q_point];
5085  return return_value;
5086 }
5087 
5088 
5089 
5090 template <int dim, int n_components_, typename Number, bool is_face>
5091 inline DEAL_II_ALWAYS_INLINE
5094  const unsigned int q_point) const
5095 {
5096  Assert(this->gradients_quad_initialized == true,
5098  AssertIndexRange(q_point, this->n_quadrature_points);
5099 
5100  Assert(jacobian != nullptr, ExcNotInitialized());
5101 
5103 
5104  // Cartesian cell
5105  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5106  {
5107  for (unsigned int comp = 0; comp < n_components; comp++)
5108  for (unsigned int d = 0; d < dim; ++d)
5109  grad_out[comp][d] =
5110  (this->gradients_quad[comp][d][q_point] * jacobian[0][d][d]);
5111  }
5112  // cell with general/affine Jacobian
5113  else
5114  {
5116  jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5117  q_point :
5118  0];
5119  for (unsigned int comp = 0; comp < n_components; comp++)
5120  for (unsigned int d = 0; d < dim; ++d)
5121  {
5122  grad_out[comp][d] =
5123  jac[d][0] * this->gradients_quad[comp][0][q_point];
5124  for (unsigned int e = 1; e < dim; ++e)
5125  grad_out[comp][d] +=
5126  jac[d][e] * this->gradients_quad[comp][e][q_point];
5127  }
5128  }
5129  return grad_out;
5130 }
5131 
5132 
5133 
5134 template <int dim, int n_components_, typename Number, bool is_face>
5135 inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArray<Number>>
5137  const unsigned int q_point) const
5138 {
5139  AssertIndexRange(q_point, this->n_quadrature_points);
5140  Assert(this->gradients_quad_initialized == true,
5142 
5143  Assert(normal_x_jacobian != nullptr, ExcNotInitialized());
5144 
5146  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5147  for (unsigned int comp = 0; comp < n_components; comp++)
5148  grad_out[comp] = this->gradients_quad[comp][dim - 1][q_point] *
5149  (this->normal_x_jacobian[0][dim - 1]);
5150  else
5151  {
5152  const unsigned int index =
5153  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5154  for (unsigned int comp = 0; comp < n_components; comp++)
5155  {
5156  grad_out[comp] = this->gradients_quad[comp][0][q_point] *
5157  this->normal_x_jacobian[index][0];
5158  for (unsigned int d = 1; d < dim; ++d)
5159  grad_out[comp] += this->gradients_quad[comp][d][q_point] *
5160  this->normal_x_jacobian[index][d];
5161  }
5162  }
5163  return grad_out;
5164 }
5165 
5166 
5167 
5168 namespace internal
5169 {
5170  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
5171  // store the lower diagonal because of symmetry
5172  template <typename Number>
5173  inline void
5174  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArray<Number>> &jac,
5175  const VectorizedArray<Number> *const hessians_quad[1],
5176  const unsigned int q_point,
5177  VectorizedArray<Number> (&tmp)[1][1])
5178  {
5179  tmp[0][0] = jac[0][0] * hessians_quad[0][q_point];
5180  }
5181 
5182  template <typename Number>
5183  inline void
5184  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArray<Number>> &jac,
5185  const VectorizedArray<Number> *const hessians_quad[3],
5186  const unsigned int q_point,
5187  VectorizedArray<Number> (&tmp)[2][2])
5188  {
5189  for (unsigned int d = 0; d < 2; ++d)
5190  {
5191  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
5192  jac[d][1] * hessians_quad[2][q_point]);
5193  tmp[1][d] = (jac[d][0] * hessians_quad[2][q_point] +
5194  jac[d][1] * hessians_quad[1][q_point]);
5195  }
5196  }
5197 
5198  template <typename Number>
5199  inline void
5200  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArray<Number>> &jac,
5201  const VectorizedArray<Number> *const hessians_quad[6],
5202  const unsigned int q_point,
5203  VectorizedArray<Number> (&tmp)[3][3])
5204  {
5205  for (unsigned int d = 0; d < 3; ++d)
5206  {
5207  tmp[0][d] = (jac[d][0] * hessians_quad[0][q_point] +
5208  jac[d][1] * hessians_quad[3][q_point] +
5209  jac[d][2] * hessians_quad[4][q_point]);
5210  tmp[1][d] = (jac[d][0] * hessians_quad[3][q_point] +
5211  jac[d][1] * hessians_quad[1][q_point] +
5212  jac[d][2] * hessians_quad[5][q_point]);
5213  tmp[2][d] = (jac[d][0] * hessians_quad[4][q_point] +
5214  jac[d][1] * hessians_quad[5][q_point] +
5215  jac[d][2] * hessians_quad[2][q_point]);
5216  }
5217  }
5218 } // namespace internal
5219 
5220 
5221 
5222 template <int dim, int n_components_, typename Number, bool is_face>
5225  const unsigned int q_point) const
5226 {
5227  Assert(!is_face, ExcNotImplemented());
5228  Assert(this->hessians_quad_initialized == true,
5230  AssertIndexRange(q_point, this->n_quadrature_points);
5231 
5232  Assert(jacobian != nullptr, ExcNotImplemented());
5234  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5235  0 :
5236  q_point];
5237 
5239 
5240  // Cartesian cell
5241  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5242  {
5243  for (unsigned int comp = 0; comp < n_components; comp++)
5244  for (unsigned int d = 0; d < dim; ++d)
5245  {
5246  hessian_out[comp][d][d] =
5247  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
5248  switch (dim)
5249  {
5250  case 1:
5251  break;
5252  case 2:
5253  hessian_out[comp][0][1] =
5254  (this->hessians_quad[comp][2][q_point] * jac[0][0] *
5255  jac[1][1]);
5256  break;
5257  case 3:
5258  hessian_out[comp][0][1] =
5259  (this->hessians_quad[comp][3][q_point] * jac[0][0] *
5260  jac[1][1]);
5261  hessian_out[comp][0][2] =
5262  (this->hessians_quad[comp][4][q_point] * jac[0][0] *
5263  jac[2][2]);
5264  hessian_out[comp][1][2] =
5265  (this->hessians_quad[comp][5][q_point] * jac[1][1] *
5266  jac[2][2]);
5267  break;
5268  default:
5269  Assert(false, ExcNotImplemented());
5270  }
5271  for (unsigned int e = d + 1; e < dim; ++e)
5272  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5273  }
5274  }
5275  // cell with general Jacobian, but constant within the cell
5276  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5277  {
5278  for (unsigned int comp = 0; comp < n_components; comp++)
5279  {
5280  // compute laplacian before the gradient because it needs to access
5281  // unscaled gradient data
5282  VectorizedArray<Number> tmp[dim][dim];
5283  internal::hessian_unit_times_jac(jac,
5284  this->hessians_quad[comp],
5285  q_point,
5286  tmp);
5287 
5288  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5289  for (unsigned int d = 0; d < dim; ++d)
5290  for (unsigned int e = d; e < dim; ++e)
5291  {
5292  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5293  for (unsigned int f = 1; f < dim; ++f)
5294  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5295  }
5296 
5297  // no J' * grad(u) part here because the Jacobian is constant
5298  // throughout the cell and hence, its derivative is zero
5299 
5300  // take symmetric part
5301  for (unsigned int d = 0; d < dim; ++d)
5302  for (unsigned int e = d + 1; e < dim; ++e)
5303  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5304  }
5305  }
5306  // cell with general Jacobian
5307  else
5308  {
5309  const Tensor<1,
5310  dim *(dim + 1) / 2,
5311  Tensor<1, dim, VectorizedArray<Number>>> &jac_grad =
5312  mapping_data->jacobian_gradients
5313  [1 - this->is_interior_face]
5314  [this->mapping_data->data_index_offsets[this->cell] + q_point];
5315  for (unsigned int comp = 0; comp < n_components; comp++)
5316  {
5317  // compute laplacian before the gradient because it needs to access
5318  // unscaled gradient data
5319  VectorizedArray<Number> tmp[dim][dim];
5320  internal::hessian_unit_times_jac(jac,
5321  this->hessians_quad[comp],
5322  q_point,
5323  tmp);
5324 
5325  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
5326  for (unsigned int d = 0; d < dim; ++d)
5327  for (unsigned int e = d; e < dim; ++e)
5328  {
5329  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
5330  for (unsigned int f = 1; f < dim; ++f)
5331  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
5332  }
5333 
5334  // add diagonal part of J' * grad(u)
5335  for (unsigned int d = 0; d < dim; ++d)
5336  for (unsigned int e = 0; e < dim; ++e)
5337  hessian_out[comp][d][d] +=
5338  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5339 
5340  // add off-diagonal part of J' * grad(u)
5341  for (unsigned int d = 0, count = dim; d < dim; ++d)
5342  for (unsigned int e = d + 1; e < dim; ++e, ++count)
5343  for (unsigned int f = 0; f < dim; ++f)
5344  hessian_out[comp][d][e] +=
5345  (jac_grad[count][f] * this->gradients_quad[comp][f][q_point]);
5346 
5347  // take symmetric part
5348  for (unsigned int d = 0; d < dim; ++d)
5349  for (unsigned int e = d + 1; e < dim; ++e)
5350  hessian_out[comp][e][d] = hessian_out[comp][d][e];
5351  }
5352  }
5354  hessian_out);
5355 }
5356 
5357 
5358 
5359 template <int dim, int n_components_, typename Number, bool is_face>
5362  const unsigned int q_point) const
5363 {
5364  Assert(!is_face, ExcNotImplemented());
5365  Assert(this->hessians_quad_initialized == true,
5367  AssertIndexRange(q_point, this->n_quadrature_points);
5368 
5369  Assert(jacobian != nullptr, ExcNotImplemented());
5371  jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
5372  0 :
5373  q_point];
5374 
5376 
5377  // Cartesian cell
5378  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5379  {
5380  for (unsigned int comp = 0; comp < n_components; comp++)
5381  for (unsigned int d = 0; d < dim; ++d)
5382  hessian_out[comp][d] =
5383  (this->hessians_quad[comp][d][q_point] * jac[d][d] * jac[d][d]);
5384  }
5385  // cell with general Jacobian, but constant within the cell
5386  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
5387  {
5388  for (unsigned int comp = 0; comp < n_components; comp++)
5389  {
5390  // compute laplacian before the gradient because it needs to access
5391  // unscaled gradient data
5392  VectorizedArray<Number> tmp[dim][dim];
5393  internal::hessian_unit_times_jac(jac,
5394  this->hessians_quad[comp],
5395  q_point,
5396  tmp);
5397 
5398  // compute only the trace part of hessian, J * tmp = J *
5399  // hess_unit(u) * J^T
5400  for (unsigned int d = 0; d < dim; ++d)
5401  {
5402  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5403  for (unsigned int f = 1; f < dim; ++f)
5404  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5405  }
5406  }
5407  }
5408  // cell with general Jacobian
5409  else
5410  {
5411  const Tensor<1,
5412  dim *(dim + 1) / 2,
5413  Tensor<1, dim, VectorizedArray<Number>>> &jac_grad =
5414  mapping_data->jacobian_gradients
5415  [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
5416  for (unsigned int comp = 0; comp < n_components; comp++)
5417  {
5418  // compute laplacian before the gradient because it needs to access
5419  // unscaled gradient data
5420  VectorizedArray<Number> tmp[dim][dim];
5421  internal::hessian_unit_times_jac(jac,
5422  this->hessians_quad[comp],
5423  q_point,
5424  tmp);
5425 
5426  // compute only the trace part of hessian, J * tmp = J *
5427  // hess_unit(u) * J^T
5428  for (unsigned int d = 0; d < dim; ++d)
5429  {
5430  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
5431  for (unsigned int f = 1; f < dim; ++f)
5432  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
5433  }
5434 
5435  for (unsigned int d = 0; d < dim; ++d)
5436  for (unsigned int e = 0; e < dim; ++e)
5437  hessian_out[comp][d] +=
5438  (jac_grad[d][e] * this->gradients_quad[comp][e][q_point]);
5439  }
5440  }
5441  return hessian_out;
5442 }
5443 
5444 
5445 
5446 template <int dim, int n_components_, typename Number, bool is_face>
5449  const unsigned int q_point) const
5450 {
5451  Assert(is_face == false, ExcNotImplemented());
5452  Assert(this->hessians_quad_initialized == true,
5454  AssertIndexRange(q_point, this->n_quadrature_points);
5455 
5458  hess_diag = get_hessian_diagonal(q_point);
5459  for (unsigned int comp = 0; comp < n_components; ++comp)
5460  {
5461  laplacian_out[comp] = hess_diag[comp][0];
5462  for (unsigned int d = 1; d < dim; ++d)
5463  laplacian_out[comp] += hess_diag[comp][d];
5464  }
5465  return laplacian_out;
5466 }
5467 
5468 
5469 
5470 template <int dim, int n_components_, typename Number, bool is_face>
5471 inline DEAL_II_ALWAYS_INLINE void
5473  const Tensor<1, n_components_, VectorizedArray<Number>> val_in,
5474  const unsigned int dof)
5475 {
5476 # ifdef DEBUG
5477  this->dof_values_initialized = true;
5478 # endif
5479  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5480  for (unsigned int comp = 0; comp < n_components; comp++)
5481  this->values_dofs[comp][dof] = val_in[comp];
5482 }
5483 
5484 
5485 
5486 template <int dim, int n_components_, typename Number, bool is_face>
5487 inline DEAL_II_ALWAYS_INLINE void
5489  const Tensor<1, n_components_, VectorizedArray<Number>> val_in,
5490  const unsigned int q_point)
5491 {
5492 # ifdef DEBUG
5494  AssertIndexRange(q_point, this->n_quadrature_points);
5495  Assert(this->J_value != nullptr, ExcNotInitialized());
5496  this->values_quad_submitted = true;
5497 # endif
5498 
5499  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5500  {
5501  const VectorizedArray<Number> JxW =
5502  J_value[0] * quadrature_weights[q_point];
5503  for (unsigned int comp = 0; comp < n_components; ++comp)
5504  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5505  }
5506  else
5507  {
5508  const VectorizedArray<Number> JxW = J_value[q_point];
5509  for (unsigned int comp = 0; comp < n_components; ++comp)
5510  this->values_quad[comp][q_point] = val_in[comp] * JxW;
5511  }
5512 }
5513 
5514 
5515 
5516 template <int dim, int n_components_, typename Number, bool is_face>
5517 inline DEAL_II_ALWAYS_INLINE void
5519  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArray<Number>>>
5520  grad_in,
5521  const unsigned int q_point)
5522 {
5523 # ifdef DEBUG
5525  AssertIndexRange(q_point, this->n_quadrature_points);
5526  this->gradients_quad_submitted = true;
5527  Assert(this->J_value != nullptr, ExcNotInitialized());
5528  Assert(this->jacobian != nullptr, ExcNotInitialized());
5529 # endif
5530 
5531  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5532  {
5533  const VectorizedArray<Number> JxW =
5534  J_value[0] * quadrature_weights[q_point];
5535  for (unsigned int comp = 0; comp < n_components; comp++)
5536  for (unsigned int d = 0; d < dim; ++d)
5537  this->gradients_quad[comp][d][q_point] =
5538  (grad_in[comp][d] * jacobian[0][d][d] * JxW);
5539  }
5540  else
5541  {
5543  this->cell_type > internal::MatrixFreeFunctions::affine ?
5544  jacobian[q_point] :
5545  jacobian[0];
5546  const VectorizedArray<Number> JxW =
5547  this->cell_type > internal::MatrixFreeFunctions::affine ?
5548  J_value[q_point] :
5549  J_value[0] * quadrature_weights[q_point];
5550  for (unsigned int comp = 0; comp < n_components; ++comp)
5551  for (unsigned int d = 0; d < dim; ++d)
5552  {
5553  VectorizedArray<Number> new_val = jac[0][d] * grad_in[comp][0];
5554  for (unsigned int e = 1; e < dim; ++e)
5555  new_val += (jac[e][d] * grad_in[comp][e]);
5556  this->gradients_quad[comp][d][q_point] = new_val * JxW;
5557  }
5558  }
5559 }
5560 
5561 
5562 
5563 template <int dim, int n_components_, typename Number, bool is_face>
5564 inline DEAL_II_ALWAYS_INLINE void
5566  const Tensor<1, n_components_, VectorizedArray<Number>> grad_in,
5567  const unsigned int q_point)
5568 {
5569 # ifdef DEBUG
5570  AssertIndexRange(q_point, this->n_quadrature_points);
5571  this->gradients_quad_submitted = true;
5572  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
5573 # endif
5574 
5575  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5576  for (unsigned int comp = 0; comp < n_components; comp++)
5577  {
5578  for (unsigned int d = 0; d < dim - 1; ++d)
5579  this->gradients_quad[comp][d][q_point] = VectorizedArray<Number>();
5580  this->gradients_quad[comp][dim - 1][q_point] =
5581  grad_in[comp] *
5582  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5583  this->quadrature_weights[q_point]);
5584  }
5585  else
5586  {
5587  const unsigned int index =
5588  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5589  for (unsigned int comp = 0; comp < n_components; comp++)
5590  {
5591  VectorizedArray<Number> factor = grad_in[comp] * this->J_value[index];
5592  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5593  factor = factor * this->quadrature_weights[q_point];
5594  for (unsigned int d = 0; d < dim; ++d)
5595  this->gradients_quad[comp][d][q_point] =
5596  factor * this->normal_x_jacobian[index][d];
5597  }
5598  }
5599 }
5600 
5601 
5602 
5603 template <int dim, int n_components_, typename Number, bool is_face>
5606 {
5607 # ifdef DEBUG
5609  Assert(this->values_quad_submitted == true,
5611 # endif
5613  for (unsigned int comp = 0; comp < n_components; ++comp)
5614  return_value[comp] = this->values_quad[comp][0];
5615  const unsigned int n_q_points = this->n_quadrature_points;
5616  for (unsigned int q = 1; q < n_q_points; ++q)
5617  for (unsigned int comp = 0; comp < n_components; ++comp)
5618  return_value[comp] += this->values_quad[comp][q];
5619  return (return_value);
5620 }
5621 
5622 
5623 
5624 /*----------------------- FEEvaluationAccess --------------------------------*/
5625 
5626 
5627 template <int dim, int n_components_, typename Number, bool is_face>
5630  const unsigned int dof_no,
5631  const unsigned int first_selected_component,
5632  const unsigned int quad_no_in,
5633  const unsigned int fe_degree,
5634  const unsigned int n_q_points,
5635  const bool is_interior_face)
5636  : FEEvaluationBase<dim, n_components_, Number, is_face>(
5637  data_in,
5638  dof_no,
5639  first_selected_component,
5640  quad_no_in,
5641  fe_degree,
5642  n_q_points,
5643  is_interior_face)
5644 {}
5645 
5646 
5647 
5648 template <int dim, int n_components_, typename Number, bool is_face>
5649 template <int n_components_other>
5652  const Mapping<dim> & mapping,
5653  const FiniteElement<dim> &fe,
5654  const Quadrature<1> & quadrature,
5655  const UpdateFlags update_flags,
5656  const unsigned int first_selected_component,
5658  : FEEvaluationBase<dim, n_components_, Number, is_face>(
5659  mapping,
5660  fe,
5661  quadrature,
5662  update_flags,
5663  first_selected_component,
5664  other)
5665 {}
5666 
5667 
5668 
5669 template <int dim, int n_components_, typename Number, bool is_face>
5673  : FEEvaluationBase<dim, n_components_, Number, is_face>(other)
5674 {}
5675 
5676 
5677 
5678 template <int dim, int n_components_, typename Number, bool is_face>
5682 {
5684  return *this;
5685 }
5686 
5687 
5688 
5689 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5690 
5691 
5692 template <int dim, typename Number, bool is_face>
5694  const MatrixFree<dim, Number> &data_in,
5695  const unsigned int dof_no,
5696  const unsigned int first_selected_component,
5697  const unsigned int quad_no_in,
5698  const unsigned int fe_degree,
5699  const unsigned int n_q_points,
5700  const bool is_interior_face)
5701  : FEEvaluationBase<dim, 1, Number, is_face>(data_in,
5702  dof_no,
5703  first_selected_component,
5704  quad_no_in,
5705  fe_degree,
5706  n_q_points,
5707  is_interior_face)
5708 {}
5709 
5710 
5711 
5712 template <int dim, typename Number, bool is_face>
5713 template <int n_components_other>
5715  const Mapping<dim> & mapping,
5716  const FiniteElement<dim> &fe,
5717  const Quadrature<1> & quadrature,
5718  const UpdateFlags update_flags,
5719  const unsigned int first_selected_component,
5721  : FEEvaluationBase<dim, 1, Number, is_face>(mapping,
5722  fe,
5723  quadrature,
5724  update_flags,
5725  first_selected_component,
5726  other)
5727 {}
5728 
5729 
5730 
5731 template <int dim, typename Number, bool is_face>
5734  : FEEvaluationBase<dim, 1, Number, is_face>(other)
5735 {}
5736 
5737 
5738 
5739 template <int dim, typename Number, bool is_face>
5743 {
5745  return *this;
5746 }
5747 
5748 
5749 
5750 template <int dim, typename Number, bool is_face>
5751 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
5753  const unsigned int dof) const
5754 {
5755  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5756  return this->values_dofs[0][dof];
5757 }
5758 
5759 
5760 
5761 template <int dim, typename Number, bool is_face>
5762 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
5764  const unsigned int q_point) const
5765 {
5766  Assert(this->values_quad_initialized == true,
5768  AssertIndexRange(q_point, this->n_quadrature_points);
5769  return this->values_quad[0][q_point];
5770 }
5771 
5772 
5773 
5774 template <int dim, typename Number, bool is_face>
5775 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
5777  const unsigned int q_point) const
5778 {
5779  return BaseClass::get_normal_derivative(q_point)[0];
5780 }
5781 
5782 
5783 
5784 template <int dim, typename Number, bool is_face>
5785 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArray<Number>>
5787  const unsigned int q_point) const
5788 {
5789  // could use the base class gradient, but that involves too many expensive
5790  // initialization operations on tensors
5791 
5792  Assert(this->gradients_quad_initialized == true,
5794  AssertIndexRange(q_point, this->n_quadrature_points);
5795 
5796  Assert(this->jacobian != nullptr, ExcNotInitialized());
5797 
5799 
5800  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5801  {
5802  for (unsigned int d = 0; d < dim; ++d)
5803  grad_out[d] =
5804  (this->gradients_quad[0][d][q_point] * this->jacobian[0][d][d]);
5805  }
5806  // cell with general/affine Jacobian
5807  else
5808  {
5810  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5811  q_point :
5812  0];
5813  for (unsigned int d = 0; d < dim; ++d)
5814  {
5815  grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point];
5816  for (unsigned int e = 1; e < dim; ++e)
5817  grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point];
5818  }
5819  }
5820  return grad_out;
5821 }
5822 
5823 
5824 
5825 template <int dim, typename Number, bool is_face>
5828  const unsigned int q_point) const
5829 {
5830  return BaseClass::get_hessian(q_point)[0];
5831 }
5832 
5833 
5834 
5835 template <int dim, typename Number, bool is_face>
5838  const unsigned int q_point) const
5839 {
5840  return BaseClass::get_hessian_diagonal(q_point)[0];
5841 }
5842 
5843 
5844 
5845 template <int dim, typename Number, bool is_face>
5848  const unsigned int q_point) const
5849 {
5850  return BaseClass::get_laplacian(q_point)[0];
5851 }
5852 
5853 
5854 
5855 template <int dim, typename Number, bool is_face>
5856 inline void DEAL_II_ALWAYS_INLINE
5858  const VectorizedArray<Number> val_in,
5859  const unsigned int dof)
5860 {
5861 # ifdef DEBUG
5862  this->dof_values_initialized = true;
5863  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5864 # endif
5865  this->values_dofs[0][dof] = val_in;
5866 }
5867 
5868 
5869 
5870 template <int dim, typename Number, bool is_face>
5871 inline void DEAL_II_ALWAYS_INLINE
5873  const VectorizedArray<Number> val_in,
5874  const unsigned int q_index)
5875 {
5876 # ifdef DEBUG
5878  AssertIndexRange(q_index, this->n_quadrature_points);
5879  Assert(this->J_value != nullptr, ExcNotInitialized());
5880  this->values_quad_submitted = true;
5881 # endif
5882  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5883  {
5885  this->J_value[0] * this->quadrature_weights[q_index];
5886  this->values_quad[0][q_index] = val_in * JxW;
5887  }
5888  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5889  {
5890  this->values_quad[0][q_index] = val_in * this->J_value[q_index];
5891  }
5892 }
5893 
5894 
5895 
5896 template <int dim, typename Number, bool is_face>
5897 inline DEAL_II_ALWAYS_INLINE void
5899  const Tensor<1, 1, VectorizedArray<Number>> val_in,
5900  const unsigned int q_point)
5901 {
5902  submit_value(val_in[0], q_point);
5903 }
5904 
5905 
5906 
5907 template <int dim, typename Number, bool is_face>
5908 inline DEAL_II_ALWAYS_INLINE void
5910  const VectorizedArray<Number> grad_in,
5911  const unsigned int q_point)
5912 {
5914  grad[0] = grad_in;
5915  BaseClass::submit_normal_derivative(grad, q_point);
5916 }
5917 
5918 
5919 
5920 template <int dim, typename Number, bool is_face>
5921 inline DEAL_II_ALWAYS_INLINE void
5923  const Tensor<1, dim, VectorizedArray<Number>> grad_in,
5924  const unsigned int q_index)
5925 {
5926 # ifdef DEBUG
5928  AssertIndexRange(q_index, this->n_quadrature_points);
5929  this->gradients_quad_submitted = true;
5930  Assert(this->J_value != nullptr, ExcNotInitialized());
5931  Assert(this->jacobian != nullptr, ExcNotInitialized());
5932 # endif
5933 
5934  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5935  {
5937  this->J_value[0] * this->quadrature_weights[q_index];
5938  for (unsigned int d = 0; d < dim; ++d)
5939  this->gradients_quad[0][d][q_index] =
5940  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5941  }
5942  // general/affine cell type
5943  else
5944  {
5946  this->cell_type > internal::MatrixFreeFunctions::affine ?
5947  this->jacobian[q_index] :
5948  this->jacobian[0];
5950  this->cell_type > internal::MatrixFreeFunctions::affine ?
5951  this->J_value[q_index] :
5952  this->J_value[0] * this->quadrature_weights[q_index];
5953  for (unsigned int d = 0; d < dim; ++d)
5954  {
5955  VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
5956  for (unsigned int e = 1; e < dim; ++e)
5957  new_val += jac[e][d] * grad_in[e];
5958  this->gradients_quad[0][d][q_index] = new_val * JxW;
5959  }
5960  }
5961 }
5962 
5963 
5964 
5965 template <int dim, typename Number, bool is_face>
5968 {
5969  return BaseClass::integrate_value()[0];
5970 }
5971 
5972 
5973 
5974 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5975 
5976 
5977 template <int dim, typename Number, bool is_face>
5979  const MatrixFree<dim, Number> &data_in,
5980  const unsigned int dof_no,
5981  const unsigned int first_selected_component,
5982  const unsigned int quad_no_in,
5983  const unsigned int fe_degree,
5984  const unsigned int n_q_points,
5985  const bool is_interior_face)
5986  : FEEvaluationBase<dim, dim, Number, is_face>(data_in,
5987  dof_no,
5988  first_selected_component,
5989  quad_no_in,
5990  fe_degree,
5991  n_q_points,
5992  is_interior_face)
5993 {}
5994 
5995 
5996 
5997 template <int dim, typename Number, bool is_face>
5998 template <int n_components_other>
6000  const Mapping<dim> & mapping,
6001  const FiniteElement<dim> &fe,
6002  const Quadrature<1> & quadrature,
6003  const UpdateFlags update_flags,
6004  const unsigned int first_selected_component,
6006  : FEEvaluationBase<dim, dim, Number, is_face>(mapping,
6007  fe,
6008  quadrature,
6009  update_flags,
6010  first_selected_component,
6011  other)
6012 {}
6013 
6014 
6015 
6016 template <int dim, typename Number, bool is_face>
6019  : FEEvaluationBase<dim, dim, Number, is_face>(other)
6020 {}
6021 
6022 
6023 
6024 template <int dim, typename Number, bool is_face>
6028 {
6030  return *this;
6031 }
6032 
6033 
6034 
6035 template <int dim, typename Number, bool is_face>
6036 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArray<Number>>
6038  const unsigned int q_point) const
6039 {
6040  return BaseClass::get_gradient(q_point);
6041 }
6042 
6043 
6044 
6045 template <int dim, typename Number, bool is_face>
6046 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
6048  const unsigned int q_point) const
6049 {
6050  Assert(this->gradients_quad_initialized == true,
6052  AssertIndexRange(q_point, this->n_quadrature_points);
6053  Assert(this->jacobian != nullptr, ExcNotInitialized());
6054 
6055  VectorizedArray<Number> divergence;
6056 
6057  // Cartesian cell
6058  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6059  {
6060  divergence =
6061  (this->gradients_quad[0][0][q_point] * this->jacobian[0][0][0]);
6062  for (unsigned int d = 1; d < dim; ++d)
6063  divergence +=
6064  (this->gradients_quad[d][d][q_point] * this->jacobian[0][d][d]);
6065  }
6066  // cell with general/constant Jacobian
6067  else
6068  {
6070  this->cell_type == internal::MatrixFreeFunctions::general ?
6071  this->jacobian[q_point] :
6072  this->jacobian[0];
6073  divergence = (jac[0][0] * this->gradients_quad[0][0][q_point]);
6074  for (unsigned int e = 1; e < dim; ++e)
6075  divergence += (jac[0][e] * this->gradients_quad[0][e][q_point]);
6076  for (unsigned int d = 1; d < dim; ++d)
6077  for (unsigned int e = 0; e < dim; ++e)
6078  divergence += (jac[d][e] * this->gradients_quad[d][e][q_point]);
6079  }
6080  return divergence;
6081 }
6082 
6083 
6084 
6085 template <int dim, typename Number, bool is_face>
6086 inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArray<Number>>
6088  const unsigned int q_point) const
6089 {
6090  // copy from generic function into dim-specialization function
6092  VectorizedArray<Number> symmetrized[(dim * dim + dim) / 2];
6093  VectorizedArray<Number> half = make_vectorized_array<Number>(0.5);
6094  for (unsigned int d = 0; d < dim; ++d)
6095  symmetrized[d] = grad[d][d];
6096  switch (dim)
6097  {
6098  case 1:
6099  break;
6100  case 2:
6101  symmetrized[2] = grad[0][1] + grad[1][0];
6102  symmetrized[2] *= half;
6103  break;
6104  case 3:
6105  symmetrized[3] = grad[0][1] + grad[1][0];
6106  symmetrized[3] *= half;
6107  symmetrized[4] = grad[0][2] + grad[2][0];
6108  symmetrized[4] *= half;
6109  symmetrized[5] = grad[1][2] + grad[2][1];
6110  symmetrized[5] *= half;
6111  break;
6112  default:
6113  Assert(false, ExcNotImplemented());
6114  }
6116 }
6117 
6118 
6119 
6120 template <int dim, typename Number, bool is_face>
6121 inline DEAL_II_ALWAYS_INLINE
6124  const unsigned int q_point) const
6125 {
6126  // copy from generic function into dim-specialization function
6129  switch (dim)
6130  {
6131  case 1:
6132  Assert(false,
6133  ExcMessage(
6134  "Computing the curl in 1d is not a useful operation"));
6135  break;
6136  case 2:
6137  curl[0] = grad[1][0] - grad[0][1];
6138  break;
6139  case 3:
6140  curl[0] = grad[2][1] - grad[1][2];
6141  curl[1] = grad[0][2] - grad[2][0];
6142  curl[2] = grad[1][0] - grad[0][1];
6143  break;
6144  default:
6145  Assert(false, ExcNotImplemented());
6146  }
6147  return curl;
6148 }
6149 
6150 
6151 
6152 template <int dim, typename Number, bool is_face>
6153 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArray<Number>>
6155  const unsigned int q_point) const
6156 {
6157  return BaseClass::get_hessian_diagonal(q_point);
6158 }
6159 
6160 
6161 
6162 template <int dim, typename Number, bool is_face>
6163 inline DEAL_II_ALWAYS_INLINE Tensor<3, dim, VectorizedArray<Number>>
6165  const unsigned int q_point) const
6166 {
6167  Assert(this->hessians_quad_initialized == true,
6169  AssertIndexRange(q_point, this->n_quadrature_points);
6170  return BaseClass::get_hessian(q_point);
6171 }
6172 
6173 
6174 
6175 template <int dim, typename Number, bool is_face>
6176 inline DEAL_II_ALWAYS_INLINE void
6178  const Tensor<2, dim, VectorizedArray<Number>> grad_in,
6179  const unsigned int q_point)
6180 {
6181  BaseClass::submit_gradient(grad_in, q_point);
6182 }
6183 
6184 
6185 
6186 template <int dim, typename Number, bool is_face>
6187 inline DEAL_II_ALWAYS_INLINE void
6189  const Tensor<1, dim, Tensor<1, dim, VectorizedArray<Number>>> grad_in,
6190  const unsigned int q_point)
6191 {
6192  BaseClass::submit_gradient(grad_in, q_point);
6193 }
6194 
6195 
6196 
6197 template <int dim, typename Number, bool is_face>
6198 inline DEAL_II_ALWAYS_INLINE void
6200  const VectorizedArray<Number> div_in,
6201  const unsigned int q_point)
6202 {
6203 # ifdef DEBUG
6205  AssertIndexRange(q_point, this->n_quadrature_points);
6206  this->gradients_quad_submitted = true;
6207  Assert(this->J_value != nullptr, ExcNotInitialized());
6208  Assert(this->jacobian != nullptr, ExcNotInitialized());
6209 # endif
6210 
6211  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6212  {
6213  const VectorizedArray<Number> fac =
6214  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6215  for (unsigned int d = 0; d < dim; ++d)
6216  {
6217  this->gradients_quad[d][d][q_point] = (fac * this->jacobian[0][d][d]);
6218  for (unsigned int e = d + 1; e < dim; ++e)
6219  {
6220  this->gradients_quad[d][e][q_point] = VectorizedArray<Number>();
6221  this->gradients_quad[e][d][q_point] = VectorizedArray<Number>();
6222  }
6223  }
6224  }
6225  else
6226  {
6228  this->cell_type == internal::MatrixFreeFunctions::general ?
6229  this->jacobian[q_point] :
6230  this->jacobian[0];
6231  const VectorizedArray<Number> fac =
6232  (this->cell_type == internal::MatrixFreeFunctions::general ?
6233  this->J_value[q_point] :
6234  this->J_value[0] * this->quadrature_weights[q_point]) *
6235  div_in;
6236  for (unsigned int d = 0; d < dim; ++d)
6237  {
6238  for (unsigned int e = 0; e < dim; ++e)
6239  this->gradients_quad[d][e][q_point] = jac[d][e] * fac;
6240  }
6241  }
6242 }
6243 
6244 
6245 
6246 template <int dim, typename Number, bool is_face>
6247 inline DEAL_II_ALWAYS_INLINE void
6249  const SymmetricTensor<2, dim, VectorizedArray<Number>> sym_grad,
6250  const unsigned int q_point)
6251 {
6252  // could have used base class operator, but that involves some overhead
6253  // which is inefficient. it is nice to have the symmetric tensor because
6254  // that saves some operations
6255 # ifdef DEBUG
6257  AssertIndexRange(q_point, this->n_quadrature_points);
6258  this->gradients_quad_submitted = true;
6259  Assert(this->J_value != nullptr, ExcNotInitialized());
6260  Assert(this->jacobian != nullptr, ExcNotInitialized());
6261 # endif
6262 
6263  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6264  {
6266  this->J_value[0] * this->quadrature_weights[q_point];
6267  for (unsigned int d = 0; d < dim; ++d)
6268  this->gradients_quad[d][d][q_point] =
6269  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6270  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6271  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6272  {
6273  const VectorizedArray<Number> value =
6274  sym_grad.access_raw_entry(counter) * JxW;
6275  this->gradients_quad[e][d][q_point] =
6276  (value * this->jacobian[0][d][d]);
6277  this->gradients_quad[d][e][q_point] =
6278  (value * this->jacobian[0][e][e]);
6279  }
6280  }
6281  // general/affine cell type
6282  else
6283  {
6285  this->cell_type == internal::MatrixFreeFunctions::general ?
6286  this->J_value[q_point] :
6287  this->J_value[0] * this->quadrature_weights[q_point];
6289  this->cell_type == internal::MatrixFreeFunctions::general ?
6290  this->jacobian[q_point] :
6291  this->jacobian[0];
6292  VectorizedArray<Number> weighted[dim][dim];
6293  for (unsigned int i = 0; i < dim; ++i)
6294  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6295  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6296  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6297  {
6298  const VectorizedArray<Number> value =
6299  sym_grad.access_raw_entry(counter) * JxW;
6300  weighted[i][j] = value;
6301  weighted[j][i] = value;
6302  }
6303  for (unsigned int comp = 0; comp < dim; ++comp)
6304  for (unsigned int d = 0; d < dim; ++d)
6305  {
6306  VectorizedArray<Number> new_val = jac[0][d] * weighted[comp][0];
6307  for (unsigned int e = 1; e < dim; ++e)
6308  new_val += jac[e][d] * weighted[comp][e];
6309  this->gradients_quad[comp][d][q_point] = new_val;
6310  }
6311  }
6312 }
6313 
6314 
6315 
6316 template <int dim, typename Number, bool is_face>
6317 inline DEAL_II_ALWAYS_INLINE void
6319  const Tensor<1, dim == 2 ? 1 : dim, VectorizedArray<Number>> curl,
6320  const unsigned int q_point)
6321 {
6323  switch (dim)
6324  {
6325  case 1:
6326  Assert(false,
6327  ExcMessage(
6328  "Testing by the curl in 1d is not a useful operation"));
6329  break;
6330  case 2:
6331  grad[1][0] = curl[0];
6332  grad[0][1] = -curl[0];
6333  break;
6334  case 3:
6335  grad[2][1] = curl[0];
6336  grad[1][2] = -curl[0];
6337  grad[0][2] = curl[1];
6338  grad[2][0] = -curl[1];
6339  grad[1][0] = curl[2];
6340  grad[0][1] = -curl[2];
6341  break;
6342  default:
6343  Assert(false, ExcNotImplemented());
6344  }
6345  submit_gradient(grad, q_point);
6346 }
6347 
6348 
6349 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6350 
6351 
6352 template <typename Number, bool is_face>
6354  const MatrixFree<1, Number> &data_in,
6355  const unsigned int dof_no,
6356  const unsigned int first_selected_component,
6357  const unsigned int quad_no_in,
6358  const unsigned int fe_degree,
6359  const unsigned int n_q_points,
6360  const bool is_interior_face)
6361  : FEEvaluationBase<1, 1, Number, is_face>(data_in,
6362  dof_no,
6363  first_selected_component,
6364  quad_no_in,
6365  fe_degree,
6366  n_q_points,
6367  is_interior_face)
6368 {}
6369 
6370 
6371 
6372 template <typename Number, bool is_face>
6373 template <int n_components_other>
6375  const Mapping<1> & mapping,
6376  const FiniteElement<1> &fe,
6377  const Quadrature<1> & quadrature,
6378  const UpdateFlags update_flags,
6379  const unsigned int first_selected_component,
6381  : FEEvaluationBase<1, 1, Number, is_face>(mapping,
6382  fe,
6383  quadrature,
6384  update_flags,
6385  first_selected_component,
6386  other)
6387 {}
6388 
6389 
6390 
6391 template <typename Number, bool is_face>
6394  : FEEvaluationBase<1, 1, Number, is_face>(other)
6395 {}
6396 
6397 
6398 
6399 template <typename Number, bool is_face>
6403 {
6405  return *this;
6406 }
6407 
6408 
6409 
6410 template <typename Number, bool is_face>
6411 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
6413  const unsigned int dof) const
6414 {
6415  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6416  return this->values_dofs[0][dof];
6417 }
6418 
6419 
6420 
6421 template <typename Number, bool is_face>
6422 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
6424  const unsigned int q_point) const
6425 {
6426  Assert(this->values_quad_initialized == true,
6428  AssertIndexRange(q_point, this->n_quadrature_points);
6429  return this->values_quad[0][q_point];
6430 }
6431 
6432 
6433 
6434 template <typename Number, bool is_face>
6435 inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArray<Number>>
6437  const unsigned int q_point) const
6438 {
6439  // could use the base class gradient, but that involves too many inefficient
6440  // initialization operations on tensors
6441 
6442  Assert(this->gradients_quad_initialized == true,
6444  AssertIndexRange(q_point, this->n_quadrature_points);
6445 
6447  this->cell_type == internal::MatrixFreeFunctions::general ?
6448  this->jacobian[q_point] :
6449  this->jacobian[0];
6450 
6452  grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];
6453 
6454  return grad_out;
6455 }
6456 
6457 
6458 
6459 template <typename Number, bool is_face>
6460 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
6462  const unsigned int q_point) const
6463 {
6464  return BaseClass::get_normal_derivative(q_point)[0];
6465 }
6466 
6467 
6468 
6469 template <typename Number, bool is_face>
6470 inline DEAL_II_ALWAYS_INLINE Tensor<2, 1, VectorizedArray<Number>>
6472  const unsigned int q_point) const
6473 {
6474  return BaseClass::get_hessian(q_point)[0];
6475 }
6476 
6477 
6478 
6479 template <typename Number, bool is_face>
6480 inline DEAL_II_ALWAYS_INLINE Tensor<1, 1, VectorizedArray<Number>>
6482  const unsigned int q_point) const
6483 {
6484  return BaseClass::get_hessian_diagonal(q_point)[0];
6485 }
6486 
6487 
6488 
6489 template <typename Number, bool is_face>
6490 inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
6492  const unsigned int q_point) const
6493 {
6494  return BaseClass::get_laplacian(q_point)[0];
6495 }
6496 
6497 
6498 
6499 template <typename Number, bool is_face>
6500 inline DEAL_II_ALWAYS_INLINE void DEAL_II_ALWAYS_INLINE
6502  const VectorizedArray<Number> val_in,
6503  const unsigned int dof)
6504 {
6505 # ifdef DEBUG
6506  this->dof_values_initialized = true;
6507  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6508 # endif
6509  this->values_dofs[0][dof] = val_in;
6510 }
6511 
6512 
6513 
6514 template <typename Number, bool is_face>
6515 inline DEAL_II_ALWAYS_INLINE void
6517  const VectorizedArray<Number> val_in,
6518  const unsigned int q_point)
6519 {
6520 # ifdef DEBUG
6522  AssertIndexRange(q_point, this->n_quadrature_points);
6523  this->values_quad_submitted = true;
6524 # endif
6525  if (this->cell_type == internal::MatrixFreeFunctions::general)
6526  {
6527  const VectorizedArray<Number> JxW = this->J_value[q_point];
6528  this->values_quad[0][q_point] = val_in * JxW;
6529  }
6530  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6531  {
6533  this->J_value[0] * this->quadrature_weights[q_point];
6534  this->values_quad[0][q_point] = val_in * JxW;
6535  }
6536 }
6537 
6538 
6539 
6540 template <typename Number, bool is_face>
6541 inline DEAL_II_ALWAYS_INLINE void
6543  const Tensor<1, 1, VectorizedArray<Number>> val_in,
6544  const unsigned int q_point)
6545 {
6546  submit_value(val_in[0], q_point);
6547 }
6548 
6549 
6550 
6551 template <typename Number, bool is_face>
6552 inline DEAL_II_ALWAYS_INLINE void
6554  const Tensor<1, 1, VectorizedArray<Number>> grad_in,
6555  const unsigned int q_point)
6556 {
6557  submit_gradient(grad_in[0], q_point);
6558 }
6559 
6560 
6561 
6562 template <typename Number, bool is_face>
6563 inline DEAL_II_ALWAYS_INLINE void
6565  const VectorizedArray<Number> grad_in,
6566  const unsigned int q_point)
6567 {
6568 # ifdef DEBUG
6570  AssertIndexRange(q_point, this->n_quadrature_points);
6571  this->gradients_quad_submitted = true;
6572 # endif
6573 
6575  this->cell_type == internal::MatrixFreeFunctions::general ?
6576  this->jacobian[q_point] :
6577  this->jacobian[0];
6579  this->cell_type == internal::MatrixFreeFunctions::general ?
6580  this->J_value[q_point] :
6581  this->J_value[0] * this->quadrature_weights[q_point];
6582 
6583  this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in * JxW;
6584 }
6585 
6586 
6587 
6588 template <typename Number, bool is_face>
6589 inline DEAL_II_ALWAYS_INLINE void
6591  const VectorizedArray<Number> grad_in,
6592  const unsigned int q_point)
6593 {
6595  grad[0] = grad_in;
6596  BaseClass::submit_normal_derivative(grad, q_point);
6597 }
6598 
6599 
6600 
6601 template <typename Number, bool is_face>
6602 inline DEAL_II_ALWAYS_INLINE void
6604  const Tensor<1, 1, VectorizedArray<Number>> grad_in,
6605  const unsigned int q_point)
6606 {
6607  BaseClass::submit_normal_derivative(grad_in, q_point);
6608 }
6609 
6610 
6611 
6612 template <typename Number, bool is_face>
6615 {
6616  return BaseClass::integrate_value()[0];
6617 }
6618 
6619 
6620 
6621 /*-------------------------- FEEvaluation -----------------------------------*/
6622 
6623 
6624 template <int dim,
6625  int fe_degree,
6626  int n_q_points_1d,
6627  int n_components_,
6628  typename Number>
6630  FEEvaluation(const MatrixFree<dim, Number> &data_in,
6631  const unsigned int fe_no,
6632  const unsigned int quad_no,
6633  const unsigned int first_selected_component)
6634  : BaseClass(data_in,
6635  fe_no,
6636  first_selected_component,
6637  quad_no,
6638  fe_degree,
6639  static_n_q_points)
6640  , dofs_per_component(this->data->dofs_per_component_on_cell)
6641  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6642  , n_q_points(this->data->n_q_points)
6643 {
6644  check_template_arguments(fe_no, 0);
6645 }
6646 
6647 
6648 
6649 template <int dim,
6650  int fe_degree,
6651  int n_q_points_1d,
6652  int n_components_,
6653  typename Number>
6655  FEEvaluation(const Mapping<dim> & mapping,
6656  const FiniteElement<dim> &fe,
6657  const Quadrature<1> & quadrature,
6658  const UpdateFlags update_flags,
6659  const unsigned int first_selected_component)
6660  : BaseClass(mapping,
6661  fe,
6662  quadrature,
6663  update_flags,
6664  first_selected_component,
6665  static_cast<FEEvaluationBase<dim, 1, Number, false> *>(nullptr))
6666  , dofs_per_component(this->data->dofs_per_component_on_cell)
6667  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6668  , n_q_points(this->data->n_q_points)
6669 {
6670  check_template_arguments(numbers::invalid_unsigned_int, 0);
6671 }
6672 
6673 
6674 
6675 template <int dim,
6676  int fe_degree,
6677  int n_q_points_1d,
6678  int n_components_,
6679  typename Number>
6682  const Quadrature<1> & quadrature,
6683  const UpdateFlags update_flags,
6684  const unsigned int first_selected_component)
6685  : BaseClass(StaticMappingQ1<dim>::mapping,
6686  fe,
6687  quadrature,
6688  update_flags,
6689  first_selected_component,
6690  static_cast<FEEvaluationBase<dim, 1, Number, false> *>(nullptr))
6691  , dofs_per_component(this->data->dofs_per_component_on_cell)
6692  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6693  , n_q_points(this->data->n_q_points)
6694 {
6695  check_template_arguments(numbers::invalid_unsigned_int, 0);
6696 }
6697 
6698 
6699 
6700 template <int dim,
6701  int fe_degree,
6702  int n_q_points_1d,
6703  int n_components_,
6704  typename Number>
6705 template <int n_components_other>
6709  const unsigned int first_selected_component)
6710  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6711  fe,
6712  other.mapped_geometry->get_quadrature(),
6713  other.mapped_geometry->get_fe_values().get_update_flags(),
6714  first_selected_component,
6715  &other)
6716  , dofs_per_component(this->data->dofs_per_component_on_cell)
6717  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6718  , n_q_points(this->data->n_q_points)
6719 {
6720  check_template_arguments(numbers::invalid_unsigned_int, 0);
6721 }
6722 
6723 
6724 
6725 template <int dim,
6726  int fe_degree,
6727  int n_q_points_1d,
6728  int n_components_,
6729  typename Number>
6731  FEEvaluation(const FEEvaluation &other)
6732  : BaseClass(other)
6733  , dofs_per_component(this->data->dofs_per_component_on_cell)
6734  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6735  , n_q_points(this->data->n_q_points)
6736 {
6737  check_template_arguments(numbers::invalid_unsigned_int, 0);
6738 }
6739 
6740 
6741 
6742 template <int dim,
6743  int fe_degree,
6744  int n_q_points_1d,
6745  int n_components_,
6746  typename Number>
6749 operator=(const FEEvaluation &other)
6750 {
6751  BaseClass::operator=(other);
6752  check_template_arguments(numbers::invalid_unsigned_int, 0);
6753  return *this;
6754 }
6755 
6756 
6757 
6758 template <int dim,
6759  int fe_degree,
6760  int n_q_points_1d,
6761  int n_components_,
6762  typename Number>
6763 inline void
6765  check_template_arguments(const unsigned int dof_no,
6766  const unsigned int first_selected_component)
6767 {
6768  (void)dof_no;
6769  (void)first_selected_component;
6770 
6771 # ifdef DEBUG
6772  // print error message when the dimensions do not match. Propose a possible
6773  // fix
6774  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6775  static_cast<unsigned int>(fe_degree) != this->data->fe_degree) ||
6776  n_q_points != this->n_quadrature_points)
6777  {
6778  std::string message =
6779  "-------------------------------------------------------\n";
6780  message += "Illegal arguments in constructor/wrong template arguments!\n";
6781  message += " Called --> FEEvaluation<dim,";
6782  message += Utilities::int_to_string(fe_degree) + ",";
6783  message += Utilities::int_to_string(n_q_points_1d);
6784  message += "," + Utilities::int_to_string(n_components);
6785  message += ",Number>(data";
6786  if (first_selected_component != numbers::invalid_unsigned_int)
6787  {
6788  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6789  message += Utilities::int_to_string(this->quad_no) + ", ";
6790  message += Utilities::int_to_string(first_selected_component);
6791  }
6792  message += ")\n";
6793 
6794  // check whether some other vector component has the correct number of
6795  // points
6796  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6797  proposed_fe_comp = numbers::invalid_unsigned_int,
6798  proposed_quad_comp = numbers::invalid_unsigned_int;
6799  if (dof_no != numbers::invalid_unsigned_int)
6800  {
6801  if (static_cast<unsigned int>(fe_degree) == this->data->fe_degree)
6802  {
6803  proposed_dof_comp = dof_no;
6804  proposed_fe_comp = first_selected_component;
6805  }
6806  else
6807  for (unsigned int no = 0; no < this->matrix_info->n_components();
6808  ++no)
6809  for (unsigned int nf = 0;
6810  nf < this->matrix_info->n_base_elements(no);
6811  ++nf)
6812  if (this->matrix_info
6813  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6814  .fe_degree == static_cast<unsigned int>(fe_degree))
6815  {
6816  proposed_dof_comp = no;
6817  proposed_fe_comp = nf;
6818  break;
6819  }
6820  if (n_q_points ==
6821  this->mapping_data->descriptor[this->active_quad_index]
6822  .n_q_points)
6823  proposed_quad_comp = this->quad_no;
6824  else
6825  for (unsigned int no = 0;
6826  no < this->matrix_info->get_mapping_info().cell_data.size();
6827  ++no)
6828  if (this->matrix_info->get_mapping_info()
6829  .cell_data[no]
6830  .descriptor[this->active_quad_index]
6831  .n_q_points == n_q_points)
6832  {
6833  proposed_quad_comp = no;
6834  break;
6835  }
6836  }
6837  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6838  proposed_quad_comp != numbers::invalid_unsigned_int)
6839  {
6840  if (proposed_dof_comp != first_selected_component)
6841  message += "Wrong vector component selection:\n";
6842  else
6843  message += "Wrong quadrature formula selection:\n";
6844  message += " Did you mean FEEvaluation<dim,";
6845  message += Utilities::int_to_string(fe_degree) + ",";
6846  message += Utilities::int_to_string(n_q_points_1d);
6847  message += "," + Utilities::int_to_string(n_components);
6848  message += ",Number>(data";
6849  if (dof_no != numbers::invalid_unsigned_int)
6850  {
6851  message +=
6852  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6853  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
6854  message += Utilities::int_to_string(proposed_fe_comp);
6855  }
6856  message += ")?\n";
6857  std::string correct_pos;
6858  if (proposed_dof_comp != dof_no)
6859  correct_pos = " ^ ";
6860  else
6861  correct_pos = " ";
6862  if (proposed_quad_comp != this->quad_no)
6863  correct_pos += " ^ ";
6864  else
6865  correct_pos += " ";
6866  if (proposed_fe_comp != first_selected_component)
6867  correct_pos += " ^\n";
6868  else
6869  correct_pos += " \n";
6870  message += " " +
6871  correct_pos;
6872  }
6873  // ok, did not find the numbers specified by the template arguments in
6874  // the given list. Suggest correct template arguments
6875  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6876  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6877  message += "Wrong template arguments:\n";
6878  message += " Did you mean FEEvaluation<dim,";
6879  message += Utilities::int_to_string(this->data->fe_degree) + ",";
6880  message += Utilities::int_to_string(proposed_n_q_points_1d);
6881  message += "," + Utilities::int_to_string(n_components);
6882  message += ",Number>(data";
6883  if (dof_no != numbers::invalid_unsigned_int)
6884  {
6885  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6886  message += Utilities::int_to_string(this->quad_no);
6887  message += ", " + Utilities::int_to_string(first_selected_component);
6888  }
6889  message += ")?\n";
6890  std::string correct_pos;
6891  if (this->data->fe_degree != static_cast<unsigned int>(fe_degree))
6892  correct_pos = " ^";
6893  else
6894  correct_pos = " ";
6895  if (proposed_n_q_points_1d != n_q_points_1d)
6896  correct_pos += " ^\n";
6897  else
6898  correct_pos += " \n";
6899  message += " " + correct_pos;
6900 
6901  Assert(static_cast<unsigned int>(fe_degree) == this->data->fe_degree &&
6902  n_q_points == this->n_quadrature_points,
6903  ExcMessage(message));
6904  }
6905  if (dof_no != numbers::invalid_unsigned_int)
6907  n_q_points,
6908  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6909 # endif
6910 }
6911 
6912 
6913 
6914 template <int dim,
6915  int fe_degree,
6916  int n_q_points_1d,
6917  int n_components_,
6918  typename Number>
6919 inline void
6921  const unsigned int cell_index)
6922 {
6923  Assert(this->mapped_geometry == nullptr,
6924  ExcMessage("FEEvaluation was initialized without a matrix-free object."
6925  " Integer indexing is not possible"));
6926  if (this->mapped_geometry != nullptr)
6927  return;
6928 
6929  Assert(this->dof_info != nullptr, ExcNotInitialized());
6930  Assert(this->mapping_data != nullptr, ExcNotInitialized());
6931  this->cell = cell_index;
6932  this->cell_type =
6933  this->matrix_info->get_mapping_info().get_cell_type(cell_index);
6934 
6935  const unsigned int offsets =
6936  this->mapping_data->data_index_offsets[cell_index];
6937  this->jacobian = &this->mapping_data->jacobians[0][offsets];
6938  this->J_value = &this->mapping_data->JxW_values[offsets];
6939 
6940 # ifdef DEBUG
6941  this->dof_values_initialized = false;
6942  this->values_quad_initialized = false;
6943  this->gradients_quad_initialized = false;
6944  this->hessians_quad_initialized = false;
6945 # endif
6946 }
6947 
6948 
6949 
6950 template <int dim,
6951  int fe_degree,
6952  int n_q_points_1d,
6953  int n_components_,
6954  typename Number>
6955 template <typename DoFHandlerType, bool level_dof_access>
6956 inline void
6959 {
6960  Assert(this->matrix_info == nullptr,
6961  ExcMessage("Cannot use initialization from cell iterator if "
6962  "initialized from MatrixFree object. Use variant for "
6963  "on the fly computation with arguments as for FEValues "
6964  "instead"));
6965  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6966  this->mapped_geometry->reinit(
6967  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6968  this->local_dof_indices.resize(cell->get_fe().dofs_per_cell);
6969  if (level_dof_access)
6970  cell->get_mg_dof_indices(this->local_dof_indices);
6971  else
6972  cell->get_dof_indices(this->local_dof_indices);
6973 }
6974 
6975 
6976 
6977 template <int dim,
6978  int fe_degree,
6979  int n_q_points_1d,
6980  int n_components_,
6981  typename Number>
6982 inline void
6984  const typename Triangulation<dim>::cell_iterator &cell)
6985 {
6986  Assert(this->matrix_info == 0,
6987  ExcMessage("Cannot use initialization from cell iterator if "
6988  "initialized from MatrixFree object. Use variant for "
6989  "on the fly computation with arguments as for FEValues "
6990  "instead"));
6991  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6992  this->mapped_geometry->reinit(cell);
6993 }
6994 
6995 
6996 
6997 template <int dim,
6998  int fe_degree,
6999  int n_q_points_1d,
7000  int n_components_,
7001  typename Number>
7004  quadrature_point(const unsigned int q) const
7005 {
7006  if (this->matrix_info == nullptr)
7007  {
7008  Assert((this->mapped_geometry->get_fe_values().get_update_flags() |
7010  ExcNotInitialized());
7011  }
7012  else
7013  {
7014  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
7015  ExcNotInitialized());
7016  }
7017 
7018  AssertIndexRange(q, n_q_points);
7019 
7020  const unsigned int n_q_points_1d_actual =
7021  fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;
7022 
7023  // Cartesian mesh: not all quadrature points are stored, only the
7024  // diagonal. Hence, need to find the tensor product index and retrieve the
7025  // value from that
7026  const Point<dim, VectorizedArray<Number>> *quadrature_points =
7027  &this->mapping_data->quadrature_points
7028  [this->mapping_data->quadrature_point_offsets[this->cell]];
7029  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
7030  {
7032  switch (dim)
7033  {
7034  case 1:
7035  return quadrature_points[q];
7036  case 2:
7037  point[0] = quadrature_points[q % n_q_points_1d_actual][0];
7038  point[1] = quadrature_points[q / n_q_points_1d_actual][1];
7039  return point;
7040  case 3:
7041  point[0] = quadrature_points[q % n_q_points_1d_actual][0];
7042  point[1] = quadrature_points[(q / n_q_points_1d_actual) %
7043  n_q_points_1d_actual][1];
7044  point[2] = quadrature_points[q / (n_q_points_1d_actual *
7045  n_q_points_1d_actual)][2];
7046  return point;
7047  default:
7048  Assert(false, ExcNotImplemented());
7049  return point;
7050  }
7051  }
7052  // all other cases: just return the respective data as it is fully stored
7053  else
7054  return quadrature_points[q];
7055 }
7056 
7057 
7058 
7059 template <int dim,
7060  int fe_degree,
7061  int n_q_points_1d,
7062  int n_components_,
7063  typename Number>
7064 inline void
7066  const bool evaluate_values,
7067  const bool evaluate_gradients,
7068  const bool evaluate_hessians)
7069 {
7070  Assert(this->dof_values_initialized == true,
7072  evaluate(this->values_dofs[0],
7073  evaluate_values,
7074  evaluate_gradients,
7075  evaluate_hessians);
7076 }
7077 
7078 
7079 
7080 template <int dim,
7081  int fe_degree,
7082  int n_q_points_1d,
7083  int n_components_,
7084  typename Number>
7085 inline void
7087  const VectorizedArray<Number> *values_array,
7088  const bool evaluate_values,
7089  const bool evaluate_gradients,
7090  const bool evaluate_hessians)
7091 {
7093  dim,
7094  fe_degree,
7095  n_q_points_1d,
7096  n_components,
7097  VectorizedArray<Number>>::evaluate(*this->data,
7098  const_cast<VectorizedArray<Number> *>(
7099  values_array),
7100  this->values_quad[0],
7101  this->gradients_quad[0][0],
7102  this->hessians_quad[0][0],
7103  this->scratch_data,
7104  evaluate_values,
7105  evaluate_gradients,
7106  evaluate_hessians);
7107 
7108 # ifdef DEBUG
7109  if (evaluate_values == true)
7110  this->values_quad_initialized = true;
7111  if (evaluate_gradients == true)
7112  this->gradients_quad_initialized = true;
7113  if (evaluate_hessians == true)
7114  this->hessians_quad_initialized = true;
7115 # endif
7116 }
7117 
7118 
7119 
7120 template <int dim,
7121  int fe_degree,
7122  int n_q_points_1d,
7123  int n_components_,
7124  typename Number>
7125 template <typename VectorType>
7126 inline void
7128  gather_evaluate(const VectorType &input_vector,
7129  const bool evaluate_values,
7130  const bool evaluate_gradients,
7131  const bool evaluate_hessians)
7132 {
7133  // If the index storage is interleaved and contiguous and the vector storage
7134  // has the correct alignment, we can directly pass the pointer into the
7135  // vector to the evaluate() call, without reading the vector entries into a
7136  // separate data field. This saves some operations.
7137  if (std::is_same<typename VectorType::value_type, Number>::value &&
7138  this->dof_info->index_storage_variants
7140  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7141  IndexStorageVariants::interleaved_contiguous &&
7142  reinterpret_cast<std::size_t>(
7143  input_vector.begin() +
7144  this->dof_info->dof_indices_contiguous
7147  sizeof(VectorizedArray<Number>) ==
7148  0)
7149  {
7150  const VectorizedArray<Number> *vec_values =
7151  reinterpret_cast<const VectorizedArray<Number> *>(
7152  input_vector.begin() +
7153  this->dof_info->dof_indices_contiguous
7156  this->dof_info
7157  ->component_dof_indices_offset[this->active_fe_index]
7158  [this->first_selected_component] *
7160 
7161  evaluate(vec_values,
7162  evaluate_values,
7163  evaluate_gradients,
7164  evaluate_hessians);
7165  }
7166  else
7167  {
7168  this->read_dof_values(input_vector);
7169  evaluate(this->begin_dof_values(),
7170  evaluate_values,
7171  evaluate_gradients,
7172  evaluate_hessians);
7173  }
7174 }
7175 
7176 
7177 
7178 template <int dim,
7179  int fe_degree,
7180  int n_q_points_1d,
7181  int n_components_,
7182  typename Number>
7183 inline void
7185  const bool integrate_values,
7186  const bool integrate_gradients)
7187 {
7188  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7189 
7190 # ifdef DEBUG
7191  this->dof_values_initialized = true;
7192 # endif
7193 }
7194 
7195 
7196 
7197 template <int dim,
7198  int fe_degree,
7199  int n_q_points_1d,
7200  int n_components_,
7201  typename Number>
7202 inline void
7204  const bool integrate_values,
7205  const bool integrate_gradients,
7206  VectorizedArray<Number> *values_array)
7207 {
7208  if (integrate_values == true)
7209  Assert(this->values_quad_submitted == true,
7211  if (integrate_gradients == true)
7212  Assert(this->gradients_quad_submitted == true,
7214  Assert(this->matrix_info != nullptr ||
7215  this->mapped_geometry->is_initialized(),
7216  ExcNotInitialized());
7217 
7218  SelectEvaluator<dim,
7219  fe_degree,
7220  n_q_points_1d,
7221  n_components,
7222  VectorizedArray<Number>>::integrate(*this->data,
7223  values_array,
7224  this->values_quad[0],
7225  this
7226  ->gradients_quad[0][0],
7227  this->scratch_data,
7228  integrate_values,
7229  integrate_gradients,
7230  false);
7231 
7232 # ifdef DEBUG
7233  this->dof_values_initialized = true;
7234 # endif
7235 }
7236 
7237 
7238 
7239 template <int dim,
7240  int fe_degree,
7241  int n_q_points_1d,
7242  int n_components_,
7243  typename Number>
7244 template <typename VectorType>
7245 inline void
7247  integrate_scatter(const bool integrate_values,
7248  const bool integrate_gradients,
7249  VectorType &destination)
7250 {
7251  // If the index storage is interleaved and contiguous and the vector storage
7252  // has the correct alignment, we can directly pass the pointer into the
7253  // vector to the integrate() call, without writing temporary results into a
7254  // separate data field that will later be added into the vector. This saves
7255  // some operations.
7256  if (std::is_same<typename VectorType::value_type, Number>::value &&
7257  this->dof_info->index_storage_variants
7259  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
7260  IndexStorageVariants::interleaved_contiguous &&
7261  reinterpret_cast<std::size_t>(
7262  destination.begin() +
7263  this->dof_info->dof_indices_contiguous
7266  sizeof(VectorizedArray<Number>) ==
7267  0)
7268  {
7269  VectorizedArray<Number> *vec_values =
7270  reinterpret_cast<VectorizedArray<Number> *>(
7271  destination.begin() +
7272  this->dof_info->dof_indices_contiguous
7275  this->dof_info
7276  ->component_dof_indices_offset[this->active_fe_index]
7277  [this->first_selected_component] *
7280  dim,
7281  fe_degree,
7282  n_q_points_1d,
7283  n_components,
7284  VectorizedArray<Number>>::integrate(*this->data,
7285  vec_values,
7286  this->values_quad[0],
7287  this->gradients_quad[0][0],
7288  this->scratch_data,
7289  integrate_values,
7290  integrate_gradients,
7291  true);
7292  }
7293  else
7294  {
7295  integrate(integrate_values,
7296  integrate_gradients,
7297  this->begin_dof_values());
7298  this->distribute_local_to_global(destination);
7299  }
7300 }
7301 
7302 
7303 
7304 /*-------------------------- FEFaceEvaluation ---------------------------*/
7305 
7306 
7307 
7308 template <int dim,
7309  int fe_degree,
7310  int n_q_points_1d,
7311  int n_components_,
7312  typename Number>
7314  FEFaceEvaluation(const MatrixFree<dim, Number> &matrix_free,
7315  const bool is_interior_face,
7316  const unsigned int dof_no,
7317  const unsigned int quad_no,
7318  const unsigned int first_selected_component)
7319  : BaseClass(matrix_free,
7320  dof_no,
7321  first_selected_component,
7322  quad_no,
7323  fe_degree,
7324  static_n_q_points,
7325  is_interior_face)
7326  , dofs_per_component(this->data->dofs_per_component_on_cell)
7327  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7328  , n_q_points(this->data->n_q_points_face)
7329 {}
7330 
7331 
7332 
7333 template <int dim,
7334  int fe_degree,
7335  int n_q_points_1d,
7336  int n_components_,
7337  typename Number>
7338 inline void
7340  const unsigned int face_index)
7341 {
7342  Assert(this->mapped_geometry == nullptr,
7343  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7344  " Integer indexing is not possible"));
7345  if (this->mapped_geometry != nullptr)
7346  return;
7347 
7348  this->cell = face_index;
7349  this->dof_access_index =
7350  this->is_interior_face ?
7353  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7354  const unsigned int n_vectors = VectorizedArray<Number>::n_array_elements;
7356  this->matrix_info->get_face_info(face_index);
7357  if (face_index >=
7358  this->matrix_info->get_task_info().face_partition_data.back() &&
7359  face_index <
7360  this->matrix_info->get_task_info().boundary_partition_data.back())
7361  Assert(this->is_interior_face,
7362  ExcMessage("Boundary faces do not have a neighbor"));
7363 
7364  this->face_no =
7365  (this->is_interior_face ? faces.interior_face_no : faces.exterior_face_no);
7366  this->subface_index = faces.subface_index;
7367  if (this->is_interior_face == true)
7368  {
7369  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7370  if (faces.face_orientation > 8)
7371  this->face_orientation = faces.face_orientation - 8;
7372  else
7373  this->face_orientation = 0;
7374  }
7375  else
7376  {
7377  if (faces.face_orientation < 8)
7378  this->face_orientation = faces.face_orientation;
7379  else
7380  this->face_orientation = 0;
7381  }
7382 
7383  this->values_quad_submitted = false;
7384 
7385  this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
7386  const unsigned int offsets =
7387  this->mapping_data->data_index_offsets[face_index];
7388  this->J_value = &this->mapping_data->JxW_values[offsets];
7389  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
7390  this->jacobian =
7391  &this->mapping_data->jacobians[!this->is_interior_face][offsets];
7392  this->normal_x_jacobian =
7393  &this->mapping_data
7394  ->normals_times_jacobians[!this->is_interior_face][offsets];
7395 
7396 # ifdef DEBUG
7397  this->dof_values_initialized = false;
7398  this->values_quad_initialized = false;
7399  this->gradients_quad_initialized = false;
7400  this->hessians_quad_initialized = false;
7401 # endif
7402 }
7403 
7404 
7405 
7406 template <int dim,
7407  int fe_degree,
7408  int n_q_points_1d,
7409  int n_components_,
7410  typename Number>
7411 inline void
7413  const unsigned int cell_index,
7414  const unsigned int face_number)
7415 {
7416  Assert(
7417  this->quad_no <
7418  this->matrix_info->get_mapping_info().face_data_by_cells.size(),
7419  ExcMessage(
7420  "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
7422  AssertIndexRange(cell_index,
7423  this->matrix_info->get_mapping_info().cell_type.size());
7424  Assert(this->mapped_geometry == nullptr,
7425  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7426  " Integer indexing is not possible"));
7427  Assert(this->is_interior_face == true,
7428  ExcMessage(
7429  "Cell-based FEFaceEvaluation::reinit only possible for the "
7430  "interior face with second argument to constructor as true"));
7431  if (this->mapped_geometry != nullptr)
7432  return;
7433  Assert(this->matrix_info != nullptr, ExcNotInitialized());
7434 
7435  this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
7436  this->cell = cell_index;
7437  this->face_orientation = 0;
7438  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
7439  this->face_no = face_number;
7440  this->dof_access_index =
7442 
7443  const unsigned int offsets =
7444  this->matrix_info->get_mapping_info()
7445  .face_data_by_cells[this->quad_no]
7446  .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
7447  face_number];
7448  AssertIndexRange(offsets,
7449  this->matrix_info->get_mapping_info()
7450  .face_data_by_cells[this->quad_no]
7451  .JxW_values.size());
7452  this->J_value = &this->matrix_info->get_mapping_info()
7453  .face_data_by_cells[this->quad_no]
7454  .JxW_values[offsets];
7455  this->normal_vectors = &this->matrix_info->get_mapping_info()
7456  .face_data_by_cells[this->quad_no]
7457  .normal_vectors[offsets];
7458  this->jacobian = &this->matrix_info->get_mapping_info()
7459  .face_data_by_cells[this->quad_no]
7460  .jacobians[0][offsets];
7461  this->normal_x_jacobian = &this->matrix_info->get_mapping_info()
7462  .face_data_by_cells[this->quad_no]
7463  .normals_times_jacobians[0][offsets];
7464 
7465 # ifdef DEBUG
7466  this->dof_values_initialized = false;
7467  this->values_quad_initialized = false;
7468  this->gradients_quad_initialized = false;
7469  this->hessians_quad_initialized = false;
7470 # endif
7471 }
7472 
7473 
7474 
7475 template <int dim,
7476  int fe_degree,
7477  int n_q_points_1d,
7478  int n_components,
7479  typename Number>
7480 inline void
7482  const bool evaluate_values,
7483  const bool evaluate_gradients)
7484 {
7485  Assert(this->dof_values_initialized, ExcNotInitialized());
7486 
7487  evaluate(this->values_dofs[0], evaluate_values, evaluate_gradients);
7488 }
7489 
7490 
7491 
7492 template <int dim,
7493  int fe_degree,
7494  int n_q_points_1d,
7495  int n_components,
7496  typename Number>
7497 inline void
7499  const VectorizedArray<Number> *values_array,
7500  const bool evaluate_values,
7501  const bool evaluate_gradients)
7502 {
7503  if (!(evaluate_values + evaluate_gradients))
7504  return;
7505 
7506  constexpr unsigned int static_dofs_per_face =
7507  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
7508  numbers::invalid_unsigned_int;
7509  const unsigned int dofs_per_face =
7510  fe_degree > -1 ? static_dofs_per_face :
7511  Utilities::pow(this->data->fe_degree + 1, dim - 1);
7512 
7513  // we allocate small amounts of data on the stack to signal the compiler
7514  // that this temporary data is only needed for the calculations but the
7515  // final results can be discarded and need not be written back to
7516  // memory. For large sizes or when the dofs per face is not a compile-time
7517  // constant, however, we want to go to the heap in the `scratch_data`
7518  // variable to not risk a stack overflow.
7519  constexpr unsigned int stack_array_size_threshold = 100;
7520 
7522  temp_data[static_dofs_per_face < stack_array_size_threshold ?
7523  n_components * 2 * static_dofs_per_face :
7524  1];
7525  VectorizedArray<Number> *temp1;
7526  if (static_dofs_per_face < stack_array_size_threshold)
7527  temp1 = &temp_data[0];
7528  else
7529  temp1 = this->scratch_data;
7530 
7531  internal::FEFaceNormalEvaluationImpl<dim,
7532  fe_degree,
7533  n_components,
7535  template interpolate<true, false>(
7536  *this->data, values_array, temp1, evaluate_gradients, this->face_no);
7537 
7538  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
7539  if (fe_degree > -1 &&
7540  this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
7541  this->data->element_type <=
7543  internal::FEFaceEvaluationImpl<
7544  true,
7545  dim,
7546  fe_degree,
7547  n_q_points_1d_actual,
7548  n_components,
7549  VectorizedArray<Number>>::evaluate_in_face(*this->data,
7550  temp1,
7551  this->begin_values(),
7552  this->begin_gradients(),
7553  this->scratch_data +
7554  2 * n_components *
7555  dofs_per_face,
7556  evaluate_values,
7557  evaluate_gradients,
7558  this->subface_index);
7559  else
7560  internal::FEFaceEvaluationImpl<
7561  false,
7562  dim,
7563  fe_degree,
7564  n_q_points_1d_actual,
7565  n_components,
7566  VectorizedArray<Number>>::evaluate_in_face(*this->data,
7567  temp1,
7568  this->begin_values(),
7569  this->begin_gradients(),
7570  this->scratch_data +
7571  2 * n_components *
7572  dofs_per_face,
7573  evaluate_values,
7574  evaluate_gradients,
7575  this->subface_index);
7576 
7577  if (this->face_orientation)
7578  adjust_for_face_orientation(false, evaluate_values, evaluate_gradients);
7579 
7580 # ifdef DEBUG
7581  if (evaluate_values == true)
7582  this->values_quad_initialized = true;
7583  if (evaluate_gradients == true)
7584  this->gradients_quad_initialized = true;
7585 # endif
7586 }
7587 
7588 
7589 
7590 template <int dim,
7591  int fe_degree,
7592  int n_q_points_1d,
7593  int n_components,
7594  typename Number>
7595 inline void
7597  integrate(const bool integrate_values, const bool integrate_gradients)
7598 {
7599  integrate(integrate_values, integrate_gradients, this->values_dofs[0]);
7600 
7601 # ifdef DEBUG
7602  this->dof_values_initialized = true;
7603 # endif
7604 }
7605 
7606 
7607 
7608 template <int dim,
7609  int fe_degree,
7610  int n_q_points_1d,
7611  int n_components,
7612  typename Number>
7613 inline void
7615  integrate(const bool integrate_values,
7616  const bool integrate_gradients,
7617  VectorizedArray<Number> *values_array)
7618 {
7619  if (!(integrate_values + integrate_gradients))
7620  return;
7621 
7622  if (this->face_orientation)
7623  adjust_for_face_orientation(true, integrate_values, integrate_gradients);
7624 
7625  constexpr unsigned int static_dofs_per_face =
7626  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
7627  numbers::invalid_unsigned_int;
7628  const unsigned int dofs_per_face =
7629  fe_degree > -1 ? static_dofs_per_face :
7630  Utilities::pow(this->data->fe_degree + 1, dim - 1);
7631 
7632  constexpr unsigned int stack_array_size_threshold = 100;
7633 
7635  temp_data[static_dofs_per_face < stack_array_size_threshold ?
7636  n_components * 2 * static_dofs_per_face :
7637  1];
7638  VectorizedArray<Number> *temp1;
7639  if (static_dofs_per_face < stack_array_size_threshold)
7640  temp1 = &temp_data[0];
7641  else
7642  temp1 = this->scratch_data;
7643 
7644  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
7645  if (fe_degree > -1 &&
7646  this->subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
7647  this->data->element_type <=
7649  internal::FEFaceEvaluationImpl<
7650  true,
7651  dim,
7652  fe_degree,
7653  n_q_points_1d_actual,
7654  n_components,
7655  VectorizedArray<Number>>::integrate_in_face(*this->data,
7656  temp1,
7657  this->begin_values(),
7658  this->begin_gradients(),
7659  this->scratch_data +
7660  2 * n_components *
7661  dofs_per_face,
7662  integrate_values,
7663  integrate_gradients,
7664  this->subface_index);
7665  else
7666  internal::FEFaceEvaluationImpl<
7667  false,
7668  dim,
7669  fe_degree,
7670  n_q_points_1d_actual,
7671  n_components,
7672  VectorizedArray<Number>>::integrate_in_face(*this->data,
7673  temp1,
7674  this->begin_values(),
7675  this->begin_gradients(),
7676  this->scratch_data +
7677  2 * n_components *
7678  dofs_per_face,
7679  integrate_values,
7680  integrate_gradients,
7681  this->subface_index);
7682 
7683  internal::FEFaceNormalEvaluationImpl<dim,
7684  fe_degree,
7685  n_components,
7687  template interpolate<false, false>(
7688  *this->data, temp1, values_array, integrate_gradients, this->face_no);
7689 }
7690 
7691 
7692 
7693 template <int dim,
7694  int fe_degree,
7695  int n_q_points_1d,
7696  int n_components_,
7697  typename Number>
7698 template <typename VectorType>
7699 inline void
7701  gather_evaluate(const VectorType &input_vector,
7702  const bool evaluate_values,
7703  const bool evaluate_gradients)
7704 {
7705  const unsigned int side = this->face_no % 2;
7706 
7707  constexpr unsigned int static_dofs_per_face =
7708  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
7709  numbers::invalid_unsigned_int;
7710  const unsigned int dofs_per_face =
7711  fe_degree > -1 ? static_dofs_per_face :
7712  Utilities::pow(this->data->fe_degree + 1, dim - 1);
7713 
7714  constexpr unsigned int stack_array_size_threshold = 100;
7715 
7717  temp_data[static_dofs_per_face < stack_array_size_threshold ?
7718  n_components_ * 2 * dofs_per_face :
7719  1];
7720  VectorizedArray<Number> *__restrict temp1;
7721  if (static_dofs_per_face < stack_array_size_threshold)
7722  temp1 = &temp_data[0];
7723  else
7724  temp1 = this->scratch_data;
7725 
7726  internal::VectorReader<Number> reader;
7727  std::integral_constant<bool,
7728  internal::is_vectorizable<VectorType, Number>::value>
7729  vector_selector;
7730 
7731  // case 1: contiguous and interleaved indices
7732  if (((evaluate_gradients == false &&
7733  this->data->nodal_at_cell_boundaries == true) ||
7734  (this->data->element_type ==
7736  fe_degree > 1)) &&
7737  this->dof_info
7738  ->index_storage_variants[this->dof_access_index][this->cell] ==
7740  interleaved_contiguous)
7741  {
7743  this->dof_info
7744  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
7746  const unsigned int dof_index =
7747  this->dof_info
7748  ->dof_indices_contiguous[this->dof_access_index]
7749  [this->cell *
7751  this->dof_info
7752  ->component_dof_indices_offset[this->active_fe_index]
7753  [this->first_selected_component] *
7755 
7756  if (fe_degree > 1 && evaluate_gradients == true)
7757  {
7758  // we know that the gradient weights for the Hermite case on the
7759  // right (side==1) are the negative from the value at the left
7760  // (side==0), so we only read out one of them.
7761  const VectorizedArray<Number> grad_weight =
7762  this->data->shape_data_on_face[0][fe_degree + 1 + side];
7763  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
7764  2 * dofs_per_face);
7765  const unsigned int *index_array =
7766  &this->data->face_to_cell_index_hermite(this->face_no, 0);
7767  for (unsigned int i = 0; i < dofs_per_face; ++i)
7768  {
7769  const unsigned int ind1 = index_array[2 * i];
7770  const unsigned int ind2 = index_array[2 * i + 1];
7771  AssertIndexRange(ind1, dofs_per_cell);
7772  AssertIndexRange(ind2, dofs_per_cell);
7773  for (unsigned int comp = 0; comp < n_components_; ++comp)
7774  {
7775  reader.process_dofs_vectorized(
7776  1,
7777  dof_index + (ind1 + comp * static_dofs_per_component) *
7779  input_vector,
7780  temp1 + i + 2 * comp * dofs_per_face,
7781  vector_selector);
7782  reader.process_dofs_vectorized(
7783  1,
7784  dof_index + (ind2 + comp * static_dofs_per_component) *
7786  input_vector,
7787  temp1 + dofs_per_face + i + 2 * comp * dofs_per_face,
7788  vector_selector);
7789  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
7790  grad_weight *
7791  (temp1[i + 2 * comp * dofs_per_face] -
7792  temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
7793  }
7794  }
7795  }
7796  else
7797  {
7798  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
7799  dofs_per_face);
7800  const unsigned int *index_array =
7801  &this->data->face_to_cell_index_nodal(this->face_no, 0);
7802  for (unsigned int i = 0; i < dofs_per_face; ++i)
7803  {
7804  const unsigned int ind = index_array[i];
7805  for (unsigned int comp = 0; comp < n_components_; ++comp)
7806  reader.process_dofs_vectorized(
7807  1,
7808  dof_index + (ind + comp * static_dofs_per_component) *
7810  input_vector,
7811  temp1 + i + 2 * comp * dofs_per_face,
7812  vector_selector);
7813  }
7814  }
7815  }
7816 
7817  // case 2: contiguous and interleaved indices with fixed stride
7818  else if (((evaluate_gradients == false &&
7819  this->data->nodal_at_cell_boundaries == true) ||
7820  (this->data->element_type ==
7822  fe_degree > 1)) &&
7823  this->dof_info
7824  ->index_storage_variants[this->dof_access_index][this->cell] ==
7826  interleaved_contiguous_strided)
7827  {
7829  this->dof_info
7830  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
7832  const unsigned int *indices =
7833  &this->dof_info
7834  ->dof_indices_contiguous[this->dof_access_index]
7835  [this->cell *
7837  if (fe_degree > 1 && evaluate_gradients == true)
7838  {
7839  // we know that the gradient weights for the Hermite case on the
7840  // right (side==1) are the negative from the value at the left
7841  // (side==0), so we only read out one of them.
7842  const VectorizedArray<Number> grad_weight =
7843  this->data->shape_data_on_face[0][fe_degree + 1 + side];
7844  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
7845  2 * dofs_per_face);
7846 
7847  const unsigned int *index_array =
7848  &this->data->face_to_cell_index_hermite(this->face_no, 0);
7849  for (unsigned int i = 0; i < dofs_per_face; ++i)
7850  {
7851  const unsigned int ind1 =
7852  index_array[2 * i] * VectorizedArray<Number>::n_array_elements;
7853  const unsigned int ind2 =
7854  index_array[2 * i + 1] *
7856  for (unsigned int comp = 0; comp < n_components_; ++comp)
7857  {
7858  reader.process_dof_gather(
7859  indices,
7860  input_vector,
7861  ind1 +
7862  comp * static_dofs_per_component *
7864  this->dof_info->component_dof_indices_offset
7865  [this->active_fe_index]
7866  [this->first_selected_component] *
7868  temp1[i + 2 * comp * dofs_per_face],
7869  vector_selector);
7871  reader.process_dof_gather(
7872  indices,
7873  input_vector,
7874  ind2 +
7875  comp * static_dofs_per_component *
7877  this->dof_info->component_dof_indices_offset
7878  [this->active_fe_index]
7879  [this->first_selected_component] *
7881  grad,
7882  vector_selector);
7883  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
7884  grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
7885  }
7886  }
7887  }
7888  else
7889  {
7890  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
7891  dofs_per_face);
7892  const unsigned int *index_array =
7893  &this->data->face_to_cell_index_nodal(this->face_no, 0);
7894  for (unsigned int i = 0; i < dofs_per_face; ++i)
7895  {
7896  const unsigned int ind =
7898  for (unsigned int comp = 0; comp < n_components_; ++comp)
7899  reader.process_dof_gather(
7900  indices,
7901  input_vector,
7902  ind +
7903  comp * static_dofs_per_component *
7905  this->dof_info->component_dof_indices_offset
7906  [this->active_fe_index]
7907  [this->first_selected_component] *
7909  temp1[i + 2 * comp * dofs_per_face],
7910  vector_selector);
7911  }
7912  }
7913  }
7914 
7915  // case 3: contiguous and interleaved indices with mixed stride
7916  else if (((evaluate_gradients == false &&
7917  this->data->nodal_at_cell_boundaries == true) ||
7918  (this->data->element_type ==
7920  fe_degree > 1)) &&
7921  this->dof_info
7922  ->index_storage_variants[this->dof_access_index][this->cell] ==
7924  interleaved_contiguous_mixed_strides)
7925  {
7926  const unsigned int *strides =
7927  &this->dof_info->dof_indices_interleave_strides
7928  [this->dof_access_index]
7930  unsigned int indices[VectorizedArray<Number>::n_array_elements];
7931  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
7932  ++v)
7933  indices[v] =
7934  this->dof_info->dof_indices_contiguous
7935  [this->dof_access_index]
7936  [this->cell * VectorizedArray<Number>::n_array_elements + v] +
7937  this->dof_info
7938  ->component_dof_indices_offset[this->active_fe_index]
7939  [this->first_selected_component] *
7940  strides[v];
7941  const unsigned int nvec =
7942  this->dof_info
7943  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
7944 
7945  if (fe_degree > 1 && evaluate_gradients == true)
7946  {
7947  // we know that the gradient weights for the Hermite case on the
7948  // right (side==1) are the negative from the value at the left
7949  // (side==0), so we only read out one of them.
7950  const VectorizedArray<Number> grad_weight =
7951  this->data->shape_data_on_face[0][fe_degree + 1 + side];
7952  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
7953  2 * dofs_per_face);
7954 
7955  const unsigned int *index_array =
7956  &this->data->face_to_cell_index_hermite(this->face_no, 0);
7958  for (unsigned int comp = 0; comp < n_components_; ++comp)
7959  for (unsigned int i = 0; i < dofs_per_face; ++i)
7960  {
7961  unsigned int ind1[VectorizedArray<Number>::n_array_elements];
7962  DEAL_II_OPENMP_SIMD_PRAGMA
7963  for (unsigned int v = 0;
7964  v < VectorizedArray<Number>::n_array_elements;
7965  ++v)
7966  ind1[v] = indices[v] + (comp * static_dofs_per_component +
7967  index_array[2 * i]) *
7968  strides[v];
7969  unsigned int ind2[VectorizedArray<Number>::n_array_elements];
7970  DEAL_II_OPENMP_SIMD_PRAGMA
7971  for (unsigned int v = 0;
7972  v < VectorizedArray<Number>::n_array_elements;
7973  ++v)
7974  ind2[v] = indices[v] + (comp * static_dofs_per_component +
7975  index_array[2 * i + 1]) *
7976  strides[v];
7977  reader.process_dof_gather(ind1,
7978  input_vector,
7979  0,
7980  temp1[i + 2 * comp * dofs_per_face],
7981  vector_selector);
7983  reader.process_dof_gather(
7984  ind2, input_vector, 0, grad, vector_selector);
7985  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
7986  grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
7987  }
7988  else
7989  {
7990  for (unsigned int i = 0; i < n_components_ * 2 * dofs_per_face;
7991  ++i)
7992  temp1[i] = VectorizedArray<Number>();
7993  for (unsigned int v = 0; v < nvec; ++v)
7994  for (unsigned int comp = 0; comp < n_components_; ++comp)
7995  for (unsigned int i = 0; i < dofs_per_face; ++i)
7996  {
7997  const unsigned int ind1 =
7998  indices[v] + (comp * static_dofs_per_component +
7999  index_array[2 * i]) *
8000  strides[v];
8001  const unsigned int ind2 =
8002  indices[v] + (comp * static_dofs_per_component +
8003  index_array[2 * i + 1]) *
8004  strides[v];
8005  reader.process_dof(
8006  ind1,
8007  const_cast<VectorType &>(input_vector),
8008  temp1[i + 2 * comp * dofs_per_face][v]);
8009  Number grad;
8010  reader.process_dof(ind2,
8011  const_cast<VectorType &>(input_vector),
8012  grad);
8013  temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v] =
8014  grad_weight[0] *
8015  (temp1[i + 2 * comp * dofs_per_face][v] - grad);
8016  }
8017  }
8018  }
8019  else
8020  {
8021  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8022  dofs_per_face);
8023  const unsigned int *index_array =
8024  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8026  for (unsigned int comp = 0; comp < n_components_; ++comp)
8027  for (unsigned int i = 0; i < dofs_per_face; ++i)
8028  {
8030  DEAL_II_OPENMP_SIMD_PRAGMA
8031  for (unsigned int v = 0;
8032  v < VectorizedArray<Number>::n_array_elements;
8033  ++v)
8034  ind[v] = indices[v] + (comp * static_dofs_per_component +
8035  index_array[i]) *
8036  strides[v];
8037  reader.process_dof_gather(ind,
8038  input_vector,
8039  0,
8040  temp1[i + 2 * comp * dofs_per_face],
8041  vector_selector);
8042  }
8043  else
8044  {
8045  for (unsigned int i = 0; i < n_components_ * dofs_per_face; ++i)
8046  temp1[i] = VectorizedArray<Number>();
8047  for (unsigned int v = 0; v < nvec; ++v)
8048  for (unsigned int comp = 0; comp < n_components_; ++comp)
8049  for (unsigned int i = 0; i < dofs_per_face; ++i)
8050  {
8051  const unsigned int ind1 =
8052  indices[v] +
8053  (comp * static_dofs_per_component + index_array[i]) *
8054  strides[v];
8055  reader.process_dof(
8056  ind1,
8057  const_cast<VectorType &>(input_vector),
8058  temp1[i + 2 * comp * dofs_per_face][v]);
8059  }
8060  }
8061  }
8062  }
8063 
8064  // case 4: contiguous indices without interleaving
8065  else if (((evaluate_gradients == false &&
8066  this->data->nodal_at_cell_boundaries == true) ||
8067  (this->data->element_type ==
8069  fe_degree > 1)) &&
8070  this->dof_info
8071  ->index_storage_variants[this->dof_access_index][this->cell] ==
8073  contiguous &&
8074  this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
8075  [this->cell] ==
8077  {
8078  const unsigned int *indices =
8079  &this->dof_info
8080  ->dof_indices_contiguous[this->dof_access_index]
8081  [this->cell *
8083  if (evaluate_gradients == true &&
8084  this->data->element_type ==
8086  {
8087  // we know that the gradient weights for the Hermite case on the
8088  // right (side==1) are the negative from the value at the left
8089  // (side==0), so we only read out one of them.
8090  const VectorizedArray<Number> grad_weight =
8091  this->data->shape_data_on_face[0][fe_degree + 1 + side];
8092  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
8093  2 * dofs_per_face);
8094 
8095  const unsigned int *index_array =
8096  &this->data->face_to_cell_index_hermite(this->face_no, 0);
8097  for (unsigned int i = 0; i < dofs_per_face; ++i)
8098  {
8099  const unsigned int ind1 = index_array[2 * i];
8100  const unsigned int ind2 = index_array[2 * i + 1];
8101  for (unsigned int comp = 0; comp < n_components_; ++comp)
8102  {
8103  reader.process_dof_gather(
8104  indices,
8105  input_vector,
8106  ind1 + comp * static_dofs_per_component +
8107  this->dof_info->component_dof_indices_offset
8108  [this->active_fe_index][this->first_selected_component],
8109  temp1[i + 2 * comp * dofs_per_face],
8110  vector_selector);
8112  reader.process_dof_gather(
8113  indices,
8114  input_vector,
8115  ind2 + comp * static_dofs_per_component +
8116  this->dof_info->component_dof_indices_offset
8117  [this->active_fe_index][this->first_selected_component],
8118  grad,
8119  vector_selector);
8120  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
8121  grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
8122  }
8123  }
8124  }
8125  else
8126  {
8127  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8128  dofs_per_face);
8129  const unsigned int *index_array =
8130  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8131  for (unsigned int i = 0; i < dofs_per_face; ++i)
8132  for (unsigned int comp = 0; comp < n_components_; ++comp)
8133  {
8134  const unsigned int ind = index_array[i];
8135  reader.process_dof_gather(
8136  indices,
8137  input_vector,
8138  ind + comp * static_dofs_per_component +
8139  this->dof_info->component_dof_indices_offset
8140  [this->active_fe_index][this->first_selected_component],
8141  temp1[i + comp * 2 * dofs_per_face],
8142  vector_selector);
8143  }
8144  }
8145  }
8146 
8147  // case 5: default vector access
8148  else
8149  {
8150  this->read_dof_values(input_vector);
8151  internal::FEFaceNormalEvaluationImpl<dim,
8152  fe_degree,
8153  n_components_,
8155  template interpolate<true, false>(*this->data,
8156  this->values_dofs[0],
8157  temp1,
8158  evaluate_gradients,
8159  this->face_no);
8160  }
8161 
8162  if (fe_degree > -1 &&
8163  this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
8164  this->data->element_type <=
8166  internal::FEFaceEvaluationImpl<
8167  true,
8168  dim,
8169  fe_degree,
8170  n_q_points_1d,
8171  n_components_,
8172  VectorizedArray<Number>>::evaluate_in_face(*this->data,
8173  temp1,
8174  this->values_quad[0],
8175  this->gradients_quad[0][0],
8176  this->scratch_data +
8177  2 * n_components_ *
8178  dofs_per_face,
8179  evaluate_values,
8180  evaluate_gradients,
8181  this->subface_index);
8182  else
8183  internal::FEFaceEvaluationImpl<
8184  false,
8185  dim,
8186  fe_degree,
8187  n_q_points_1d,
8188  n_components_,
8189  VectorizedArray<Number>>::evaluate_in_face(*this->data,
8190  temp1,
8191  this->values_quad[0],
8192  this->gradients_quad[0][0],
8193  this->scratch_data +
8194  2 * n_components_ *
8195  dofs_per_face,
8196  evaluate_values,
8197  evaluate_gradients,
8198  this->subface_index);
8199 
8200  if (this->face_orientation)
8201  adjust_for_face_orientation(false, evaluate_values, evaluate_gradients);
8202 
8203 # ifdef DEBUG
8204  if (evaluate_values == true)
8205  this->values_quad_initialized = true;
8206  if (evaluate_gradients == true)
8207  this->gradients_quad_initialized = true;
8208 # endif
8209 }
8210 
8211 
8212 
8213 template <int dim,
8214  int fe_degree,
8215  int n_q_points_1d,
8216  int n_components_,
8217  typename Number>
8218 template <typename VectorType>
8219 inline void
8221  integrate_scatter(const bool integrate_values,
8222  const bool integrate_gradients,
8223  VectorType &destination)
8224 {
8225  const unsigned int side = this->face_no % 2;
8226  const unsigned int dofs_per_face =
8227  fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
8228  Utilities::pow(this->data->fe_degree + 1, dim - 1);
8229 
8230  constexpr unsigned int stack_array_size_threshold = 100;
8231 
8232  VectorizedArray<Number> temp_data[dofs_per_face < stack_array_size_threshold ?
8233  n_components_ * 2 * dofs_per_face :
8234  1];
8235  VectorizedArray<Number> *__restrict temp1;
8236  if (dofs_per_face < stack_array_size_threshold)
8237  temp1 = &temp_data[0];
8238  else
8239  temp1 = this->scratch_data;
8240 
8241  if (this->face_orientation)
8242  adjust_for_face_orientation(true, integrate_values, integrate_gradients);
8243  if (fe_degree > -1 &&
8244  this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
8245  this->data->element_type <=
8247  internal::FEFaceEvaluationImpl<
8248  true,
8249  dim,
8250  fe_degree,
8251  n_q_points_1d,
8252  n_components_,
8253  VectorizedArray<Number>>::integrate_in_face(*this->data,
8254  temp1,
8255  this->values_quad[0],
8256  this->gradients_quad[0][0],
8257  this->scratch_data +
8258  2 * n_components_ *
8259  dofs_per_face,
8260  integrate_values,
8261  integrate_gradients,
8262  this->subface_index);
8263  else
8264  internal::FEFaceEvaluationImpl<
8265  false,
8266  dim,
8267  fe_degree,
8268  n_q_points_1d,
8269  n_components_,
8270  VectorizedArray<Number>>::integrate_in_face(*this->data,
8271  temp1,
8272  this->values_quad[0],
8273  this->gradients_quad[0][0],
8274  this->scratch_data +
8275  2 * n_components_ *
8276  dofs_per_face,
8277  integrate_values,
8278  integrate_gradients,
8279  this->subface_index);
8280 
8281 # ifdef DEBUG
8282  this->dof_values_initialized = true;
8283 # endif
8284 
8285  internal::VectorDistributorLocalToGlobal<Number> writer;
8286  std::integral_constant<bool,
8287  internal::is_vectorizable<VectorType, Number>::value>
8288  vector_selector;
8289 
8290  // case 1: contiguous and interleaved indices
8291  if (((integrate_gradients == false &&
8292  this->data->nodal_at_cell_boundaries == true) ||
8293  (this->data->element_type ==
8295  fe_degree > 1)) &&
8296  this->dof_info
8297  ->index_storage_variants[this->dof_access_index][this->cell] ==
8299  interleaved_contiguous)
8300  {
8302  this->dof_info
8303  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
8305  const unsigned int dof_index =
8306  this->dof_info
8307  ->dof_indices_contiguous[this->dof_access_index]
8308  [this->cell *
8310  this->dof_info
8311  ->component_dof_indices_offset[this->active_fe_index]
8312  [this->first_selected_component] *
8314 
8315  if (fe_degree > 1 && integrate_gradients == true)
8316  {
8317  // we know that the gradient weights for the Hermite case on the
8318  // right (side==1) are the negative from the value at the left
8319  // (side==0), so we only read out one of them.
8320  const VectorizedArray<Number> grad_weight =
8321  this->data->shape_data_on_face[0][fe_degree + 2 - side];
8322  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
8323  2 * dofs_per_face);
8324  const unsigned int *index_array =
8325  &this->data->face_to_cell_index_hermite(this->face_no, 0);
8326  for (unsigned int i = 0; i < dofs_per_face; ++i)
8327  {
8328  const unsigned int ind1 = index_array[2 * i];
8329  const unsigned int ind2 = index_array[2 * i + 1];
8330  AssertIndexRange(ind1, dofs_per_cell);
8331  AssertIndexRange(ind2, dofs_per_cell);
8332  for (unsigned int comp = 0; comp < n_components_; ++comp)
8333  {
8335  temp1[i + 2 * comp * dofs_per_face] -
8336  grad_weight *
8337  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8339  grad_weight *
8340  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8341  writer.process_dofs_vectorized(
8342  1,
8343  dof_index + (ind1 + comp * static_dofs_per_component) *
8345  destination,
8346  &val,
8347  vector_selector);
8348  writer.process_dofs_vectorized(
8349  1,
8350  dof_index + (ind2 + comp * static_dofs_per_component) *
8352  destination,
8353  &grad,
8354  vector_selector);
8355  }
8356  }
8357  }
8358  else
8359  {
8360  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8361  dofs_per_face);
8362  const unsigned int *index_array =
8363  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8364  for (unsigned int i = 0; i < dofs_per_face; ++i)
8365  {
8366  const unsigned int ind = index_array[i];
8367  for (unsigned int comp = 0; comp < n_components_; ++comp)
8368  writer.process_dofs_vectorized(
8369  1,
8370  dof_index + (ind + comp * static_dofs_per_component) *
8372  destination,
8373  temp1 + i + 2 * comp * dofs_per_face,
8374  vector_selector);
8375  }
8376  }
8377  }
8378 
8379  // case 2: contiguous and interleaved indices with fixed stride
8380  else if (((integrate_gradients == false &&
8381  this->data->nodal_at_cell_boundaries == true) ||
8382  (this->data->element_type ==
8384  fe_degree > 1)) &&
8385  this->dof_info
8386  ->index_storage_variants[this->dof_access_index][this->cell] ==
8388  interleaved_contiguous_strided)
8389  {
8391  this->dof_info
8392  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
8394  const unsigned int *indices =
8395  &this->dof_info
8396  ->dof_indices_contiguous[this->dof_access_index]
8397  [this->cell *
8399  if (fe_degree > 1 && integrate_gradients == true)
8400  {
8401  // we know that the gradient weights for the Hermite case on the
8402  // right (side==1) are the negative from the value at the left
8403  // (side==0), so we only read out one of them.
8404  const VectorizedArray<Number> grad_weight =
8405  this->data->shape_data_on_face[0][fe_degree + 2 - side];
8406  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
8407  2 * dofs_per_face);
8408 
8409  const unsigned int *index_array =
8410  &this->data->face_to_cell_index_hermite(this->face_no, 0);
8411  for (unsigned int i = 0; i < dofs_per_face; ++i)
8412  {
8413  const unsigned int ind1 =
8414  index_array[2 * i] * VectorizedArray<Number>::n_array_elements;
8415  const unsigned int ind2 =
8416  index_array[2 * i + 1] *
8418  for (unsigned int comp = 0; comp < n_components_; ++comp)
8419  {
8421  temp1[i + 2 * comp * dofs_per_face] -
8422  grad_weight *
8423  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8425  grad_weight *
8426  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8427  writer.process_dof_gather(
8428  indices,
8429  destination,
8430  ind1 +
8431  comp * static_dofs_per_component *
8433  this->dof_info->component_dof_indices_offset
8434  [this->active_fe_index]
8435  [this->first_selected_component] *
8437  val,
8438  vector_selector);
8439  writer.process_dof_gather(
8440  indices,
8441  destination,
8442  ind2 +
8443  comp * static_dofs_per_component *
8445  this->dof_info->component_dof_indices_offset
8446  [this->active_fe_index]
8447  [this->first_selected_component] *
8449  grad,
8450  vector_selector);
8451  }
8452  }
8453  }
8454  else
8455  {
8456  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8457  dofs_per_face);
8458  const unsigned int *index_array =
8459  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8460  for (unsigned int i = 0; i < dofs_per_face; ++i)
8461  {
8462  const unsigned int ind =
8464  for (unsigned int comp = 0; comp < n_components_; ++comp)
8465  writer.process_dof_gather(
8466  indices,
8467  destination,
8468  ind +
8469  comp * static_dofs_per_component *
8471  this->dof_info->component_dof_indices_offset
8472  [this->active_fe_index]
8473  [this->first_selected_component] *
8475  temp1[i + 2 * comp * dofs_per_face],
8476  vector_selector);
8477  }
8478  }
8479  }
8480 
8481  // case 3: contiguous and interleaved indices with mixed stride
8482  else if (((integrate_gradients == false &&
8483  this->data->nodal_at_cell_boundaries == true) ||
8484  (this->data->element_type ==
8486  fe_degree > 1)) &&
8487  this->dof_info
8488  ->index_storage_variants[this->dof_access_index][this->cell] ==
8490  interleaved_contiguous_mixed_strides)
8491  {
8492  const unsigned int *strides =
8493  &this->dof_info->dof_indices_interleave_strides
8494  [this->dof_access_index]
8496  unsigned int indices[VectorizedArray<Number>::n_array_elements];
8497  for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
8498  ++v)
8499  indices[v] =
8500  this->dof_info->dof_indices_contiguous
8501  [this->dof_access_index]
8502  [this->cell * VectorizedArray<Number>::n_array_elements + v] +
8503  this->dof_info
8504  ->component_dof_indices_offset[this->active_fe_index]
8505  [this->first_selected_component] *
8506  strides[v];
8507  const unsigned int nvec =
8508  this->dof_info
8509  ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
8510 
8511  if (fe_degree > 1 && integrate_gradients == true)
8512  {
8513  // we know that the gradient weights for the Hermite case on the
8514  // right (side==1) are the negative from the value at the left
8515  // (side==0), so we only read out one of them.
8516  const VectorizedArray<Number> grad_weight =
8517  this->data->shape_data_on_face[0][fe_degree + 2 - side];
8518  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
8519  2 * dofs_per_face);
8520 
8521  const unsigned int *index_array =
8522  &this->data->face_to_cell_index_hermite(this->face_no, 0);
8524  for (unsigned int comp = 0; comp < n_components_; ++comp)
8525  for (unsigned int i = 0; i < dofs_per_face; ++i)
8526  {
8527  unsigned int ind1[VectorizedArray<Number>::n_array_elements];
8528  DEAL_II_OPENMP_SIMD_PRAGMA
8529  for (unsigned int v = 0;
8530  v < VectorizedArray<Number>::n_array_elements;
8531  ++v)
8532  ind1[v] = indices[v] + (comp * static_dofs_per_component +
8533  index_array[2 * i]) *
8534  strides[v];
8535  unsigned int ind2[VectorizedArray<Number>::n_array_elements];
8536  DEAL_II_OPENMP_SIMD_PRAGMA
8537  for (unsigned int v = 0;
8538  v < VectorizedArray<Number>::n_array_elements;
8539  ++v)
8540  ind2[v] = indices[v] + (comp * static_dofs_per_component +
8541  index_array[2 * i + 1]) *
8542  strides[v];
8544  temp1[i + 2 * comp * dofs_per_face] -
8545  grad_weight *
8546  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8548  grad_weight *
8549  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8550  writer.process_dof_gather(
8551  ind1, destination, 0, val, vector_selector);
8552  writer.process_dof_gather(
8553  ind2, destination, 0, grad, vector_selector);
8554  }
8555  else
8556  {
8557  for (unsigned int v = 0; v < nvec; ++v)
8558  for (unsigned int comp = 0; comp < n_components_; ++comp)
8559  for (unsigned int i = 0; i < dofs_per_face; ++i)
8560  {
8561  const unsigned int ind1 =
8562  indices[v] + (comp * static_dofs_per_component +
8563  index_array[2 * i]) *
8564  strides[v];
8565  const unsigned int ind2 =
8566  indices[v] + (comp * static_dofs_per_component +
8567  index_array[2 * i + 1]) *
8568  strides[v];
8569  Number val =
8570  temp1[i + 2 * comp * dofs_per_face][v] -
8571  grad_weight[0] * temp1[i + dofs_per_face +
8572  2 * comp * dofs_per_face][v];
8573  Number grad =
8574  grad_weight[0] *
8575  temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v];
8576  writer.process_dof(ind1, destination, val);
8577  writer.process_dof(ind2, destination, grad);
8578  }
8579  }
8580  }
8581  else
8582  {
8583  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8584  dofs_per_face);
8585  const unsigned int *index_array =
8586  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8588  for (unsigned int comp = 0; comp < n_components_; ++comp)
8589  for (unsigned int i = 0; i < dofs_per_face; ++i)
8590  {
8592  DEAL_II_OPENMP_SIMD_PRAGMA
8593  for (unsigned int v = 0;
8594  v < VectorizedArray<Number>::n_array_elements;
8595  ++v)
8596  ind[v] = indices[v] + (comp * static_dofs_per_component +
8597  index_array[i]) *
8598  strides[v];
8599  writer.process_dof_gather(ind,
8600  destination,
8601  0,
8602  temp1[i + 2 * comp * dofs_per_face],
8603  vector_selector);
8604  }
8605  else
8606  {
8607  for (unsigned int v = 0; v < nvec; ++v)
8608  for (unsigned int comp = 0; comp < n_components_; ++comp)
8609  for (unsigned int i = 0; i < dofs_per_face; ++i)
8610  {
8611  const unsigned int ind1 =
8612  indices[v] +
8613  (comp * static_dofs_per_component + index_array[i]) *
8614  strides[v];
8615  writer.process_dof(
8616  ind1,
8617  destination,
8618  temp1[i + 2 * comp * dofs_per_face][v]);
8619  }
8620  }
8621  }
8622  }
8623 
8624  // case 4: contiguous indices without interleaving
8625  else if (((integrate_gradients == false &&
8626  this->data->nodal_at_cell_boundaries == true) ||
8627  (this->data->element_type ==
8629  fe_degree > 1)) &&
8630  this->dof_info
8631  ->index_storage_variants[this->dof_access_index][this->cell] ==
8633  contiguous &&
8634  this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
8635  [this->cell] ==
8637  {
8638  const unsigned int *indices =
8639  &this->dof_info
8640  ->dof_indices_contiguous[this->dof_access_index]
8641  [this->cell *
8643 
8644  if (integrate_gradients == true &&
8645  this->data->element_type ==
8647  {
8648  // we know that the gradient weights for the Hermite case on the
8649  // right (side==1) are the negative from the value at the left
8650  // (side==0), so we only read out one of them.
8651  const VectorizedArray<Number> grad_weight =
8652  this->data->shape_data_on_face[0][fe_degree + 2 - side];
8653  AssertDimension(this->data->face_to_cell_index_hermite.size(1),
8654  2 * dofs_per_face);
8655  const unsigned int *index_array =
8656  &this->data->face_to_cell_index_hermite(this->face_no, 0);
8657  for (unsigned int i = 0; i < dofs_per_face; ++i)
8658  {
8659  const unsigned int ind1 = index_array[2 * i];
8660  const unsigned int ind2 = index_array[2 * i + 1];
8661  for (unsigned int comp = 0; comp < n_components_; ++comp)
8662  {
8664  temp1[i + 2 * comp * dofs_per_face] -
8665  grad_weight *
8666  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8668  grad_weight *
8669  temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
8670  writer.process_dof_gather(
8671  indices,
8672  destination,
8673  comp * static_dofs_per_component + ind1 +
8674  this->dof_info->component_dof_indices_offset
8675  [this->active_fe_index][this->first_selected_component],
8676  val,
8677  vector_selector);
8678  writer.process_dof_gather(
8679  indices,
8680  destination,
8681  comp * static_dofs_per_component + ind2 +
8682  this->dof_info->component_dof_indices_offset
8683  [this->active_fe_index][this->first_selected_component],
8684  grad,
8685  vector_selector);
8686  }
8687  }
8688  }
8689  else
8690  {
8691  AssertDimension(this->data->face_to_cell_index_nodal.size(1),
8692  dofs_per_face);
8693  const unsigned int *index_array =
8694  &this->data->face_to_cell_index_nodal(this->face_no, 0);
8695  for (unsigned int i = 0; i < dofs_per_face; ++i)
8696  {
8697  const unsigned int ind = index_array[i];
8698  for (unsigned int comp = 0; comp < n_components_; ++comp)
8699  writer.process_dof_gather(
8700  indices,
8701  destination,
8702  comp * static_dofs_per_component + ind +
8703  this->dof_info->component_dof_indices_offset
8704  [this->active_fe_index][this->first_selected_component],
8705  temp1[i + 2 * comp * dofs_per_face],
8706  vector_selector);
8707  }
8708  }
8709  }
8710 
8711  // case 5: default vector access
8712  else
8713  {
8714  internal::FEFaceNormalEvaluationImpl<dim,
8715  fe_degree,
8716  n_components_,
8718  template interpolate<false, false>(*this->data,
8719  temp1,
8720  this->values_dofs[0],
8721  integrate_gradients,
8722  this->face_no);
8723  this->distribute_local_to_global(destination);
8724  }
8725 }
8726 
8727 
8728 
8729 template <int dim,
8730  int fe_degree,
8731  int n_q_points_1d,
8732  int n_components,
8733  typename Number>
8734 inline void
8736  adjust_for_face_orientation(const bool integrate,
8737  const bool values,
8738  const bool gradients)
8739 {
8740  VectorizedArray<Number> *tmp_values = this->scratch_data;
8741  const unsigned int * orientations =
8742  &this->mapping_data->descriptor[this->active_fe_index]
8743  .face_orientations[this->face_orientation][0];
8744  for (unsigned int c = 0; c < n_components; ++c)
8745  {
8746  if (values == true)
8747  {
8748  if (integrate)
8749  for (unsigned int q = 0; q < n_q_points; ++q)
8750  tmp_values[q] = this->values_quad[c][orientations[q]];
8751  else
8752  for (unsigned int q = 0; q < n_q_points; ++q)
8753  tmp_values[orientations[q]] = this->values_quad[c][q];
8754  for (unsigned int q = 0; q < n_q_points; ++q)
8755  this->values_quad[c][q] = tmp_values[q];
8756  }
8757  if (gradients == true)
8758  for (unsigned int d = 0; d < dim; ++d)
8759  {
8760  if (integrate)
8761  for (unsigned int q = 0; q < n_q_points; ++q)
8762  tmp_values[q] = this->gradients_quad[c][d][orientations[q]];
8763  else
8764  for (unsigned int q = 0; q < n_q_points; ++q)
8765  tmp_values[orientations[q]] = this->gradients_quad[c][d][q];
8766  for (unsigned int q = 0; q < n_q_points; ++q)
8767  this->gradients_quad[c][d][q] = tmp_values[q];
8768  }
8769  }
8770 }
8771 
8772 
8773 
8774 template <int dim,
8775  int fe_degree,
8776  int n_q_points_1d,
8777  int n_components_,
8778  typename Number>
8781  quadrature_point(const unsigned int q) const
8782 {
8783  AssertIndexRange(q, n_q_points);
8784  if (this->dof_access_index < 2)
8785  {
8786  Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
8787  ExcNotImplemented());
8788  AssertIndexRange(this->cell,
8789  this->mapping_data->quadrature_point_offsets.size());
8790  return this->mapping_data->quadrature_points
8791  [this->mapping_data->quadrature_point_offsets[this->cell] + q];
8792  }
8793  else
8794  {
8795  Assert(this->matrix_info->get_mapping_info()
8796  .face_data_by_cells[this->quad_no]
8797  .quadrature_point_offsets.empty() == false,
8798  ExcNotImplemented());
8799  const unsigned int index =
8800  this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
8801  AssertIndexRange(index,
8802  this->matrix_info->get_mapping_info()
8803  .face_data_by_cells[this->quad_no]
8804  .quadrature_point_offsets.size());
8805  return this->matrix_info->get_mapping_info()
8806  .face_data_by_cells[this->quad_no]
8807  .quadrature_points[this->matrix_info->get_mapping_info()
8808  .face_data_by_cells[this->quad_no]
8809  .quadrature_point_offsets[index] +
8810  q];
8811  }
8812 }
8813 
8814 
8815 
8816 /*------------------------- end FEFaceEvaluation ------------------------- */
8817 
8818 
8819 #endif // ifndef DOXYGEN
8820 
8821 
8822 DEAL_II_NAMESPACE_CLOSE
8823 
8824 #endif
gradient_type get_gradient(const unsigned int q_point) const
const unsigned int active_quad_index
std::vector< unsigned int > plain_dof_indices
Definition: dof_info.h:505
static const unsigned int invalid_unsigned_int
Definition: types.h:173
const Number * constraint_pool_begin(const unsigned int pool_index) const
void store(Number *ptr) const
const VectorizedArray< Number > * begin_hessians() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
void submit_divergence(const VectorizedArray< Number > div_in, const unsigned int q_point)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
FEEvaluation & operator=(const FEEvaluation &other)
internal::MatrixFreeFunctions::GeometryType get_cell_type() const
internal::MatrixFreeFunctions::GeometryType cell_type
static ::ExceptionBase & ExcAccessToUninitializedField()
std::vector< unsigned int > component_to_base_index
Definition: dof_info.h:535
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< IndexStorageVariants > index_storage_variants[3]
Definition: dof_info.h:388
const unsigned int dofs_per_cell
unsigned int n_components() const
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
void integrate(const bool integrate_values, const bool integrate_gradients)
void read_write_operation(const VectorOperation &operation, VectorType *vectors[], const std::bitset< VectorizedArray< Number >::n_array_elements > &mask, const bool apply_constraints=true) const
Point< dim, VectorizedArray< Number > > quadrature_point(const unsigned int q_point) const
std::vector< unsigned int > dof_indices_interleave_strides[3]
Definition: dof_info.h:450
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
static constexpr unsigned int static_dofs_per_cell
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
std::shared_ptr< const Utilities::MPI::Partitioner > vector_partitioner
Definition: dof_info.h:469
ArrayView< VectorizedArray< Number > > get_scratch_data() const
FEEvaluationBase & operator=(const FEEvaluationBase &other)
const internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 :dim), dim, Number > * mapping_data
const unsigned int active_fe_index
const VectorizedArray< Number > * begin_dof_values() const
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > * data
STL namespace.
Transformed quadrature points.
static ::ExceptionBase & ExcNotInitialized()
std::vector< unsigned int > cell_partition_data
Definition: task_info.h:465
void read_dof_values(const VectorType &src, const unsigned int first_index=0)
static constexpr unsigned int n_components
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
AlignedVector< VectorizedArray< Number > > * scratch_data_array
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void scatter(const unsigned int *offsets, Number *base_ptr) const
std::vector< unsigned int > dof_indices
Definition: dof_info.h:413
SymmetricTensor< 2, dim, VectorizedArray< Number > > get_symmetric_gradient(const unsigned int q_point) const
value_type integrate_value() const
Definition: point.h:110
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:392
const unsigned int dofs_per_component
value_type get_value(const unsigned int q_point) const
void submit_value(const value_type val_in, const unsigned int q_point)
std::vector< QuadratureDescriptor > descriptor
Definition: mapping_info.h:173
VectorizedArray< Number > JxW(const unsigned int q_index) const
std::vector< unsigned char > n_vectorization_lanes_filled[3]
Definition: dof_info.h:461
AlignedVector< VectorizedArray< Number > > * acquire_scratch_data() const
void set_data_pointers()
const Tensor< 1, dim, VectorizedArray< Number > > * normal_vectors
T fixed_power(const T t)
Definition: utilities.h:1031
const Number * quadrature_weights
static constexpr unsigned int tensor_dofs_per_cell
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArray< Number >::n_array_elements > &mask=std::bitset< VectorizedArray< Number >::n_array_elements >().flip()) const
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0)
const unsigned int dofs_per_cell
internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index
static ::ExceptionBase & ExcMessage(std::string arg1)
const unsigned int n_q_points
gradient_type get_hessian_diagonal(const unsigned int q_point) const
value_type get_laplacian(const unsigned int q_point) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArray< Number >> grad_in, const unsigned int q_point)
const Tensor< 2, dim, VectorizedArray< Number > > * jacobian
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArray< Number > > & get_shape_info() const
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
Definition: fe.h:3235
T sum(const T &t, const MPI_Comm &mpi_communicator)
VectorizedArray< Number > get_divergence(const unsigned int q_point) const
std::vector< unsigned int > dof_indices_interleaved
Definition: dof_info.h:430
const unsigned int quad_no
VectorizedArray< Number > * scratch_data
#define Assert(cond, exc)
Definition: exceptions.h:1407
unsigned int element_multiplicity(const unsigned int index) const
Definition: fe.h:3134
UpdateFlags
Point< dim, VectorizedArray< Number > > quadrature_point(const unsigned int q_point) const
const internal::MatrixFreeFunctions::DoFInfo * dof_info
const VectorizedArray< Number > * begin_values() const
unsigned int subface_index
std::vector< unsigned int > row_starts_plain_indices
Definition: dof_info.h:495
void read_write_operation_contiguous(const VectorOperation &operation, VectorType *vectors[], const std::bitset< VectorizedArray< Number >::n_array_elements > &mask) const
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:428
const MatrixFree< dim, Number > * matrix_info
#define DeclException0(Exception0)
Definition: exceptions.h:473
std::vector< unsigned int > boundary_partition_data
Definition: task_info.h:483
const Number * constraint_pool_end(const unsigned int pool_index) const
Tensor< 2, dim, VectorizedArray< Number > > inverse_jacobian(const unsigned int q_index) const
void read_write_operation_global(const VectorOperation &operation, VectorType *vectors[]) const
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition: dof_info.h:396
unsigned int get_mapping_data_index_offset() const
void reinit(const unsigned int face_batch_number)
unsigned int face_no
VectorizedArray< Number > * gradients_quad[n_components][dim]
VectorizedArray< Number > pow(const ::VectorizedArray< Number > &x, const Number p)
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArray< Number >> curl_in, const unsigned int q_point)
static constexpr unsigned int static_dofs_per_cell
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
friend class FEEvaluationBase
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition: dof_info.h:425
const VectorizedArray< Number > * J_value
bool hessians_quad_initialized
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:383
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:180
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArray< Number > > > get_hessian(const unsigned int q_point) const
void reinit(const unsigned int cell_batch_index)
void adjust_for_face_orientation(const bool integrate, const bool values, const bool gradients)
unsigned int cell
void integrate(const bool integrate_values, const bool integrate_gradients)
static constexpr unsigned int n_components
std::vector< unsigned int > face_partition_data
Definition: task_info.h:474
unsigned int global_dof_index
Definition: types.h:89
static constexpr unsigned int static_n_q_points_cell
FEFaceEvaluation(const MatrixFree< dim, Number > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
Definition: cuda.h:31
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArray< Number >::n_array_elements > &mask=std::bitset< VectorizedArray< Number >::n_array_elements >().flip()) const
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
VectorizedArray< Number > * hessians_quad[n_components][(dim *(dim+1))/2]
value_type get_dof_value(const unsigned int dof) const
Definition: mpi.h:90
static constexpr unsigned int static_n_q_points
std::vector< unsigned int > dof_indices_contiguous[3]
Definition: dof_info.h:440
void load(const Number *ptr)
const unsigned int n_q_points
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition: dof_info.h:548
VectorizedArray< Number > * values_dofs[n_components]
const VectorizedArray< Number > * begin_gradients() const
Tensor< 1, dim, VectorizedArray< Number > > get_normal_vector(const unsigned int q_point) const
static constexpr unsigned int tensor_dofs_per_cell
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArray< Number >::n_array_elements > & get_face_info(const unsigned int face_batch_number) const
VectorizedArray< Number > read_cell_data(const AlignedVector< VectorizedArray< Number >> &array) const
unsigned int get_cell_data_number() const
static constexpr unsigned int dimension
unsigned int n_base_elements(const unsigned int dof_handler_index) const
static constexpr unsigned int static_dofs_per_component
const unsigned int n_quadrature_points
void gather(const Number *base_ptr, const unsigned int *offsets)
const unsigned int dofs_per_component
static ::ExceptionBase & ExcNotImplemented()
const Tensor< 1, dim, VectorizedArray< Number > > * normal_x_jacobian
std::vector< types::global_dof_index > local_dof_indices
value_type get_normal_derivative(const unsigned int q_point) const
void fill_JxW_values(AlignedVector< VectorizedArray< Number >> &JxW_values) const
bool indices_initialized() const
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
VectorizedArray< Number > * values_quad[n_components]
void submit_dof_value(const value_type val_in, const unsigned int dof)
bool gradients_quad_initialized
std::vector< unsigned int > start_components
Definition: dof_info.h:529
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
std::vector< unsigned int > lexicographic_numbering
Definition: shape_info.h:242
FEEvaluation(const MatrixFree< dim, Number > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
static constexpr unsigned int static_dofs_per_component
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
FEEvaluationAccess(const MatrixFree< dim, Number > &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)
unsigned int face_orientation
static constexpr unsigned int static_n_q_points
const std::vector< unsigned int > & get_internal_dof_numbering() const
const unsigned int n_fe_components
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArray< Number > > get_curl(const unsigned int q_point) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
const unsigned int first_selected_component
static ::ExceptionBase & ExcInternalError()
void release_scratch_data(const AlignedVector< VectorizedArray< Number >> *memory) const