Reference documentation for deal.II version GIT d6cf33b98c 2023-09-22 19:50:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_matrix_free_fe_evaluation_h
18 #define dealii_matrix_free_fe_evaluation_h
19 
20 
21 #include <deal.II/base/config.h>
22 
29 
31 
44 
45 #include <type_traits>
46 
47 
49 
50 
51 
89 template <int dim,
90  int n_components_,
91  typename Number,
92  bool is_face,
93  typename VectorizedArrayType>
95  : public FEEvaluationData<dim, VectorizedArrayType, is_face>
96 {
97 public:
98  using number_type = Number;
102  using hessian_type =
104  static constexpr unsigned int dimension = dim;
105  static constexpr unsigned int n_components = n_components_;
106  static constexpr unsigned int n_lanes = VectorizedArrayType::size();
107 
144  template <typename VectorType>
145  void
147  const VectorType &src,
148  const unsigned int first_index = 0,
149  const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
150 
179  template <typename VectorType>
180  void
182  const VectorType &src,
183  const unsigned int first_index = 0,
184  const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
185 
217  template <typename VectorType>
218  void
220  VectorType &dst,
221  const unsigned int first_index = 0,
222  const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
223 
262  template <typename VectorType>
263  void
265  VectorType &dst,
266  const unsigned int first_index = 0,
267  const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
268 
272  template <typename VectorType>
273  void
275  VectorType &dst,
276  const unsigned int first_index = 0,
277  const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
278 
301  value_type
302  get_dof_value(const unsigned int dof) const;
303 
314  void
315  submit_dof_value(const value_type val_in, const unsigned int dof);
316 
329  value_type
330  get_value(const unsigned int q_point) const;
331 
344  void
345  submit_value(const value_type val_in, const unsigned int q_point);
346 
358  get_gradient(const unsigned int q_point) const;
359 
374  value_type
375  get_normal_derivative(const unsigned int q_point) const;
376 
389  void
390  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
391 
410  void
412  const unsigned int q_point);
413 
426  void
427  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
428 
441  get_hessian(const unsigned int q_point) const;
442 
453  get_hessian_diagonal(const unsigned int q_point) const;
454 
466  value_type
467  get_laplacian(const unsigned int q_point) const;
468 
469 #ifdef DOXYGEN
470  // doxygen does not anyhow mention functions coming from partial template
471  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
472  // For now, hack in those functions manually only to fix documentation:
473 
480  VectorizedArrayType
481  get_divergence(const unsigned int q_point) const;
482 
492  get_symmetric_gradient(const unsigned int q_point) const;
493 
500  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
501  get_curl(const unsigned int q_point) const;
502 
518  void
519  submit_divergence(const VectorizedArrayType div_in,
520  const unsigned int q_point);
521 
538  void
541  const unsigned int q_point);
542 
555  void
557  const unsigned int q_point);
558 
559 #endif
560 
577  value_type
579 
587 
588 protected:
599  const unsigned int dof_no,
600  const unsigned int first_selected_component,
601  const unsigned int quad_no,
602  const unsigned int fe_degree,
603  const unsigned int n_q_points,
604  const bool is_interior_face,
605  const unsigned int active_fe_index,
606  const unsigned int active_quad_index,
607  const unsigned int face_type);
608 
646  const Mapping<dim> &mapping,
647  const FiniteElement<dim> &fe,
648  const Quadrature<1> &quadrature,
649  const UpdateFlags update_flags,
650  const unsigned int first_selected_component,
652 
660 
669 
674 
681  template <typename VectorType, typename VectorOperation>
682  void
684  const VectorOperation &operation,
685  const std::array<VectorType *, n_components_> &vectors,
686  const std::array<
688  n_components_> &vectors_sm,
689  const std::bitset<n_lanes> &mask,
690  const bool apply_constraints = true) const;
691 
699  template <typename VectorType, typename VectorOperation>
700  void
702  const VectorOperation &operation,
703  const std::array<VectorType *, n_components_> &vectors,
704  const std::array<
706  n_components_> &vectors_sm,
707  const std::bitset<n_lanes> &mask) const;
708 
716  template <typename VectorType, typename VectorOperation>
717  void
719  const VectorOperation &operation,
720  const std::array<VectorType *, n_components_> &vectors) const;
721 
725  void
727 
732 
737 
742  mutable std::vector<types::global_dof_index> local_dof_indices;
743 };
744 
745 
746 
754 template <int dim,
755  int n_components_,
756  typename Number,
757  bool is_face,
758  typename VectorizedArrayType = VectorizedArray<Number>>
760  n_components_,
761  Number,
762  is_face,
763  VectorizedArrayType>
764 {
765  static_assert(
766  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
767  "Type of Number and of VectorizedArrayType do not match.");
768 
769 public:
770  using number_type = Number;
774  static constexpr unsigned int dimension = dim;
775  static constexpr unsigned int n_components = n_components_;
776  using BaseClass =
778 
779 protected:
789  const unsigned int dof_no,
790  const unsigned int first_selected_component,
791  const unsigned int quad_no,
792  const unsigned int fe_degree,
793  const unsigned int n_q_points,
794  const bool is_interior_face = true,
795  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
797  const unsigned int face_type = numbers::invalid_unsigned_int);
798 
804  const Mapping<dim> &mapping,
805  const FiniteElement<dim> &fe,
806  const Quadrature<1> &quadrature,
807  const UpdateFlags update_flags,
808  const unsigned int first_selected_component,
810 
815 
821 };
822 
823 
824 
833 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
834 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
835  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
836 {
837  static_assert(
838  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
839  "Type of Number and of VectorizedArrayType do not match.");
840 
841 public:
842  using number_type = Number;
843  using value_type = VectorizedArrayType;
846  static constexpr unsigned int dimension = dim;
847  using BaseClass =
849 
853  value_type
854  get_dof_value(const unsigned int dof) const;
855 
859  void
860  submit_dof_value(const value_type val_in, const unsigned int dof);
861 
865  value_type
866  get_value(const unsigned int q_point) const;
867 
871  void
872  submit_value(const value_type val_in, const unsigned int q_point);
873 
877  void
879  const unsigned int q_point);
880 
885  get_gradient(const unsigned int q_point) const;
886 
890  value_type
891  get_normal_derivative(const unsigned int q_point) const;
892 
896  void
897  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
898 
902  void
904  const unsigned int q_point);
905 
910  get_hessian(unsigned int q_point) const;
911 
916  get_hessian_diagonal(const unsigned int q_point) const;
917 
921  void
922  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
923 
927  value_type
928  get_laplacian(const unsigned int q_point) const;
929 
933  value_type
935 
936 protected:
946  const unsigned int dof_no,
947  const unsigned int first_selected_component,
948  const unsigned int quad_no,
949  const unsigned int fe_degree,
950  const unsigned int n_q_points,
951  const bool is_interior_face = true,
952  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
954  const unsigned int face_type = numbers::invalid_unsigned_int);
955 
961  const Mapping<dim> &mapping,
962  const FiniteElement<dim> &fe,
963  const Quadrature<1> &quadrature,
964  const UpdateFlags update_flags,
965  const unsigned int first_selected_component,
967 
972 
978 };
979 
980 
981 
991 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
992 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
993  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
994 {
995  static_assert(
996  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
997  "Type of Number and of VectorizedArrayType do not match.");
998 
999 public:
1000  using number_type = Number;
1003  static constexpr unsigned int dimension = dim;
1004  static constexpr unsigned int n_components = dim;
1005  using BaseClass =
1007 
1011  value_type
1012  get_value(const unsigned int q_point) const;
1013 
1018  get_gradient(const unsigned int q_point) const;
1019 
1024  VectorizedArrayType
1025  get_divergence(const unsigned int q_point) const;
1026 
1034  get_symmetric_gradient(const unsigned int q_point) const;
1035 
1040  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1041  get_curl(const unsigned int q_point) const;
1042 
1047  get_hessian(const unsigned int q_point) const;
1048 
1053  get_hessian_diagonal(const unsigned int q_point) const;
1054 
1058  void
1060  const unsigned int q_point);
1061 
1065  void
1066  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1067 
1076  void
1078  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1079  const unsigned int q_point);
1080 
1089  void
1090  submit_divergence(const VectorizedArrayType div_in,
1091  const unsigned int q_point);
1092 
1101  void
1104  const unsigned int q_point);
1105 
1110  void
1112  const unsigned int q_point);
1113 
1114 protected:
1124  const unsigned int dof_no,
1125  const unsigned int first_selected_component,
1126  const unsigned int quad_no,
1127  const unsigned int dofs_per_cell,
1128  const unsigned int n_q_points,
1129  const bool is_interior_face = true,
1130  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1131  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1132  const unsigned int face_type = numbers::invalid_unsigned_int);
1133 
1139  const Mapping<dim> &mapping,
1140  const FiniteElement<dim> &fe,
1141  const Quadrature<1> &quadrature,
1142  const UpdateFlags update_flags,
1143  const unsigned int first_selected_component,
1145 
1150 
1156 };
1157 
1158 
1167 template <typename Number, bool is_face, typename VectorizedArrayType>
1168 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1169  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1170 {
1171  static_assert(
1172  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1173  "Type of Number and of VectorizedArrayType do not match.");
1174 
1175 public:
1176  using number_type = Number;
1177  using value_type = VectorizedArrayType;
1180  static constexpr unsigned int dimension = 1;
1181  using BaseClass =
1183 
1187  value_type
1188  get_dof_value(const unsigned int dof) const;
1189 
1193  void
1194  submit_dof_value(const value_type val_in, const unsigned int dof);
1195 
1199  value_type
1200  get_value(const unsigned int q_point) const;
1201 
1205  void
1206  submit_value(const value_type val_in, const unsigned int q_point);
1207 
1211  void
1212  submit_value(const gradient_type val_in, const unsigned int q_point);
1213 
1218  get_gradient(const unsigned int q_point) const;
1219 
1223  value_type
1224  get_divergence(const unsigned int q_point) const;
1225 
1229  value_type
1230  get_normal_derivative(const unsigned int q_point) const;
1231 
1235  void
1236  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1237 
1241  void
1242  submit_gradient(const value_type grad_in, const unsigned int q_point);
1243 
1247  void
1249  const unsigned int q_point);
1250 
1254  void
1256  const unsigned int q_point);
1257 
1261  void
1263  const unsigned int q_point);
1264 
1268  hessian_type
1269  get_hessian(unsigned int q_point) const;
1270 
1275  get_hessian_diagonal(const unsigned int q_point) const;
1276 
1280  void
1281  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1282 
1286  value_type
1287  get_laplacian(const unsigned int q_point) const;
1288 
1292  value_type
1294 
1295 protected:
1305  const unsigned int dof_no,
1306  const unsigned int first_selected_component,
1307  const unsigned int quad_no,
1308  const unsigned int fe_degree,
1309  const unsigned int n_q_points,
1310  const bool is_interior_face = true,
1311  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1312  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1313  const unsigned int face_type = numbers::invalid_unsigned_int);
1314 
1320  const Mapping<1> &mapping,
1321  const FiniteElement<1> &fe,
1322  const Quadrature<1> &quadrature,
1323  const UpdateFlags update_flags,
1324  const unsigned int first_selected_component,
1326 
1331 
1337 };
1338 
1339 
1340 
1904 template <int dim,
1905  int fe_degree,
1906  int n_q_points_1d,
1907  int n_components_,
1908  typename Number,
1909  typename VectorizedArrayType>
1911  n_components_,
1912  Number,
1913  false,
1914  VectorizedArrayType>
1915 {
1916  static_assert(
1917  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1918  "Type of Number and of VectorizedArrayType do not match.");
1919 
1920 public:
1924  using BaseClass =
1926 
1930  using number_type = Number;
1931 
1938 
1945 
1949  static constexpr unsigned int dimension = dim;
1950 
1955  static constexpr unsigned int n_components = n_components_;
1956 
1960  static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1961 
1970  static constexpr unsigned int static_n_q_points =
1971  Utilities::pow(n_q_points_1d, dim);
1972 
1982  static constexpr unsigned int static_dofs_per_component =
1983  Utilities::pow(fe_degree + 1, dim);
1984 
1994  static constexpr unsigned int tensor_dofs_per_cell =
1996 
2006  static constexpr unsigned int static_dofs_per_cell =
2008 
2045  const unsigned int dof_no = 0,
2046  const unsigned int quad_no = 0,
2047  const unsigned int first_selected_component = 0,
2048  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2049  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2050 
2059  const std::pair<unsigned int, unsigned int> &range,
2060  const unsigned int dof_no = 0,
2061  const unsigned int quad_no = 0,
2062  const unsigned int first_selected_component = 0);
2063 
2092  FEEvaluation(const Mapping<dim> &mapping,
2093  const FiniteElement<dim> &fe,
2094  const Quadrature<1> &quadrature,
2095  const UpdateFlags update_flags,
2096  const unsigned int first_selected_component = 0);
2097 
2104  const Quadrature<1> &quadrature,
2105  const UpdateFlags update_flags,
2106  const unsigned int first_selected_component = 0);
2107 
2120  const unsigned int first_selected_component = 0);
2121 
2129 
2136  FEEvaluation &
2137  operator=(const FEEvaluation &other);
2138 
2147  void
2148  reinit(const unsigned int cell_batch_index);
2149 
2156  void
2157  reinit(const std::array<unsigned int, n_lanes> &cell_ids);
2158 
2171  template <bool level_dof_access>
2172  void
2174 
2185  void
2187 
2191  static bool
2192  fast_evaluation_supported(const unsigned int given_degree,
2193  const unsigned int given_n_q_points_1d);
2194 
2204  void
2206 
2219  void
2220  evaluate(const VectorizedArrayType *values_array,
2221  const EvaluationFlags::EvaluationFlags evaluation_flag);
2222 
2236  template <typename VectorType>
2237  void
2238  gather_evaluate(const VectorType &input_vector,
2239  const EvaluationFlags::EvaluationFlags evaluation_flag);
2240 
2251  void
2253 
2265  void
2267  VectorizedArrayType *values_array,
2268  const bool sum_into_values = false);
2269 
2283  template <typename VectorType>
2284  void
2286  VectorType &output_vector);
2287 
2294  dof_indices() const;
2295 
2302  const unsigned int dofs_per_component;
2303 
2310  const unsigned int dofs_per_cell;
2311 
2319  const unsigned int n_q_points;
2320 
2321 private:
2326  void
2327  check_template_arguments(const unsigned int fe_no,
2328  const unsigned int first_selected_component);
2329 };
2330 
2331 
2332 
2368 template <int dim,
2369  int fe_degree,
2370  int n_q_points_1d = fe_degree + 1,
2371  int n_components_ = 1,
2372  typename Number = double,
2373  typename VectorizedArrayType = VectorizedArray<Number>>
2375  n_components_,
2376  Number,
2377  true,
2378  VectorizedArrayType>
2379 {
2380  static_assert(
2381  std::is_same_v<Number, typename VectorizedArrayType::value_type>,
2382  "Type of Number and of VectorizedArrayType do not match.");
2383 
2384 public:
2388  using BaseClass =
2390 
2394  using number_type = Number;
2395 
2402 
2409 
2413  static constexpr unsigned int dimension = dim;
2414 
2419  static constexpr unsigned int n_components = n_components_;
2420 
2424  static constexpr unsigned int n_lanes = VectorizedArrayType::size();
2425 
2435  static constexpr unsigned int static_n_q_points =
2436  Utilities::pow(n_q_points_1d, dim - 1);
2437 
2446  static constexpr unsigned int static_n_q_points_cell =
2447  Utilities::pow(n_q_points_1d, dim);
2448 
2457  static constexpr unsigned int static_dofs_per_component =
2458  Utilities::pow(fe_degree + 1, dim);
2459 
2468  static constexpr unsigned int tensor_dofs_per_cell =
2470 
2479  static constexpr unsigned int static_dofs_per_cell =
2481 
2525  const bool is_interior_face = true,
2526  const unsigned int dof_no = 0,
2527  const unsigned int quad_no = 0,
2528  const unsigned int first_selected_component = 0,
2529  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2530  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
2531  const unsigned int face_type = numbers::invalid_unsigned_int);
2532 
2542  const std::pair<unsigned int, unsigned int> &range,
2543  const bool is_interior_face = true,
2544  const unsigned int dof_no = 0,
2545  const unsigned int quad_no = 0,
2546  const unsigned int first_selected_component = 0);
2547 
2558  void
2559  reinit(const unsigned int face_batch_number);
2560 
2568  void
2569  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2570 
2574  static bool
2575  fast_evaluation_supported(const unsigned int given_degree,
2576  const unsigned int given_n_q_points_1d);
2577 
2588  void
2590 
2603  void
2604  evaluate(const VectorizedArrayType *values_array,
2605  const EvaluationFlags::EvaluationFlags evaluation_flag);
2606 
2618  template <typename VectorType>
2619  void
2620  gather_evaluate(const VectorType &input_vector,
2621  const EvaluationFlags::EvaluationFlags evaluation_flag);
2622 
2632  void
2634 
2643  void
2645  VectorizedArrayType *values_array);
2646 
2658  template <typename VectorType>
2659  void
2661  VectorType &output_vector);
2662 
2666  template <typename VectorType>
2667  void
2668  integrate_scatter(const bool integrate_values,
2669  const bool integrate_gradients,
2670  VectorType &output_vector);
2671 
2678  dof_indices() const;
2679 
2684  bool
2685  at_boundary() const;
2686 
2700  boundary_id() const;
2701 
2708  const unsigned int dofs_per_component;
2709 
2716  const unsigned int dofs_per_cell;
2717 
2725  const unsigned int n_q_points;
2726 };
2727 
2728 
2729 
2730 namespace internal
2731 {
2732  namespace MatrixFreeFunctions
2733  {
2734  // a helper function to compute the number of DoFs of a DGP element at
2735  // compile time, depending on the degree
2736  template <int dim, int degree>
2738  {
2739  // this division is always without remainder
2740  static constexpr unsigned int value =
2741  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2742  };
2743 
2744  // base specialization: 1d elements have 'degree+1' degrees of freedom
2745  template <int degree>
2746  struct DGP_dofs_per_component<1, degree>
2747  {
2748  static constexpr unsigned int value = degree + 1;
2749  };
2750  } // namespace MatrixFreeFunctions
2751 } // namespace internal
2752 
2753 
2754 /*----------------------- Inline functions ----------------------------------*/
2755 
2756 #ifndef DOXYGEN
2757 
2758 
2759 namespace internal
2760 {
2761  // Extract all internal data pointers and indices in a single function that
2762  // get passed on to the constructor of FEEvaluationData, avoiding to look
2763  // things up multiple times
2764  template <bool is_face,
2765  int dim,
2766  typename Number,
2767  typename VectorizedArrayType>
2769  InitializationData
2770  extract_initialization_data(
2772  const unsigned int dof_no,
2773  const unsigned int first_selected_component,
2774  const unsigned int quad_no,
2775  const unsigned int fe_degree,
2776  const unsigned int n_q_points,
2777  const unsigned int active_fe_index_given,
2778  const unsigned int active_quad_index_given,
2779  const unsigned int face_type)
2780  {
2782  InitializationData init_data;
2783 
2784  init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2785  init_data.mapping_data =
2786  &internal::MatrixFreeFunctions::
2787  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2788  matrix_free.get_mapping_info(), quad_no);
2789 
2790  init_data.active_fe_index =
2791  fe_degree != numbers::invalid_unsigned_int ?
2792  init_data.dof_info->fe_index_from_degree(first_selected_component,
2793  fe_degree) :
2794  (active_fe_index_given != numbers::invalid_unsigned_int ?
2795  active_fe_index_given :
2796  0);
2797  init_data.active_quad_index =
2798  fe_degree == numbers::invalid_unsigned_int ?
2799  (active_quad_index_given != numbers::invalid_unsigned_int ?
2800  active_quad_index_given :
2801  std::min<unsigned int>(init_data.active_fe_index,
2802  init_data.mapping_data->descriptor.size() -
2803  1)) :
2804  init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2805 
2806  init_data.shape_info = &matrix_free.get_shape_info(
2807  dof_no,
2808  quad_no,
2809  init_data.dof_info->component_to_base_index[first_selected_component],
2810  init_data.active_fe_index,
2811  init_data.active_quad_index);
2812  init_data.descriptor =
2813  &init_data.mapping_data->descriptor
2814  [is_face ?
2815  (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2816  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2817  init_data.active_quad_index];
2818 
2819  return init_data;
2820  }
2821 } // namespace internal
2822 
2823 
2824 
2825 /*----------------------- FEEvaluationBase ----------------------------------*/
2826 
2827 template <int dim,
2828  int n_components_,
2829  typename Number,
2830  bool is_face,
2831  typename VectorizedArrayType>
2832 inline FEEvaluationBase<dim,
2833  n_components_,
2834  Number,
2835  is_face,
2836  VectorizedArrayType>::
2837  FEEvaluationBase(
2839  const unsigned int dof_no,
2840  const unsigned int first_selected_component,
2841  const unsigned int quad_no,
2842  const unsigned int fe_degree,
2843  const unsigned int n_q_points,
2844  const bool is_interior_face,
2845  const unsigned int active_fe_index,
2846  const unsigned int active_quad_index,
2847  const unsigned int face_type)
2848  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2849  internal::extract_initialization_data<is_face>(matrix_free,
2850  dof_no,
2851  first_selected_component,
2852  quad_no,
2853  fe_degree,
2854  n_q_points,
2855  active_fe_index,
2856  active_quad_index,
2857  face_type),
2858  is_interior_face,
2859  quad_no,
2860  first_selected_component)
2861  , scratch_data_array(matrix_free.acquire_scratch_data())
2862  , matrix_free(&matrix_free)
2863 {
2864  this->set_data_pointers(scratch_data_array, n_components_);
2865  Assert(
2866  this->dof_info->start_components.back() == 1 ||
2867  static_cast<int>(n_components_) <=
2868  static_cast<int>(
2869  this->dof_info->start_components
2870  [this->dof_info->component_to_base_index[first_selected_component] +
2871  1]) -
2872  first_selected_component,
2873  ExcMessage(
2874  "You tried to construct a vector-valued evaluator with " +
2875  std::to_string(n_components) +
2876  " components. However, "
2877  "the current base element has only " +
2879  this->dof_info->start_components
2880  [this->dof_info->component_to_base_index[first_selected_component] +
2881  1] -
2882  first_selected_component) +
2883  " components left when starting from local element index " +
2885  first_selected_component -
2886  this->dof_info->start_components
2887  [this->dof_info->component_to_base_index[first_selected_component]]) +
2888  " (global index " + std::to_string(first_selected_component) + ")"));
2889 
2890  // do not check for correct dimensions of data fields here, should be done
2891  // in derived classes
2892 }
2893 
2894 
2895 
2896 template <int dim,
2897  int n_components_,
2898  typename Number,
2899  bool is_face,
2900  typename VectorizedArrayType>
2901 inline FEEvaluationBase<dim,
2902  n_components_,
2903  Number,
2904  is_face,
2905  VectorizedArrayType>::
2906  FEEvaluationBase(
2907  const Mapping<dim> &mapping,
2908  const FiniteElement<dim> &fe,
2909  const Quadrature<1> &quadrature,
2910  const UpdateFlags update_flags,
2911  const unsigned int first_selected_component,
2913  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2914  other != nullptr &&
2915  other->mapped_geometry->get_quadrature() == quadrature ?
2916  other->mapped_geometry :
2917  std::make_shared<internal::MatrixFreeFunctions::
2918  MappingDataOnTheFly<dim, VectorizedArrayType>>(
2919  mapping,
2920  quadrature,
2921  update_flags),
2922  n_components_,
2923  first_selected_component)
2924  , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2925  , matrix_free(nullptr)
2926 {
2927  const unsigned int base_element_number =
2928  fe.component_to_base_index(first_selected_component).first;
2929  Assert(fe.element_multiplicity(base_element_number) == 1 ||
2930  fe.element_multiplicity(base_element_number) -
2931  first_selected_component >=
2932  n_components_,
2933  ExcMessage("The underlying element must at least contain as many "
2934  "components as requested by this class"));
2935  (void)base_element_number;
2936 
2937  Assert(this->data == nullptr, ExcInternalError());
2938  this->data =
2940  Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2941  fe,
2942  fe.component_to_base_index(first_selected_component).first);
2943 
2944  this->set_data_pointers(scratch_data_array, n_components_);
2945 }
2946 
2947 
2948 
2949 template <int dim,
2950  int n_components_,
2951  typename Number,
2952  bool is_face,
2953  typename VectorizedArrayType>
2954 inline FEEvaluationBase<dim,
2955  n_components_,
2956  Number,
2957  is_face,
2958  VectorizedArrayType>::
2959  FEEvaluationBase(const FEEvaluationBase<dim,
2960  n_components_,
2961  Number,
2962  is_face,
2963  VectorizedArrayType> &other)
2964  : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
2965  , scratch_data_array(other.matrix_free == nullptr ?
2966  new AlignedVector<VectorizedArrayType>() :
2967  other.matrix_free->acquire_scratch_data())
2968  , matrix_free(other.matrix_free)
2969 {
2970  if (other.matrix_free == nullptr)
2971  {
2972  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
2973  this->data =
2975  *other.data);
2976 
2977  // Create deep copy of mapped geometry for use in parallel
2978  this->mapped_geometry =
2979  std::make_shared<internal::MatrixFreeFunctions::
2980  MappingDataOnTheFly<dim, VectorizedArrayType>>(
2981  other.mapped_geometry->get_fe_values().get_mapping(),
2982  other.mapped_geometry->get_quadrature(),
2983  other.mapped_geometry->get_fe_values().get_update_flags());
2984  this->mapping_data = &this->mapped_geometry->get_data_storage();
2985  this->cell = 0;
2986 
2987  this->jacobian =
2988  this->mapped_geometry->get_data_storage().jacobians[0].begin();
2989  this->J_value =
2990  this->mapped_geometry->get_data_storage().JxW_values.begin();
2991  this->jacobian_gradients =
2992  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
2993  this->jacobian_gradients_non_inverse =
2994  this->mapped_geometry->get_data_storage()
2995  .jacobian_gradients_non_inverse[0]
2996  .begin();
2997  this->quadrature_points =
2998  this->mapped_geometry->get_data_storage().quadrature_points.begin();
2999  }
3000 
3001  this->set_data_pointers(scratch_data_array, n_components_);
3002 }
3003 
3004 
3005 
3006 template <int dim,
3007  int n_components_,
3008  typename Number,
3009  bool is_face,
3010  typename VectorizedArrayType>
3011 inline FEEvaluationBase<dim,
3012  n_components_,
3013  Number,
3014  is_face,
3015  VectorizedArrayType> &
3017 operator=(const FEEvaluationBase<dim,
3018  n_components_,
3019  Number,
3020  is_face,
3021  VectorizedArrayType> &other)
3022 {
3023  // release old memory
3024  if (matrix_free == nullptr)
3025  {
3026  delete this->data;
3027  delete scratch_data_array;
3028  }
3029  else
3030  {
3031  matrix_free->release_scratch_data(scratch_data_array);
3032  }
3033 
3035 
3036  matrix_free = other.matrix_free;
3037 
3038  if (other.matrix_free == nullptr)
3039  {
3040  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3041  this->data =
3043  *other.data);
3044  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3045 
3046  // Create deep copy of mapped geometry for use in parallel
3047  this->mapped_geometry =
3048  std::make_shared<internal::MatrixFreeFunctions::
3049  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3050  other.mapped_geometry->get_fe_values().get_mapping(),
3051  other.mapped_geometry->get_quadrature(),
3052  other.mapped_geometry->get_fe_values().get_update_flags());
3053  this->cell = 0;
3054  this->mapping_data = &this->mapped_geometry->get_data_storage();
3055  this->jacobian =
3056  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3057  this->J_value =
3058  this->mapped_geometry->get_data_storage().JxW_values.begin();
3059  this->jacobian_gradients =
3060  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3061  this->jacobian_gradients_non_inverse =
3062  this->mapped_geometry->get_data_storage()
3063  .jacobian_gradients_non_inverse[0]
3064  .begin();
3065  this->quadrature_points =
3066  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3067  }
3068  else
3069  {
3070  scratch_data_array = matrix_free->acquire_scratch_data();
3071  }
3072 
3073  this->set_data_pointers(scratch_data_array, n_components_);
3074 
3075  return *this;
3076 }
3077 
3078 
3079 
3080 template <int dim,
3081  int n_components_,
3082  typename Number,
3083  bool is_face,
3084  typename VectorizedArrayType>
3085 inline FEEvaluationBase<dim,
3086  n_components_,
3087  Number,
3088  is_face,
3089  VectorizedArrayType>::~FEEvaluationBase()
3090 {
3091  if (matrix_free != nullptr)
3092  {
3093  try
3094  {
3095  matrix_free->release_scratch_data(scratch_data_array);
3096  }
3097  catch (...)
3098  {}
3099  }
3100  else
3101  {
3102  delete scratch_data_array;
3103  delete this->data;
3104  }
3105 }
3106 
3107 
3108 
3109 template <int dim,
3110  int n_components_,
3111  typename Number,
3112  bool is_face,
3113  typename VectorizedArrayType>
3116  get_matrix_free() const
3117 {
3118  Assert(matrix_free != nullptr,
3119  ExcMessage(
3120  "FEEvaluation was not initialized with a MatrixFree object!"));
3121  return *matrix_free;
3122 }
3123 
3124 
3125 
3126 namespace internal
3127 {
3128  // given a block vector return the underlying vector type
3129  // including constness (specified by bool)
3130  template <typename VectorType, bool>
3131  struct ConstBlockVectorSelector;
3132 
3133  template <typename VectorType>
3134  struct ConstBlockVectorSelector<VectorType, true>
3135  {
3136  using BaseVectorType = const typename VectorType::BlockType;
3137  };
3138 
3139  template <typename VectorType>
3140  struct ConstBlockVectorSelector<VectorType, false>
3141  {
3142  using BaseVectorType = typename VectorType::BlockType;
3143  };
3144 
3145  // allows to select between block vectors and non-block vectors, which
3146  // allows to use a unified interface for extracting blocks on block vectors
3147  // and doing nothing on usual vectors
3148  template <typename VectorType, bool>
3149  struct BlockVectorSelector;
3150 
3151  template <typename VectorType>
3152  struct BlockVectorSelector<VectorType, true>
3153  {
3154  using BaseVectorType = typename ConstBlockVectorSelector<
3155  VectorType,
3156  std::is_const_v<VectorType>>::BaseVectorType;
3157 
3158  static BaseVectorType *
3159  get_vector_component(VectorType &vec, const unsigned int component)
3160  {
3161  AssertIndexRange(component, vec.n_blocks());
3162  return &vec.block(component);
3163  }
3164  };
3165 
3166  template <typename VectorType>
3167  struct BlockVectorSelector<VectorType, false>
3168  {
3169  using BaseVectorType = VectorType;
3170 
3171  static BaseVectorType *
3172  get_vector_component(VectorType &vec, const unsigned int component)
3173  {
3174  // FEEvaluation allows to combine several vectors from a scalar
3175  // FiniteElement into a "vector-valued" FEEvaluation object with
3176  // multiple components. These components can be extracted with the other
3177  // get_vector_component functions. If we do not get a vector of vectors
3178  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3179  // must make sure that we do not duplicate the components in input
3180  // and/or duplicate the resulting integrals. In such a case, we should
3181  // only get the zeroth component in the vector contained set nullptr for
3182  // the others which allows us to catch unintended use in
3183  // read_write_operation.
3184  if (component == 0)
3185  return &vec;
3186  else
3187  return nullptr;
3188  }
3189  };
3190 
3191  template <typename VectorType>
3192  struct BlockVectorSelector<std::vector<VectorType>, false>
3193  {
3194  using BaseVectorType = VectorType;
3195 
3196  static BaseVectorType *
3197  get_vector_component(std::vector<VectorType> &vec,
3198  const unsigned int component)
3199  {
3200  AssertIndexRange(component, vec.size());
3201  return &vec[component];
3202  }
3203  };
3204 
3205  template <typename VectorType>
3206  struct BlockVectorSelector<const std::vector<VectorType>, false>
3207  {
3208  using BaseVectorType = const VectorType;
3209 
3210  static const BaseVectorType *
3211  get_vector_component(const std::vector<VectorType> &vec,
3212  const unsigned int component)
3213  {
3214  AssertIndexRange(component, vec.size());
3215  return &vec[component];
3216  }
3217  };
3218 
3219  template <typename VectorType>
3220  struct BlockVectorSelector<std::vector<VectorType *>, false>
3221  {
3222  using BaseVectorType = VectorType;
3223 
3224  static BaseVectorType *
3225  get_vector_component(std::vector<VectorType *> &vec,
3226  const unsigned int component)
3227  {
3228  AssertIndexRange(component, vec.size());
3229  return vec[component];
3230  }
3231  };
3232 
3233  template <typename VectorType>
3234  struct BlockVectorSelector<const std::vector<VectorType *>, false>
3235  {
3236  using BaseVectorType = const VectorType;
3237 
3238  static const BaseVectorType *
3239  get_vector_component(const std::vector<VectorType *> &vec,
3240  const unsigned int component)
3241  {
3242  AssertIndexRange(component, vec.size());
3243  return vec[component];
3244  }
3245  };
3246 } // namespace internal
3247 
3248 
3249 
3250 template <int dim,
3251  int n_components_,
3252  typename Number,
3253  bool is_face,
3254  typename VectorizedArrayType>
3255 template <typename VectorType, typename VectorOperation>
3256 inline void
3259  const VectorOperation &operation,
3260  const std::array<VectorType *, n_components_> &src,
3261  const std::array<
3263  n_components_> &src_sm,
3264  const std::bitset<n_lanes> &mask,
3265  const bool apply_constraints) const
3266 {
3267  // Case 1: No MatrixFree object given, simple case because we do not need to
3268  // process constraints and need not care about vectorization -> go to
3269  // separate function
3270  if (this->matrix_free == nullptr)
3271  {
3272  read_write_operation_global(operation, src);
3273  return;
3274  }
3275 
3276  Assert(this->dof_info != nullptr, ExcNotInitialized());
3277  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3278  Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3279  if (this->n_fe_components == 1)
3280  for (unsigned int comp = 0; comp < n_components; ++comp)
3281  {
3282  Assert(src[comp] != nullptr,
3283  ExcMessage("The finite element underlying this FEEvaluation "
3284  "object is scalar, but you requested " +
3285  std::to_string(n_components) +
3286  " components via the template argument in "
3287  "FEEvaluation. In that case, you must pass an "
3288  "std::vector<VectorType> or a BlockVector to " +
3289  "read_dof_values and distribute_local_to_global."));
3291  *this->matrix_free,
3292  *this->dof_info);
3293  }
3294  else
3295  {
3297  *this->matrix_free,
3298  *this->dof_info);
3299  }
3300 
3301  const bool accesses_exterior_dofs =
3302  this->dof_access_index ==
3304  this->is_interior_face() == false;
3305 
3306  // Case 2: contiguous indices which use reduced storage of indices and can
3307  // use vectorized load/store operations -> go to separate function
3308  if (this->cell != numbers::invalid_unsigned_int)
3309  {
3311  this->cell,
3312  dof_info.index_storage_variants[this->dof_access_index].size());
3313 
3314  bool is_contiguous = true;
3315  // check if exterior cells are not contiguous (ECL case)
3316  if (accesses_exterior_dofs)
3317  {
3318  const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3319  const unsigned int n_filled_lanes =
3322  [this->cell];
3323  // we have to check all filled lanes which are active in the mask
3324  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3325  if (mask[v] == true &&
3326  dof_info.index_storage_variants
3328  [cells[v] / n_lanes] <
3330  contiguous)
3331  is_contiguous = false;
3332  } // or if cell/face batch is not contiguous
3333  else if (dof_info.index_storage_variants
3334  [is_face ?
3335  this->dof_access_index :
3338  IndexStorageVariants::contiguous)
3339  {
3340  is_contiguous = false;
3341  }
3342 
3343  if (is_contiguous)
3344  {
3345  read_write_operation_contiguous(operation, src, src_sm, mask);
3346  return;
3347  }
3348  }
3349 
3350  // Case 3: standard operation with one index per degree of freedom -> go on
3351  // here
3352  std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
3353 
3354  const bool masking_is_active = mask.count() < n_lanes;
3355  if (masking_is_active)
3356  for (unsigned int v = 0; v < n_lanes; ++v)
3357  if (mask[v] == false)
3358  cells[v] = numbers::invalid_unsigned_int;
3359 
3360  bool has_hn_constraints = false;
3361 
3362  if (is_face == false)
3363  {
3364  if (!dof_info.hanging_node_constraint_masks.empty() &&
3365  !dof_info.hanging_node_constraint_masks_comp.empty() &&
3366  dof_info
3367  .hanging_node_constraint_masks_comp[this->active_fe_index]
3368  [this->first_selected_component])
3369  for (unsigned int v = 0; v < n_lanes; ++v)
3370  if (cells[v] != numbers::invalid_unsigned_int &&
3371  dof_info.hanging_node_constraint_masks[cells[v]] !=
3374  has_hn_constraints = true;
3375  }
3376 
3377  std::integral_constant<bool,
3378  internal::is_vectorizable<VectorType, Number>::value>
3379  vector_selector;
3380 
3381  const bool use_vectorized_path =
3382  !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
3383 
3384  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3385  std::array<VectorizedArrayType *, n_components> values_dofs;
3386  for (unsigned int c = 0; c < n_components; ++c)
3387  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3388  c * dofs_per_component;
3389 
3390  if (this->cell != numbers::invalid_unsigned_int &&
3391  dof_info.index_storage_variants
3392  [is_face ? this->dof_access_index :
3394  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3395  IndexStorageVariants::interleaved &&
3396  use_vectorized_path)
3397  {
3398  const unsigned int *dof_indices =
3399  dof_info.dof_indices_interleaved.data() +
3400  dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
3401  .first +
3402  this->dof_info
3403  ->component_dof_indices_offset[this->active_fe_index]
3404  [this->first_selected_component] *
3405  n_lanes;
3406 
3407  std::array<typename VectorType::value_type *, n_components> src_ptrs;
3408  if (n_components == 1 || this->n_fe_components == 1)
3409  for (unsigned int comp = 0; comp < n_components; ++comp)
3410  src_ptrs[comp] =
3411  const_cast<typename VectorType::value_type *>(src[comp]->begin());
3412  else
3413  src_ptrs[0] =
3414  const_cast<typename VectorType::value_type *>(src[0]->begin());
3415 
3416  if (n_components == 1 || this->n_fe_components == 1)
3417  for (unsigned int i = 0; i < dofs_per_component;
3418  ++i, dof_indices += n_lanes)
3419  for (unsigned int comp = 0; comp < n_components; ++comp)
3420  operation.process_dof_gather(dof_indices,
3421  *src[comp],
3422  0,
3423  src_ptrs[comp],
3424  values_dofs[comp][i],
3425  vector_selector);
3426  else
3427  for (unsigned int comp = 0; comp < n_components; ++comp)
3428  for (unsigned int i = 0; i < dofs_per_component;
3429  ++i, dof_indices += n_lanes)
3430  operation.process_dof_gather(dof_indices,
3431  *src[0],
3432  0,
3433  src_ptrs[0],
3434  values_dofs[comp][i],
3435  vector_selector);
3436  return;
3437  }
3438 
3439  // Allocate pointers, then initialize all of them to nullptrs and
3440  // below overwrite the ones we actually use:
3441  std::array<const unsigned int *, n_lanes> dof_indices;
3442  dof_indices.fill(nullptr);
3443 
3444  // Assign the appropriate cell ids for face/cell case and get the pointers
3445  // to the dof indices of the cells on all lanes
3446 
3447  bool has_constraints = false;
3448  const unsigned int n_components_read =
3449  this->n_fe_components > 1 ? n_components : 1;
3450 
3451  if (is_face)
3452  {
3453  for (unsigned int v = 0; v < n_lanes; ++v)
3454  {
3455  if (cells[v] == numbers::invalid_unsigned_int)
3456  continue;
3457 
3458  Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
3459  const std::pair<unsigned int, unsigned int> *my_index_start =
3460  &dof_info.row_starts[cells[v] * this->n_fe_components +
3461  this->first_selected_component];
3462 
3463  // check whether any of the SIMD lanes has constraints, i.e., the
3464  // constraint indicator which is the second entry of row_starts
3465  // increments on this cell
3466  if (my_index_start[n_components_read].second !=
3467  my_index_start[0].second)
3468  has_constraints = true;
3469 
3470  dof_indices[v] =
3471  dof_info.dof_indices.data() + my_index_start[0].first;
3472  }
3473  }
3474  else
3475  {
3476  for (unsigned int v = 0; v < n_lanes; ++v)
3477  {
3478  if (cells[v] == numbers::invalid_unsigned_int)
3479  continue;
3480 
3481  const std::pair<unsigned int, unsigned int> *my_index_start =
3482  &dof_info.row_starts[cells[v] * this->n_fe_components +
3483  this->first_selected_component];
3484  if (my_index_start[n_components_read].second !=
3485  my_index_start[0].second)
3486  has_constraints = true;
3487 
3488  if (dof_info.hanging_node_constraint_masks.size() > 0 &&
3489  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3490  dof_info.hanging_node_constraint_masks[cells[v]] !=
3494  [this->active_fe_index][this->first_selected_component])
3495  has_hn_constraints = true;
3496 
3497  Assert(my_index_start[n_components_read].first ==
3498  my_index_start[0].first ||
3499  my_index_start[0].first < dof_info.dof_indices.size(),
3500  ExcIndexRange(0,
3501  my_index_start[0].first,
3502  dof_info.dof_indices.size()));
3503  dof_indices[v] =
3504  dof_info.dof_indices.data() + my_index_start[0].first;
3505  }
3506  }
3507 
3508  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3509  return i != numbers::invalid_unsigned_int;
3510  }) < n_lanes)
3511  for (unsigned int comp = 0; comp < n_components; ++comp)
3512  for (unsigned int i = 0; i < dofs_per_component; ++i)
3513  operation.process_empty(values_dofs[comp][i]);
3514 
3515  // Case where we have no constraints throughout the whole cell: Can go
3516  // through the list of DoFs directly
3517  if (!has_constraints && apply_constraints)
3518  {
3519  if (n_components == 1 || this->n_fe_components == 1)
3520  {
3521  for (unsigned int v = 0; v < n_lanes; ++v)
3522  {
3523  if (cells[v] == numbers::invalid_unsigned_int)
3524  continue;
3525 
3526  for (unsigned int i = 0; i < dofs_per_component; ++i)
3527  for (unsigned int comp = 0; comp < n_components; ++comp)
3528  operation.process_dof(dof_indices[v][i],
3529  *src[comp],
3530  values_dofs[comp][i][v]);
3531  }
3532  }
3533  else
3534  {
3535  for (unsigned int comp = 0; comp < n_components; ++comp)
3536  for (unsigned int v = 0; v < n_lanes; ++v)
3537  {
3538  if (cells[v] == numbers::invalid_unsigned_int)
3539  continue;
3540 
3541  for (unsigned int i = 0; i < dofs_per_component; ++i)
3542  operation.process_dof(
3543  dof_indices[v][comp * dofs_per_component + i],
3544  *src[0],
3545  values_dofs[comp][i][v]);
3546  }
3547  }
3548  return;
3549  }
3550 
3551  // In the case where there are some constraints to be resolved, loop over
3552  // all vector components that are filled and then over local dofs. ind_local
3553  // holds local number on cell, index iterates over the elements of
3554  // index_local_to_global and dof_indices points to the global indices stored
3555  // in index_local_to_global
3556 
3557  for (unsigned int v = 0; v < n_lanes; ++v)
3558  {
3559  if (cells[v] == numbers::invalid_unsigned_int)
3560  continue;
3561 
3562  const unsigned int cell_index = cells[v];
3563  const unsigned int cell_dof_index =
3564  cell_index * this->n_fe_components + this->first_selected_component;
3565  const unsigned int n_components_read =
3566  this->n_fe_components > 1 ? n_components : 1;
3567  unsigned int index_indicators =
3568  dof_info.row_starts[cell_dof_index].second;
3569  unsigned int next_index_indicators =
3570  dof_info.row_starts[cell_dof_index + 1].second;
3571 
3572  // For read_dof_values_plain, redirect the dof_indices field to the
3573  // unconstrained indices
3574  if (apply_constraints == false &&
3575  (dof_info.row_starts[cell_dof_index].second !=
3576  dof_info.row_starts[cell_dof_index + n_components_read].second ||
3577  ((dof_info.hanging_node_constraint_masks.size() > 0 &&
3578  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3583  [this->active_fe_index][this->first_selected_component])))
3584  {
3587  ExcNotInitialized());
3588  dof_indices[v] =
3589  dof_info.plain_dof_indices.data() +
3590  this->dof_info
3591  ->component_dof_indices_offset[this->active_fe_index]
3592  [this->first_selected_component] +
3594  next_index_indicators = index_indicators;
3595  }
3596 
3597  if (n_components == 1 || this->n_fe_components == 1)
3598  {
3599  unsigned int ind_local = 0;
3600  for (; index_indicators != next_index_indicators; ++index_indicators)
3601  {
3602  const std::pair<unsigned short, unsigned short> indicator =
3603  dof_info.constraint_indicator[index_indicators];
3604  // run through values up to next constraint
3605  for (unsigned int j = 0; j < indicator.first; ++j)
3606  for (unsigned int comp = 0; comp < n_components; ++comp)
3607  operation.process_dof(dof_indices[v][j],
3608  *src[comp],
3609  values_dofs[comp][ind_local + j][v]);
3610 
3611  ind_local += indicator.first;
3612  dof_indices[v] += indicator.first;
3613 
3614  // constrained case: build the local value as a linear
3615  // combination of the global value according to constraints
3616  Number value[n_components];
3617  for (unsigned int comp = 0; comp < n_components; ++comp)
3618  operation.pre_constraints(values_dofs[comp][ind_local][v],
3619  value[comp]);
3620 
3621  const Number *data_val =
3622  this->matrix_free->constraint_pool_begin(indicator.second);
3623  const Number *end_pool =
3624  this->matrix_free->constraint_pool_end(indicator.second);
3625  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3626  for (unsigned int comp = 0; comp < n_components; ++comp)
3627  operation.process_constraint(*dof_indices[v],
3628  *data_val,
3629  *src[comp],
3630  value[comp]);
3631 
3632  for (unsigned int comp = 0; comp < n_components; ++comp)
3633  operation.post_constraints(value[comp],
3634  values_dofs[comp][ind_local][v]);
3635  ind_local++;
3636  }
3637 
3638  AssertIndexRange(ind_local, dofs_per_component + 1);
3639 
3640  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3641  for (unsigned int comp = 0; comp < n_components; ++comp)
3642  operation.process_dof(*dof_indices[v],
3643  *src[comp],
3644  values_dofs[comp][ind_local][v]);
3645  }
3646  else
3647  {
3648  // case with vector-valued finite elements where all components are
3649  // included in one single vector. Assumption: first come all entries
3650  // to the first component, then all entries to the second one, and
3651  // so on. This is ensured by the way MatrixFree reads out the
3652  // indices.
3653  for (unsigned int comp = 0; comp < n_components; ++comp)
3654  {
3655  unsigned int ind_local = 0;
3656 
3657  // check whether there is any constraint on the current cell
3658  for (; index_indicators != next_index_indicators;
3659  ++index_indicators)
3660  {
3661  const std::pair<unsigned short, unsigned short> indicator =
3662  dof_info.constraint_indicator[index_indicators];
3663 
3664  // run through values up to next constraint
3665  for (unsigned int j = 0; j < indicator.first; ++j)
3666  operation.process_dof(dof_indices[v][j],
3667  *src[0],
3668  values_dofs[comp][ind_local + j][v]);
3669  ind_local += indicator.first;
3670  dof_indices[v] += indicator.first;
3671 
3672  // constrained case: build the local value as a linear
3673  // combination of the global value according to constraints
3674  Number value;
3675  operation.pre_constraints(values_dofs[comp][ind_local][v],
3676  value);
3677 
3678  const Number *data_val =
3679  this->matrix_free->constraint_pool_begin(indicator.second);
3680  const Number *end_pool =
3681  this->matrix_free->constraint_pool_end(indicator.second);
3682 
3683  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3684  operation.process_constraint(*dof_indices[v],
3685  *data_val,
3686  *src[0],
3687  value);
3688 
3689  operation.post_constraints(value,
3690  values_dofs[comp][ind_local][v]);
3691  ind_local++;
3692  }
3693 
3694  AssertIndexRange(ind_local, dofs_per_component + 1);
3695 
3696  // get the dof values past the last constraint
3697  for (; ind_local < dofs_per_component;
3698  ++dof_indices[v], ++ind_local)
3699  {
3700  AssertIndexRange(*dof_indices[v], src[0]->size());
3701  operation.process_dof(*dof_indices[v],
3702  *src[0],
3703  values_dofs[comp][ind_local][v]);
3704  }
3705 
3706  if (apply_constraints == true && comp + 1 < n_components)
3707  next_index_indicators =
3708  dof_info.row_starts[cell_dof_index + comp + 2].second;
3709  }
3710  }
3711  }
3712 }
3713 
3714 
3715 
3716 template <int dim,
3717  int n_components_,
3718  typename Number,
3719  bool is_face,
3720  typename VectorizedArrayType>
3721 template <typename VectorType, typename VectorOperation>
3722 inline void
3725  const VectorOperation &operation,
3726  const std::array<VectorType *, n_components_> &src) const
3727 {
3728  Assert(!local_dof_indices.empty(), ExcNotInitialized());
3729 
3730  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3731  unsigned int index = this->first_selected_component * dofs_per_component;
3732  for (unsigned int comp = 0; comp < n_components; ++comp)
3733  {
3734  for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3735  {
3736  operation.process_empty(
3737  this->values_dofs[comp * dofs_per_component + i]);
3738  operation.process_dof_global(
3739  local_dof_indices[this->data->lexicographic_numbering[index]],
3740  *src[0],
3741  this->values_dofs[comp * dofs_per_component + i][0]);
3742  }
3743  }
3744 }
3745 
3746 
3747 
3748 template <int dim,
3749  int n_components_,
3750  typename Number,
3751  bool is_face,
3752  typename VectorizedArrayType>
3753 template <typename VectorType, typename VectorOperation>
3754 inline void
3757  const VectorOperation &operation,
3758  const std::array<VectorType *, n_components_> &src,
3759  const std::array<
3761  n_components_> &vectors_sm,
3762  const std::bitset<n_lanes> &mask) const
3763 {
3764  // This functions processes the functions read_dof_values,
3765  // distribute_local_to_global, and set_dof_values with the same code for
3766  // contiguous cell indices (DG case). The distinction between these three
3767  // cases is made by the input VectorOperation that either reads values from
3768  // a vector and puts the data into the local data field or write local data
3769  // into the vector. Certain operations are no-ops for the given use case.
3770 
3771  std::integral_constant<bool,
3772  internal::is_vectorizable<VectorType, Number>::value>
3773  vector_selector;
3775  is_face ? this->dof_access_index :
3777  const unsigned int n_active_lanes = mask.count();
3778 
3779  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3780  const std::vector<unsigned int> &dof_indices_cont =
3781  dof_info.dof_indices_contiguous[ind];
3782 
3783  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3784  std::array<VectorizedArrayType *, n_components> values_dofs;
3785  for (unsigned int c = 0; c < n_components; ++c)
3786  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3787  c * dofs_per_component;
3788 
3790 
3791  const bool accesses_exterior_dofs =
3792  this->dof_access_index ==
3794  this->is_interior_face() == false;
3795 
3796  // Simple case: We have contiguous storage, so we can simply copy out the
3797  // data
3798  if (dof_info.index_storage_variants[ind][this->cell] ==
3800  interleaved_contiguous &&
3801  n_active_lanes == n_lanes && !accesses_exterior_dofs)
3802  {
3803  const unsigned int dof_index =
3804  dof_indices_cont[this->cell * n_lanes] +
3805  this->dof_info
3806  ->component_dof_indices_offset[this->active_fe_index]
3807  [this->first_selected_component] *
3808  n_lanes;
3809  if (n_components == 1 || this->n_fe_components == 1)
3810  for (unsigned int comp = 0; comp < n_components; ++comp)
3811  operation.process_dofs_vectorized(dofs_per_component,
3812  dof_index,
3813  *src[comp],
3814  values_dofs[comp],
3815  vector_selector);
3816  else
3817  operation.process_dofs_vectorized(dofs_per_component * n_components,
3818  dof_index,
3819  *src[0],
3820  values_dofs[0],
3821  vector_selector);
3822  return;
3823  }
3824 
3825  const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3826 
3827  // More general case: Must go through the components one by one and apply
3828  // some transformations
3829  const unsigned int n_filled_lanes =
3830  dof_info.n_vectorization_lanes_filled[ind][this->cell];
3831 
3832  const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3833  n_active_lanes == n_lanes &&
3834  !accesses_exterior_dofs;
3835 
3836  if (vectors_sm[0] != nullptr)
3837  {
3838  const auto compute_vector_ptrs = [&](const unsigned int comp) {
3839  std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3840 
3841  const auto upper_bound =
3842  std::min<unsigned int>(n_filled_lanes, n_lanes);
3843  for (unsigned int v = 0; v < upper_bound; ++v)
3844  {
3845  if (mask[v] == false)
3846  {
3847  vector_ptrs[v] = nullptr;
3848  continue;
3849  }
3850 
3852  ExcNotImplemented());
3853  Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3854  ExcIndexRange(ind,
3855  0,
3856  dof_info.dof_indices_contiguous_sm.size()));
3857  Assert(
3858  cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3859  ExcIndexRange(cells[v],
3860  0,
3861  dof_info.dof_indices_contiguous_sm[ind].size()));
3862 
3863  const auto &temp =
3864  dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3865 
3866  if (temp.first != numbers::invalid_unsigned_int)
3867  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3868  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3870  [this->active_fe_index][this->first_selected_component]);
3871  else
3872  vector_ptrs[v] = nullptr;
3873  }
3874  for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3875  vector_ptrs[v] = nullptr;
3876 
3877  return vector_ptrs;
3878  };
3879 
3880  if (use_vectorized_path)
3881  {
3882  if (n_components == 1 || this->n_fe_components == 1)
3883  {
3884  for (unsigned int comp = 0; comp < n_components; ++comp)
3885  {
3886  auto vector_ptrs = compute_vector_ptrs(comp);
3887  operation.process_dofs_vectorized_transpose(
3888  dofs_per_component,
3889  vector_ptrs,
3890  values_dofs[comp],
3891  vector_selector);
3892  }
3893  }
3894  else
3895  {
3896  auto vector_ptrs = compute_vector_ptrs(0);
3897  operation.process_dofs_vectorized_transpose(dofs_per_component *
3898  n_components,
3899  vector_ptrs,
3900  &values_dofs[0][0],
3901  vector_selector);
3902  }
3903  }
3904  else
3905  for (unsigned int comp = 0; comp < n_components; ++comp)
3906  {
3907  auto vector_ptrs = compute_vector_ptrs(
3908  (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3909 
3910  for (unsigned int i = 0; i < dofs_per_component; ++i)
3911  operation.process_empty(values_dofs[comp][i]);
3912 
3913  if (n_components == 1 || this->n_fe_components == 1)
3914  {
3915  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3916  if (mask[v] == true)
3917  for (unsigned int i = 0; i < dofs_per_component; ++i)
3918  operation.process_dof(vector_ptrs[v][i],
3919  values_dofs[comp][i][v]);
3920  }
3921  else
3922  {
3923  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3924  if (mask[v] == true)
3925  for (unsigned int i = 0; i < dofs_per_component; ++i)
3926  operation.process_dof(
3927  vector_ptrs[v][i + comp * dofs_per_component],
3928  values_dofs[comp][i][v]);
3929  }
3930  }
3931  return;
3932  }
3933 
3934  std::array<unsigned int, n_lanes> dof_indices;
3935  std::fill(dof_indices.begin(),
3936  dof_indices.end(),
3938 
3939  Assert(n_filled_lanes <= n_lanes, ExcInternalError());
3940  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3941  {
3942  Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
3943  ExcNotImplemented());
3944  if (mask[v] == true)
3945  dof_indices[v] =
3946  dof_indices_cont[cells[v]] +
3947  this->dof_info
3948  ->component_dof_indices_offset[this->active_fe_index]
3949  [this->first_selected_component] *
3950  dof_info.dof_indices_interleave_strides[ind][cells[v]];
3951  }
3952 
3953  // In the case with contiguous cell indices, we know that there are no
3954  // constraints and that the indices within each element are contiguous
3955  if (use_vectorized_path)
3956  {
3957  if (dof_info.index_storage_variants[ind][this->cell] ==
3959  contiguous)
3960  {
3961  if (n_components == 1 || this->n_fe_components == 1)
3962  for (unsigned int comp = 0; comp < n_components; ++comp)
3963  operation.process_dofs_vectorized_transpose(dofs_per_component,
3964  dof_indices.data(),
3965  *src[comp],
3966  values_dofs[comp],
3967  vector_selector);
3968  else
3969  operation.process_dofs_vectorized_transpose(dofs_per_component *
3970  n_components,
3971  dof_indices.data(),
3972  *src[0],
3973  &values_dofs[0][0],
3974  vector_selector);
3975  }
3976  else if (dof_info.index_storage_variants[ind][this->cell] ==
3978  interleaved_contiguous_strided)
3979  {
3980  std::array<typename VectorType::value_type *, n_components> src_ptrs;
3981  if (n_components == 1 || this->n_fe_components == 1)
3982  for (unsigned int comp = 0; comp < n_components; ++comp)
3983  src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
3984  src[comp]->begin());
3985  else
3986  src_ptrs[0] =
3987  const_cast<typename VectorType::value_type *>(src[0]->begin());
3988 
3989  if (n_components == 1 || this->n_fe_components == 1)
3990  for (unsigned int i = 0; i < dofs_per_component; ++i)
3991  {
3992  for (unsigned int comp = 0; comp < n_components; ++comp)
3993  operation.process_dof_gather(dof_indices.data(),
3994  *src[comp],
3995  i * n_lanes,
3996  src_ptrs[comp] + i * n_lanes,
3997  values_dofs[comp][i],
3998  vector_selector);
3999  }
4000  else
4001  for (unsigned int comp = 0; comp < n_components; ++comp)
4002  for (unsigned int i = 0; i < dofs_per_component; ++i)
4003  {
4004  operation.process_dof_gather(
4005  dof_indices.data(),
4006  *src[0],
4007  (comp * dofs_per_component + i) * n_lanes,
4008  src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
4009  values_dofs[comp][i],
4010  vector_selector);
4011  }
4012  }
4013  else
4014  {
4015  Assert(dof_info.index_storage_variants[ind][this->cell] ==
4017  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4018  ExcNotImplemented());
4019  std::array<typename VectorType::value_type *, n_components> src_ptrs;
4020  if (n_components == 1 || this->n_fe_components == 1)
4021  for (unsigned int comp = 0; comp < n_components; ++comp)
4022  src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
4023  src[comp]->begin());
4024  else
4025  src_ptrs[0] =
4026  const_cast<typename VectorType::value_type *>(src[0]->begin());
4027 
4028  const unsigned int *offsets =
4029  &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
4030  if (n_components == 1 || this->n_fe_components == 1)
4031  for (unsigned int i = 0; i < dofs_per_component; ++i)
4032  {
4033  for (unsigned int comp = 0; comp < n_components; ++comp)
4034  operation.process_dof_gather(dof_indices.data(),
4035  *src[comp],
4036  0,
4037  src_ptrs[comp],
4038  values_dofs[comp][i],
4039  vector_selector);
4041  for (unsigned int v = 0; v < n_lanes; ++v)
4042  dof_indices[v] += offsets[v];
4043  }
4044  else
4045  for (unsigned int comp = 0; comp < n_components; ++comp)
4046  for (unsigned int i = 0; i < dofs_per_component; ++i)
4047  {
4048  operation.process_dof_gather(dof_indices.data(),
4049  *src[0],
4050  0,
4051  src_ptrs[0],
4052  values_dofs[comp][i],
4053  vector_selector);
4055  for (unsigned int v = 0; v < n_lanes; ++v)
4056  dof_indices[v] += offsets[v];
4057  }
4058  }
4059  }
4060  else
4061  for (unsigned int comp = 0; comp < n_components; ++comp)
4062  {
4063  for (unsigned int i = 0; i < dofs_per_component; ++i)
4064  operation.process_empty(values_dofs[comp][i]);
4065  if (accesses_exterior_dofs)
4066  {
4067  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4068  if (mask[v] == true)
4069  {
4070  if (dof_info.index_storage_variants
4071  [ind][cells[v] / VectorizedArrayType::size()] ==
4074  {
4075  if (n_components == 1 || this->n_fe_components == 1)
4076  {
4077  for (unsigned int i = 0; i < dofs_per_component; ++i)
4078  operation.process_dof(dof_indices[v] + i,
4079  *src[comp],
4080  values_dofs[comp][i][v]);
4081  }
4082  else
4083  {
4084  for (unsigned int i = 0; i < dofs_per_component; ++i)
4085  operation.process_dof(dof_indices[v] + i +
4086  comp * dofs_per_component,
4087  *src[0],
4088  values_dofs[comp][i][v]);
4089  }
4090  }
4091  else
4092  {
4093  const unsigned int offset =
4094  dof_info.dof_indices_interleave_strides[ind][cells[v]];
4095  AssertIndexRange(offset, VectorizedArrayType::size() + 1);
4096  if (n_components == 1 || this->n_fe_components == 1)
4097  {
4098  for (unsigned int i = 0; i < dofs_per_component; ++i)
4099  operation.process_dof(dof_indices[v] + i * offset,
4100  *src[comp],
4101  values_dofs[comp][i][v]);
4102  }
4103  else
4104  {
4105  for (unsigned int i = 0; i < dofs_per_component; ++i)
4106  operation.process_dof(
4107  dof_indices[v] +
4108  (i + comp * dofs_per_component) * offset,
4109  *src[0],
4110  values_dofs[comp][i][v]);
4111  }
4112  }
4113  }
4114  }
4115  else
4116  {
4117  if (dof_info.index_storage_variants[ind][this->cell] ==
4119  contiguous)
4120  {
4121  if (n_components == 1 || this->n_fe_components == 1)
4122  {
4123  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4124  if (mask[v] == true)
4125  for (unsigned int i = 0; i < dofs_per_component; ++i)
4126  operation.process_dof(dof_indices[v] + i,
4127  *src[comp],
4128  values_dofs[comp][i][v]);
4129  }
4130  else
4131  {
4132  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4133  if (mask[v] == true)
4134  for (unsigned int i = 0; i < dofs_per_component; ++i)
4135  operation.process_dof(dof_indices[v] + i +
4136  comp * dofs_per_component,
4137  *src[0],
4138  values_dofs[comp][i][v]);
4139  }
4140  }
4141  else
4142  {
4143  const unsigned int *offsets =
4145  [ind][VectorizedArrayType::size() * this->cell];
4146  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4147  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4148  if (n_components == 1 || this->n_fe_components == 1)
4149  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4150  {
4151  if (mask[v] == true)
4152  for (unsigned int i = 0; i < dofs_per_component; ++i)
4153  operation.process_dof(dof_indices[v] + i * offsets[v],
4154  *src[comp],
4155  values_dofs[comp][i][v]);
4156  }
4157  else
4158  {
4159  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4160  if (mask[v] == true)
4161  for (unsigned int i = 0; i < dofs_per_component; ++i)
4162  operation.process_dof(
4163  dof_indices[v] +
4164  (i + comp * dofs_per_component) * offsets[v],
4165  *src[0],
4166  values_dofs[comp][i][v]);
4167  }
4168  }
4169  }
4170  }
4171 }
4172 
4173 namespace internal
4174 {
4175  template <
4176  typename Number,
4177  typename VectorType,
4178  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4179  decltype(std::declval<VectorType>().begin())
4180  get_beginning(VectorType &vec)
4181  {
4182  return vec.begin();
4183  }
4184 
4185  template <
4186  typename Number,
4187  typename VectorType,
4188  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4189  typename VectorType::value_type *
4190  get_beginning(VectorType &)
4191  {
4192  return nullptr;
4193  }
4194 
4195  template <typename VectorType,
4196  std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4197  nullptr>
4198  const std::vector<ArrayView<const typename VectorType::value_type>> *
4199  get_shared_vector_data(VectorType *vec,
4200  const bool is_valid_mode_for_sm,
4201  const unsigned int active_fe_index,
4203  {
4204  // note: no hp is supported
4205  if (is_valid_mode_for_sm &&
4206  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4207  .size() > 0 &&
4208  active_fe_index == 0)
4209  return &vec->shared_vector_data();
4210  else
4211  return nullptr;
4212  }
4213 
4214  template <typename VectorType,
4215  std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4216  * = nullptr>
4217  const std::vector<ArrayView<const typename VectorType::value_type>> *
4218  get_shared_vector_data(VectorType *,
4219  const bool,
4220  const unsigned int,
4222  {
4223  return nullptr;
4224  }
4225 
4226  template <int n_components, typename VectorType>
4227  std::pair<
4228  std::array<typename internal::BlockVectorSelector<
4229  VectorType,
4230  IsBlockVector<VectorType>::value>::BaseVectorType *,
4231  n_components>,
4232  std::array<
4233  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4234  VectorType,
4235  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4236  n_components>>
4237  get_vector_data(VectorType &src,
4238  const unsigned int first_index,
4239  const bool is_valid_mode_for_sm,
4240  const unsigned int active_fe_index,
4242  {
4243  // select between block vectors and non-block vectors. Note that the number
4244  // of components is checked in the internal data
4245  std::pair<
4246  std::array<typename internal::BlockVectorSelector<
4247  VectorType,
4248  IsBlockVector<VectorType>::value>::BaseVectorType *,
4249  n_components>,
4250  std::array<
4251  const std::vector<
4252  ArrayView<const typename internal::BlockVectorSelector<
4253  VectorType,
4254  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4255  n_components>>
4256  src_data;
4257 
4258  for (unsigned int d = 0; d < n_components; ++d)
4259  src_data.first[d] = internal::BlockVectorSelector<
4260  VectorType,
4261  IsBlockVector<VectorType>::value>::get_vector_component(src,
4262  d +
4263  first_index);
4264 
4265  for (unsigned int d = 0; d < n_components; ++d)
4266  src_data.second[d] = get_shared_vector_data(
4267  const_cast<typename internal::BlockVectorSelector<
4268  std::remove_const_t<VectorType>,
4269  IsBlockVector<std::remove_const_t<VectorType>>::value>::BaseVectorType
4270  *>(src_data.first[d]),
4271  is_valid_mode_for_sm,
4272  active_fe_index,
4273  dof_info);
4274 
4275  return src_data;
4276  }
4277 } // namespace internal
4278 
4279 
4280 
4281 template <int dim,
4282  int n_components_,
4283  typename Number,
4284  bool is_face,
4285  typename VectorizedArrayType>
4286 inline void
4289 {
4290  if (this->dof_info == nullptr ||
4291  this->dof_info->hanging_node_constraint_masks.empty() ||
4292  this->dof_info->hanging_node_constraint_masks_comp.empty() ||
4293  this->dof_info->hanging_node_constraint_masks_comp
4294  [this->active_fe_index][this->first_selected_component] == false)
4295  return; // nothing to do with faces
4296 
4297  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4298  constraint_mask;
4299 
4300  bool hn_available = false;
4301 
4302  const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
4303 
4304  for (unsigned int v = 0; v < n_lanes; ++v)
4305  {
4306  if (cells[v] == numbers::invalid_unsigned_int)
4307  {
4308  constraint_mask[v] = internal::MatrixFreeFunctions::
4310  continue;
4311  }
4312 
4313  const unsigned int cell_index = cells[v];
4314  const auto mask =
4315  this->dof_info->hanging_node_constraint_masks[cell_index];
4316  constraint_mask[v] = mask;
4317 
4318  hn_available |= (mask != internal::MatrixFreeFunctions::
4320  }
4321 
4322  if (hn_available == false)
4323  return; // no hanging node on cell batch -> nothing to do
4324 
4326  apply(n_components,
4327  this->data->data.front().fe_degree,
4328  this->get_shape_info(),
4329  transpose,
4330  constraint_mask,
4331  this->values_dofs);
4332 }
4333 
4334 
4335 
4336 template <int dim,
4337  int n_components_,
4338  typename Number,
4339  bool is_face,
4340  typename VectorizedArrayType>
4341 template <typename VectorType>
4342 inline void
4344  read_dof_values(const VectorType &src,
4345  const unsigned int first_index,
4346  const std::bitset<n_lanes> &mask)
4347 {
4348  const auto src_data = internal::get_vector_data<n_components_>(
4349  src,
4350  first_index,
4351  this->dof_access_index ==
4353  this->active_fe_index,
4354  this->dof_info);
4355 
4357  read_write_operation(reader, src_data.first, src_data.second, mask, true);
4358 
4359  apply_hanging_node_constraints(false);
4360 
4361 # ifdef DEBUG
4362  this->dof_values_initialized = true;
4363 # endif
4364 }
4365 
4366 
4367 
4368 template <int dim,
4369  int n_components_,
4370  typename Number,
4371  bool is_face,
4372  typename VectorizedArrayType>
4373 template <typename VectorType>
4374 inline void
4376  read_dof_values_plain(const VectorType &src,
4377  const unsigned int first_index,
4378  const std::bitset<n_lanes> &mask)
4379 {
4380  const auto src_data = internal::get_vector_data<n_components_>(
4381  src,
4382  first_index,
4383  this->dof_access_index ==
4385  this->active_fe_index,
4386  this->dof_info);
4387 
4389  read_write_operation(reader, src_data.first, src_data.second, mask, false);
4390 
4391 # ifdef DEBUG
4392  this->dof_values_initialized = true;
4393 # endif
4394 }
4395 
4396 
4397 
4398 template <int dim,
4399  int n_components_,
4400  typename Number,
4401  bool is_face,
4402  typename VectorizedArrayType>
4403 template <typename VectorType>
4404 inline void
4406  distribute_local_to_global(VectorType &dst,
4407  const unsigned int first_index,
4408  const std::bitset<n_lanes> &mask) const
4409 {
4410 # ifdef DEBUG
4411  Assert(this->dof_values_initialized == true,
4413 # endif
4414 
4415  apply_hanging_node_constraints(true);
4416 
4417  const auto dst_data = internal::get_vector_data<n_components_>(
4418  dst,
4419  first_index,
4420  this->dof_access_index ==
4422  this->active_fe_index,
4423  this->dof_info);
4424 
4426  distributor;
4427  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4428 }
4429 
4430 
4431 
4432 template <int dim,
4433  int n_components_,
4434  typename Number,
4435  bool is_face,
4436  typename VectorizedArrayType>
4437 template <typename VectorType>
4438 inline void
4440  set_dof_values(VectorType &dst,
4441  const unsigned int first_index,
4442  const std::bitset<n_lanes> &mask) const
4443 {
4444 # ifdef DEBUG
4445  Assert(this->dof_values_initialized == true,
4447 # endif
4448 
4449  const auto dst_data = internal::get_vector_data<n_components_>(
4450  dst,
4451  first_index,
4452  this->dof_access_index ==
4454  this->active_fe_index,
4455  this->dof_info);
4456 
4458  read_write_operation(setter, dst_data.first, dst_data.second, mask);
4459 }
4460 
4461 
4462 
4463 template <int dim,
4464  int n_components_,
4465  typename Number,
4466  bool is_face,
4467  typename VectorizedArrayType>
4468 template <typename VectorType>
4469 inline void
4471  set_dof_values_plain(VectorType &dst,
4472  const unsigned int first_index,
4473  const std::bitset<n_lanes> &mask) const
4474 {
4475 # ifdef DEBUG
4476  Assert(this->dof_values_initialized == true,
4478 # endif
4479 
4480  const auto dst_data = internal::get_vector_data<n_components_>(
4481  dst,
4482  first_index,
4483  this->dof_access_index ==
4485  this->active_fe_index,
4486  this->dof_info);
4487 
4489  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4490 }
4491 
4492 
4493 
4494 /*------------------------------ access to data fields ----------------------*/
4495 
4496 
4497 
4498 template <int dim,
4499  int n_components_,
4500  typename Number,
4501  bool is_face,
4502  typename VectorizedArrayType>
4505  get_dof_value(const unsigned int dof) const
4506 {
4507  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4508  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4510  for (unsigned int comp = 0; comp < n_components; ++comp)
4511  return_value[comp] = this->values_dofs[comp * dofs + dof];
4512  return return_value;
4513 }
4514 
4515 
4516 
4517 template <int dim,
4518  int n_components_,
4519  typename Number,
4520  bool is_face,
4521  typename VectorizedArrayType>
4524  get_value(const unsigned int q_point) const
4525 {
4526 # ifdef DEBUG
4527  Assert(this->values_quad_initialized == true,
4529 # endif
4530 
4531  AssertIndexRange(q_point, this->n_quadrature_points);
4532  const std::size_t nqp = this->n_quadrature_points;
4534  for (unsigned int comp = 0; comp < n_components; ++comp)
4535  return_value[comp] = this->values_quad[comp * nqp + q_point];
4536  return return_value;
4537 }
4538 
4539 
4540 
4541 template <int dim,
4542  int n_components_,
4543  typename Number,
4544  bool is_face,
4545  typename VectorizedArrayType>
4546 inline DEAL_II_ALWAYS_INLINE
4549  get_gradient(const unsigned int q_point) const
4550 {
4551 # ifdef DEBUG
4552  Assert(this->gradients_quad_initialized == true,
4554 # endif
4555 
4556  AssertIndexRange(q_point, this->n_quadrature_points);
4557  Assert(this->jacobian != nullptr,
4559  "update_gradients"));
4560  const std::size_t nqp = this->n_quadrature_points;
4562 
4563  // Cartesian cell
4564  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4565  {
4566  for (unsigned int comp = 0; comp < n_components; ++comp)
4567  for (unsigned int d = 0; d < dim; ++d)
4568  grad_out[comp][d] =
4569  this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4570  this->jacobian[0][d][d];
4571  }
4572  // cell with general/affine Jacobian
4573  else
4574  {
4576  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4577  q_point :
4578  0];
4579  for (unsigned int comp = 0; comp < n_components; ++comp)
4580  for (unsigned int d = 0; d < dim; ++d)
4581  {
4582  grad_out[comp][d] =
4583  jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4584  for (unsigned int e = 1; e < dim; ++e)
4585  grad_out[comp][d] +=
4586  jac[d][e] *
4587  this->gradients_quad[(comp * nqp + q_point) * dim + e];
4588  }
4589  }
4590  return grad_out;
4591 }
4592 
4593 
4594 
4595 template <int dim,
4596  int n_components_,
4597  typename Number,
4598  bool is_face,
4599  typename VectorizedArrayType>
4602  get_normal_derivative(const unsigned int q_point) const
4603 {
4604  AssertIndexRange(q_point, this->n_quadrature_points);
4605 # ifdef DEBUG
4606  Assert(this->gradients_quad_initialized == true,
4608 # endif
4609 
4610  Assert(this->normal_x_jacobian != nullptr,
4612  "update_gradients"));
4613 
4614  const std::size_t nqp = this->n_quadrature_points;
4616 
4617  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4618  for (unsigned int comp = 0; comp < n_components; ++comp)
4619  grad_out[comp] =
4620  this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4621  (this->normal_x_jacobian[0][dim - 1]);
4622  else
4623  {
4624  const std::size_t index =
4625  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4626  for (unsigned int comp = 0; comp < n_components; ++comp)
4627  {
4628  grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4629  this->normal_x_jacobian[index][0];
4630  for (unsigned int d = 1; d < dim; ++d)
4631  grad_out[comp] +=
4632  this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4633  this->normal_x_jacobian[index][d];
4634  }
4635  }
4636  return grad_out;
4637 }
4638 
4639 
4640 
4641 namespace internal
4642 {
4643  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4644  // store the lower diagonal because of symmetry
4645  template <typename VectorizedArrayType>
4646  inline void
4647  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4648  const VectorizedArrayType *const hessians,
4649  const unsigned int,
4650  VectorizedArrayType (&tmp)[1][1])
4651  {
4652  tmp[0][0] = jac[0][0] * hessians[0];
4653  }
4654 
4655  template <typename VectorizedArrayType>
4656  inline void
4657  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4658  const VectorizedArrayType *const hessians,
4659  const unsigned int nqp,
4660  VectorizedArrayType (&tmp)[2][2])
4661  {
4662  for (unsigned int d = 0; d < 2; ++d)
4663  {
4664  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4665  tmp[1][d] =
4666  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4667  }
4668  }
4669 
4670  template <typename VectorizedArrayType>
4671  inline void
4672  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4673  const VectorizedArrayType *const hessians,
4674  const unsigned int nqp,
4675  VectorizedArrayType (&tmp)[3][3])
4676  {
4677  for (unsigned int d = 0; d < 3; ++d)
4678  {
4679  tmp[0][d] =
4680  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4681  jac[d][2] * hessians[4 * nqp]);
4682  tmp[1][d] =
4683  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4684  jac[d][2] * hessians[5 * nqp]);
4685  tmp[2][d] =
4686  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4687  jac[d][2] * hessians[2 * nqp]);
4688  }
4689  }
4690 } // namespace internal
4691 
4692 
4693 
4694 template <int dim,
4695  int n_components_,
4696  typename Number,
4697  bool is_face,
4698  typename VectorizedArrayType>
4701  get_hessian(const unsigned int q_point) const
4702 {
4703 # ifdef DEBUG
4704  Assert(this->hessians_quad_initialized == true,
4706 # endif
4707  AssertIndexRange(q_point, this->n_quadrature_points);
4708 
4709  Assert(this->jacobian != nullptr,
4711  "update_hessian"));
4713  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4714  0 :
4715  q_point];
4716 
4718 
4719  const std::size_t nqp = this->n_quadrature_points;
4720  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4721 
4722  // Cartesian cell
4723  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4724  {
4725  for (unsigned int comp = 0; comp < n_components; ++comp)
4726  {
4727  for (unsigned int d = 0; d < dim; ++d)
4728  hessian_out[comp][d][d] =
4729  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4730  (jac[d][d] * jac[d][d]);
4731  switch (dim)
4732  {
4733  case 1:
4734  break;
4735  case 2:
4736  hessian_out[comp][0][1] =
4737  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4738  (jac[0][0] * jac[1][1]);
4739  break;
4740  case 3:
4741  hessian_out[comp][0][1] =
4742  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4743  (jac[0][0] * jac[1][1]);
4744  hessian_out[comp][0][2] =
4745  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4746  (jac[0][0] * jac[2][2]);
4747  hessian_out[comp][1][2] =
4748  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4749  (jac[1][1] * jac[2][2]);
4750  break;
4751  default:
4752  Assert(false, ExcNotImplemented());
4753  }
4754  for (unsigned int d = 0; d < dim; ++d)
4755  for (unsigned int e = d + 1; e < dim; ++e)
4756  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4757  }
4758  }
4759  // cell with general Jacobian, but constant within the cell
4760  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4761  {
4762  for (unsigned int comp = 0; comp < n_components; ++comp)
4763  {
4764  VectorizedArrayType tmp[dim][dim];
4765  internal::hessian_unit_times_jac(
4766  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4767 
4768  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4769  for (unsigned int d = 0; d < dim; ++d)
4770  for (unsigned int e = d; e < dim; ++e)
4771  {
4772  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4773  for (unsigned int f = 1; f < dim; ++f)
4774  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4775  }
4776 
4777  // no J' * grad(u) part here because the Jacobian is constant
4778  // throughout the cell and hence, its derivative is zero
4779 
4780  // take symmetric part
4781  for (unsigned int d = 0; d < dim; ++d)
4782  for (unsigned int e = d + 1; e < dim; ++e)
4783  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4784  }
4785  }
4786  // cell with general Jacobian
4787  else
4788  {
4789  const auto &jac_grad = this->jacobian_gradients[q_point];
4790  for (unsigned int comp = 0; comp < n_components; ++comp)
4791  {
4792  VectorizedArrayType tmp[dim][dim];
4793  internal::hessian_unit_times_jac(
4794  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4795 
4796  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4797  for (unsigned int d = 0; d < dim; ++d)
4798  for (unsigned int e = d; e < dim; ++e)
4799  {
4800  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4801  for (unsigned int f = 1; f < dim; ++f)
4802  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4803  }
4804 
4805  // add diagonal part of J' * grad(u)
4806  for (unsigned int d = 0; d < dim; ++d)
4807  for (unsigned int e = 0; e < dim; ++e)
4808  hessian_out[comp][d][d] +=
4809  jac_grad[d][e] *
4810  this->gradients_quad[(comp * nqp + q_point) * dim + e];
4811 
4812  // add off-diagonal part of J' * grad(u)
4813  for (unsigned int d = 0, count = dim; d < dim; ++d)
4814  for (unsigned int e = d + 1; e < dim; ++e, ++count)
4815  for (unsigned int f = 0; f < dim; ++f)
4816  hessian_out[comp][d][e] +=
4817  jac_grad[count][f] *
4818  this->gradients_quad[(comp * nqp + q_point) * dim + f];
4819 
4820  // take symmetric part
4821  for (unsigned int d = 0; d < dim; ++d)
4822  for (unsigned int e = d + 1; e < dim; ++e)
4823  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4824  }
4825  }
4826  return hessian_out;
4827 }
4828 
4829 
4830 
4831 template <int dim,
4832  int n_components_,
4833  typename Number,
4834  bool is_face,
4835  typename VectorizedArrayType>
4838  get_hessian_diagonal(const unsigned int q_point) const
4839 {
4840  Assert(!is_face, ExcNotImplemented());
4841 # ifdef DEBUG
4842  Assert(this->hessians_quad_initialized == true,
4844 # endif
4845  AssertIndexRange(q_point, this->n_quadrature_points);
4846 
4847  Assert(this->jacobian != nullptr, ExcNotImplemented());
4849  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4850  0 :
4851  q_point];
4852 
4853  const std::size_t nqp = this->n_quadrature_points;
4854  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4856 
4857  // Cartesian cell
4858  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4859  {
4860  for (unsigned int comp = 0; comp < n_components; ++comp)
4861  for (unsigned int d = 0; d < dim; ++d)
4862  hessian_out[comp][d] =
4863  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4864  (jac[d][d] * jac[d][d]);
4865  }
4866  // cell with general Jacobian, but constant within the cell
4867  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4868  {
4869  for (unsigned int comp = 0; comp < n_components; ++comp)
4870  {
4871  // compute laplacian before the gradient because it needs to access
4872  // unscaled gradient data
4873  VectorizedArrayType tmp[dim][dim];
4874  internal::hessian_unit_times_jac(
4875  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4876 
4877  // compute only the trace part of hessian, J * tmp = J *
4878  // hess_unit(u) * J^T
4879  for (unsigned int d = 0; d < dim; ++d)
4880  {
4881  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4882  for (unsigned int f = 1; f < dim; ++f)
4883  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4884  }
4885  }
4886  }
4887  // cell with general Jacobian
4888  else
4889  {
4890  const auto &jac_grad = this->jacobian_gradients[q_point];
4891  for (unsigned int comp = 0; comp < n_components; ++comp)
4892  {
4893  // compute laplacian before the gradient because it needs to access
4894  // unscaled gradient data
4895  VectorizedArrayType tmp[dim][dim];
4896  internal::hessian_unit_times_jac(
4897  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4898 
4899  // compute only the trace part of hessian, J * tmp = J *
4900  // hess_unit(u) * J^T
4901  for (unsigned int d = 0; d < dim; ++d)
4902  {
4903  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4904  for (unsigned int f = 1; f < dim; ++f)
4905  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4906  }
4907 
4908  for (unsigned int d = 0; d < dim; ++d)
4909  for (unsigned int e = 0; e < dim; ++e)
4910  hessian_out[comp][d] +=
4911  jac_grad[d][e] *
4912  this->gradients_quad[(comp * nqp + q_point) * dim + e];
4913  }
4914  }
4915  return hessian_out;
4916 }
4917 
4918 
4919 
4920 template <int dim,
4921  int n_components_,
4922  typename Number,
4923  bool is_face,
4924  typename VectorizedArrayType>
4927  get_laplacian(const unsigned int q_point) const
4928 {
4929  Assert(is_face == false, ExcNotImplemented());
4930 # ifdef DEBUG
4931  Assert(this->hessians_quad_initialized == true,
4933 # endif
4934  AssertIndexRange(q_point, this->n_quadrature_points);
4935 
4937  const auto hess_diag = get_hessian_diagonal(q_point);
4938  for (unsigned int comp = 0; comp < n_components; ++comp)
4939  {
4940  laplacian_out[comp] = hess_diag[comp][0];
4941  for (unsigned int d = 1; d < dim; ++d)
4942  laplacian_out[comp] += hess_diag[comp][d];
4943  }
4944  return laplacian_out;
4945 }
4946 
4947 
4948 
4949 template <int dim,
4950  int n_components_,
4951  typename Number,
4952  bool is_face,
4953  typename VectorizedArrayType>
4954 inline DEAL_II_ALWAYS_INLINE void
4957  const unsigned int dof)
4958 {
4959 # ifdef DEBUG
4960  this->dof_values_initialized = true;
4961 # endif
4962  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4963  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4964  for (unsigned int comp = 0; comp < n_components; ++comp)
4965  this->values_dofs[comp * dofs + dof] = val_in[comp];
4966 }
4967 
4968 
4969 
4970 template <int dim,
4971  int n_components_,
4972  typename Number,
4973  bool is_face,
4974  typename VectorizedArrayType>
4975 inline DEAL_II_ALWAYS_INLINE void
4978  const unsigned int q_point)
4979 {
4980 # ifdef DEBUG
4981  Assert(this->is_reinitialized, ExcNotInitialized());
4982 # endif
4983  AssertIndexRange(q_point, this->n_quadrature_points);
4984  Assert(this->J_value != nullptr,
4986  "update_values"));
4987 # ifdef DEBUG
4988  this->values_quad_submitted = true;
4989 # endif
4990 
4991  const std::size_t nqp = this->n_quadrature_points;
4992  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4993  {
4994  const VectorizedArrayType JxW =
4995  this->J_value[0] * this->quadrature_weights[q_point];
4996  for (unsigned int comp = 0; comp < n_components; ++comp)
4997  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4998  }
4999  else
5000  {
5001  const VectorizedArrayType JxW = this->J_value[q_point];
5002  for (unsigned int comp = 0; comp < n_components; ++comp)
5003  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5004  }
5005 }
5006 
5007 
5008 
5009 template <int dim,
5010  int n_components_,
5011  typename Number,
5012  bool is_face,
5013  typename VectorizedArrayType>
5014 inline DEAL_II_ALWAYS_INLINE void
5017  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5018  const unsigned int q_point)
5019 {
5020 # ifdef DEBUG
5021  Assert(this->is_reinitialized, ExcNotInitialized());
5022 # endif
5023  AssertIndexRange(q_point, this->n_quadrature_points);
5024  Assert(this->J_value != nullptr,
5026  "update_gradients"));
5027  Assert(this->jacobian != nullptr,
5029  "update_gradients"));
5030 # ifdef DEBUG
5031  this->gradients_quad_submitted = true;
5032 # endif
5033 
5034  const std::size_t nqp_d = this->n_quadrature_points * dim;
5035  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5036 
5037  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5038  {
5039  const VectorizedArrayType JxW =
5040  this->J_value[0] * this->quadrature_weights[q_point];
5041  std::array<VectorizedArrayType, dim> jac;
5042  for (unsigned int d = 0; d < dim; ++d)
5043  jac[d] = this->jacobian[0][d][d];
5044  for (unsigned int d = 0; d < dim; ++d)
5045  {
5046  const VectorizedArrayType factor = jac[d] * JxW;
5047  for (unsigned int comp = 0; comp < n_components; ++comp)
5048  gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5049  }
5050  }
5051  else
5052  {
5054  this->cell_type > internal::MatrixFreeFunctions::affine ?
5055  this->jacobian[q_point] :
5056  this->jacobian[0];
5057  const VectorizedArrayType JxW =
5058  this->cell_type > internal::MatrixFreeFunctions::affine ?
5059  this->J_value[q_point] :
5060  this->J_value[0] * this->quadrature_weights[q_point];
5061  for (unsigned int comp = 0; comp < n_components; ++comp)
5062  for (unsigned int d = 0; d < dim; ++d)
5063  {
5064  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5065  for (unsigned int e = 1; e < dim; ++e)
5066  new_val += (jac[e][d] * grad_in[comp][e]);
5067  gradients[comp * nqp_d + d] = new_val * JxW;
5068  }
5069  }
5070 }
5071 
5072 
5073 
5074 template <int dim,
5075  int n_components_,
5076  typename Number,
5077  bool is_face,
5078  typename VectorizedArrayType>
5079 inline DEAL_II_ALWAYS_INLINE void
5083  const unsigned int q_point)
5084 {
5085  AssertIndexRange(q_point, this->n_quadrature_points);
5086  Assert(this->normal_x_jacobian != nullptr,
5088  "update_gradients"));
5089 # ifdef DEBUG
5090  this->gradients_quad_submitted = true;
5091 # endif
5092 
5093  const std::size_t nqp_d = this->n_quadrature_points * dim;
5094  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5095 
5096  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5097  {
5098  const VectorizedArrayType JxW_jac = this->J_value[0] *
5099  this->quadrature_weights[q_point] *
5100  this->normal_x_jacobian[0][dim - 1];
5101  for (unsigned int comp = 0; comp < n_components; ++comp)
5102  {
5103  for (unsigned int d = 0; d < dim - 1; ++d)
5104  gradients[comp * nqp_d + d] = VectorizedArrayType();
5105  gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5106  }
5107  }
5108  else
5109  {
5110  const unsigned int index =
5111  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5113  this->normal_x_jacobian[index];
5114  const VectorizedArrayType JxW =
5115  (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5116  this->J_value[index] * this->quadrature_weights[q_point] :
5117  this->J_value[index];
5118  for (unsigned int comp = 0; comp < n_components; ++comp)
5119  {
5120  for (unsigned int d = 0; d < dim; ++d)
5121  gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
5122  }
5123  }
5124 }
5125 
5126 
5127 
5128 template <int dim,
5129  int n_components_,
5130  typename Number,
5131  bool is_face,
5132  typename VectorizedArrayType>
5133 inline DEAL_II_ALWAYS_INLINE void
5136  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5137  hessian_in,
5138  const unsigned int q_point)
5139 {
5140 # ifdef DEBUG
5141  Assert(this->is_reinitialized, ExcNotInitialized());
5142 # endif
5143  AssertIndexRange(q_point, this->n_quadrature_points);
5144  Assert(this->J_value != nullptr,
5146  "update_hessians"));
5147  Assert(this->jacobian != nullptr,
5149  "update_hessians"));
5150 # ifdef DEBUG
5151  this->hessians_quad_submitted = true;
5152 # endif
5153 
5154  // compute hessian_unit = J^T * hessian_in(u) * J
5155  const std::size_t nqp = this->n_quadrature_points;
5156  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5157  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5158  {
5159  const VectorizedArrayType JxW =
5160  this->J_value[0] * this->quadrature_weights[q_point];
5161 
5162  // diagonal part
5163  for (unsigned int d = 0; d < dim; ++d)
5164  {
5165  const auto jac_d = this->jacobian[0][d][d];
5166  const VectorizedArrayType factor = jac_d * jac_d * JxW;
5167  for (unsigned int comp = 0; comp < n_components; ++comp)
5168  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5169  hessian_in[comp][d][d] * factor;
5170  }
5171 
5172  // off diagonal part
5173  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5174  for (unsigned int e = 0; e < d; ++e, ++off_dia)
5175  {
5176  const auto jac_d = this->jacobian[0][d][d];
5177  const auto jac_e = this->jacobian[0][e][e];
5178  const VectorizedArrayType factor = jac_d * jac_e * JxW;
5179  for (unsigned int comp = 0; comp < n_components; ++comp)
5180  this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5181  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5182  }
5183  }
5184  // cell with general Jacobian, but constant within the cell
5185  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5186  {
5187  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5188  const VectorizedArrayType JxW =
5189  this->J_value[0] * this->quadrature_weights[q_point];
5190  for (unsigned int comp = 0; comp < n_components; ++comp)
5191  {
5192  // 1. tmp = hessian_in(u) * J
5193  VectorizedArrayType tmp[dim][dim];
5194  for (unsigned int i = 0; i < dim; ++i)
5195  for (unsigned int j = 0; j < dim; ++j)
5196  {
5197  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5198  for (unsigned int k = 1; k < dim; ++k)
5199  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5200  }
5201 
5202  // 2. hessian_unit = J^T * tmp
5203  VectorizedArrayType tmp2[dim][dim];
5204  for (unsigned int i = 0; i < dim; ++i)
5205  for (unsigned int j = 0; j < dim; ++j)
5206  {
5207  tmp2[i][j] = jac[0][i] * tmp[0][j];
5208  for (unsigned int k = 1; k < dim; ++k)
5209  tmp2[i][j] += jac[k][i] * tmp[k][j];
5210  }
5211 
5212  // diagonal part
5213  for (unsigned int d = 0; d < dim; ++d)
5214  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5215  tmp2[d][d] * JxW;
5216 
5217  // off diagonal part
5218  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5219  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5220  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5221  (tmp2[d][e] + tmp2[e][d]) * JxW;
5222  }
5223  }
5224  else
5225  {
5226  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5227  const VectorizedArrayType JxW = this->J_value[q_point];
5228  const auto &jac_grad = this->jacobian_gradients[q_point];
5229  for (unsigned int comp = 0; comp < n_components; ++comp)
5230  {
5231  // 1. tmp = hessian_in(u) * J
5232  VectorizedArrayType tmp[dim][dim];
5233  for (unsigned int i = 0; i < dim; ++i)
5234  for (unsigned int j = 0; j < dim; ++j)
5235  {
5236  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5237  for (unsigned int k = 1; k < dim; ++k)
5238  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5239  }
5240 
5241  // 2. hessian_unit = J^T * tmp
5242  VectorizedArrayType tmp2[dim][dim];
5243  for (unsigned int i = 0; i < dim; ++i)
5244  for (unsigned int j = 0; j < dim; ++j)
5245  {
5246  tmp2[i][j] = jac[0][i] * tmp[0][j];
5247  for (unsigned int k = 1; k < dim; ++k)
5248  tmp2[i][j] += jac[k][i] * tmp[k][j];
5249  }
5250 
5251  // diagonal part
5252  for (unsigned int d = 0; d < dim; ++d)
5253  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5254  tmp2[d][d] * JxW;
5255 
5256  // off diagonal part
5257  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5258  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5259  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5260  (tmp2[d][e] + tmp2[e][d]) * JxW;
5261 
5262  // 3. gradient_unit = J' ** hessian_in
5263  for (unsigned int d = 0; d < dim; ++d)
5264  {
5265  VectorizedArrayType sum = 0;
5266  for (unsigned int e = 0; e < dim; ++e)
5267  sum += hessian_in[comp][e][e] * jac_grad[e][d];
5268  for (unsigned int e = 0, count = dim; e < dim; ++e)
5269  for (unsigned int f = e + 1; f < dim; ++f, ++count)
5270  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5271  jac_grad[count][d];
5272  this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5273  d] = sum * JxW;
5274  }
5275  }
5276  }
5277 }
5278 
5279 
5280 
5281 template <int dim,
5282  int n_components_,
5283  typename Number,
5284  bool is_face,
5285  typename VectorizedArrayType>
5288  integrate_value() const
5289 {
5290 # ifdef DEBUG
5291  Assert(this->is_reinitialized, ExcNotInitialized());
5292  Assert(this->values_quad_submitted == true,
5294 # endif
5295 
5297  const std::size_t nqp = this->n_quadrature_points;
5298  for (unsigned int q = 0; q < nqp; ++q)
5299  for (unsigned int comp = 0; comp < n_components; ++comp)
5300  return_value[comp] += this->values_quad[comp * nqp + q];
5301  return (return_value);
5302 }
5303 
5304 
5305 
5306 /*----------------------- FEEvaluationAccess --------------------------------*/
5307 
5308 
5309 template <int dim,
5310  int n_components_,
5311  typename Number,
5312  bool is_face,
5313  typename VectorizedArrayType>
5314 inline FEEvaluationAccess<dim,
5315  n_components_,
5316  Number,
5317  is_face,
5318  VectorizedArrayType>::
5319  FEEvaluationAccess(
5321  const unsigned int dof_no,
5322  const unsigned int first_selected_component,
5323  const unsigned int quad_no,
5324  const unsigned int fe_degree,
5325  const unsigned int n_q_points,
5326  const bool is_interior_face,
5327  const unsigned int active_fe_index,
5328  const unsigned int active_quad_index,
5329  const unsigned int face_type)
5330  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5331  matrix_free,
5332  dof_no,
5333  first_selected_component,
5334  quad_no,
5335  fe_degree,
5336  n_q_points,
5337  is_interior_face,
5338  active_fe_index,
5339  active_quad_index,
5340  face_type)
5341 {}
5342 
5343 
5344 
5345 template <int dim,
5346  int n_components_,
5347  typename Number,
5348  bool is_face,
5349  typename VectorizedArrayType>
5350 inline FEEvaluationAccess<dim,
5351  n_components_,
5352  Number,
5353  is_face,
5354  VectorizedArrayType>::
5355  FEEvaluationAccess(
5356  const Mapping<dim> &mapping,
5357  const FiniteElement<dim> &fe,
5358  const Quadrature<1> &quadrature,
5359  const UpdateFlags update_flags,
5360  const unsigned int first_selected_component,
5362  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5363  mapping,
5364  fe,
5365  quadrature,
5366  update_flags,
5367  first_selected_component,
5368  other)
5369 {}
5370 
5371 
5372 
5373 template <int dim,
5374  int n_components_,
5375  typename Number,
5376  bool is_face,
5377  typename VectorizedArrayType>
5378 inline FEEvaluationAccess<dim,
5379  n_components_,
5380  Number,
5381  is_face,
5382  VectorizedArrayType>::
5383  FEEvaluationAccess(const FEEvaluationAccess<dim,
5384  n_components_,
5385  Number,
5386  is_face,
5387  VectorizedArrayType> &other)
5388  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5389  other)
5390 {}
5391 
5392 
5393 
5394 template <int dim,
5395  int n_components_,
5396  typename Number,
5397  bool is_face,
5398  typename VectorizedArrayType>
5399 inline FEEvaluationAccess<dim,
5400  n_components_,
5401  Number,
5402  is_face,
5403  VectorizedArrayType> &
5405 operator=(const FEEvaluationAccess<dim,
5406  n_components_,
5407  Number,
5408  is_face,
5409  VectorizedArrayType> &other)
5410 {
5411  this->FEEvaluationBase<dim,
5412  n_components_,
5413  Number,
5414  is_face,
5415  VectorizedArrayType>::operator=(other);
5416  return *this;
5417 }
5418 
5419 
5420 
5421 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5422 
5423 
5424 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5428  const unsigned int dof_no,
5429  const unsigned int first_selected_component,
5430  const unsigned int quad_no,
5431  const unsigned int fe_degree,
5432  const unsigned int n_q_points,
5433  const bool is_interior_face,
5434  const unsigned int active_fe_index,
5435  const unsigned int active_quad_index,
5436  const unsigned int face_type)
5437  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5438  matrix_free,
5439  dof_no,
5440  first_selected_component,
5441  quad_no,
5442  fe_degree,
5443  n_q_points,
5444  is_interior_face,
5445  active_fe_index,
5446  active_quad_index,
5447  face_type)
5448 {}
5449 
5450 
5451 
5452 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5455  const Mapping<dim> &mapping,
5456  const FiniteElement<dim> &fe,
5457  const Quadrature<1> &quadrature,
5458  const UpdateFlags update_flags,
5459  const unsigned int first_selected_component,
5461  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5462  mapping,
5463  fe,
5464  quadrature,
5465  update_flags,
5466  first_selected_component,
5467  other)
5468 {}
5469 
5470 
5471 
5472 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5476  &other)
5477  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5478 {}
5479 
5480 
5481 
5482 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5486 {
5487  this
5488  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5489  other);
5490  return *this;
5491 }
5492 
5493 
5494 
5495 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5496 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5498  const unsigned int dof) const
5499 {
5500  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5501  return this->values_dofs[dof];
5502 }
5503 
5504 
5505 
5506 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5507 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5509  const unsigned int q_point) const
5510 {
5511 # ifdef DEBUG
5512  Assert(this->values_quad_initialized == true,
5514 # endif
5515  AssertIndexRange(q_point, this->n_quadrature_points);
5516  return this->values_quad[q_point];
5517 }
5518 
5519 
5520 
5521 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5522 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5524  get_normal_derivative(const unsigned int q_point) const
5525 {
5526  return BaseClass::get_normal_derivative(q_point)[0];
5527 }
5528 
5529 
5530 
5531 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5534  const unsigned int q_point_in) const
5535 {
5536  // could use the base class gradient, but that involves too many expensive
5537  // initialization operations on tensors
5538 
5539 # ifdef DEBUG
5540  Assert(this->gradients_quad_initialized == true,
5542 # endif
5543  AssertIndexRange(q_point_in, this->n_quadrature_points);
5544 
5545  Assert(this->jacobian != nullptr,
5547  "update_gradients"));
5548 
5550 
5551  const std::size_t q_point = q_point_in;
5552  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5553  {
5554  for (unsigned int d = 0; d < dim; ++d)
5555  grad_out[d] =
5556  this->gradients_quad[dim * q_point + d] * this->jacobian[0][d][d];
5557  }
5558  // cell with general/affine Jacobian
5559  else
5560  {
5562  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5563  q_point :
5564  0];
5565  for (unsigned int d = 0; d < dim; ++d)
5566  {
5567  grad_out[d] = jac[d][0] * this->gradients_quad[dim * q_point];
5568  for (unsigned int e = 1; e < dim; ++e)
5569  grad_out[d] += jac[d][e] * this->gradients_quad[dim * q_point + e];
5570  }
5571  }
5572  return grad_out;
5573 }
5574 
5575 
5576 
5577 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5580  const unsigned int q_point) const
5581 {
5582  return BaseClass::get_hessian(q_point)[0];
5583 }
5584 
5585 
5586 
5587 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5590  get_hessian_diagonal(const unsigned int q_point) const
5591 {
5592  return BaseClass::get_hessian_diagonal(q_point)[0];
5593 }
5594 
5595 
5596 
5597 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5598 inline VectorizedArrayType
5600  const unsigned int q_point) const
5601 {
5602  return BaseClass::get_laplacian(q_point)[0];
5603 }
5604 
5605 
5606 
5607 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5608 inline void DEAL_II_ALWAYS_INLINE
5610  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5611 {
5612 # ifdef DEBUG
5613  this->dof_values_initialized = true;
5614  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5615 # endif
5616  this->values_dofs[dof] = val_in;
5617 }
5618 
5619 
5620 
5621 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5622 inline void DEAL_II_ALWAYS_INLINE
5624  const VectorizedArrayType val_in,
5625  const unsigned int q_point)
5626 {
5627 # ifdef DEBUG
5628  Assert(this->is_reinitialized, ExcNotInitialized());
5629 # endif
5630  AssertIndexRange(q_point, this->n_quadrature_points);
5631  Assert(this->J_value != nullptr,
5633  "update_value"));
5634 # ifdef DEBUG
5635  this->values_quad_submitted = true;
5636 # endif
5637 
5638  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5639  {
5640  const VectorizedArrayType JxW =
5641  this->J_value[0] * this->quadrature_weights[q_point];
5642  this->values_quad[q_point] = val_in * JxW;
5643  }
5644  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5645  {
5646  this->values_quad[q_point] = val_in * this->J_value[q_point];
5647  }
5648 }
5649 
5650 
5651 
5652 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5653 inline DEAL_II_ALWAYS_INLINE void
5655  const Tensor<1, 1, VectorizedArrayType> val_in,
5656  const unsigned int q_point)
5657 {
5658  submit_value(val_in[0], q_point);
5659 }
5660 
5661 
5662 
5663 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5664 inline DEAL_II_ALWAYS_INLINE void
5666  submit_normal_derivative(const VectorizedArrayType grad_in,
5667  const unsigned int q_point)
5668 {
5670  grad[0] = grad_in;
5671  BaseClass::submit_normal_derivative(grad, q_point);
5672 }
5673 
5674 
5675 
5676 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5677 inline DEAL_II_ALWAYS_INLINE void
5680  const unsigned int q_point_in)
5681 {
5682 # ifdef DEBUG
5683  Assert(this->is_reinitialized, ExcNotInitialized());
5684 # endif
5685  AssertIndexRange(q_point_in, this->n_quadrature_points);
5686  Assert(this->J_value != nullptr,
5688  "update_gradients"));
5689  Assert(this->jacobian != nullptr,
5691  "update_gradients"));
5692 # ifdef DEBUG
5693  this->gradients_quad_submitted = true;
5694 # endif
5695 
5696  const std::size_t q_point = q_point_in;
5697  VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point;
5698  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5699  {
5700  const VectorizedArrayType JxW =
5701  this->J_value[0] * this->quadrature_weights[q_point];
5702 
5703  // Make sure the compiler does not think 'jacobian' is aliased with
5704  // 'gradients_quad'
5705  std::array<VectorizedArrayType, dim> jac;
5706  for (unsigned int d = 0; d < dim; ++d)
5707  jac[d] = this->jacobian[0][d][d];
5708 
5709  for (unsigned int d = 0; d < dim; ++d)
5710  grad_ptr[d] = grad_in[d] * jac[d] * JxW;
5711  }
5712  // general/affine cell type
5713  else
5714  {
5716  this->cell_type > internal::MatrixFreeFunctions::affine ?
5717  this->jacobian[q_point] :
5718  this->jacobian[0];
5719  const VectorizedArrayType JxW =
5720  this->cell_type > internal::MatrixFreeFunctions::affine ?
5721  this->J_value[q_point] :
5722  this->J_value[0] * this->quadrature_weights[q_point];
5723  for (unsigned int d = 0; d < dim; ++d)
5724  {
5725  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5726  for (unsigned int e = 1; e < dim; ++e)
5727  new_val += jac[e][d] * grad_in[e];
5728  grad_ptr[d] = new_val * JxW;
5729  }
5730  }
5731 }
5732 
5733 
5734 
5735 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5736 inline DEAL_II_ALWAYS_INLINE void
5739  const unsigned int q_point)
5740 {
5742  hessian[0] = hessian_in;
5743  BaseClass::submit_hessian(hessian, q_point);
5744 }
5745 
5746 
5747 
5748 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5749 inline VectorizedArrayType
5751  integrate_value() const
5752 {
5753  return BaseClass::integrate_value()[0];
5754 }
5755 
5756 
5757 
5758 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5759 
5760 
5761 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5765  const unsigned int dof_no,
5766  const unsigned int first_selected_component,
5767  const unsigned int quad_no,
5768  const unsigned int fe_degree,
5769  const unsigned int n_q_points,
5770  const bool is_interior_face,
5771  const unsigned int active_fe_index,
5772  const unsigned int active_quad_index,
5773  const unsigned int face_type)
5774  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5775  matrix_free,
5776  dof_no,
5777  first_selected_component,
5778  quad_no,
5779  fe_degree,
5780  n_q_points,
5781  is_interior_face,
5782  active_fe_index,
5783  active_quad_index,
5784  face_type)
5785 {}
5786 
5787 
5788 
5789 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5792  const Mapping<dim> &mapping,
5793  const FiniteElement<dim> &fe,
5794  const Quadrature<1> &quadrature,
5795  const UpdateFlags update_flags,
5796  const unsigned int first_selected_component,
5798  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5799  mapping,
5800  fe,
5801  quadrature,
5802  update_flags,
5803  first_selected_component,
5804  other)
5805 {}
5806 
5807 
5808 
5809 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5813  &other)
5814  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5815 {}
5816 
5817 
5818 
5819 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5823  &other)
5824 {
5826  operator=(other);
5827  return *this;
5828 }
5829 
5830 
5831 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5834  const unsigned int q_point) const
5835 {
5836  if (this->data->element_type ==
5838  {
5839  // Piola transform is required
5840 # ifdef DEBUG
5841  Assert(this->values_quad_initialized == true,
5843 # endif
5844 
5845  AssertIndexRange(q_point, this->n_quadrature_points);
5846  Assert(this->J_value != nullptr,
5848  "update_values"));
5849  const std::size_t nqp = this->n_quadrature_points;
5851 
5852  if (!is_face &&
5853  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5854  {
5855  // Cartesian cell
5856  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5857  const VectorizedArrayType inv_det =
5858  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5859  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5860  this->jacobian[0][2][2];
5861 
5862  // J * u * det(J^-1)
5863  for (unsigned int comp = 0; comp < n_components; ++comp)
5864  value_out[comp] = this->values_quad[comp * nqp + q_point] *
5865  jac[comp][comp] * inv_det;
5866  }
5867  else
5868  {
5869  // Affine or general cell
5870  const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5871  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5872  this->jacobian[q_point] :
5873  this->jacobian[0];
5875  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5876  transpose(invert(inv_t_jac)) :
5877  this->jacobian[1];
5878 
5879  // Derivatives are reordered for faces. Need to take this into account
5880  const VectorizedArrayType inv_det =
5881  (is_face && dim == 2 && this->get_face_no() < 2) ?
5882  -determinant(inv_t_jac) :
5883  determinant(inv_t_jac);
5884  // J * u * det(J^-1)
5885  for (unsigned int comp = 0; comp < n_components; ++comp)
5886  {
5887  value_out[comp] = this->values_quad[q_point] * jac[comp][0];
5888  for (unsigned int e = 1; e < dim; ++e)
5889  value_out[comp] +=
5890  this->values_quad[e * nqp + q_point] * jac[comp][e];
5891  value_out[comp] *= inv_det;
5892  }
5893  }
5894  return value_out;
5895  }
5896  else
5897  {
5898  // No Piola needed
5899  return BaseClass::get_value(q_point);
5900  }
5901 }
5902 
5903 
5904 
5905 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5908  get_gradient(const unsigned int q_point) const
5909 {
5910  if (this->data->element_type ==
5912  {
5913  // Piola transform is required
5914 # ifdef DEBUG
5915  Assert(this->gradients_quad_initialized == true,
5917 # endif
5918 
5919  AssertIndexRange(q_point, this->n_quadrature_points);
5920  Assert(this->jacobian != nullptr,
5922  "update_gradients"));
5923  const std::size_t nqp = this->n_quadrature_points;
5924  const std::size_t nqp_d = nqp * dim;
5926  const VectorizedArrayType *gradients =
5927  this->gradients_quad + q_point * dim;
5928 
5929 
5930  if (!is_face &&
5931  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5932  {
5933  // Cartesian cell
5934  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5935  this->jacobian[0];
5936  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5937  const VectorizedArrayType inv_det =
5938  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5939  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5940  this->jacobian[0][2][2];
5941 
5942  // J * grad_quad * J^-1 * det(J^-1)
5943  for (unsigned int d = 0; d < dim; ++d)
5944  for (unsigned int comp = 0; comp < n_components; ++comp)
5945  grad_out[comp][d] = gradients[comp * nqp_d + d] *
5946  inv_t_jac[d][d] * (jac[comp][comp] * inv_det);
5947  }
5948  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5949  {
5950  // Affine cell
5951  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5952  this->jacobian[0];
5953  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5954 
5955  // Derivatives are reordered for faces. Need to take this into account
5956  const VectorizedArrayType inv_det =
5957  (is_face && dim == 2 && this->get_face_no() < 2) ?
5958  -determinant(inv_t_jac) :
5959  determinant(inv_t_jac);
5960 
5961  VectorizedArrayType tmp[dim][dim];
5962  // J * grad_quad * J^-1 * det(J^-1)
5963  for (unsigned int d = 0; d < dim; ++d)
5964  for (unsigned int e = 0; e < dim; ++e)
5965  {
5966  tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
5967  for (unsigned int f = 1; f < dim; ++f)
5968  tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
5969  }
5970  for (unsigned int comp = 0; comp < n_components; ++comp)
5971  for (unsigned int d = 0; d < dim; ++d)
5972  {
5973  VectorizedArrayType res = jac[comp][0] * tmp[d][0];
5974  for (unsigned int f = 1; f < dim; ++f)
5975  res += jac[comp][f] * tmp[d][f];
5976 
5977  grad_out[comp][d] = res * inv_det;
5978  }
5979  }
5980  else
5981  {
5982  // General cell
5983 
5984  // This assert could be removed if we make sure that this is updated
5985  // even though update_hessians or update_jacobian_grads is not passed,
5986  // i.e make the necessary changes in
5987  // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
5988  Assert(this->jacobian_gradients_non_inverse != nullptr,
5990  "update_hessians"));
5991 
5992  const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
5993  const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5994  this->jacobian[q_point];
5995 
5996  // Derivatives are reordered for faces. Need to take this into account
5997  const VectorizedArrayType inv_det =
5998  (is_face && dim == 2 && this->get_face_no() < 2) ?
5999  -determinant(inv_t_jac) :
6000  determinant(inv_t_jac);
6001  const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac);
6002 
6003  // (J * grad_quad) * J^-1 * det(J^-1), part in braces
6004  VectorizedArrayType tmp[dim][dim];
6005  for (unsigned int d = 0; d < dim; ++d)
6006  for (unsigned int e = 0; e < dim; ++e)
6007  {
6008  tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
6009  for (unsigned int f = 1; f < dim; ++f)
6010  tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
6011  }
6012 
6013  // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
6014  // outside braces with gradient part from above
6015  for (unsigned int d = 0; d < dim; ++d)
6016  {
6017  for (unsigned int e = 0; e < dim; ++e)
6018  tmp[e][d] +=
6019  jac_grad[e][d] * this->values_quad[e * nqp + q_point];
6020  for (unsigned int f = 0, r = dim; f < dim; ++f)
6021  for (unsigned int k = f + 1; k < dim; ++k, ++r)
6022  {
6023  tmp[k][d] +=
6024  jac_grad[r][d] * this->values_quad[f * nqp + q_point];
6025  tmp[f][d] +=
6026  jac_grad[r][d] * this->values_quad[k * nqp + q_point];
6027  }
6028  }
6029 
6030  // Apply J^{-1} appearing in both terms outside braces above
6031  for (unsigned int d = 0; d < dim; ++d)
6032  for (unsigned int e = 0; e < dim; ++e)
6033  {
6034  VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
6035  for (unsigned int f = 1; f < dim; ++f)
6036  res += tmp[f][d] * inv_t_jac[e][f];
6037  grad_out[d][e] = res;
6038  }
6039 
6040  // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
6041  // which can be expressed as a rank-1 update tmp[d] * tmp4[e], where
6042  // tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
6043  VectorizedArrayType tmp3[dim], tmp4[dim];
6044  for (unsigned int d = 0; d < dim; ++d)
6045  {
6046  tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
6047  for (unsigned int e = 1; e < dim; ++e)
6048  tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
6049  }
6050  for (unsigned int e = 0, k = dim; e < dim; ++e)
6051  for (unsigned int f = e + 1; f < dim; ++k, ++f)
6052  for (unsigned int d = 0; d < dim; ++d)
6053  {
6054  tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
6055  tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
6056  }
6057  for (unsigned int d = 0; d < dim; ++d)
6058  {
6059  tmp4[d] = tmp3[0] * inv_t_jac[d][0];
6060  for (unsigned int e = 1; e < dim; ++e)
6061  tmp4[d] += tmp3[e] * inv_t_jac[d][e];
6062  }
6063 
6064  VectorizedArrayType tmp2[dim];
6065  for (unsigned int d = 0; d < dim; ++d)
6066  {
6067  tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
6068  for (unsigned e = 1; e < dim; ++e)
6069  tmp2[d] += t_jac[e][d] * this->values_quad[e * nqp + q_point];
6070  }
6071 
6072  for (unsigned int d = 0; d < dim; ++d)
6073  for (unsigned int e = 0; e < dim; ++e)
6074  {
6075  grad_out[d][e] -= tmp4[e] * tmp2[d];
6076 
6077  // finally multiply by det(J^{-1}) necessary in all
6078  // contributions above
6079  grad_out[d][e] *= inv_det;
6080  }
6081  }
6082  return grad_out;
6083  }
6084  else
6085  {
6086  return BaseClass::get_gradient(q_point);
6087  }
6088 }
6089 
6090 
6091 
6092 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6093 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6095  get_divergence(const unsigned int q_point) const
6096 {
6097 # ifdef DEBUG
6098  Assert(this->gradients_quad_initialized == true,
6100 # endif
6101  AssertIndexRange(q_point, this->n_quadrature_points);
6102  Assert(this->jacobian != nullptr,
6104  "update_gradients"));
6105 
6106  VectorizedArrayType divergence;
6107  const std::size_t nqp = this->n_quadrature_points;
6108 
6109  if (this->data->element_type ==
6111  {
6112  VectorizedArrayType inv_det =
6113  (!is_face &&
6114  this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
6115  this->jacobian[0][0][0] *
6116  ((dim == 2) ? this->jacobian[0][1][1] :
6117  this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
6118  determinant(this->jacobian[this->cell_type >
6120  q_point :
6121  0]);
6122 
6123  // on faces in 2d, the determinant has the wrong sign due to ordering of
6124  // derivatives
6125  if (is_face && dim == 2 && this->get_face_no() < 2)
6126  inv_det = -inv_det;
6127 
6128  // div * det(J^-1)
6129  divergence = this->gradients_quad[q_point * dim];
6130  for (unsigned int d = 1; d < dim; ++d)
6131  divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
6132  divergence *= inv_det;
6133  }
6134  else
6135  {
6136  if (!is_face &&
6137  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6138  {
6139  // Cartesian cell
6140  divergence =
6141  this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
6142  for (unsigned int d = 1; d < dim; ++d)
6143  divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
6144  this->jacobian[0][d][d];
6145  }
6146  else
6147  {
6148  // cell with general/constant Jacobian
6150  this->cell_type == internal::MatrixFreeFunctions::general ?
6151  this->jacobian[q_point] :
6152  this->jacobian[0];
6153  divergence = jac[0][0] * this->gradients_quad[q_point * dim];
6154  for (unsigned int e = 1; e < dim; ++e)
6155  divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
6156  for (unsigned int d = 1; d < dim; ++d)
6157  for (unsigned int e = 0; e < dim; ++e)
6158  divergence +=
6159  jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
6160  }
6161  }
6162  return divergence;
6163 }
6164 
6165 
6166 
6167 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6170  get_symmetric_gradient(const unsigned int q_point) const
6171 {
6172  // copy from generic function into dim-specialization function
6173  const auto grad = get_gradient(q_point);
6174  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6175  VectorizedArrayType half = Number(0.5);
6176  for (unsigned int d = 0; d < dim; ++d)
6177  symmetrized[d] = grad[d][d];
6178  switch (dim)
6179  {
6180  case 1:
6181  break;
6182  case 2:
6183  symmetrized[2] = grad[0][1] + grad[1][0];
6184  symmetrized[2] *= half;
6185  break;
6186  case 3:
6187  symmetrized[3] = grad[0][1] + grad[1][0];
6188  symmetrized[3] *= half;
6189  symmetrized[4] = grad[0][2] + grad[2][0];
6190  symmetrized[4] *= half;
6191  symmetrized[5] = grad[1][2] + grad[2][1];
6192  symmetrized[5] *= half;
6193  break;
6194  default:
6195  Assert(false, ExcNotImplemented());
6196  }
6198 }
6199 
6200 
6201 
6202 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6203 inline DEAL_II_ALWAYS_INLINE
6204  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6206  const unsigned int q_point) const
6207 {
6208  // copy from generic function into dim-specialization function
6209  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6210  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6211  switch (dim)
6212  {
6213  case 1:
6214  Assert(false,
6215  ExcMessage(
6216  "Computing the curl in 1d is not a useful operation"));
6217  break;
6218  case 2:
6219  curl[0] = grad[1][0] - grad[0][1];
6220  break;
6221  case 3:
6222  curl[0] = grad[2][1] - grad[1][2];
6223  curl[1] = grad[0][2] - grad[2][0];
6224  curl[2] = grad[1][0] - grad[0][1];
6225  break;
6226  default:
6227  Assert(false, ExcNotImplemented());
6228  }
6229  return curl;
6230 }
6231 
6232 
6233 
6234 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6237  get_hessian_diagonal(const unsigned int q_point) const
6238 {
6239  return BaseClass::get_hessian_diagonal(q_point);
6240 }
6241 
6242 
6243 
6244 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6247  const unsigned int q_point) const
6248 {
6249 # ifdef DEBUG
6250  Assert(this->hessians_quad_initialized == true,
6252 # endif
6253  AssertIndexRange(q_point, this->n_quadrature_points);
6254  return BaseClass::get_hessian(q_point);
6255 }
6256 
6257 
6258 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6259 inline DEAL_II_ALWAYS_INLINE void
6262  const unsigned int q_point)
6263 {
6264  if (this->data->element_type ==
6266  {
6267  // Piola transform is required
6268  AssertIndexRange(q_point, this->n_quadrature_points);
6269  Assert(this->J_value != nullptr,
6271  "update_value"));
6272 # ifdef DEBUG
6273  Assert(this->is_reinitialized, ExcNotInitialized());
6274  this->values_quad_submitted = true;
6275 # endif
6276 
6277  VectorizedArrayType *values = this->values_quad + q_point;
6278  const std::size_t nqp = this->n_quadrature_points;
6279 
6280  if (!is_face &&
6281  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6282  {
6283  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6284  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6285 
6286  for (unsigned int comp = 0; comp < n_components; ++comp)
6287  values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
6288  }
6289  else
6290  {
6291  // Affine or general cell
6292  const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
6293  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6294  this->jacobian[q_point] :
6295  this->jacobian[0];
6296 
6297  // Derivatives are reordered for faces. Need to take this into account
6298  // and 1/inv_det != J_value for faces
6299  const VectorizedArrayType fac =
6300  (!is_face) ?
6301  this->quadrature_weights[q_point] :
6302  (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6303  this->J_value[q_point] :
6304  this->J_value[0] * this->quadrature_weights[q_point]) *
6305  ((dim == 2 && this->get_face_no() < 2) ?
6306  -determinant(inv_t_jac) :
6307  determinant(inv_t_jac)));
6309  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6310  transpose(invert(inv_t_jac)) :
6311  this->jacobian[1];
6312 
6313  // J^T * u * factor
6314  for (unsigned int comp = 0; comp < n_components; ++comp)
6315  {
6316  values[comp * nqp] = val_in[0] * jac[0][comp];
6317  for (unsigned int e = 1; e < dim; ++e)
6318  values[comp * nqp] += val_in[e] * jac[e][comp];
6319  values[comp * nqp] *= fac;
6320  }
6321  }
6322  }
6323  else
6324  {
6325  // No Piola transform
6326  BaseClass::submit_value(val_in, q_point);
6327  }
6328 }
6329 
6330 
6331 
6332 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6333 inline DEAL_II_ALWAYS_INLINE void
6336  const unsigned int q_point)
6337 {
6338  if (this->data->element_type ==
6340  {
6341  // Piola transform is required
6342 
6343 # ifdef DEBUG
6344  Assert(this->is_reinitialized, ExcNotInitialized());
6345 # endif
6346  AssertIndexRange(q_point, this->n_quadrature_points);
6347  Assert(this->J_value != nullptr,
6349  "update_gradients"));
6350  Assert(this->jacobian != nullptr,
6352  "update_gradients"));
6353 # ifdef DEBUG
6354  this->gradients_quad_submitted = true;
6355 # endif
6356 
6357  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
6358  VectorizedArrayType *values = this->values_from_gradients_quad + q_point;
6359  const std::size_t nqp = this->n_quadrature_points;
6360  const std::size_t nqp_d = nqp * dim;
6361 
6362  if (!is_face &&
6363  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6364  {
6365  // Cartesian cell
6366  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6367  this->jacobian[0];
6368  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6369  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6370  for (unsigned int d = 0; d < dim; ++d)
6371  for (unsigned int comp = 0; comp < n_components; ++comp)
6372  gradients[comp * nqp_d + d] =
6373  grad_in[comp][d] * inv_t_jac[d][d] * (jac[comp][comp] * weight);
6374  }
6375  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6376  {
6377  // Affine cell
6378  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6379  this->jacobian[0];
6380  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6381 
6382  // Derivatives are reordered for faces. Need to take this into account
6383  // and 1/inv_det != J_value for faces
6384  const VectorizedArrayType fac =
6385  (!is_face) ? this->quadrature_weights[q_point] :
6386  this->J_value[0] * this->quadrature_weights[q_point] *
6387  ((dim == 2 && this->get_face_no() < 2) ?
6388  -determinant(inv_t_jac) :
6389  determinant(inv_t_jac));
6390 
6391  // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6392  VectorizedArrayType tmp[dim][dim];
6393  for (unsigned int d = 0; d < dim; ++d)
6394  for (unsigned int e = 0; e < dim; ++e)
6395  {
6396  tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
6397  for (unsigned int f = 1; f < dim; ++f)
6398  tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
6399  }
6400  for (unsigned int comp = 0; comp < n_components; ++comp)
6401  for (unsigned int d = 0; d < dim; ++d)
6402  {
6403  VectorizedArrayType res = jac[0][comp] * tmp[d][0];
6404  for (unsigned int f = 1; f < dim; ++f)
6405  res += jac[f][comp] * tmp[d][f];
6406 
6407  gradients[comp * nqp_d + d] = res * fac;
6408  }
6409  }
6410  else
6411  {
6412  // General cell
6413 
6414  const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
6415  const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
6416  this->jacobian[q_point];
6417 
6418  // Derivatives are reordered for faces. Need to take this into account
6419  // and 1/inv_det != J_value for faces
6420  const VectorizedArrayType fac =
6421  (!is_face) ?
6422  this->quadrature_weights[q_point] :
6423  this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6424  -determinant(inv_t_jac) :
6425  determinant(inv_t_jac));
6426  const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac);
6427 
6428  // Start evaluation for values part below to enable the compiler to
6429  // possibly re-use the same computation in get_gradient() without
6430  // interfering with stores to 'gradients'
6431  VectorizedArrayType tmp3[dim], tmp4[dim];
6432  for (unsigned int d = 0; d < dim; ++d)
6433  {
6434  tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
6435  for (unsigned int e = 1; e < dim; ++e)
6436  tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
6437  }
6438  for (unsigned int e = 0, k = dim; e < dim; ++e)
6439  for (unsigned int f = e + 1; f < dim; ++k, ++f)
6440  for (unsigned int d = 0; d < dim; ++d)
6441  {
6442  tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
6443  tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
6444  }
6445  for (unsigned int d = 0; d < dim; ++d)
6446  {
6447  tmp4[d] = tmp3[0] * inv_t_jac[d][0];
6448  for (unsigned int e = 1; e < dim; ++e)
6449  tmp4[d] += tmp3[e] * inv_t_jac[d][e];
6450  }
6451 
6452  const Tensor<2, dim, VectorizedArrayType> grad_in_scaled =
6453  fac * grad_in;
6454 
6455  VectorizedArrayType tmp[dim][dim];
6456 
6457  // J * (J^{-1} * (grad_in * factor))
6458  for (unsigned int d = 0; d < dim; ++d)
6459  for (unsigned int e = 0; e < dim; ++e)
6460  {
6461  tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
6462  for (unsigned int f = 1; f < dim; ++f)
6463  tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
6464  }
6465 
6466  for (unsigned int d = 0; d < dim; ++d)
6467  for (unsigned int e = 0; e < dim; ++e)
6468  {
6469  VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
6470  for (unsigned int f = 1; f < dim; ++f)
6471  res += t_jac[d][f] * tmp[e][f];
6472 
6473  gradients[d * nqp_d + e] = res;
6474  }
6475 
6476  // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
6477  // as 'tmp' from above
6478  VectorizedArrayType value[dim];
6479  for (unsigned int d = 0; d < dim; ++d)
6480  {
6481  value[d] = tmp[d][0] * jac_grad[d][0];
6482  for (unsigned int e = 1; e < dim; ++e)
6483  value[d] += tmp[d][e] * jac_grad[d][e];
6484  }
6485  for (unsigned int e = 0, k = dim; e < dim; ++e)
6486  for (unsigned int f = e + 1; f < dim; ++k, ++f)
6487  for (unsigned int d = 0; d < dim; ++d)
6488  {
6489  value[e] += tmp[f][d] * jac_grad[k][d];
6490  value[f] += tmp[e][d] * jac_grad[k][d];
6491  }
6492 
6493  // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
6494  // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
6495  for (unsigned int d = 0; d < dim; ++d)
6496  {
6497  VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
6498  for (unsigned int e = 1; e < dim; ++e)
6499  tmp2 += grad_in_scaled[d][e] * tmp4[e];
6500  for (unsigned int e = 0; e < dim; ++e)
6501  value[e] -= t_jac[e][d] * tmp2;
6502  }
6503 
6504  for (unsigned int d = 0; d < dim; ++d)
6505  values[d * nqp] = value[d];
6506  }
6507  }
6508  else
6509  {
6510  BaseClass::submit_gradient(grad_in, q_point);
6511  }
6512 }
6513 
6514 
6515 
6516 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6517 inline DEAL_II_ALWAYS_INLINE void
6520  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6521  const unsigned int q_point)
6522 {
6523  if (this->data->element_type ==
6525  {
6526  // Piola transform is required
6527  const Tensor<2, dim, VectorizedArrayType> grad = grad_in;
6529  submit_gradient(grad, q_point);
6530  }
6531  else
6532  {
6533  BaseClass::submit_gradient(grad_in, q_point);
6534  }
6535 }
6536 
6537 
6538 
6539 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6540 inline DEAL_II_ALWAYS_INLINE void
6542  submit_divergence(const VectorizedArrayType div_in,
6543  const unsigned int q_point)
6544 {
6545 # ifdef DEBUG
6546  Assert(this->is_reinitialized, ExcNotInitialized());
6547 # endif
6548  AssertIndexRange(q_point, this->n_quadrature_points);
6549  Assert(this->J_value != nullptr,
6551  "update_gradients"));
6552  Assert(this->jacobian != nullptr,
6554  "update_gradients"));
6555 # ifdef DEBUG
6556  this->gradients_quad_submitted = true;
6557 # endif
6558 
6559  const std::size_t nqp_d = this->n_quadrature_points * dim;
6560  VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
6561 
6562  if (this->data->element_type ==
6564  {
6565  // General cell
6566 
6567  // Derivatives are reordered for faces. Need to take this into account
6568  // and 1/inv_det != J_value for faces
6569  const VectorizedArrayType fac =
6570  (!is_face) ?
6571  this->quadrature_weights[q_point] * div_in :
6572  (this->cell_type > internal::MatrixFreeFunctions::affine ?
6573  this->J_value[q_point] :
6574  this->J_value[0] * this->quadrature_weights[q_point]) *
6575  div_in *
6576  determinant(
6577  this->jacobian[this->cell_type >
6579  q_point :
6580  0]) *
6581  Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6582 
6583  for (unsigned int d = 0; d < dim; ++d)
6584  {
6585  for (unsigned int e = 0; e < dim; ++e)
6586  gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
6587  }
6588  this->divergence_is_requested = true;
6589  }
6590  else
6591  {
6592  if (!is_face &&
6593  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6594  {
6595  const VectorizedArrayType fac =
6596  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6597  for (unsigned int d = 0; d < dim; ++d)
6598  {
6599  const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
6600  for (unsigned int e = 0; e < dim; ++e)
6601  gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
6602  }
6603  }
6604  else
6605  {
6607  this->cell_type == internal::MatrixFreeFunctions::general ?
6608  this->jacobian[q_point] :
6609  this->jacobian[0];
6610  const VectorizedArrayType fac =
6611  (this->cell_type == internal::MatrixFreeFunctions::general ?
6612  this->J_value[q_point] :
6613  this->J_value[0] * this->quadrature_weights[q_point]) *
6614  div_in;
6615  for (unsigned int d = 0; d < dim; ++d)
6616  {
6617  for (unsigned int e = 0; e < dim; ++e)
6618  gradients[d * nqp_d + e] = jac[d][e] * fac;
6619  }
6620  }
6621  }
6622 }
6623 
6624 
6625 
6626 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6627 inline DEAL_II_ALWAYS_INLINE void
6631  const unsigned int q_point)
6632 {
6633  AssertThrow(
6634  this->data->element_type !=
6636  ExcNotImplemented());
6637 
6638  // could have used base class operator, but that involves some overhead
6639  // which is inefficient. it is nice to have the symmetric tensor because
6640  // that saves some operations
6641 # ifdef DEBUG
6642  Assert(this->is_reinitialized, ExcNotInitialized());
6643 # endif
6644  AssertIndexRange(q_point, this->n_quadrature_points);
6645  Assert(this->J_value != nullptr,
6647  "update_gradients"));
6648  Assert(this->jacobian != nullptr,
6650  "update_gradients"));
6651 # ifdef DEBUG
6652  this->gradients_quad_submitted = true;
6653 # endif
6654 
6655  const std::size_t nqp_d = this->n_quadrature_points * dim;
6656  VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
6657  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6658  {
6659  const VectorizedArrayType JxW =
6660  this->J_value[0] * this->quadrature_weights[q_point];
6661  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
6662  for (unsigned int d = 0; d < dim; ++d)
6663  gradients[d * nqp_d + d] =
6664  (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
6665  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6666  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6667  {
6668  const VectorizedArrayType value =
6669  sym_grad.access_raw_entry(counter) * JxW;
6670  gradients[e * nqp_d + d] = value * jac[d][d];
6671  gradients[d * nqp_d + e] = value * jac[e][e];
6672  }
6673  }
6674  // general/affine cell type
6675  else
6676  {
6677  const VectorizedArrayType JxW =
6678  this->cell_type == internal::MatrixFreeFunctions::general ?
6679  this->J_value[q_point] :
6680  this->J_value[0] * this->quadrature_weights[q_point];
6682  this->cell_type == internal::MatrixFreeFunctions::general ?
6683  this->jacobian[q_point] :
6684  this->jacobian[0];
6685  VectorizedArrayType weighted[dim][dim];
6686  for (unsigned int i = 0; i < dim; ++i)
6687  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6688  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6689  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6690  {
6691  const VectorizedArrayType value =
6692  sym_grad.access_raw_entry(counter) * JxW;
6693  weighted[i][j] = value;
6694  weighted[j][i] = value;
6695  }
6696  for (unsigned int comp = 0; comp < dim; ++comp)
6697  for (unsigned int d = 0; d < dim; ++d)
6698  {
6699  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6700  for (unsigned int e = 1; e < dim; ++e)
6701  new_val += jac[e][d] * weighted[comp][e];
6702  gradients[comp * nqp_d + d] = new_val;
6703  }
6704  }
6705 }
6706 
6707 
6708 
6709 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6710 inline DEAL_II_ALWAYS_INLINE void
6713  const unsigned int q_point)
6714 {
6716  switch (dim)
6717  {
6718  case 1:
6719  Assert(false,
6720  ExcMessage(
6721  "Testing by the curl in 1d is not a useful operation"));
6722  break;
6723  case 2:
6724  grad[1][0] = curl[0];
6725  grad[0][1] = -curl[0];
6726  break;
6727  case 3:
6728  grad[2][1] = curl[0];
6729  grad[1][2] = -curl[0];
6730  grad[0][2] = curl[1];
6731  grad[2][0] = -curl[1];
6732  grad[1][0] = curl[2];
6733  grad[0][1] = -curl[2];
6734  break;
6735  default:
6736  Assert(false, ExcNotImplemented());
6737  }
6738  submit_gradient(grad, q_point);
6739 }
6740 
6741 
6742 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6743 
6744 
6745 template <typename Number, bool is_face, typename VectorizedArrayType>
6748  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
6749  const unsigned int dof_no,
6750  const unsigned int first_selected_component,
6751  const unsigned int quad_no,
6752  const unsigned int fe_degree,
6753  const unsigned int n_q_points,
6754  const bool is_interior_face,
6755  const unsigned int active_fe_index,
6756  const unsigned int active_quad_index,
6757  const unsigned int face_type)
6758  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6759  matrix_free,
6760  dof_no,
6761  first_selected_component,
6762  quad_no,
6763  fe_degree,
6764  n_q_points,
6765  is_interior_face,
6766  active_fe_index,
6767  active_quad_index,
6768  face_type)
6769 {}
6770 
6771 
6772 
6773 template <typename Number, bool is_face, typename VectorizedArrayType>
6776  const Mapping<1> &mapping,
6777  const FiniteElement<1> &fe,
6778  const Quadrature<1> &quadrature,
6779  const UpdateFlags update_flags,
6780  const unsigned int first_selected_component,
6782  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6783  mapping,
6784  fe,
6785  quadrature,
6786  update_flags,
6787  first_selected_component,
6788  other)
6789 {}
6790 
6791 
6792 
6793 template <typename Number, bool is_face, typename VectorizedArrayType>
6797  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6798 {}
6799 
6800 
6801 
6802 template <typename Number, bool is_face, typename VectorizedArrayType>
6806 {
6808  other);
6809  return *this;
6810 }
6811 
6812 
6813 
6814 template <typename Number, bool is_face, typename VectorizedArrayType>
6815 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6817  const unsigned int dof) const
6818 {
6819  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6820  return this->values_dofs[dof];
6821 }
6822 
6823 
6824 
6825 template <typename Number, bool is_face, typename VectorizedArrayType>
6826 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6828  const unsigned int q_point) const
6829 {
6830 # ifdef DEBUG
6831  Assert(this->values_quad_initialized == true,
6833 # endif
6834  AssertIndexRange(q_point, this->n_quadrature_points);
6835  return this->values_quad[q_point];
6836 }
6837 
6838 
6839 
6840 template <typename Number, bool is_face, typename VectorizedArrayType>
6843  const unsigned int q_point) const
6844 {
6845  // could use the base class gradient, but that involves too many inefficient
6846  // initialization operations on tensors
6847 
6848 # ifdef DEBUG
6849  Assert(this->gradients_quad_initialized == true,
6851 # endif
6852  AssertIndexRange(q_point, this->n_quadrature_points);
6853 
6855  this->cell_type == internal::MatrixFreeFunctions::general ?
6856  this->jacobian[q_point] :
6857  this->jacobian[0];
6858 
6860  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6861 
6862  return grad_out;
6863 }
6864 
6865 
6866 
6867 template <typename Number, bool is_face, typename VectorizedArrayType>
6868 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6870  const unsigned int q_point) const
6871 {
6872  return get_gradient(q_point)[0];
6873 }
6874 
6875 
6876 
6877 template <typename Number, bool is_face, typename VectorizedArrayType>
6878 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6880  get_normal_derivative(const unsigned int q_point) const
6881 {
6882  return BaseClass::get_normal_derivative(q_point)[0];
6883 }
6884 
6885 
6886 
6887 template <typename Number, bool is_face, typename VectorizedArrayType>
6890  const unsigned int q_point) const
6891 {
6892  return BaseClass::get_hessian(q_point)[0];
6893 }
6894 
6895 
6896 
6897 template <typename Number, bool is_face, typename VectorizedArrayType>
6900  get_hessian_diagonal(const unsigned int q_point) const
6901 {
6902  return BaseClass::get_hessian_diagonal(q_point)[0];
6903 }
6904 
6905 
6906 
6907 template <typename Number, bool is_face, typename VectorizedArrayType>
6908 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6910  const unsigned int q_point) const
6911 {
6912  return BaseClass::get_laplacian(q_point)[0];
6913 }
6914 
6915 
6916 
6917 template <typename Number, bool is_face, typename VectorizedArrayType>
6920  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6921 {
6922 # ifdef DEBUG
6923  this->dof_values_initialized = true;
6924  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6925 # endif
6926  this->values_dofs[dof] = val_in;
6927 }
6928 
6929 
6930 
6931 template <typename Number, bool is_face, typename VectorizedArrayType>
6932 inline DEAL_II_ALWAYS_INLINE void
6934  const VectorizedArrayType val_in,
6935  const unsigned int q_point)
6936 {
6937 # ifdef DEBUG
6938  Assert(this->is_reinitialized, ExcNotInitialized());
6939 # endif
6940  AssertIndexRange(q_point, this->n_quadrature_points);
6941 # ifdef DEBUG
6942  this->values_quad_submitted = true;
6943 # endif
6944 
6945  if (this->cell_type == internal::MatrixFreeFunctions::general)
6946  {
6947  const VectorizedArrayType JxW = this->J_value[q_point];
6948  this->values_quad[q_point] = val_in * JxW;
6949  }
6950  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6951  {
6952  const VectorizedArrayType JxW =
6953  this->J_value[0] * this->quadrature_weights[q_point];
6954  this->values_quad[q_point] = val_in * JxW;
6955  }
6956 }
6957 
6958 
6959 
6960 template <typename Number, bool is_face, typename VectorizedArrayType>
6961 inline DEAL_II_ALWAYS_INLINE void
6963  const Tensor<1, 1, VectorizedArrayType> val_in,
6964  const unsigned int q_point)
6965 {
6966  submit_value(val_in[0], q_point);
6967 }
6968 
6969 
6970 
6971 template <typename Number, bool is_face, typename VectorizedArrayType>
6972 inline DEAL_II_ALWAYS_INLINE void
6974  const Tensor<1, 1, VectorizedArrayType> grad_in,
6975  const unsigned int q_point)
6976 {
6977  submit_gradient(grad_in[0], q_point);
6978 }
6979 
6980 
6981 
6982 template <typename Number, bool is_face, typename VectorizedArrayType>
6983 inline DEAL_II_ALWAYS_INLINE void
6985  const VectorizedArrayType grad_in,
6986  const unsigned int q_point)
6987 {
6988 # ifdef DEBUG
6989  Assert(this->is_reinitialized, ExcNotInitialized());
6990 # endif
6991  AssertIndexRange(q_point, this->n_quadrature_points);
6992 # ifdef DEBUG
6993  this->gradients_quad_submitted = true;
6994 # endif
6995 
6997  this->cell_type == internal::MatrixFreeFunctions::general ?
6998  this->jacobian[q_point] :
6999  this->jacobian[0];
7000  const VectorizedArrayType JxW =
7001  this->cell_type == internal::MatrixFreeFunctions::general ?
7002  this->J_value[q_point] :
7003  this->J_value[0] * this->quadrature_weights[q_point];
7004 
7005  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7006 }
7007 
7008 
7009 
7010 template <typename Number, bool is_face, typename VectorizedArrayType>
7011 inline DEAL_II_ALWAYS_INLINE void
7013  const Tensor<2, 1, VectorizedArrayType> grad_in,
7014  const unsigned int q_point)
7015 {
7016  submit_gradient(grad_in[0][0], q_point);
7017 }
7018 
7019 
7020 
7021 template <typename Number, bool is_face, typename VectorizedArrayType>
7022 inline DEAL_II_ALWAYS_INLINE void
7024  submit_normal_derivative(const VectorizedArrayType grad_in,
7025  const unsigned int q_point)
7026 {
7028  grad[0] = grad_in;
7029  BaseClass::submit_normal_derivative(grad, q_point);
7030 }
7031 
7032 
7033 
7034 template <typename Number, bool is_face, typename VectorizedArrayType>
7035 inline DEAL_II_ALWAYS_INLINE void
7038  const unsigned int q_point)
7039 {
7040  BaseClass::submit_normal_derivative(grad_in, q_point);
7041 }
7042 
7043 
7044 template <typename Number, bool is_face, typename VectorizedArrayType>
7045 inline DEAL_II_ALWAYS_INLINE void
7047  const Tensor<2, 1, VectorizedArrayType> hessian_in,
7048  const unsigned int q_point)
7049 {
7051  hessian[0] = hessian_in;
7052  BaseClass::submit_hessian(hessian, q_point);
7053 }
7054 
7055 
7056 template <typename Number, bool is_face, typename VectorizedArrayType>
7057 inline VectorizedArrayType
7059  integrate_value() const
7060 {
7061  return BaseClass::integrate_value()[0];
7062 }
7063 
7064 
7065 
7066 /*-------------------------- FEEvaluation -----------------------------------*/
7067 
7068 
7069 template <int dim,
7070  int fe_degree,
7071  int n_q_points_1d,
7072  int n_components_,
7073  typename Number,
7074  typename VectorizedArrayType>
7075 inline FEEvaluation<dim,
7076  fe_degree,
7077  n_q_points_1d,
7078  n_components_,
7079  Number,
7080  VectorizedArrayType>::
7081  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7082  const unsigned int fe_no,
7083  const unsigned int quad_no,
7084  const unsigned int first_selected_component,
7085  const unsigned int active_fe_index,
7086  const unsigned int active_quad_index)
7087  : BaseClass(matrix_free,
7088  fe_no,
7089  first_selected_component,
7090  quad_no,
7091  fe_degree,
7092  static_n_q_points,
7093  true /*note: this is not a face*/,
7094  active_fe_index,
7095  active_quad_index)
7096  , dofs_per_component(this->data->dofs_per_component_on_cell)
7097  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7098  , n_q_points(this->data->n_q_points)
7099 {
7100  check_template_arguments(fe_no, 0);
7101 }
7102 
7103 
7104 
7105 template <int dim,
7106  int fe_degree,
7107  int n_q_points_1d,
7108  int n_components_,
7109  typename Number,
7110  typename VectorizedArrayType>
7111 inline FEEvaluation<dim,
7112  fe_degree,
7113  n_q_points_1d,
7114  n_components_,
7115  Number,
7116  VectorizedArrayType>::
7117  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7118  const std::pair<unsigned int, unsigned int> &range,
7119  const unsigned int dof_no,
7120  const unsigned int quad_no,
7121  const unsigned int first_selected_component)
7122  : FEEvaluation(matrix_free,
7123  dof_no,
7124  quad_no,
7125  first_selected_component,
7126  matrix_free.get_cell_active_fe_index(range))
7127 {}
7128 
7129 
7130 
7131 template <int dim,
7132  int fe_degree,
7133  int n_q_points_1d,
7134  int n_components_,
7135  typename Number,
7136  typename VectorizedArrayType>
7137 inline FEEvaluation<dim,
7138  fe_degree,
7139  n_q_points_1d,
7140  n_components_,
7141  Number,
7142  VectorizedArrayType>::
7143  FEEvaluation(const Mapping<dim> &mapping,
7144  const FiniteElement<dim> &fe,
7145  const Quadrature<1> &quadrature,
7146  const UpdateFlags update_flags,
7147  const unsigned int first_selected_component)
7148  : BaseClass(mapping,
7149  fe,
7150  quadrature,
7151  update_flags,
7152  first_selected_component,
7153  nullptr)
7154  , dofs_per_component(this->data->dofs_per_component_on_cell)
7155  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7156  , n_q_points(this->data->n_q_points)
7157 {
7158  check_template_arguments(numbers::invalid_unsigned_int, 0);
7159 }
7160 
7161 
7162 
7163 template <int dim,
7164  int fe_degree,
7165  int n_q_points_1d,
7166  int n_components_,
7167  typename Number,
7168  typename VectorizedArrayType>
7169 inline FEEvaluation<dim,
7170  fe_degree,
7171  n_q_points_1d,
7172  n_components_,
7173  Number,
7174  VectorizedArrayType>::
7175  FEEvaluation(const FiniteElement<dim> &fe,
7176  const Quadrature<1> &quadrature,
7177  const UpdateFlags update_flags,
7178  const unsigned int first_selected_component)
7179  : BaseClass(StaticMappingQ1<dim>::mapping,
7180  fe,
7181  quadrature,
7182  update_flags,
7183  first_selected_component,
7184  nullptr)
7185  , dofs_per_component(this->data->dofs_per_component_on_cell)
7186  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7187  , n_q_points(this->data->n_q_points)
7188 {
7189  check_template_arguments(numbers::invalid_unsigned_int, 0);
7190 }
7191 
7192 
7193 
7194 template <int dim,
7195  int fe_degree,
7196  int n_q_points_1d,
7197  int n_components_,
7198  typename Number,
7199  typename VectorizedArrayType>
7200 inline FEEvaluation<dim,
7201  fe_degree,
7202  n_q_points_1d,
7203  n_components_,
7204  Number,
7205  VectorizedArrayType>::
7206  FEEvaluation(const FiniteElement<dim> &fe,
7208  const unsigned int first_selected_component)
7209  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7210  fe,
7211  other.mapped_geometry->get_quadrature(),
7212  other.mapped_geometry->get_fe_values().get_update_flags(),
7213  first_selected_component,
7214  &other)
7215  , dofs_per_component(this->data->dofs_per_component_on_cell)
7216  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7217  , n_q_points(this->data->n_q_points)
7218 {
7219  check_template_arguments(numbers::invalid_unsigned_int, 0);
7220 }
7221 
7222 
7223 
7224 template <int dim,
7225  int fe_degree,
7226  int n_q_points_1d,
7227  int n_components_,
7228  typename Number,
7229  typename VectorizedArrayType>
7230 inline FEEvaluation<dim,
7231  fe_degree,
7232  n_q_points_1d,
7233  n_components_,
7234  Number,
7235  VectorizedArrayType>::FEEvaluation(const FEEvaluation
7236  &other)
7237  : BaseClass(other)
7238  , dofs_per_component(this->data->dofs_per_component_on_cell)
7239  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7240  , n_q_points(this->data->n_q_points)
7241 {
7242  check_template_arguments(numbers::invalid_unsigned_int, 0);
7243 }
7244 
7245 
7246 
7247 template <int dim,
7248  int fe_degree,
7249  int n_q_points_1d,
7250  int n_components_,
7251  typename Number,
7252  typename VectorizedArrayType>
7253 inline FEEvaluation<dim,
7254  fe_degree,
7255  n_q_points_1d,
7256  n_components_,
7257  Number,
7258  VectorizedArrayType> &
7259 FEEvaluation<dim,
7260  fe_degree,
7261  n_q_points_1d,
7262  n_components_,
7263  Number,
7264  VectorizedArrayType>::operator=(const FEEvaluation &other)
7265 {
7266  BaseClass::operator=(other);
7267  check_template_arguments(numbers::invalid_unsigned_int, 0);
7268  return *this;
7269 }
7270 
7271 
7272 
7273 template <int dim,
7274  int fe_degree,
7275  int n_q_points_1d,
7276  int n_components_,
7277  typename Number,
7278  typename VectorizedArrayType>
7279 inline void
7280 FEEvaluation<dim,
7281  fe_degree,
7282  n_q_points_1d,
7283  n_components_,
7284  Number,
7285  VectorizedArrayType>::
7286  check_template_arguments(const unsigned int dof_no,
7287  const unsigned int first_selected_component)
7288 {
7289  (void)dof_no;
7290  (void)first_selected_component;
7291 
7292  Assert(
7293  this->data->dofs_per_component_on_cell > 0,
7294  ExcMessage(
7295  "There is nothing useful you can do with an FEEvaluation object with "
7296  "FE_Nothing, i.e., without DoFs! If you have passed to "
7297  "MatrixFree::reinit() a collection of finite elements also containing "
7298  "FE_Nothing, please check - before creating FEEvaluation - the category "
7299  "of the current range by calling either "
7300  "MatrixFree::get_cell_range_category(range) or "
7301  "MatrixFree::get_face_range_category(range). The returned category "
7302  "is the index of the active FE, which you can use to exclude "
7303  "FE_Nothing."));
7304 
7305 # ifdef DEBUG
7306  // print error message when the dimensions do not match. Propose a possible
7307  // fix
7308  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7309  static_cast<unsigned int>(fe_degree) !=
7310  this->data->data.front().fe_degree) ||
7311  n_q_points != this->n_quadrature_points)
7312  {
7313  std::string message =
7314  "-------------------------------------------------------\n";
7315  message += "Illegal arguments in constructor/wrong template arguments!\n";
7316  message += " Called --> FEEvaluation<dim,";
7317  message += Utilities::int_to_string(fe_degree) + ",";
7318  message += Utilities::int_to_string(n_q_points_1d);
7319  message += "," + Utilities::int_to_string(n_components);
7320  message += ",Number>(data";
7321  if (first_selected_component != numbers::invalid_unsigned_int)
7322  {
7323  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7324  message += Utilities::int_to_string(this->quad_no) + ", ";
7325  message += Utilities::int_to_string(first_selected_component);
7326  }
7327  message += ")\n";
7328 
7329  // check whether some other vector component has the correct number of
7330  // points
7331  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7332  proposed_fe_comp = numbers::invalid_unsigned_int,
7333  proposed_quad_comp = numbers::invalid_unsigned_int;
7334  if (dof_no != numbers::invalid_unsigned_int)
7335  {
7336  if (static_cast<unsigned int>(fe_degree) ==
7337  this->data->data.front().fe_degree)
7338  {
7339  proposed_dof_comp = dof_no;
7340  proposed_fe_comp = first_selected_component;
7341  }
7342  else
7343  for (unsigned int no = 0; no < this->matrix_free->n_components();
7344  ++no)
7345  for (unsigned int nf = 0;
7346  nf < this->matrix_free->n_base_elements(no);
7347  ++nf)
7348  if (this->matrix_free
7349  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7350  .data.front()
7351  .fe_degree == static_cast<unsigned int>(fe_degree))
7352  {
7353  proposed_dof_comp = no;
7354  proposed_fe_comp = nf;
7355  break;
7356  }
7357  if (n_q_points ==
7358  this->mapping_data->descriptor[this->active_quad_index]
7359  .n_q_points)
7360  proposed_quad_comp = this->quad_no;
7361  else
7362  for (unsigned int no = 0;
7363  no < this->matrix_free->get_mapping_info().cell_data.size();
7364  ++no)
7365  if (this->matrix_free->get_mapping_info()
7366  .cell_data[no]
7367  .descriptor[this->active_quad_index]
7368  .n_q_points == n_q_points)
7369  {
7370  proposed_quad_comp = no;
7371  break;
7372  }
7373  }
7374  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7375  proposed_quad_comp != numbers::invalid_unsigned_int)
7376  {
7377  if (proposed_dof_comp != first_selected_component)
7378  message += "Wrong vector component selection:\n";
7379  else
7380  message += "Wrong quadrature formula selection:\n";
7381  message += " Did you mean FEEvaluation<dim,";
7382  message += Utilities::int_to_string(fe_degree) + ",";
7383  message += Utilities::int_to_string(n_q_points_1d);
7384  message += "," + Utilities::int_to_string(n_components);
7385  message += ",Number>(data";
7386  if (dof_no != numbers::invalid_unsigned_int)
7387  {
7388  message +=
7389  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7390  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7391  message += Utilities::int_to_string(proposed_fe_comp);
7392  }
7393  message += ")?\n";
7394  std::string correct_pos;
7395  if (proposed_dof_comp != dof_no)
7396  correct_pos = " ^ ";
7397  else
7398  correct_pos = " ";
7399  if (proposed_quad_comp != this->quad_no)
7400  correct_pos += " ^ ";
7401  else
7402  correct_pos += " ";
7403  if (proposed_fe_comp != first_selected_component)
7404  correct_pos += " ^\n";
7405  else
7406  correct_pos += "