Reference documentation for deal.II version GIT 9c182271f7 2023-03-28 14:30:01+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 - 2022 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 
143  template <typename VectorType>
144  void
145  read_dof_values(const VectorType & src,
146  const unsigned int first_index = 0,
147  const std::bitset<VectorizedArrayType::size()> &mask =
148  std::bitset<VectorizedArrayType::size()>().flip());
149 
178  template <typename VectorType>
179  void
180  read_dof_values_plain(const VectorType & src,
181  const unsigned int first_index = 0,
182  const std::bitset<VectorizedArrayType::size()> &mask =
183  std::bitset<VectorizedArrayType::size()>().flip());
184 
216  template <typename VectorType>
217  void
219  VectorType & dst,
220  const unsigned int first_index = 0,
221  const std::bitset<VectorizedArrayType::size()> &mask =
222  std::bitset<VectorizedArrayType::size()>().flip()) const;
223 
262  template <typename VectorType>
263  void
264  set_dof_values(VectorType & dst,
265  const unsigned int first_index = 0,
266  const std::bitset<VectorizedArrayType::size()> &mask =
267  std::bitset<VectorizedArrayType::size()>().flip()) const;
268 
272  template <typename VectorType>
273  void
275  VectorType & dst,
276  const unsigned int first_index = 0,
277  const std::bitset<VectorizedArrayType::size()> &mask =
278  std::bitset<VectorizedArrayType::size()>().flip()) const;
279 
302  value_type
303  get_dof_value(const unsigned int dof) const;
304 
315  void
316  submit_dof_value(const value_type val_in, const unsigned int dof);
317 
330  value_type
331  get_value(const unsigned int q_point) const;
332 
345  void
346  submit_value(const value_type val_in, const unsigned int q_point);
347 
359  get_gradient(const unsigned int q_point) const;
360 
375  value_type
376  get_normal_derivative(const unsigned int q_point) const;
377 
390  void
391  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
392 
411  void
413  const unsigned int q_point);
414 
427  void
428  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
429 
442  get_hessian(const unsigned int q_point) const;
443 
454  get_hessian_diagonal(const unsigned int q_point) const;
455 
467  value_type
468  get_laplacian(const unsigned int q_point) const;
469 
470 #ifdef DOXYGEN
471  // doxygen does not anyhow mention functions coming from partial template
472  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
473  // For now, hack in those functions manually only to fix documentation:
474 
481  VectorizedArrayType
482  get_divergence(const unsigned int q_point) const;
483 
493  get_symmetric_gradient(const unsigned int q_point) const;
494 
501  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
502  get_curl(const unsigned int q_point) const;
503 
519  void
520  submit_divergence(const VectorizedArrayType div_in,
521  const unsigned int q_point);
522 
539  void
542  const unsigned int q_point);
543 
556  void
558  const unsigned int q_point);
559 
560 #endif
561 
578  value_type
580 
588 
589 protected:
600  const unsigned int dof_no,
601  const unsigned int first_selected_component,
602  const unsigned int quad_no,
603  const unsigned int fe_degree,
604  const unsigned int n_q_points,
605  const bool is_interior_face,
606  const unsigned int active_fe_index,
607  const unsigned int active_quad_index,
608  const unsigned int face_type);
609 
647  const Mapping<dim> & mapping,
648  const FiniteElement<dim> &fe,
649  const Quadrature<1> & quadrature,
650  const UpdateFlags update_flags,
651  const unsigned int first_selected_component,
653 
661 
670 
675 
682  template <typename VectorType, typename VectorOperation>
683  void
685  const VectorOperation & operation,
686  const std::array<VectorType *, n_components_> &vectors,
687  const std::array<
689  n_components_> & vectors_sm,
690  const std::bitset<VectorizedArrayType::size()> &mask,
691  const bool apply_constraints = true) const;
692 
700  template <typename VectorType, typename VectorOperation>
701  void
703  const VectorOperation & operation,
704  const std::array<VectorType *, n_components_> &vectors,
705  const std::array<
707  n_components_> & vectors_sm,
708  const std::bitset<VectorizedArrayType::size()> &mask) const;
709 
717  template <typename VectorType, typename VectorOperation>
718  void
720  const VectorOperation & operation,
721  const std::array<VectorType *, n_components_> &vectors) const;
722 
726  void
728 
733 
738 
743  mutable std::vector<types::global_dof_index> local_dof_indices;
744 };
745 
746 
747 
755 template <int dim,
756  int n_components_,
757  typename Number,
758  bool is_face,
759  typename VectorizedArrayType = VectorizedArray<Number>>
761  n_components_,
762  Number,
763  is_face,
764  VectorizedArrayType>
765 {
766  static_assert(
767  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
768  "Type of Number and of VectorizedArrayType do not match.");
769 
770 public:
771  using number_type = Number;
775  static constexpr unsigned int dimension = dim;
776  static constexpr unsigned int n_components = n_components_;
777  using BaseClass =
779 
780 protected:
790  const unsigned int dof_no,
791  const unsigned int first_selected_component,
792  const unsigned int quad_no,
793  const unsigned int fe_degree,
794  const unsigned int n_q_points,
795  const bool is_interior_face = true,
796  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
798  const unsigned int face_type = numbers::invalid_unsigned_int);
799 
805  const Mapping<dim> & mapping,
806  const FiniteElement<dim> &fe,
807  const Quadrature<1> & quadrature,
808  const UpdateFlags update_flags,
809  const unsigned int first_selected_component,
811 
816 
822 };
823 
824 
825 
834 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
835 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
836  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
837 {
838  static_assert(
839  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
840  "Type of Number and of VectorizedArrayType do not match.");
841 
842 public:
843  using number_type = Number;
844  using value_type = VectorizedArrayType;
847  static constexpr unsigned int dimension = dim;
848  using BaseClass =
850 
854  value_type
855  get_dof_value(const unsigned int dof) const;
856 
860  void
861  submit_dof_value(const value_type val_in, const unsigned int dof);
862 
866  value_type
867  get_value(const unsigned int q_point) const;
868 
872  void
873  submit_value(const value_type val_in, const unsigned int q_point);
874 
878  void
880  const unsigned int q_point);
881 
886  get_gradient(const unsigned int q_point) const;
887 
891  value_type
892  get_normal_derivative(const unsigned int q_point) const;
893 
897  void
898  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
899 
903  void
905  const unsigned int q_point);
906 
911  get_hessian(unsigned int q_point) const;
912 
917  get_hessian_diagonal(const unsigned int q_point) const;
918 
922  void
923  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
924 
928  value_type
929  get_laplacian(const unsigned int q_point) const;
930 
934  value_type
936 
937 protected:
947  const unsigned int dof_no,
948  const unsigned int first_selected_component,
949  const unsigned int quad_no,
950  const unsigned int fe_degree,
951  const unsigned int n_q_points,
952  const bool is_interior_face = true,
953  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
955  const unsigned int face_type = numbers::invalid_unsigned_int);
956 
962  const Mapping<dim> & mapping,
963  const FiniteElement<dim> &fe,
964  const Quadrature<1> & quadrature,
965  const UpdateFlags update_flags,
966  const unsigned int first_selected_component,
968 
973 
979 };
980 
981 
982 
992 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
993 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
994  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
995 {
996  static_assert(
997  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
998  "Type of Number and of VectorizedArrayType do not match.");
999 
1000 public:
1001  using number_type = Number;
1004  static constexpr unsigned int dimension = dim;
1005  static constexpr unsigned int n_components = dim;
1006  using BaseClass =
1008 
1012  value_type
1013  get_value(const unsigned int q_point) const;
1014 
1019  get_gradient(const unsigned int q_point) const;
1020 
1025  VectorizedArrayType
1026  get_divergence(const unsigned int q_point) const;
1027 
1035  get_symmetric_gradient(const unsigned int q_point) const;
1036 
1041  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1042  get_curl(const unsigned int q_point) const;
1043 
1048  get_hessian(const unsigned int q_point) const;
1049 
1054  get_hessian_diagonal(const unsigned int q_point) const;
1055 
1059  void
1061  const unsigned int q_point);
1062 
1066  void
1067  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1068 
1077  void
1079  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1080  const unsigned int q_point);
1081 
1090  void
1091  submit_divergence(const VectorizedArrayType div_in,
1092  const unsigned int q_point);
1093 
1102  void
1105  const unsigned int q_point);
1106 
1111  void
1113  const unsigned int q_point);
1114 
1115 protected:
1125  const unsigned int dof_no,
1126  const unsigned int first_selected_component,
1127  const unsigned int quad_no,
1128  const unsigned int dofs_per_cell,
1129  const unsigned int n_q_points,
1130  const bool is_interior_face = true,
1131  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1132  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1133  const unsigned int face_type = numbers::invalid_unsigned_int);
1134 
1140  const Mapping<dim> & mapping,
1141  const FiniteElement<dim> &fe,
1142  const Quadrature<1> & quadrature,
1143  const UpdateFlags update_flags,
1144  const unsigned int first_selected_component,
1146 
1151 
1157 };
1158 
1159 
1168 template <typename Number, bool is_face, typename VectorizedArrayType>
1169 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1170  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1171 {
1172  static_assert(
1173  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1174  "Type of Number and of VectorizedArrayType do not match.");
1175 
1176 public:
1177  using number_type = Number;
1178  using value_type = VectorizedArrayType;
1181  static constexpr unsigned int dimension = 1;
1182  using BaseClass =
1184 
1188  value_type
1189  get_dof_value(const unsigned int dof) const;
1190 
1194  void
1195  submit_dof_value(const value_type val_in, const unsigned int dof);
1196 
1200  value_type
1201  get_value(const unsigned int q_point) const;
1202 
1206  void
1207  submit_value(const value_type val_in, const unsigned int q_point);
1208 
1212  void
1213  submit_value(const gradient_type val_in, const unsigned int q_point);
1214 
1219  get_gradient(const unsigned int q_point) const;
1220 
1224  value_type
1225  get_divergence(const unsigned int q_point) const;
1226 
1230  value_type
1231  get_normal_derivative(const unsigned int q_point) const;
1232 
1236  void
1237  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1238 
1242  void
1243  submit_gradient(const value_type grad_in, const unsigned int q_point);
1244 
1248  void
1250  const unsigned int q_point);
1251 
1255  void
1257  const unsigned int q_point);
1258 
1262  void
1264  const unsigned int q_point);
1265 
1269  hessian_type
1270  get_hessian(unsigned int q_point) const;
1271 
1276  get_hessian_diagonal(const unsigned int q_point) const;
1277 
1281  void
1282  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1283 
1287  value_type
1288  get_laplacian(const unsigned int q_point) const;
1289 
1293  value_type
1295 
1296 protected:
1306  const unsigned int dof_no,
1307  const unsigned int first_selected_component,
1308  const unsigned int quad_no,
1309  const unsigned int fe_degree,
1310  const unsigned int n_q_points,
1311  const bool is_interior_face = true,
1312  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1313  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1314  const unsigned int face_type = numbers::invalid_unsigned_int);
1315 
1321  const Mapping<1> & mapping,
1322  const FiniteElement<1> &fe,
1323  const Quadrature<1> & quadrature,
1324  const UpdateFlags update_flags,
1325  const unsigned int first_selected_component,
1327 
1332 
1338 };
1339 
1340 
1341 
1897 template <int dim,
1898  int fe_degree,
1899  int n_q_points_1d,
1900  int n_components_,
1901  typename Number,
1902  typename VectorizedArrayType>
1904  n_components_,
1905  Number,
1906  false,
1907  VectorizedArrayType>
1908 {
1909  static_assert(
1910  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1911  "Type of Number and of VectorizedArrayType do not match.");
1912 
1913 public:
1917  using BaseClass =
1919 
1923  using number_type = Number;
1924 
1931 
1938 
1942  static constexpr unsigned int dimension = dim;
1943 
1948  static constexpr unsigned int n_components = n_components_;
1949 
1958  static constexpr unsigned int static_n_q_points =
1959  Utilities::pow(n_q_points_1d, dim);
1960 
1970  static constexpr unsigned int static_dofs_per_component =
1971  Utilities::pow(fe_degree + 1, dim);
1972 
1982  static constexpr unsigned int tensor_dofs_per_cell =
1984 
1994  static constexpr unsigned int static_dofs_per_cell =
1996 
2033  const unsigned int dof_no = 0,
2034  const unsigned int quad_no = 0,
2035  const unsigned int first_selected_component = 0,
2036  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2037  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2038 
2047  const std::pair<unsigned int, unsigned int> & range,
2048  const unsigned int dof_no = 0,
2049  const unsigned int quad_no = 0,
2050  const unsigned int first_selected_component = 0);
2051 
2080  FEEvaluation(const Mapping<dim> & mapping,
2081  const FiniteElement<dim> &fe,
2082  const Quadrature<1> & quadrature,
2083  const UpdateFlags update_flags,
2084  const unsigned int first_selected_component = 0);
2085 
2092  const Quadrature<1> & quadrature,
2093  const UpdateFlags update_flags,
2094  const unsigned int first_selected_component = 0);
2095 
2108  const unsigned int first_selected_component = 0);
2109 
2117 
2124  FEEvaluation &
2125  operator=(const FEEvaluation &other);
2126 
2135  void
2136  reinit(const unsigned int cell_batch_index);
2137 
2144  void
2145  reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids);
2146 
2159  template <bool level_dof_access>
2160  void
2162 
2173  void
2175 
2179  static bool
2180  fast_evaluation_supported(const unsigned int given_degree,
2181  const unsigned int given_n_q_points_1d);
2182 
2192  void
2194 
2199  DEAL_II_DEPRECATED void
2200  evaluate(const bool evaluate_values,
2201  const bool evaluate_gradients,
2202  const bool evaluate_hessians = false);
2203 
2216  void
2217  evaluate(const VectorizedArrayType * values_array,
2218  const EvaluationFlags::EvaluationFlags evaluation_flag);
2219 
2224  DEAL_II_DEPRECATED void
2225  evaluate(const VectorizedArrayType *values_array,
2226  const bool evaluate_values,
2227  const bool evaluate_gradients,
2228  const bool evaluate_hessians = false);
2229 
2243  template <typename VectorType>
2244  void
2245  gather_evaluate(const VectorType & input_vector,
2246  const EvaluationFlags::EvaluationFlags evaluation_flag);
2247 
2251  template <typename VectorType>
2252  DEAL_II_DEPRECATED void
2253  gather_evaluate(const VectorType &input_vector,
2254  const bool evaluate_values,
2255  const bool evaluate_gradients,
2256  const bool evaluate_hessians = false);
2257 
2268  void
2270 
2274  DEAL_II_DEPRECATED void
2275  integrate(const bool integrate_values, const bool integrate_gradients);
2276 
2288  void
2290  VectorizedArrayType * values_array,
2291  const bool sum_into_values = false);
2292 
2296  DEAL_II_DEPRECATED void
2297  integrate(const bool integrate_values,
2298  const bool integrate_gradients,
2299  VectorizedArrayType *values_array);
2300 
2314  template <typename VectorType>
2315  void
2317  VectorType & output_vector);
2318 
2322  template <typename VectorType>
2323  DEAL_II_DEPRECATED void
2324  integrate_scatter(const bool integrate_values,
2325  const bool integrate_gradients,
2326  VectorType &output_vector);
2327 
2334  dof_indices() const;
2335 
2342  const unsigned int dofs_per_component;
2343 
2350  const unsigned int dofs_per_cell;
2351 
2359  const unsigned int n_q_points;
2360 
2361 private:
2366  void
2367  check_template_arguments(const unsigned int fe_no,
2368  const unsigned int first_selected_component);
2369 };
2370 
2371 
2372 
2408 template <int dim,
2409  int fe_degree,
2410  int n_q_points_1d = fe_degree + 1,
2411  int n_components_ = 1,
2412  typename Number = double,
2413  typename VectorizedArrayType = VectorizedArray<Number>>
2415  n_components_,
2416  Number,
2417  true,
2418  VectorizedArrayType>
2419 {
2420  static_assert(
2421  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2422  "Type of Number and of VectorizedArrayType do not match.");
2423 
2424 public:
2428  using BaseClass =
2430 
2434  using number_type = Number;
2435 
2442 
2449 
2453  static constexpr unsigned int dimension = dim;
2454 
2459  static constexpr unsigned int n_components = n_components_;
2460 
2470  static constexpr unsigned int static_n_q_points =
2471  Utilities::pow(n_q_points_1d, dim - 1);
2472 
2481  static constexpr unsigned int static_n_q_points_cell =
2482  Utilities::pow(n_q_points_1d, dim);
2483 
2492  static constexpr unsigned int static_dofs_per_component =
2493  Utilities::pow(fe_degree + 1, dim);
2494 
2503  static constexpr unsigned int tensor_dofs_per_cell =
2505 
2514  static constexpr unsigned int static_dofs_per_cell =
2516 
2560  const bool is_interior_face = true,
2561  const unsigned int dof_no = 0,
2562  const unsigned int quad_no = 0,
2563  const unsigned int first_selected_component = 0,
2564  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2565  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
2566  const unsigned int face_type = numbers::invalid_unsigned_int);
2567 
2577  const std::pair<unsigned int, unsigned int> & range,
2578  const bool is_interior_face = true,
2579  const unsigned int dof_no = 0,
2580  const unsigned int quad_no = 0,
2581  const unsigned int first_selected_component = 0);
2582 
2593  void
2594  reinit(const unsigned int face_batch_number);
2595 
2603  void
2604  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2605 
2609  static bool
2610  fast_evaluation_supported(const unsigned int given_degree,
2611  const unsigned int given_n_q_points_1d);
2612 
2623  void
2625 
2629  DEAL_II_DEPRECATED void
2630  evaluate(const bool evaluate_values, const bool evaluate_gradients);
2631 
2644  void
2645  evaluate(const VectorizedArrayType * values_array,
2646  const EvaluationFlags::EvaluationFlags evaluation_flag);
2647 
2651  DEAL_II_DEPRECATED void
2652  evaluate(const VectorizedArrayType *values_array,
2653  const bool evaluate_values,
2654  const bool evaluate_gradients);
2655 
2667  template <typename VectorType>
2668  void
2669  gather_evaluate(const VectorType & input_vector,
2670  const EvaluationFlags::EvaluationFlags evaluation_flag);
2671 
2675  template <typename VectorType>
2676  DEAL_II_DEPRECATED void
2677  gather_evaluate(const VectorType &input_vector,
2678  const bool evaluate_values,
2679  const bool evaluate_gradients);
2680 
2690  void
2692 
2696  DEAL_II_DEPRECATED void
2697  integrate(const bool integrate_values, const bool integrate_gradients);
2698 
2707  void
2709  VectorizedArrayType * values_array);
2710 
2714  DEAL_II_DEPRECATED void
2715  integrate(const bool integrate_values,
2716  const bool integrate_gradients,
2717  VectorizedArrayType *values_array);
2718 
2730  template <typename VectorType>
2731  void
2733  VectorType & output_vector);
2734 
2738  template <typename VectorType>
2739  void
2740  integrate_scatter(const bool integrate_values,
2741  const bool integrate_gradients,
2742  VectorType &output_vector);
2743 
2750  dof_indices() const;
2751 
2756  bool
2757  at_boundary() const;
2758 
2772  boundary_id() const;
2773 
2780  const unsigned int dofs_per_component;
2781 
2788  const unsigned int dofs_per_cell;
2789 
2797  const unsigned int n_q_points;
2798 };
2799 
2800 
2801 
2802 namespace internal
2803 {
2804  namespace MatrixFreeFunctions
2805  {
2806  // a helper function to compute the number of DoFs of a DGP element at
2807  // compile time, depending on the degree
2808  template <int dim, int degree>
2810  {
2811  // this division is always without remainder
2812  static constexpr unsigned int value =
2813  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2814  };
2815 
2816  // base specialization: 1d elements have 'degree+1' degrees of freedom
2817  template <int degree>
2818  struct DGP_dofs_per_component<1, degree>
2819  {
2820  static constexpr unsigned int value = degree + 1;
2821  };
2822  } // namespace MatrixFreeFunctions
2823 } // namespace internal
2824 
2825 
2826 /*----------------------- Inline functions ----------------------------------*/
2827 
2828 #ifndef DOXYGEN
2829 
2830 
2831 namespace internal
2832 {
2833  // Extract all internal data pointers and indices in a single function that
2834  // get passed on to the constructor of FEEvaluationData, avoiding to look
2835  // things up multiple times
2836  template <bool is_face,
2837  int dim,
2838  typename Number,
2839  typename VectorizedArrayType>
2841  InitializationData
2842  extract_initialization_data(
2844  const unsigned int dof_no,
2845  const unsigned int first_selected_component,
2846  const unsigned int quad_no,
2847  const unsigned int fe_degree,
2848  const unsigned int n_q_points,
2849  const unsigned int active_fe_index_given,
2850  const unsigned int active_quad_index_given,
2851  const unsigned int face_type)
2852  {
2854  InitializationData init_data;
2855 
2856  init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2857  init_data.mapping_data =
2858  &internal::MatrixFreeFunctions::
2859  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2860  matrix_free.get_mapping_info(), quad_no);
2861 
2862  init_data.active_fe_index =
2863  fe_degree != numbers::invalid_unsigned_int ?
2864  init_data.dof_info->fe_index_from_degree(first_selected_component,
2865  fe_degree) :
2866  (active_fe_index_given != numbers::invalid_unsigned_int ?
2867  active_fe_index_given :
2868  0);
2869  init_data.active_quad_index =
2870  fe_degree == numbers::invalid_unsigned_int ?
2871  (active_quad_index_given != numbers::invalid_unsigned_int ?
2872  active_quad_index_given :
2873  std::min<unsigned int>(init_data.active_fe_index,
2874  init_data.mapping_data->descriptor.size() -
2875  1)) :
2876  init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2877 
2878  init_data.shape_info = &matrix_free.get_shape_info(
2879  dof_no,
2880  quad_no,
2881  init_data.dof_info->component_to_base_index[first_selected_component],
2882  init_data.active_fe_index,
2883  init_data.active_quad_index);
2884  init_data.descriptor =
2885  &init_data.mapping_data->descriptor
2886  [is_face ?
2887  (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2888  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2889  init_data.active_quad_index];
2890 
2891  return init_data;
2892  }
2893 } // namespace internal
2894 
2895 
2896 
2897 /*----------------------- FEEvaluationBase ----------------------------------*/
2898 
2899 template <int dim,
2900  int n_components_,
2901  typename Number,
2902  bool is_face,
2903  typename VectorizedArrayType>
2904 inline FEEvaluationBase<dim,
2905  n_components_,
2906  Number,
2907  is_face,
2908  VectorizedArrayType>::
2909  FEEvaluationBase(
2911  const unsigned int dof_no,
2912  const unsigned int first_selected_component,
2913  const unsigned int quad_no,
2914  const unsigned int fe_degree,
2915  const unsigned int n_q_points,
2916  const bool is_interior_face,
2917  const unsigned int active_fe_index,
2918  const unsigned int active_quad_index,
2919  const unsigned int face_type)
2920  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2921  internal::extract_initialization_data<is_face>(matrix_free,
2922  dof_no,
2923  first_selected_component,
2924  quad_no,
2925  fe_degree,
2926  n_q_points,
2927  active_fe_index,
2928  active_quad_index,
2929  face_type),
2930  is_interior_face,
2931  quad_no,
2932  first_selected_component)
2933  , scratch_data_array(matrix_free.acquire_scratch_data())
2934  , matrix_free(&matrix_free)
2935 {
2936  this->set_data_pointers(scratch_data_array, n_components_);
2937  Assert(
2938  this->dof_info->start_components.back() == 1 ||
2939  static_cast<int>(n_components_) <=
2940  static_cast<int>(
2941  this->dof_info->start_components
2942  [this->dof_info->component_to_base_index[first_selected_component] +
2943  1]) -
2944  first_selected_component,
2945  ExcMessage(
2946  "You tried to construct a vector-valued evaluator with " +
2947  std::to_string(n_components) +
2948  " components. However, "
2949  "the current base element has only " +
2951  this->dof_info->start_components
2952  [this->dof_info->component_to_base_index[first_selected_component] +
2953  1] -
2954  first_selected_component) +
2955  " components left when starting from local element index " +
2957  first_selected_component -
2958  this->dof_info->start_components
2959  [this->dof_info->component_to_base_index[first_selected_component]]) +
2960  " (global index " + std::to_string(first_selected_component) + ")"));
2961 
2962  // do not check for correct dimensions of data fields here, should be done
2963  // in derived classes
2964 }
2965 
2966 
2967 
2968 template <int dim,
2969  int n_components_,
2970  typename Number,
2971  bool is_face,
2972  typename VectorizedArrayType>
2973 inline FEEvaluationBase<dim,
2974  n_components_,
2975  Number,
2976  is_face,
2977  VectorizedArrayType>::
2978  FEEvaluationBase(
2979  const Mapping<dim> & mapping,
2980  const FiniteElement<dim> &fe,
2981  const Quadrature<1> & quadrature,
2982  const UpdateFlags update_flags,
2983  const unsigned int first_selected_component,
2985  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2986  other != nullptr &&
2987  other->mapped_geometry->get_quadrature() == quadrature ?
2988  other->mapped_geometry :
2989  std::make_shared<internal::MatrixFreeFunctions::
2990  MappingDataOnTheFly<dim, VectorizedArrayType>>(
2991  mapping,
2992  quadrature,
2993  update_flags),
2994  n_components_,
2995  first_selected_component)
2996  , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2997  , matrix_free(nullptr)
2998 {
2999  const unsigned int base_element_number =
3000  fe.component_to_base_index(first_selected_component).first;
3001  Assert(fe.element_multiplicity(base_element_number) == 1 ||
3002  fe.element_multiplicity(base_element_number) -
3003  first_selected_component >=
3004  n_components_,
3005  ExcMessage("The underlying element must at least contain as many "
3006  "components as requested by this class"));
3007  (void)base_element_number;
3008 
3009  Assert(this->data == nullptr, ExcInternalError());
3010  this->data =
3012  Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
3013  fe,
3014  fe.component_to_base_index(first_selected_component).first);
3015 
3016  this->set_data_pointers(scratch_data_array, n_components_);
3017 }
3018 
3019 
3020 
3021 template <int dim,
3022  int n_components_,
3023  typename Number,
3024  bool is_face,
3025  typename VectorizedArrayType>
3026 inline FEEvaluationBase<dim,
3027  n_components_,
3028  Number,
3029  is_face,
3030  VectorizedArrayType>::
3031  FEEvaluationBase(const FEEvaluationBase<dim,
3032  n_components_,
3033  Number,
3034  is_face,
3035  VectorizedArrayType> &other)
3036  : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
3037  , scratch_data_array(other.matrix_free == nullptr ?
3038  new AlignedVector<VectorizedArrayType>() :
3039  other.matrix_free->acquire_scratch_data())
3040  , matrix_free(other.matrix_free)
3041 {
3042  if (other.matrix_free == nullptr)
3043  {
3044  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3045  this->data =
3047  *other.data);
3048 
3049  // Create deep copy of mapped geometry for use in parallel
3050  this->mapped_geometry =
3051  std::make_shared<internal::MatrixFreeFunctions::
3052  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3053  other.mapped_geometry->get_fe_values().get_mapping(),
3054  other.mapped_geometry->get_quadrature(),
3055  other.mapped_geometry->get_fe_values().get_update_flags());
3056  this->mapping_data = &this->mapped_geometry->get_data_storage();
3057  this->cell = 0;
3058 
3059  this->jacobian =
3060  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3061  this->J_value =
3062  this->mapped_geometry->get_data_storage().JxW_values.begin();
3063  this->jacobian_gradients =
3064  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3065  this->jacobian_gradients_non_inverse =
3066  this->mapped_geometry->get_data_storage()
3067  .jacobian_gradients_non_inverse[0]
3068  .begin();
3069  this->quadrature_points =
3070  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3071  }
3072 
3073  this->set_data_pointers(scratch_data_array, n_components_);
3074 }
3075 
3076 
3077 
3078 template <int dim,
3079  int n_components_,
3080  typename Number,
3081  bool is_face,
3082  typename VectorizedArrayType>
3083 inline FEEvaluationBase<dim,
3084  n_components_,
3085  Number,
3086  is_face,
3087  VectorizedArrayType> &
3089 operator=(const FEEvaluationBase<dim,
3090  n_components_,
3091  Number,
3092  is_face,
3093  VectorizedArrayType> &other)
3094 {
3095  // release old memory
3096  if (matrix_free == nullptr)
3097  {
3098  delete this->data;
3099  delete scratch_data_array;
3100  }
3101  else
3102  {
3103  matrix_free->release_scratch_data(scratch_data_array);
3104  }
3105 
3107 
3108  matrix_free = other.matrix_free;
3109 
3110  if (other.matrix_free == nullptr)
3111  {
3112  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3113  this->data =
3115  *other.data);
3116  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3117 
3118  // Create deep copy of mapped geometry for use in parallel
3119  this->mapped_geometry =
3120  std::make_shared<internal::MatrixFreeFunctions::
3121  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3122  other.mapped_geometry->get_fe_values().get_mapping(),
3123  other.mapped_geometry->get_quadrature(),
3124  other.mapped_geometry->get_fe_values().get_update_flags());
3125  this->cell = 0;
3126  this->mapping_data = &this->mapped_geometry->get_data_storage();
3127  this->jacobian =
3128  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3129  this->J_value =
3130  this->mapped_geometry->get_data_storage().JxW_values.begin();
3131  this->jacobian_gradients =
3132  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3133  this->jacobian_gradients_non_inverse =
3134  this->mapped_geometry->get_data_storage()
3135  .jacobian_gradients_non_inverse[0]
3136  .begin();
3137  this->quadrature_points =
3138  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3139  }
3140  else
3141  {
3142  scratch_data_array = matrix_free->acquire_scratch_data();
3143  }
3144 
3145  this->set_data_pointers(scratch_data_array, n_components_);
3146 
3147  return *this;
3148 }
3149 
3150 
3151 
3152 template <int dim,
3153  int n_components_,
3154  typename Number,
3155  bool is_face,
3156  typename VectorizedArrayType>
3157 inline FEEvaluationBase<dim,
3158  n_components_,
3159  Number,
3160  is_face,
3161  VectorizedArrayType>::~FEEvaluationBase()
3162 {
3163  if (matrix_free != nullptr)
3164  {
3165  try
3166  {
3167  matrix_free->release_scratch_data(scratch_data_array);
3168  }
3169  catch (...)
3170  {}
3171  }
3172  else
3173  {
3174  delete scratch_data_array;
3175  delete this->data;
3176  }
3177 }
3178 
3179 
3180 
3181 template <int dim,
3182  int n_components_,
3183  typename Number,
3184  bool is_face,
3185  typename VectorizedArrayType>
3188  get_matrix_free() const
3189 {
3190  Assert(matrix_free != nullptr,
3191  ExcMessage(
3192  "FEEvaluation was not initialized with a MatrixFree object!"));
3193  return *matrix_free;
3194 }
3195 
3196 
3197 
3198 namespace internal
3199 {
3200  // given a block vector return the underlying vector type
3201  // including constness (specified by bool)
3202  template <typename VectorType, bool>
3203  struct ConstBlockVectorSelector;
3204 
3205  template <typename VectorType>
3206  struct ConstBlockVectorSelector<VectorType, true>
3207  {
3208  using BaseVectorType = const typename VectorType::BlockType;
3209  };
3210 
3211  template <typename VectorType>
3212  struct ConstBlockVectorSelector<VectorType, false>
3213  {
3214  using BaseVectorType = typename VectorType::BlockType;
3215  };
3216 
3217  // allows to select between block vectors and non-block vectors, which
3218  // allows to use a unified interface for extracting blocks on block vectors
3219  // and doing nothing on usual vectors
3220  template <typename VectorType, bool>
3221  struct BlockVectorSelector;
3222 
3223  template <typename VectorType>
3224  struct BlockVectorSelector<VectorType, true>
3225  {
3226  using BaseVectorType = typename ConstBlockVectorSelector<
3227  VectorType,
3228  std::is_const<VectorType>::value>::BaseVectorType;
3229 
3230  static BaseVectorType *
3231  get_vector_component(VectorType &vec, const unsigned int component)
3232  {
3233  AssertIndexRange(component, vec.n_blocks());
3234  return &vec.block(component);
3235  }
3236  };
3237 
3238  template <typename VectorType>
3239  struct BlockVectorSelector<VectorType, false>
3240  {
3241  using BaseVectorType = VectorType;
3242 
3243  static BaseVectorType *
3244  get_vector_component(VectorType &vec, const unsigned int component)
3245  {
3246  // FEEvaluation allows to combine several vectors from a scalar
3247  // FiniteElement into a "vector-valued" FEEvaluation object with
3248  // multiple components. These components can be extracted with the other
3249  // get_vector_component functions. If we do not get a vector of vectors
3250  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3251  // must make sure that we do not duplicate the components in input
3252  // and/or duplicate the resulting integrals. In such a case, we should
3253  // only get the zeroth component in the vector contained set nullptr for
3254  // the others which allows us to catch unintended use in
3255  // read_write_operation.
3256  if (component == 0)
3257  return &vec;
3258  else
3259  return nullptr;
3260  }
3261  };
3262 
3263  template <typename VectorType>
3264  struct BlockVectorSelector<std::vector<VectorType>, false>
3265  {
3266  using BaseVectorType = VectorType;
3267 
3268  static BaseVectorType *
3269  get_vector_component(std::vector<VectorType> &vec,
3270  const unsigned int component)
3271  {
3272  AssertIndexRange(component, vec.size());
3273  return &vec[component];
3274  }
3275  };
3276 
3277  template <typename VectorType>
3278  struct BlockVectorSelector<const std::vector<VectorType>, false>
3279  {
3280  using BaseVectorType = const VectorType;
3281 
3282  static const BaseVectorType *
3283  get_vector_component(const std::vector<VectorType> &vec,
3284  const unsigned int component)
3285  {
3286  AssertIndexRange(component, vec.size());
3287  return &vec[component];
3288  }
3289  };
3290 
3291  template <typename VectorType>
3292  struct BlockVectorSelector<std::vector<VectorType *>, false>
3293  {
3294  using BaseVectorType = VectorType;
3295 
3296  static BaseVectorType *
3297  get_vector_component(std::vector<VectorType *> &vec,
3298  const unsigned int component)
3299  {
3300  AssertIndexRange(component, vec.size());
3301  return vec[component];
3302  }
3303  };
3304 
3305  template <typename VectorType>
3306  struct BlockVectorSelector<const std::vector<VectorType *>, false>
3307  {
3308  using BaseVectorType = const VectorType;
3309 
3310  static const BaseVectorType *
3311  get_vector_component(const std::vector<VectorType *> &vec,
3312  const unsigned int component)
3313  {
3314  AssertIndexRange(component, vec.size());
3315  return vec[component];
3316  }
3317  };
3318 } // namespace internal
3319 
3320 
3321 
3322 template <int dim,
3323  int n_components_,
3324  typename Number,
3325  bool is_face,
3326  typename VectorizedArrayType>
3327 template <typename VectorType, typename VectorOperation>
3328 inline void
3331  const VectorOperation & operation,
3332  const std::array<VectorType *, n_components_> &src,
3333  const std::array<
3335  n_components_> & src_sm,
3336  const std::bitset<VectorizedArrayType::size()> &mask,
3337  const bool apply_constraints) const
3338 {
3339  // Case 1: No MatrixFree object given, simple case because we do not need to
3340  // process constraints and need not care about vectorization -> go to
3341  // separate function
3342  if (this->matrix_free == nullptr)
3343  {
3344  read_write_operation_global(operation, src);
3345  return;
3346  }
3347 
3348  Assert(this->dof_info != nullptr, ExcNotInitialized());
3349  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3350  Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3351  if (this->n_fe_components == 1)
3352  for (unsigned int comp = 0; comp < n_components; ++comp)
3353  {
3354  Assert(src[comp] != nullptr,
3355  ExcMessage("The finite element underlying this FEEvaluation "
3356  "object is scalar, but you requested " +
3357  std::to_string(n_components) +
3358  " components via the template argument in "
3359  "FEEvaluation. In that case, you must pass an "
3360  "std::vector<VectorType> or a BlockVector to " +
3361  "read_dof_values and distribute_local_to_global."));
3363  *this->matrix_free,
3364  *this->dof_info);
3365  }
3366  else
3367  {
3369  *this->matrix_free,
3370  *this->dof_info);
3371  }
3372 
3373  // Case 2: contiguous indices which use reduced storage of indices and can
3374  // use vectorized load/store operations -> go to separate function
3375  if (this->cell != numbers::invalid_unsigned_int)
3376  {
3378  this->cell,
3379  dof_info.index_storage_variants[this->dof_access_index].size());
3380  if (dof_info.index_storage_variants
3381  [is_face ? this->dof_access_index :
3383  [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::
3384  IndexStorageVariants::contiguous)
3385  {
3386  read_write_operation_contiguous(operation, src, src_sm, mask);
3387  return;
3388  }
3389  }
3390 
3391  // Case 3: standard operation with one index per degree of freedom -> go on
3392  // here
3393  constexpr unsigned int n_lanes = VectorizedArrayType::size();
3394 
3395  std::array<unsigned int, VectorizedArrayType::size()> cells =
3396  this->get_cell_ids();
3397 
3398  const bool masking_is_active = mask.count() < n_lanes;
3399  if (masking_is_active)
3400  for (unsigned int v = 0; v < n_lanes; ++v)
3401  if (mask[v] == false)
3402  cells[v] = numbers::invalid_unsigned_int;
3403 
3404  bool has_hn_constraints = false;
3405 
3406  if (is_face == false)
3407  {
3408  if (!dof_info.hanging_node_constraint_masks.empty() &&
3409  !dof_info.hanging_node_constraint_masks_comp.empty() &&
3410  dof_info
3411  .hanging_node_constraint_masks_comp[this->active_fe_index]
3412  [this->first_selected_component])
3413  for (unsigned int v = 0; v < n_lanes; ++v)
3414  if (cells[v] != numbers::invalid_unsigned_int &&
3415  dof_info.hanging_node_constraint_masks[cells[v]] !=
3418  has_hn_constraints = true;
3419  }
3420 
3421  std::integral_constant<bool,
3422  internal::is_vectorizable<VectorType, Number>::value>
3423  vector_selector;
3424 
3425  const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
3426 
3427  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3428  std::array<VectorizedArrayType *, n_components> values_dofs;
3429  for (unsigned int c = 0; c < n_components; ++c)
3430  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3431  c * dofs_per_component;
3432 
3433  if (this->cell != numbers::invalid_unsigned_int &&
3434  dof_info.index_storage_variants
3435  [is_face ? this->dof_access_index :
3437  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3438  IndexStorageVariants::interleaved &&
3439  use_vectorized_path)
3440  {
3441  const unsigned int *dof_indices =
3442  dof_info.dof_indices_interleaved.data() +
3443  dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
3444  .first +
3445  this->dof_info
3446  ->component_dof_indices_offset[this->active_fe_index]
3447  [this->first_selected_component] *
3448  n_lanes;
3449  if (n_components == 1 || this->n_fe_components == 1)
3450  for (unsigned int i = 0; i < dofs_per_component;
3451  ++i, dof_indices += n_lanes)
3452  for (unsigned int comp = 0; comp < n_components; ++comp)
3453  operation.process_dof_gather(dof_indices,
3454  *src[comp],
3455  0,
3456  values_dofs[comp][i],
3457  vector_selector);
3458  else
3459  for (unsigned int comp = 0; comp < n_components; ++comp)
3460  for (unsigned int i = 0; i < dofs_per_component;
3461  ++i, dof_indices += n_lanes)
3462  operation.process_dof_gather(
3463  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3464  return;
3465  }
3466 
3467  // Allocate pointers, then initialize all of them to nullptrs and
3468  // below overwrite the ones we actually use:
3469  std::array<const unsigned int *, n_lanes> dof_indices;
3470  dof_indices.fill(nullptr);
3471 
3472  // Assign the appropriate cell ids for face/cell case and get the pointers
3473  // to the dof indices of the cells on all lanes
3474 
3475  bool has_constraints = false;
3476  const unsigned int n_components_read =
3477  this->n_fe_components > 1 ? n_components : 1;
3478 
3479  if (is_face)
3480  {
3481  for (unsigned int v = 0; v < n_lanes; ++v)
3482  {
3483  if (cells[v] == numbers::invalid_unsigned_int)
3484  continue;
3485 
3486  Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
3487  const std::pair<unsigned int, unsigned int> *my_index_start =
3488  &dof_info.row_starts[cells[v] * this->n_fe_components +
3489  this->first_selected_component];
3490 
3491  // check whether any of the SIMD lanes has constraints, i.e., the
3492  // constraint indicator which is the second entry of row_starts
3493  // increments on this cell
3494  if (my_index_start[n_components_read].second !=
3495  my_index_start[0].second)
3496  has_constraints = true;
3497 
3498  dof_indices[v] =
3499  dof_info.dof_indices.data() + my_index_start[0].first;
3500  }
3501  }
3502  else
3503  {
3504  for (unsigned int v = 0; v < n_lanes; ++v)
3505  {
3506  if (cells[v] == numbers::invalid_unsigned_int)
3507  continue;
3508 
3509  const std::pair<unsigned int, unsigned int> *my_index_start =
3510  &dof_info.row_starts[cells[v] * this->n_fe_components +
3511  this->first_selected_component];
3512  if (my_index_start[n_components_read].second !=
3513  my_index_start[0].second)
3514  has_constraints = true;
3515 
3516  if (dof_info.hanging_node_constraint_masks.size() > 0 &&
3517  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3518  dof_info.hanging_node_constraint_masks[cells[v]] !=
3522  [this->active_fe_index][this->first_selected_component])
3523  has_hn_constraints = true;
3524 
3525  Assert(my_index_start[n_components_read].first ==
3526  my_index_start[0].first ||
3527  my_index_start[0].first < dof_info.dof_indices.size(),
3528  ExcIndexRange(0,
3529  my_index_start[0].first,
3530  dof_info.dof_indices.size()));
3531  dof_indices[v] =
3532  dof_info.dof_indices.data() + my_index_start[0].first;
3533  }
3534  }
3535 
3536  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3537  return i != numbers::invalid_unsigned_int;
3538  }) < n_lanes)
3539  for (unsigned int comp = 0; comp < n_components; ++comp)
3540  for (unsigned int i = 0; i < dofs_per_component; ++i)
3541  operation.process_empty(values_dofs[comp][i]);
3542 
3543  // Case where we have no constraints throughout the whole cell: Can go
3544  // through the list of DoFs directly
3545  if (!has_constraints && apply_constraints)
3546  {
3547  if (n_components == 1 || this->n_fe_components == 1)
3548  {
3549  for (unsigned int v = 0; v < n_lanes; ++v)
3550  {
3551  if (cells[v] == numbers::invalid_unsigned_int)
3552  continue;
3553 
3554  for (unsigned int i = 0; i < dofs_per_component; ++i)
3555  for (unsigned int comp = 0; comp < n_components; ++comp)
3556  operation.process_dof(dof_indices[v][i],
3557  *src[comp],
3558  values_dofs[comp][i][v]);
3559  }
3560  }
3561  else
3562  {
3563  for (unsigned int comp = 0; comp < n_components; ++comp)
3564  for (unsigned int v = 0; v < n_lanes; ++v)
3565  {
3566  if (cells[v] == numbers::invalid_unsigned_int)
3567  continue;
3568 
3569  for (unsigned int i = 0; i < dofs_per_component; ++i)
3570  operation.process_dof(
3571  dof_indices[v][comp * dofs_per_component + i],
3572  *src[0],
3573  values_dofs[comp][i][v]);
3574  }
3575  }
3576  return;
3577  }
3578 
3579  // In the case where there are some constraints to be resolved, loop over
3580  // all vector components that are filled and then over local dofs. ind_local
3581  // holds local number on cell, index iterates over the elements of
3582  // index_local_to_global and dof_indices points to the global indices stored
3583  // in index_local_to_global
3584 
3585  for (unsigned int v = 0; v < n_lanes; ++v)
3586  {
3587  if (cells[v] == numbers::invalid_unsigned_int)
3588  continue;
3589 
3590  const unsigned int cell_index = cells[v];
3591  const unsigned int cell_dof_index =
3592  cell_index * this->n_fe_components + this->first_selected_component;
3593  const unsigned int n_components_read =
3594  this->n_fe_components > 1 ? n_components : 1;
3595  unsigned int index_indicators =
3596  dof_info.row_starts[cell_dof_index].second;
3597  unsigned int next_index_indicators =
3598  dof_info.row_starts[cell_dof_index + 1].second;
3599 
3600  // For read_dof_values_plain, redirect the dof_indices field to the
3601  // unconstrained indices
3602  if (apply_constraints == false &&
3603  (dof_info.row_starts[cell_dof_index].second !=
3604  dof_info.row_starts[cell_dof_index + n_components_read].second ||
3605  ((dof_info.hanging_node_constraint_masks.size() > 0 &&
3606  dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3611  [this->active_fe_index][this->first_selected_component])))
3612  {
3615  ExcNotInitialized());
3616  dof_indices[v] =
3617  dof_info.plain_dof_indices.data() +
3618  this->dof_info
3619  ->component_dof_indices_offset[this->active_fe_index]
3620  [this->first_selected_component] +
3622  next_index_indicators = index_indicators;
3623  }
3624 
3625  if (n_components == 1 || this->n_fe_components == 1)
3626  {
3627  unsigned int ind_local = 0;
3628  for (; index_indicators != next_index_indicators; ++index_indicators)
3629  {
3630  const std::pair<unsigned short, unsigned short> indicator =
3631  dof_info.constraint_indicator[index_indicators];
3632  // run through values up to next constraint
3633  for (unsigned int j = 0; j < indicator.first; ++j)
3634  for (unsigned int comp = 0; comp < n_components; ++comp)
3635  operation.process_dof(dof_indices[v][j],
3636  *src[comp],
3637  values_dofs[comp][ind_local + j][v]);
3638 
3639  ind_local += indicator.first;
3640  dof_indices[v] += indicator.first;
3641 
3642  // constrained case: build the local value as a linear
3643  // combination of the global value according to constraints
3644  Number value[n_components];
3645  for (unsigned int comp = 0; comp < n_components; ++comp)
3646  operation.pre_constraints(values_dofs[comp][ind_local][v],
3647  value[comp]);
3648 
3649  const Number *data_val =
3650  this->matrix_free->constraint_pool_begin(indicator.second);
3651  const Number *end_pool =
3652  this->matrix_free->constraint_pool_end(indicator.second);
3653  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3654  for (unsigned int comp = 0; comp < n_components; ++comp)
3655  operation.process_constraint(*dof_indices[v],
3656  *data_val,
3657  *src[comp],
3658  value[comp]);
3659 
3660  for (unsigned int comp = 0; comp < n_components; ++comp)
3661  operation.post_constraints(value[comp],
3662  values_dofs[comp][ind_local][v]);
3663  ind_local++;
3664  }
3665 
3666  AssertIndexRange(ind_local, dofs_per_component + 1);
3667 
3668  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3669  for (unsigned int comp = 0; comp < n_components; ++comp)
3670  operation.process_dof(*dof_indices[v],
3671  *src[comp],
3672  values_dofs[comp][ind_local][v]);
3673  }
3674  else
3675  {
3676  // case with vector-valued finite elements where all components are
3677  // included in one single vector. Assumption: first come all entries
3678  // to the first component, then all entries to the second one, and
3679  // so on. This is ensured by the way MatrixFree reads out the
3680  // indices.
3681  for (unsigned int comp = 0; comp < n_components; ++comp)
3682  {
3683  unsigned int ind_local = 0;
3684 
3685  // check whether there is any constraint on the current cell
3686  for (; index_indicators != next_index_indicators;
3687  ++index_indicators)
3688  {
3689  const std::pair<unsigned short, unsigned short> indicator =
3690  dof_info.constraint_indicator[index_indicators];
3691 
3692  // run through values up to next constraint
3693  for (unsigned int j = 0; j < indicator.first; ++j)
3694  operation.process_dof(dof_indices[v][j],
3695  *src[0],
3696  values_dofs[comp][ind_local + j][v]);
3697  ind_local += indicator.first;
3698  dof_indices[v] += indicator.first;
3699 
3700  // constrained case: build the local value as a linear
3701  // combination of the global value according to constraints
3702  Number value;
3703  operation.pre_constraints(values_dofs[comp][ind_local][v],
3704  value);
3705 
3706  const Number *data_val =
3707  this->matrix_free->constraint_pool_begin(indicator.second);
3708  const Number *end_pool =
3709  this->matrix_free->constraint_pool_end(indicator.second);
3710 
3711  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3712  operation.process_constraint(*dof_indices[v],
3713  *data_val,
3714  *src[0],
3715  value);
3716 
3717  operation.post_constraints(value,
3718  values_dofs[comp][ind_local][v]);
3719  ind_local++;
3720  }
3721 
3722  AssertIndexRange(ind_local, dofs_per_component + 1);
3723 
3724  // get the dof values past the last constraint
3725  for (; ind_local < dofs_per_component;
3726  ++dof_indices[v], ++ind_local)
3727  {
3728  AssertIndexRange(*dof_indices[v], src[0]->size());
3729  operation.process_dof(*dof_indices[v],
3730  *src[0],
3731  values_dofs[comp][ind_local][v]);
3732  }
3733 
3734  if (apply_constraints == true && comp + 1 < n_components)
3735  next_index_indicators =
3736  dof_info.row_starts[cell_dof_index + comp + 2].second;
3737  }
3738  }
3739  }
3740 }
3741 
3742 
3743 
3744 template <int dim,
3745  int n_components_,
3746  typename Number,
3747  bool is_face,
3748  typename VectorizedArrayType>
3749 template <typename VectorType, typename VectorOperation>
3750 inline void
3753  const VectorOperation & operation,
3754  const std::array<VectorType *, n_components_> &src) const
3755 {
3756  Assert(!local_dof_indices.empty(), ExcNotInitialized());
3757 
3758  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3759  unsigned int index = this->first_selected_component * dofs_per_component;
3760  for (unsigned int comp = 0; comp < n_components; ++comp)
3761  {
3762  for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3763  {
3764  operation.process_empty(
3765  this->values_dofs[comp * dofs_per_component + i]);
3766  operation.process_dof_global(
3767  local_dof_indices[this->data->lexicographic_numbering[index]],
3768  *src[0],
3769  this->values_dofs[comp * dofs_per_component + i][0]);
3770  }
3771  }
3772 }
3773 
3774 
3775 
3776 template <int dim,
3777  int n_components_,
3778  typename Number,
3779  bool is_face,
3780  typename VectorizedArrayType>
3781 template <typename VectorType, typename VectorOperation>
3782 inline void
3785  const VectorOperation & operation,
3786  const std::array<VectorType *, n_components_> &src,
3787  const std::array<
3789  n_components_> & vectors_sm,
3790  const std::bitset<VectorizedArrayType::size()> &mask) const
3791 {
3792  // This functions processes the functions read_dof_values,
3793  // distribute_local_to_global, and set_dof_values with the same code for
3794  // contiguous cell indices (DG case). The distinction between these three
3795  // cases is made by the input VectorOperation that either reads values from
3796  // a vector and puts the data into the local data field or write local data
3797  // into the vector. Certain operations are no-ops for the given use case.
3798 
3799  std::integral_constant<bool,
3800  internal::is_vectorizable<VectorType, Number>::value>
3801  vector_selector;
3803  is_face ? this->dof_access_index :
3805  const unsigned int n_lanes = mask.count();
3806 
3807  const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3808  const std::vector<unsigned int> & dof_indices_cont =
3809  dof_info.dof_indices_contiguous[ind];
3810 
3811  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3812  std::array<VectorizedArrayType *, n_components> values_dofs;
3813  for (unsigned int c = 0; c < n_components; ++c)
3814  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3815  c * dofs_per_component;
3816 
3818 
3819  // Simple case: We have contiguous storage, so we can simply copy out the
3820  // data
3821  if ((dof_info.index_storage_variants[ind][this->cell] ==
3823  interleaved_contiguous &&
3824  n_lanes == VectorizedArrayType::size()) &&
3825  !(is_face &&
3826  this->dof_access_index ==
3828  this->is_interior_face() == false) &&
3829  !(!is_face && !this->is_interior_face()))
3830  {
3831  const unsigned int dof_index =
3832  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3833  this->dof_info
3834  ->component_dof_indices_offset[this->active_fe_index]
3835  [this->first_selected_component] *
3836  VectorizedArrayType::size();
3837  if (n_components == 1 || this->n_fe_components == 1)
3838  for (unsigned int comp = 0; comp < n_components; ++comp)
3839  operation.process_dofs_vectorized(dofs_per_component,
3840  dof_index,
3841  *src[comp],
3842  values_dofs[comp],
3843  vector_selector);
3844  else
3845  operation.process_dofs_vectorized(dofs_per_component * n_components,
3846  dof_index,
3847  *src[0],
3848  values_dofs[0],
3849  vector_selector);
3850  return;
3851  }
3852 
3853  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3854  this->get_cell_or_face_ids();
3855 
3856  // More general case: Must go through the components one by one and apply
3857  // some transformations
3858  const unsigned int n_filled_lanes =
3859  dof_info.n_vectorization_lanes_filled[ind][this->cell];
3860 
3861  const bool is_ecl =
3862  (this->dof_access_index ==
3864  this->is_interior_face() == false) ||
3865  (!is_face && !this->is_interior_face());
3866 
3867  if (vectors_sm[0] != nullptr)
3868  {
3869  const auto compute_vector_ptrs = [&](const unsigned int comp) {
3870  std::array<typename VectorType::value_type *,
3871  VectorizedArrayType::size()>
3872  vector_ptrs = {};
3873 
3874  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3875  {
3876  if (mask[v] == false)
3877  {
3878  vector_ptrs[v] = nullptr;
3879  continue;
3880  }
3881 
3883  ExcNotImplemented());
3884  Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3885  ExcIndexRange(ind,
3886  0,
3887  dof_info.dof_indices_contiguous_sm.size()));
3888  Assert(
3889  cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3890  ExcIndexRange(cells[v],
3891  0,
3892  dof_info.dof_indices_contiguous_sm[ind].size()));
3893 
3894  const auto &temp =
3895  dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3896 
3897  if (temp.first != numbers::invalid_unsigned_int)
3898  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3899  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3901  [this->active_fe_index][this->first_selected_component]);
3902  else
3903  vector_ptrs[v] = nullptr;
3904  }
3905  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3906  ++v)
3907  vector_ptrs[v] = nullptr;
3908 
3909  return vector_ptrs;
3910  };
3911 
3912  if (n_filled_lanes == VectorizedArrayType::size() &&
3913  n_lanes == VectorizedArrayType::size() && !is_ecl)
3914  {
3915  if (n_components == 1 || this->n_fe_components == 1)
3916  {
3917  for (unsigned int comp = 0; comp < n_components; ++comp)
3918  {
3919  auto vector_ptrs = compute_vector_ptrs(comp);
3920  operation.process_dofs_vectorized_transpose(
3921  dofs_per_component,
3922  vector_ptrs,
3923  values_dofs[comp],
3924  vector_selector);
3925  }
3926  }
3927  else
3928  {
3929  auto vector_ptrs = compute_vector_ptrs(0);
3930  operation.process_dofs_vectorized_transpose(dofs_per_component *
3931  n_components,
3932  vector_ptrs,
3933  &values_dofs[0][0],
3934  vector_selector);
3935  }
3936  }
3937  else
3938  for (unsigned int comp = 0; comp < n_components; ++comp)
3939  {
3940  auto vector_ptrs = compute_vector_ptrs(
3941  (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3942 
3943  for (unsigned int i = 0; i < dofs_per_component; ++i)
3944  operation.process_empty(values_dofs[comp][i]);
3945 
3946  if (n_components == 1 || this->n_fe_components == 1)
3947  {
3948  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3949  if (mask[v] == true)
3950  for (unsigned int i = 0; i < dofs_per_component; ++i)
3951  operation.process_dof(vector_ptrs[v][i],
3952  values_dofs[comp][i][v]);
3953  }
3954  else
3955  {
3956  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3957  if (mask[v] == true)
3958  for (unsigned int i = 0; i < dofs_per_component; ++i)
3959  operation.process_dof(
3960  vector_ptrs[v][i + comp * dofs_per_component],
3961  values_dofs[comp][i][v]);
3962  }
3963  }
3964  return;
3965  }
3966 
3967  unsigned int dof_indices[VectorizedArrayType::size()];
3968 
3969  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3970  {
3971  Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
3972  ExcNotImplemented());
3973  if (mask[v] == true)
3974  dof_indices[v] =
3975  dof_indices_cont[cells[v]] +
3976  this->dof_info
3977  ->component_dof_indices_offset[this->active_fe_index]
3978  [this->first_selected_component] *
3979  dof_info.dof_indices_interleave_strides[ind][cells[v]];
3980  else
3981  dof_indices[v] = numbers::invalid_unsigned_int;
3982  }
3983 
3984  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3985  dof_indices[v] = numbers::invalid_unsigned_int;
3986 
3987  // In the case with contiguous cell indices, we know that there are no
3988  // constraints and that the indices within each element are contiguous
3989  if (n_filled_lanes == VectorizedArrayType::size() &&
3990  n_lanes == VectorizedArrayType::size() && !is_ecl)
3991  {
3992  if (dof_info.index_storage_variants[ind][this->cell] ==
3994  contiguous)
3995  {
3996  if (n_components == 1 || this->n_fe_components == 1)
3997  for (unsigned int comp = 0; comp < n_components; ++comp)
3998  operation.process_dofs_vectorized_transpose(dofs_per_component,
3999  dof_indices,
4000  *src[comp],
4001  values_dofs[comp],
4002  vector_selector);
4003  else
4004  operation.process_dofs_vectorized_transpose(dofs_per_component *
4005  n_components,
4006  dof_indices,
4007  *src[0],
4008  &values_dofs[0][0],
4009  vector_selector);
4010  }
4011  else if (dof_info.index_storage_variants[ind][this->cell] ==
4013  interleaved_contiguous_strided)
4014  {
4015  if (n_components == 1 || this->n_fe_components == 1)
4016  for (unsigned int i = 0; i < dofs_per_component; ++i)
4017  {
4018  for (unsigned int comp = 0; comp < n_components; ++comp)
4019  operation.process_dof_gather(dof_indices,
4020  *src[comp],
4021  i * VectorizedArrayType::size(),
4022  values_dofs[comp][i],
4023  vector_selector);
4024  }
4025  else
4026  for (unsigned int comp = 0; comp < n_components; ++comp)
4027  for (unsigned int i = 0; i < dofs_per_component; ++i)
4028  {
4029  operation.process_dof_gather(dof_indices,
4030  *src[0],
4031  (comp * dofs_per_component + i) *
4032  VectorizedArrayType::size(),
4033  values_dofs[comp][i],
4034  vector_selector);
4035  }
4036  }
4037  else
4038  {
4039  Assert(dof_info.index_storage_variants[ind][this->cell] ==
4041  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4042  ExcNotImplemented());
4043  const unsigned int *offsets =
4045  [ind][VectorizedArrayType::size() * this->cell];
4046  if (n_components == 1 || this->n_fe_components == 1)
4047  for (unsigned int i = 0; i < dofs_per_component; ++i)
4048  {
4049  for (unsigned int comp = 0; comp < n_components; ++comp)
4050  operation.process_dof_gather(dof_indices,
4051  *src[comp],
4052  0,
4053  values_dofs[comp][i],
4054  vector_selector);
4056  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4057  dof_indices[v] += offsets[v];
4058  }
4059  else
4060  for (unsigned int comp = 0; comp < n_components; ++comp)
4061  for (unsigned int i = 0; i < dofs_per_component; ++i)
4062  {
4063  operation.process_dof_gather(dof_indices,
4064  *src[0],
4065  0,
4066  values_dofs[comp][i],
4067  vector_selector);
4069  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4070  dof_indices[v] += offsets[v];
4071  }
4072  }
4073  }
4074  else
4075  for (unsigned int comp = 0; comp < n_components; ++comp)
4076  {
4077  for (unsigned int i = 0; i < dofs_per_component; ++i)
4078  operation.process_empty(values_dofs[comp][i]);
4079  if (dof_info.index_storage_variants[ind][this->cell] ==
4081  contiguous)
4082  {
4083  if (n_components == 1 || this->n_fe_components == 1)
4084  {
4085  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4086  if (mask[v] == true)
4087  for (unsigned int i = 0; i < dofs_per_component; ++i)
4088  operation.process_dof(dof_indices[v] + i,
4089  *src[comp],
4090  values_dofs[comp][i][v]);
4091  }
4092  else
4093  {
4094  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4095  if (mask[v] == true)
4096  for (unsigned int i = 0; i < dofs_per_component; ++i)
4097  operation.process_dof(dof_indices[v] + i +
4098  comp * dofs_per_component,
4099  *src[0],
4100  values_dofs[comp][i][v]);
4101  }
4102  }
4103  else
4104  {
4105  const unsigned int *offsets =
4107  [ind][VectorizedArrayType::size() * this->cell];
4108  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4109  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4110  if (n_components == 1 || this->n_fe_components == 1)
4111  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4112  {
4113  if (mask[v] == true)
4114  for (unsigned int i = 0; i < dofs_per_component; ++i)
4115  operation.process_dof(dof_indices[v] + i * offsets[v],
4116  *src[comp],
4117  values_dofs[comp][i][v]);
4118  }
4119  else
4120  {
4121  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4122  if (mask[v] == true)
4123  for (unsigned int i = 0; i < dofs_per_component; ++i)
4124  operation.process_dof(dof_indices[v] +
4125  (i + comp * dofs_per_component) *
4126  offsets[v],
4127  *src[0],
4128  values_dofs[comp][i][v]);
4129  }
4130  }
4131  }
4132 }
4133 
4134 namespace internal
4135 {
4136  template <
4137  typename Number,
4138  typename VectorType,
4139  std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4140  decltype(std::declval<VectorType>().begin())
4141  get_beginning(VectorType &vec)
4142  {
4143  return vec.begin();
4144  }
4145 
4146  template <
4147  typename Number,
4148  typename VectorType,
4149  std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4150  typename VectorType::value_type *
4151  get_beginning(VectorType &)
4152  {
4153  return nullptr;
4154  }
4155 
4156  template <typename VectorType,
4157  std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4158  nullptr>
4159  const std::vector<ArrayView<const typename VectorType::value_type>> *
4160  get_shared_vector_data(VectorType * vec,
4161  const bool is_valid_mode_for_sm,
4162  const unsigned int active_fe_index,
4164  {
4165  // note: no hp is supported
4166  if (is_valid_mode_for_sm &&
4167  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4168  .size() > 0 &&
4169  active_fe_index == 0)
4170  return &vec->shared_vector_data();
4171  else
4172  return nullptr;
4173  }
4174 
4175  template <typename VectorType,
4176  std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4177  * = nullptr>
4178  const std::vector<ArrayView<const typename VectorType::value_type>> *
4179  get_shared_vector_data(VectorType *,
4180  const bool,
4181  const unsigned int,
4183  {
4184  return nullptr;
4185  }
4186 
4187  template <int n_components, typename VectorType>
4188  std::pair<
4189  std::array<typename internal::BlockVectorSelector<
4190  VectorType,
4191  IsBlockVector<VectorType>::value>::BaseVectorType *,
4192  n_components>,
4193  std::array<
4194  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4195  VectorType,
4196  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4197  n_components>>
4198  get_vector_data(VectorType & src,
4199  const unsigned int first_index,
4200  const bool is_valid_mode_for_sm,
4201  const unsigned int active_fe_index,
4203  {
4204  // select between block vectors and non-block vectors. Note that the number
4205  // of components is checked in the internal data
4206  std::pair<
4207  std::array<typename internal::BlockVectorSelector<
4208  VectorType,
4209  IsBlockVector<VectorType>::value>::BaseVectorType *,
4210  n_components>,
4211  std::array<
4212  const std::vector<
4213  ArrayView<const typename internal::BlockVectorSelector<
4214  VectorType,
4215  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4216  n_components>>
4217  src_data;
4218 
4219  for (unsigned int d = 0; d < n_components; ++d)
4220  src_data.first[d] = internal::BlockVectorSelector<
4221  VectorType,
4222  IsBlockVector<VectorType>::value>::get_vector_component(src,
4223  d +
4224  first_index);
4225 
4226  for (unsigned int d = 0; d < n_components; ++d)
4227  src_data.second[d] = get_shared_vector_data(
4228  const_cast<typename internal::BlockVectorSelector<
4229  typename std::remove_const<VectorType>::type,
4231  BaseVectorType *>(src_data.first[d]),
4232  is_valid_mode_for_sm,
4233  active_fe_index,
4234  dof_info);
4235 
4236  return src_data;
4237  }
4238 } // namespace internal
4239 
4240 
4241 
4242 template <int dim,
4243  int n_components_,
4244  typename Number,
4245  bool is_face,
4246  typename VectorizedArrayType>
4247 inline void
4250 {
4251  if (this->dof_info == nullptr ||
4252  this->dof_info->hanging_node_constraint_masks.size() == 0 ||
4253  this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4254  this->dof_info->hanging_node_constraint_masks_comp
4255  [this->active_fe_index][this->first_selected_component] == false)
4256  return; // nothing to do with faces
4257 
4258  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4259  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4260  constraint_mask;
4261 
4262  bool hn_available = false;
4263 
4264  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
4265  this->get_cell_ids();
4266 
4267  for (unsigned int v = 0; v < n_lanes; ++v)
4268  {
4269  if (cells[v] == numbers::invalid_unsigned_int)
4270  {
4271  constraint_mask[v] = internal::MatrixFreeFunctions::
4273  continue;
4274  }
4275 
4276  const unsigned int cell_index = cells[v];
4277  const auto mask =
4278  this->dof_info->hanging_node_constraint_masks[cell_index];
4279  constraint_mask[v] = mask;
4280 
4281  hn_available |= (mask != internal::MatrixFreeFunctions::
4283  }
4284 
4285  if (hn_available == false)
4286  return; // no hanging node on cell batch -> nothing to do
4287 
4289  apply(n_components,
4290  this->data->data.front().fe_degree,
4291  this->get_shape_info(),
4292  transpose,
4293  constraint_mask,
4294  this->values_dofs);
4295 }
4296 
4297 
4298 
4299 template <int dim,
4300  int n_components_,
4301  typename Number,
4302  bool is_face,
4303  typename VectorizedArrayType>
4304 template <typename VectorType>
4305 inline void
4307  read_dof_values(const VectorType & src,
4308  const unsigned int first_index,
4309  const std::bitset<VectorizedArrayType::size()> &mask)
4310 {
4311  const auto src_data = internal::get_vector_data<n_components_>(
4312  src,
4313  first_index,
4314  this->dof_access_index ==
4316  this->active_fe_index,
4317  this->dof_info);
4318 
4320  read_write_operation(reader, src_data.first, src_data.second, mask, true);
4321 
4322  apply_hanging_node_constraints(false);
4323 
4324 # ifdef DEBUG
4325  this->dof_values_initialized = true;
4326 # endif
4327 }
4328 
4329 
4330 
4331 template <int dim,
4332  int n_components_,
4333  typename Number,
4334  bool is_face,
4335  typename VectorizedArrayType>
4336 template <typename VectorType>
4337 inline void
4339  read_dof_values_plain(const VectorType & src,
4340  const unsigned int first_index,
4341  const std::bitset<VectorizedArrayType::size()> &mask)
4342 {
4343  const auto src_data = internal::get_vector_data<n_components_>(
4344  src,
4345  first_index,
4346  this->dof_access_index ==
4348  this->active_fe_index,
4349  this->dof_info);
4350 
4352  read_write_operation(reader, src_data.first, src_data.second, mask, false);
4353 
4354 # ifdef DEBUG
4355  this->dof_values_initialized = true;
4356 # endif
4357 }
4358 
4359 
4360 
4361 template <int dim,
4362  int n_components_,
4363  typename Number,
4364  bool is_face,
4365  typename VectorizedArrayType>
4366 template <typename VectorType>
4367 inline void
4370  VectorType & dst,
4371  const unsigned int first_index,
4372  const std::bitset<VectorizedArrayType::size()> &mask) const
4373 {
4374 # ifdef DEBUG
4375  Assert(this->dof_values_initialized == true,
4377 # endif
4378 
4379  apply_hanging_node_constraints(true);
4380 
4381  const auto dst_data = internal::get_vector_data<n_components_>(
4382  dst,
4383  first_index,
4384  this->dof_access_index ==
4386  this->active_fe_index,
4387  this->dof_info);
4388 
4390  distributor;
4391  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4392 }
4393 
4394 
4395 
4396 template <int dim,
4397  int n_components_,
4398  typename Number,
4399  bool is_face,
4400  typename VectorizedArrayType>
4401 template <typename VectorType>
4402 inline void
4404  set_dof_values(VectorType & dst,
4405  const unsigned int first_index,
4406  const std::bitset<VectorizedArrayType::size()> &mask) const
4407 {
4408 # ifdef DEBUG
4409  Assert(this->dof_values_initialized == true,
4411 # endif
4412 
4413  const auto dst_data = internal::get_vector_data<n_components_>(
4414  dst,
4415  first_index,
4416  this->dof_access_index ==
4418  this->active_fe_index,
4419  this->dof_info);
4420 
4422  read_write_operation(setter, dst_data.first, dst_data.second, mask);
4423 }
4424 
4425 
4426 
4427 template <int dim,
4428  int n_components_,
4429  typename Number,
4430  bool is_face,
4431  typename VectorizedArrayType>
4432 template <typename VectorType>
4433 inline void
4436  VectorType & dst,
4437  const unsigned int first_index,
4438  const std::bitset<VectorizedArrayType::size()> &mask) const
4439 {
4440 # ifdef DEBUG
4441  Assert(this->dof_values_initialized == true,
4443 # endif
4444 
4445  const auto dst_data = internal::get_vector_data<n_components_>(
4446  dst,
4447  first_index,
4448  this->dof_access_index ==
4450  this->active_fe_index,
4451  this->dof_info);
4452 
4454  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4455 }
4456 
4457 
4458 
4459 /*------------------------------ access to data fields ----------------------*/
4460 
4461 
4462 
4463 template <int dim,
4464  int n_components_,
4465  typename Number,
4466  bool is_face,
4467  typename VectorizedArrayType>
4470  get_dof_value(const unsigned int dof) const
4471 {
4472  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4473  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4475  for (unsigned int comp = 0; comp < n_components; ++comp)
4476  return_value[comp] = this->values_dofs[comp * dofs + dof];
4477  return return_value;
4478 }
4479 
4480 
4481 
4482 template <int dim,
4483  int n_components_,
4484  typename Number,
4485  bool is_face,
4486  typename VectorizedArrayType>
4489  get_value(const unsigned int q_point) const
4490 {
4491 # ifdef DEBUG
4492  Assert(this->values_quad_initialized == true,
4494 # endif
4495 
4496  AssertIndexRange(q_point, this->n_quadrature_points);
4497  const std::size_t nqp = this->n_quadrature_points;
4499  for (unsigned int comp = 0; comp < n_components; ++comp)
4500  return_value[comp] = this->values_quad[comp * nqp + q_point];
4501  return return_value;
4502 }
4503 
4504 
4505 
4506 template <int dim,
4507  int n_components_,
4508  typename Number,
4509  bool is_face,
4510  typename VectorizedArrayType>
4511 inline DEAL_II_ALWAYS_INLINE
4514  get_gradient(const unsigned int q_point) const
4515 {
4516 # ifdef DEBUG
4517  Assert(this->gradients_quad_initialized == true,
4519 # endif
4520 
4521  AssertIndexRange(q_point, this->n_quadrature_points);
4522  Assert(this->jacobian != nullptr,
4524  "update_gradients"));
4525  const std::size_t nqp = this->n_quadrature_points;
4527 
4528  // Cartesian cell
4529  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4530  {
4531  for (unsigned int d = 0; d < dim; ++d)
4532  for (unsigned int comp = 0; comp < n_components; ++comp)
4533  grad_out[comp][d] =
4534  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4535  this->jacobian[0][d][d];
4536  }
4537  // cell with general/affine Jacobian
4538  else
4539  {
4541  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4542  q_point :
4543  0];
4544  for (unsigned int comp = 0; comp < n_components; ++comp)
4545  for (unsigned int d = 0; d < dim; ++d)
4546  {
4547  grad_out[comp][d] =
4548  jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4549  for (unsigned int e = 1; e < dim; ++e)
4550  grad_out[comp][d] +=
4551  jac[d][e] *
4552  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4553  }
4554  }
4555  return grad_out;
4556 }
4557 
4558 
4559 
4560 template <int dim,
4561  int n_components_,
4562  typename Number,
4563  bool is_face,
4564  typename VectorizedArrayType>
4567  get_normal_derivative(const unsigned int q_point) const
4568 {
4569  AssertIndexRange(q_point, this->n_quadrature_points);
4570 # ifdef DEBUG
4571  Assert(this->gradients_quad_initialized == true,
4573 # endif
4574 
4575  Assert(this->normal_x_jacobian != nullptr,
4577  "update_gradients"));
4578 
4579  const std::size_t nqp = this->n_quadrature_points;
4581 
4582  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4583  for (unsigned int comp = 0; comp < n_components; ++comp)
4584  grad_out[comp] =
4585  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4586  (this->normal_x_jacobian[0][dim - 1]);
4587  else
4588  {
4589  const std::size_t index =
4590  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4591  for (unsigned int comp = 0; comp < n_components; ++comp)
4592  {
4593  grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4594  this->normal_x_jacobian[index][0];
4595  for (unsigned int d = 1; d < dim; ++d)
4596  grad_out[comp] +=
4597  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4598  this->normal_x_jacobian[index][d];
4599  }
4600  }
4601  return grad_out;
4602 }
4603 
4604 
4605 
4606 namespace internal
4607 {
4608  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4609  // store the lower diagonal because of symmetry
4610  template <typename VectorizedArrayType>
4611  inline void
4612  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4613  const VectorizedArrayType *const hessians,
4614  const unsigned int,
4615  VectorizedArrayType (&tmp)[1][1])
4616  {
4617  tmp[0][0] = jac[0][0] * hessians[0];
4618  }
4619 
4620  template <typename VectorizedArrayType>
4621  inline void
4622  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4623  const VectorizedArrayType *const hessians,
4624  const unsigned int nqp,
4625  VectorizedArrayType (&tmp)[2][2])
4626  {
4627  for (unsigned int d = 0; d < 2; ++d)
4628  {
4629  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4630  tmp[1][d] =
4631  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4632  }
4633  }
4634 
4635  template <typename VectorizedArrayType>
4636  inline void
4637  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4638  const VectorizedArrayType *const hessians,
4639  const unsigned int nqp,
4640  VectorizedArrayType (&tmp)[3][3])
4641  {
4642  for (unsigned int d = 0; d < 3; ++d)
4643  {
4644  tmp[0][d] =
4645  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4646  jac[d][2] * hessians[4 * nqp]);
4647  tmp[1][d] =
4648  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4649  jac[d][2] * hessians[5 * nqp]);
4650  tmp[2][d] =
4651  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4652  jac[d][2] * hessians[2 * nqp]);
4653  }
4654  }
4655 } // namespace internal
4656 
4657 
4658 
4659 template <int dim,
4660  int n_components_,
4661  typename Number,
4662  bool is_face,
4663  typename VectorizedArrayType>
4666  get_hessian(const unsigned int q_point) const
4667 {
4668 # ifdef DEBUG
4669  Assert(this->hessians_quad_initialized == true,
4671 # endif
4672  AssertIndexRange(q_point, this->n_quadrature_points);
4673 
4674  Assert(this->jacobian != nullptr,
4676  "update_hessian"));
4678  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4679  0 :
4680  q_point];
4681 
4683 
4684  const std::size_t nqp = this->n_quadrature_points;
4685  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4686 
4687  // Cartesian cell
4688  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4689  {
4690  for (unsigned int comp = 0; comp < n_components; ++comp)
4691  {
4692  for (unsigned int d = 0; d < dim; ++d)
4693  hessian_out[comp][d][d] =
4694  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4695  (jac[d][d] * jac[d][d]);
4696  switch (dim)
4697  {
4698  case 1:
4699  break;
4700  case 2:
4701  hessian_out[comp][0][1] =
4702  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4703  (jac[0][0] * jac[1][1]);
4704  break;
4705  case 3:
4706  hessian_out[comp][0][1] =
4707  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4708  (jac[0][0] * jac[1][1]);
4709  hessian_out[comp][0][2] =
4710  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4711  (jac[0][0] * jac[2][2]);
4712  hessian_out[comp][1][2] =
4713  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4714  (jac[1][1] * jac[2][2]);
4715  break;
4716  default:
4717  Assert(false, ExcNotImplemented());
4718  }
4719  for (unsigned int d = 0; d < dim; ++d)
4720  for (unsigned int e = d + 1; e < dim; ++e)
4721  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4722  }
4723  }
4724  // cell with general Jacobian, but constant within the cell
4725  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4726  {
4727  for (unsigned int comp = 0; comp < n_components; ++comp)
4728  {
4729  VectorizedArrayType tmp[dim][dim];
4730  internal::hessian_unit_times_jac(
4731  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4732 
4733  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4734  for (unsigned int d = 0; d < dim; ++d)
4735  for (unsigned int e = d; e < dim; ++e)
4736  {
4737  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4738  for (unsigned int f = 1; f < dim; ++f)
4739  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4740  }
4741 
4742  // no J' * grad(u) part here because the Jacobian is constant
4743  // throughout the cell and hence, its derivative is zero
4744 
4745  // take symmetric part
4746  for (unsigned int d = 0; d < dim; ++d)
4747  for (unsigned int e = d + 1; e < dim; ++e)
4748  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4749  }
4750  }
4751  // cell with general Jacobian
4752  else
4753  {
4754  const auto &jac_grad = this->jacobian_gradients[q_point];
4755  for (unsigned int comp = 0; comp < n_components; ++comp)
4756  {
4757  VectorizedArrayType tmp[dim][dim];
4758  internal::hessian_unit_times_jac(
4759  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4760 
4761  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4762  for (unsigned int d = 0; d < dim; ++d)
4763  for (unsigned int e = d; e < dim; ++e)
4764  {
4765  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4766  for (unsigned int f = 1; f < dim; ++f)
4767  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4768  }
4769 
4770  // add diagonal part of J' * grad(u)
4771  for (unsigned int d = 0; d < dim; ++d)
4772  for (unsigned int e = 0; e < dim; ++e)
4773  hessian_out[comp][d][d] +=
4774  jac_grad[d][e] *
4775  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4776 
4777  // add off-diagonal part of J' * grad(u)
4778  for (unsigned int d = 0, count = dim; d < dim; ++d)
4779  for (unsigned int e = d + 1; e < dim; ++e, ++count)
4780  for (unsigned int f = 0; f < dim; ++f)
4781  hessian_out[comp][d][e] +=
4782  jac_grad[count][f] *
4783  this->gradients_quad[(comp * dim + f) * nqp + q_point];
4784 
4785  // take symmetric part
4786  for (unsigned int d = 0; d < dim; ++d)
4787  for (unsigned int e = d + 1; e < dim; ++e)
4788  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4789  }
4790  }
4791  return hessian_out;
4792 }
4793 
4794 
4795 
4796 template <int dim,
4797  int n_components_,
4798  typename Number,
4799  bool is_face,
4800  typename VectorizedArrayType>
4803  get_hessian_diagonal(const unsigned int q_point) const
4804 {
4805  Assert(!is_face, ExcNotImplemented());
4806 # ifdef DEBUG
4807  Assert(this->hessians_quad_initialized == true,
4809 # endif
4810  AssertIndexRange(q_point, this->n_quadrature_points);
4811 
4812  Assert(this->jacobian != nullptr, ExcNotImplemented());
4814  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4815  0 :
4816  q_point];
4817 
4818  const std::size_t nqp = this->n_quadrature_points;
4819  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4821 
4822  // Cartesian cell
4823  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4824  {
4825  for (unsigned int comp = 0; comp < n_components; ++comp)
4826  for (unsigned int d = 0; d < dim; ++d)
4827  hessian_out[comp][d] =
4828  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4829  (jac[d][d] * jac[d][d]);
4830  }
4831  // cell with general Jacobian, but constant within the cell
4832  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4833  {
4834  for (unsigned int comp = 0; comp < n_components; ++comp)
4835  {
4836  // compute laplacian before the gradient because it needs to access
4837  // unscaled gradient data
4838  VectorizedArrayType tmp[dim][dim];
4839  internal::hessian_unit_times_jac(
4840  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4841 
4842  // compute only the trace part of hessian, J * tmp = J *
4843  // hess_unit(u) * J^T
4844  for (unsigned int d = 0; d < dim; ++d)
4845  {
4846  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4847  for (unsigned int f = 1; f < dim; ++f)
4848  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4849  }
4850  }
4851  }
4852  // cell with general Jacobian
4853  else
4854  {
4855  const auto &jac_grad = this->jacobian_gradients[q_point];
4856  for (unsigned int comp = 0; comp < n_components; ++comp)
4857  {
4858  // compute laplacian before the gradient because it needs to access
4859  // unscaled gradient data
4860  VectorizedArrayType tmp[dim][dim];
4861  internal::hessian_unit_times_jac(
4862  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4863 
4864  // compute only the trace part of hessian, J * tmp = J *
4865  // hess_unit(u) * J^T
4866  for (unsigned int d = 0; d < dim; ++d)
4867  {
4868  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4869  for (unsigned int f = 1; f < dim; ++f)
4870  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4871  }
4872 
4873  for (unsigned int d = 0; d < dim; ++d)
4874  for (unsigned int e = 0; e < dim; ++e)
4875  hessian_out[comp][d] +=
4876  jac_grad[d][e] *
4877  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4878  }
4879  }
4880  return hessian_out;
4881 }
4882 
4883 
4884 
4885 template <int dim,
4886  int n_components_,
4887  typename Number,
4888  bool is_face,
4889  typename VectorizedArrayType>
4892  get_laplacian(const unsigned int q_point) const
4893 {
4894  Assert(is_face == false, ExcNotImplemented());
4895 # ifdef DEBUG
4896  Assert(this->hessians_quad_initialized == true,
4898 # endif
4899  AssertIndexRange(q_point, this->n_quadrature_points);
4900 
4902  const auto hess_diag = get_hessian_diagonal(q_point);
4903  for (unsigned int comp = 0; comp < n_components; ++comp)
4904  {
4905  laplacian_out[comp] = hess_diag[comp][0];
4906  for (unsigned int d = 1; d < dim; ++d)
4907  laplacian_out[comp] += hess_diag[comp][d];
4908  }
4909  return laplacian_out;
4910 }
4911 
4912 
4913 
4914 template <int dim,
4915  int n_components_,
4916  typename Number,
4917  bool is_face,
4918  typename VectorizedArrayType>
4919 inline DEAL_II_ALWAYS_INLINE void
4922  const unsigned int dof)
4923 {
4924 # ifdef DEBUG
4925  this->dof_values_initialized = true;
4926 # endif
4927  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4928  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4929  for (unsigned int comp = 0; comp < n_components; ++comp)
4930  this->values_dofs[comp * dofs + dof] = val_in[comp];
4931 }
4932 
4933 
4934 
4935 template <int dim,
4936  int n_components_,
4937  typename Number,
4938  bool is_face,
4939  typename VectorizedArrayType>
4940 inline DEAL_II_ALWAYS_INLINE void
4943  const unsigned int q_point)
4944 {
4945 # ifdef DEBUG
4946  Assert(this->is_reinitialized, ExcNotInitialized());
4947 # endif
4948  AssertIndexRange(q_point, this->n_quadrature_points);
4949  Assert(this->J_value != nullptr,
4951  "update_values"));
4952 # ifdef DEBUG
4953  this->values_quad_submitted = true;
4954 # endif
4955 
4956  const std::size_t nqp = this->n_quadrature_points;
4957  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4958  {
4959  const VectorizedArrayType JxW =
4960  this->J_value[0] * this->quadrature_weights[q_point];
4961  for (unsigned int comp = 0; comp < n_components; ++comp)
4962  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4963  }
4964  else
4965  {
4966  const VectorizedArrayType JxW = this->J_value[q_point];
4967  for (unsigned int comp = 0; comp < n_components; ++comp)
4968  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4969  }
4970 }
4971 
4972 
4973 
4974 template <int dim,
4975  int n_components_,
4976  typename Number,
4977  bool is_face,
4978  typename VectorizedArrayType>
4979 inline DEAL_II_ALWAYS_INLINE void
4982  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
4983  const unsigned int q_point)
4984 {
4985 # ifdef DEBUG
4986  Assert(this->is_reinitialized, ExcNotInitialized());
4987 # endif
4988  AssertIndexRange(q_point, this->n_quadrature_points);
4989  Assert(this->J_value != nullptr,
4991  "update_gradients"));
4992  Assert(this->jacobian != nullptr,
4994  "update_gradients"));
4995 # ifdef DEBUG
4996  this->gradients_quad_submitted = true;
4997 # endif
4998 
4999  const std::size_t nqp = this->n_quadrature_points;
5000  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5001  {
5002  const VectorizedArrayType JxW =
5003  this->J_value[0] * this->quadrature_weights[q_point];
5004  std::array<VectorizedArrayType, dim> jac;
5005  for (unsigned int d = 0; d < dim; ++d)
5006  jac[d] = this->jacobian[0][d][d];
5007  for (unsigned int d = 0; d < dim; ++d)
5008  {
5009  const VectorizedArrayType factor = jac[d] * JxW;
5010  for (unsigned int comp = 0; comp < n_components; ++comp)
5011  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5012  grad_in[comp][d] * factor;
5013  }
5014  }
5015  else
5016  {
5018  this->cell_type > internal::MatrixFreeFunctions::affine ?
5019  this->jacobian[q_point] :
5020  this->jacobian[0];
5021  const VectorizedArrayType JxW =
5022  this->cell_type > internal::MatrixFreeFunctions::affine ?
5023  this->J_value[q_point] :
5024  this->J_value[0] * this->quadrature_weights[q_point];
5025  for (unsigned int comp = 0; comp < n_components; ++comp)
5026  for (unsigned int d = 0; d < dim; ++d)
5027  {
5028  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5029  for (unsigned int e = 1; e < dim; ++e)
5030  new_val += (jac[e][d] * grad_in[comp][e]);
5031  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5032  new_val * JxW;
5033  }
5034  }
5035 }
5036 
5037 
5038 
5039 template <int dim,
5040  int n_components_,
5041  typename Number,
5042  bool is_face,
5043  typename VectorizedArrayType>
5044 inline DEAL_II_ALWAYS_INLINE void
5048  const unsigned int q_point)
5049 {
5050  AssertIndexRange(q_point, this->n_quadrature_points);
5051  Assert(this->normal_x_jacobian != nullptr,
5053  "update_gradients"));
5054 # ifdef DEBUG
5055  this->gradients_quad_submitted = true;
5056 # endif
5057 
5058  const std::size_t nqp = this->n_quadrature_points;
5059  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5060  {
5061  const VectorizedArrayType JxW_jac = this->J_value[0] *
5062  this->quadrature_weights[q_point] *
5063  this->normal_x_jacobian[0][dim - 1];
5064  for (unsigned int comp = 0; comp < n_components; ++comp)
5065  {
5066  for (unsigned int d = 0; d < dim - 1; ++d)
5067  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5068  VectorizedArrayType();
5069  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5070  grad_in[comp] * JxW_jac;
5071  }
5072  }
5073  else
5074  {
5075  const unsigned int index =
5076  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5078  this->normal_x_jacobian[index];
5079  const VectorizedArrayType JxW =
5080  (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5081  this->J_value[index] * this->quadrature_weights[q_point] :
5082  this->J_value[index];
5083  for (unsigned int comp = 0; comp < n_components; ++comp)
5084  {
5085  for (unsigned int d = 0; d < dim; ++d)
5086  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5087  (grad_in[comp] * JxW) * jac[d];
5088  }
5089  }
5090 }
5091 
5092 
5093 
5094 template <int dim,
5095  int n_components_,
5096  typename Number,
5097  bool is_face,
5098  typename VectorizedArrayType>
5099 inline DEAL_II_ALWAYS_INLINE void
5102  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5103  hessian_in,
5104  const unsigned int q_point)
5105 {
5106 # ifdef DEBUG
5107  Assert(this->is_reinitialized, ExcNotInitialized());
5108 # endif
5109  AssertIndexRange(q_point, this->n_quadrature_points);
5110  Assert(this->J_value != nullptr,
5112  "update_hessians"));
5113  Assert(this->jacobian != nullptr,
5115  "update_hessians"));
5116 # ifdef DEBUG
5117  this->hessians_quad_submitted = true;
5118 # endif
5119 
5120  // compute hessian_unit = J^T * hessian_in(u) * J
5121  const std::size_t nqp = this->n_quadrature_points;
5122  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5123  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5124  {
5125  const VectorizedArrayType JxW =
5126  this->J_value[0] * this->quadrature_weights[q_point];
5127 
5128  // diagonal part
5129  for (unsigned int d = 0; d < dim; ++d)
5130  {
5131  const auto jac_d = this->jacobian[0][d][d];
5132  const VectorizedArrayType factor = jac_d * jac_d * JxW;
5133  for (unsigned int comp = 0; comp < n_components; ++comp)
5134  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5135  hessian_in[comp][d][d] * factor;
5136  }
5137 
5138  // off diagonal part
5139  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5140  for (unsigned int e = 0; e < d; ++e, ++off_dia)
5141  {
5142  const auto jac_d = this->jacobian[0][d][d];
5143  const auto jac_e = this->jacobian[0][e][e];
5144  const VectorizedArrayType factor = jac_d * jac_e * JxW;
5145  for (unsigned int comp = 0; comp < n_components; ++comp)
5146  this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5147  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5148  }
5149  }
5150  // cell with general Jacobian, but constant within the cell
5151  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5152  {
5153  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5154  const VectorizedArrayType JxW =
5155  this->J_value[0] * this->quadrature_weights[q_point];
5156  for (unsigned int comp = 0; comp < n_components; ++comp)
5157  {
5158  // 1. tmp = hessian_in(u) * J
5159  VectorizedArrayType tmp[dim][dim];
5160  for (unsigned int i = 0; i < dim; ++i)
5161  for (unsigned int j = 0; j < dim; ++j)
5162  {
5163  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5164  for (unsigned int k = 1; k < dim; ++k)
5165  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5166  }
5167 
5168  // 2. hessian_unit = J^T * tmp
5169  VectorizedArrayType tmp2[dim][dim];
5170  for (unsigned int i = 0; i < dim; ++i)
5171  for (unsigned int j = 0; j < dim; ++j)
5172  {
5173  tmp2[i][j] = jac[0][i] * tmp[0][j];
5174  for (unsigned int k = 1; k < dim; ++k)
5175  tmp2[i][j] += jac[k][i] * tmp[k][j];
5176  }
5177 
5178  // diagonal part
5179  for (unsigned int d = 0; d < dim; ++d)
5180  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5181  tmp2[d][d] * JxW;
5182 
5183  // off diagonal part
5184  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5185  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5186  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5187  (tmp2[d][e] + tmp2[e][d]) * JxW;
5188  }
5189  }
5190  else
5191  {
5192  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5193  const VectorizedArrayType JxW = this->J_value[q_point];
5194  const auto &jac_grad = this->jacobian_gradients[q_point];
5195  for (unsigned int comp = 0; comp < n_components; ++comp)
5196  {
5197  // 1. tmp = hessian_in(u) * J
5198  VectorizedArrayType tmp[dim][dim];
5199  for (unsigned int i = 0; i < dim; ++i)
5200  for (unsigned int j = 0; j < dim; ++j)
5201  {
5202  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5203  for (unsigned int k = 1; k < dim; ++k)
5204  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5205  }
5206 
5207  // 2. hessian_unit = J^T * tmp
5208  VectorizedArrayType tmp2[dim][dim];
5209  for (unsigned int i = 0; i < dim; ++i)
5210  for (unsigned int j = 0; j < dim; ++j)
5211  {
5212  tmp2[i][j] = jac[0][i] * tmp[0][j];
5213  for (unsigned int k = 1; k < dim; ++k)
5214  tmp2[i][j] += jac[k][i] * tmp[k][j];
5215  }
5216 
5217  // diagonal part
5218  for (unsigned int d = 0; d < dim; ++d)
5219  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5220  tmp2[d][d] * JxW;
5221 
5222  // off diagonal part
5223  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5224  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5225  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5226  (tmp2[d][e] + tmp2[e][d]) * JxW;
5227 
5228  // 3. gradient_unit = J' ** hessian_in
5229  for (unsigned int d = 0; d < dim; ++d)
5230  {
5231  VectorizedArrayType sum = 0;
5232  for (unsigned int e = 0; e < dim; ++e)
5233  sum += hessian_in[comp][e][e] * jac_grad[e][d];
5234  for (unsigned int e = 0, count = dim; e < dim; ++e)
5235  for (unsigned int f = e + 1; f < dim; ++f, ++count)
5236  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5237  jac_grad[count][d];
5238  this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
5239  q_point] = sum * JxW;
5240  }
5241  }
5242  }
5243 }
5244 
5245 
5246 
5247 template <int dim,
5248  int n_components_,
5249  typename Number,
5250  bool is_face,
5251  typename VectorizedArrayType>
5254  integrate_value() const
5255 {
5256 # ifdef DEBUG
5257  Assert(this->is_reinitialized, ExcNotInitialized());
5258  Assert(this->values_quad_submitted == true,
5260 # endif
5261 
5263  const std::size_t nqp = this->n_quadrature_points;
5264  for (unsigned int q = 0; q < nqp; ++q)
5265  for (unsigned int comp = 0; comp < n_components; ++comp)
5266  return_value[comp] += this->values_quad[comp * nqp + q];
5267  return (return_value);
5268 }
5269 
5270 
5271 
5272 /*----------------------- FEEvaluationAccess --------------------------------*/
5273 
5274 
5275 template <int dim,
5276  int n_components_,
5277  typename Number,
5278  bool is_face,
5279  typename VectorizedArrayType>
5280 inline FEEvaluationAccess<dim,
5281  n_components_,
5282  Number,
5283  is_face,
5284  VectorizedArrayType>::
5285  FEEvaluationAccess(
5287  const unsigned int dof_no,
5288  const unsigned int first_selected_component,
5289  const unsigned int quad_no,
5290  const unsigned int fe_degree,
5291  const unsigned int n_q_points,
5292  const bool is_interior_face,
5293  const unsigned int active_fe_index,
5294  const unsigned int active_quad_index,
5295  const unsigned int face_type)
5296  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5297  matrix_free,
5298  dof_no,
5299  first_selected_component,
5300  quad_no,
5301  fe_degree,
5302  n_q_points,
5303  is_interior_face,
5304  active_fe_index,
5305  active_quad_index,
5306  face_type)
5307 {}
5308 
5309 
5310 
5311 template <int dim,
5312  int n_components_,
5313  typename Number,
5314  bool is_face,
5315  typename VectorizedArrayType>
5316 inline FEEvaluationAccess<dim,
5317  n_components_,
5318  Number,
5319  is_face,
5320  VectorizedArrayType>::
5321  FEEvaluationAccess(
5322  const Mapping<dim> & mapping,
5323  const FiniteElement<dim> &fe,
5324  const Quadrature<1> & quadrature,
5325  const UpdateFlags update_flags,
5326  const unsigned int first_selected_component,
5328  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5329  mapping,
5330  fe,
5331  quadrature,
5332  update_flags,
5333  first_selected_component,
5334  other)
5335 {}
5336 
5337 
5338 
5339 template <int dim,
5340  int n_components_,
5341  typename Number,
5342  bool is_face,
5343  typename VectorizedArrayType>
5344 inline FEEvaluationAccess<dim,
5345  n_components_,
5346  Number,
5347  is_face,
5348  VectorizedArrayType>::
5349  FEEvaluationAccess(const FEEvaluationAccess<dim,
5350  n_components_,
5351  Number,
5352  is_face,
5353  VectorizedArrayType> &other)
5354  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5355  other)
5356 {}
5357 
5358 
5359 
5360 template <int dim,
5361  int n_components_,
5362  typename Number,
5363  bool is_face,
5364  typename VectorizedArrayType>
5365 inline FEEvaluationAccess<dim,
5366  n_components_,
5367  Number,
5368  is_face,
5369  VectorizedArrayType> &
5371 operator=(const FEEvaluationAccess<dim,
5372  n_components_,
5373  Number,
5374  is_face,
5375  VectorizedArrayType> &other)
5376 {
5377  this->FEEvaluationBase<dim,
5378  n_components_,
5379  Number,
5380  is_face,
5381  VectorizedArrayType>::operator=(other);
5382  return *this;
5383 }
5384 
5385 
5386 
5387 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5388 
5389 
5390 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5394  const unsigned int dof_no,
5395  const unsigned int first_selected_component,
5396  const unsigned int quad_no,
5397  const unsigned int fe_degree,
5398  const unsigned int n_q_points,
5399  const bool is_interior_face,
5400  const unsigned int active_fe_index,
5401  const unsigned int active_quad_index,
5402  const unsigned int face_type)
5403  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5404  matrix_free,
5405  dof_no,
5406  first_selected_component,
5407  quad_no,
5408  fe_degree,
5409  n_q_points,
5410  is_interior_face,
5411  active_fe_index,
5412  active_quad_index,
5413  face_type)
5414 {}
5415 
5416 
5417 
5418 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5421  const Mapping<dim> & mapping,
5422  const FiniteElement<dim> &fe,
5423  const Quadrature<1> & quadrature,
5424  const UpdateFlags update_flags,
5425  const unsigned int first_selected_component,
5427  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5428  mapping,
5429  fe,
5430  quadrature,
5431  update_flags,
5432  first_selected_component,
5433  other)
5434 {}
5435 
5436 
5437 
5438 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5442  &other)
5443  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5444 {}
5445 
5446 
5447 
5448 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5452 {
5453  this
5454  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5455  other);
5456  return *this;
5457 }
5458 
5459 
5460 
5461 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5462 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5464  const unsigned int dof) const
5465 {
5466  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5467  return this->values_dofs[dof];
5468 }
5469 
5470 
5471 
5472 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5473 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5475  const unsigned int q_point) const
5476 {
5477 # ifdef DEBUG
5478  Assert(this->values_quad_initialized == true,
5480 # endif
5481  AssertIndexRange(q_point, this->n_quadrature_points);
5482  return this->values_quad[q_point];
5483 }
5484 
5485 
5486 
5487 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5488 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5490  get_normal_derivative(const unsigned int q_point) const
5491 {
5492  return BaseClass::get_normal_derivative(q_point)[0];
5493 }
5494 
5495 
5496 
5497 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5500  const unsigned int q_point) const
5501 {
5502  // could use the base class gradient, but that involves too many expensive
5503  // initialization operations on tensors
5504 
5505 # ifdef DEBUG
5506  Assert(this->gradients_quad_initialized == true,
5508 # endif
5509  AssertIndexRange(q_point, this->n_quadrature_points);
5510 
5511  Assert(this->jacobian != nullptr,
5513  "update_gradients"));
5514 
5516 
5517  const std::size_t nqp = this->n_quadrature_points;
5518  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5519  {
5520  for (unsigned int d = 0; d < dim; ++d)
5521  grad_out[d] =
5522  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5523  }
5524  // cell with general/affine Jacobian
5525  else
5526  {
5528  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5529  q_point :
5530  0];
5531  for (unsigned int d = 0; d < dim; ++d)
5532  {
5533  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5534  for (unsigned int e = 1; e < dim; ++e)
5535  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5536  }
5537  }
5538  return grad_out;
5539 }
5540 
5541 
5542 
5543 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5546  const unsigned int q_point) const
5547 {
5548  return BaseClass::get_hessian(q_point)[0];
5549 }
5550 
5551 
5552 
5553 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5556  get_hessian_diagonal(const unsigned int q_point) const
5557 {
5558  return BaseClass::get_hessian_diagonal(q_point)[0];
5559 }
5560 
5561 
5562 
5563 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5564 inline VectorizedArrayType
5566  const unsigned int q_point) const
5567 {
5568  return BaseClass::get_laplacian(q_point)[0];
5569 }
5570 
5571 
5572 
5573 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5574 inline void DEAL_II_ALWAYS_INLINE
5576  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5577 {
5578 # ifdef DEBUG
5579  this->dof_values_initialized = true;
5580  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5581 # endif
5582  this->values_dofs[dof] = val_in;
5583 }
5584 
5585 
5586 
5587 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5588 inline void DEAL_II_ALWAYS_INLINE
5590  const VectorizedArrayType val_in,
5591  const unsigned int q_point)
5592 {
5593 # ifdef DEBUG
5594  Assert(this->is_reinitialized, ExcNotInitialized());
5595 # endif
5596  AssertIndexRange(q_point, this->n_quadrature_points);
5597  Assert(this->J_value != nullptr,
5599  "update_value"));
5600 # ifdef DEBUG
5601  this->values_quad_submitted = true;
5602 # endif
5603 
5604  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5605  {
5606  const VectorizedArrayType JxW =
5607  this->J_value[0] * this->quadrature_weights[q_point];
5608  this->values_quad[q_point] = val_in * JxW;
5609  }
5610  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5611  {
5612  this->values_quad[q_point] = val_in * this->J_value[q_point];
5613  }
5614 }
5615 
5616 
5617 
5618 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5619 inline DEAL_II_ALWAYS_INLINE void
5621  const Tensor<1, 1, VectorizedArrayType> val_in,
5622  const unsigned int q_point)
5623 {
5624  submit_value(val_in[0], q_point);
5625 }
5626 
5627 
5628 
5629 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5630 inline DEAL_II_ALWAYS_INLINE void
5632  submit_normal_derivative(const VectorizedArrayType grad_in,
5633  const unsigned int q_point)
5634 {
5636  grad[0] = grad_in;
5637  BaseClass::submit_normal_derivative(grad, q_point);
5638 }
5639 
5640 
5641 
5642 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5643 inline DEAL_II_ALWAYS_INLINE void
5646  const unsigned int q_point)
5647 {
5648 # ifdef DEBUG
5649  Assert(this->is_reinitialized, ExcNotInitialized());
5650 # endif
5651  AssertIndexRange(q_point, this->n_quadrature_points);
5652  Assert(this->J_value != nullptr,
5654  "update_gradients"));
5655  Assert(this->jacobian != nullptr,
5657  "update_gradients"));
5658 # ifdef DEBUG
5659  this->gradients_quad_submitted = true;
5660 # endif
5661 
5662  const std::size_t nqp = this->n_quadrature_points;
5663  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5664  {
5665  const VectorizedArrayType JxW =
5666  this->J_value[0] * this->quadrature_weights[q_point];
5667 
5668  // Make sure the compiler does not think 'jacobian' is aliased with
5669  // 'gradients_quad'
5670  std::array<VectorizedArrayType, dim> jac;
5671  for (unsigned int d = 0; d < dim; ++d)
5672  jac[d] = this->jacobian[0][d][d];
5673 
5674  for (unsigned int d = 0; d < dim; ++d)
5675  this->gradients_quad[d * nqp + q_point] = grad_in[d] * jac[d] * JxW;
5676  }
5677  // general/affine cell type
5678  else
5679  {
5681  this->cell_type > internal::MatrixFreeFunctions::affine ?
5682  this->jacobian[q_point] :
5683  this->jacobian[0];
5684  const VectorizedArrayType JxW =
5685  this->cell_type > internal::MatrixFreeFunctions::affine ?
5686  this->J_value[q_point] :
5687  this->J_value[0] * this->quadrature_weights[q_point];
5688  for (unsigned int d = 0; d < dim; ++d)
5689  {
5690  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5691  for (unsigned int e = 1; e < dim; ++e)
5692  new_val += jac[e][d] * grad_in[e];
5693  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
5694  }
5695  }
5696 }
5697 
5698 
5699 
5700 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5701 inline DEAL_II_ALWAYS_INLINE void
5704  const unsigned int q_point)
5705 {
5707  hessian[0] = hessian_in;
5708  BaseClass::submit_hessian(hessian, q_point);
5709 }
5710 
5711 
5712 
5713 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5714 inline VectorizedArrayType
5716  integrate_value() const
5717 {
5718  return BaseClass::integrate_value()[0];
5719 }
5720 
5721 
5722 
5723 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5724 
5725 
5726 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5730  const unsigned int dof_no,
5731  const unsigned int first_selected_component,
5732  const unsigned int quad_no,
5733  const unsigned int fe_degree,
5734  const unsigned int n_q_points,
5735  const bool is_interior_face,
5736  const unsigned int active_fe_index,
5737  const unsigned int active_quad_index,
5738  const unsigned int face_type)
5739  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5740  matrix_free,
5741  dof_no,
5742  first_selected_component,
5743  quad_no,
5744  fe_degree,
5745  n_q_points,
5746  is_interior_face,
5747  active_fe_index,
5748  active_quad_index,
5749  face_type)
5750 {}
5751 
5752 
5753 
5754 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5757  const Mapping<dim> & mapping,
5758  const FiniteElement<dim> &fe,
5759  const Quadrature<1> & quadrature,
5760  const UpdateFlags update_flags,
5761  const unsigned int first_selected_component,
5763  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5764  mapping,
5765  fe,
5766  quadrature,
5767  update_flags,
5768  first_selected_component,
5769  other)
5770 {}
5771 
5772 
5773 
5774 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5778  &other)
5779  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5780 {}
5781 
5782 
5783 
5784 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5788  &other)
5789 {
5791  operator=(other);
5792  return *this;
5793 }
5794 
5795 
5796 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5799  const unsigned int q_point) const
5800 {
5801  if (this->data->element_type ==
5803  {
5804  // Piola transform is required
5805 # ifdef DEBUG
5806  Assert(this->values_quad_initialized == true,
5808 # endif
5809 
5810  AssertIndexRange(q_point, this->n_quadrature_points);
5811  Assert(this->J_value != nullptr,
5813  "update_values"));
5814  const std::size_t nqp = this->n_quadrature_points;
5816 
5817  if (!is_face &&
5818  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5819  {
5820  // Cartesian cell
5821  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5822  const VectorizedArrayType inv_det =
5823  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5824  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5825  this->jacobian[0][2][2];
5826 
5827  // J * u * det(J^-1)
5828  for (unsigned int comp = 0; comp < n_components; ++comp)
5829  value_out[comp] = this->values_quad[comp * nqp + q_point] *
5830  jac[comp][comp] * inv_det;
5831  }
5832  else
5833  {
5834  // Affine or general cell
5835  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5836  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5837  this->jacobian[q_point] :
5838  this->jacobian[0];
5840  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5841  transpose(invert(inv_t_jac)) :
5842  this->jacobian[1];
5843 
5844  // Derivatives are reordered for faces. Need to take this into account
5845  const VectorizedArrayType inv_det =
5846  (is_face && dim == 2 && this->get_face_no() < 2) ?
5847  -determinant(inv_t_jac) :
5848  determinant(inv_t_jac);
5849  // J * u * det(J^-1)
5850  for (unsigned int comp = 0; comp < n_components; ++comp)
5851  {
5852  value_out[comp] =
5853  this->values_quad[q_point] * jac[comp][0] * inv_det;
5854  for (unsigned int e = 1; e < dim; ++e)
5855  value_out[comp] +=
5856  this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
5857  }
5858  }
5859  return value_out;
5860  }
5861  else
5862  {
5863  // No Piola needed
5864  return BaseClass::get_value(q_point);
5865  }
5866 }
5867 
5868 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5871  get_gradient(const unsigned int q_point) const
5872 {
5873  if (this->data->element_type ==
5875  {
5876  // Piola transform is required
5877 # ifdef DEBUG
5878  Assert(this->gradients_quad_initialized == true,
5880 # endif
5881 
5882  AssertIndexRange(q_point, this->n_quadrature_points);
5883  Assert(this->jacobian != nullptr,
5885  "update_gradients"));
5886  const std::size_t nqp = this->n_quadrature_points;
5888 
5889  if (!is_face &&
5890  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5891  {
5892  // Cartesian cell
5893  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5894  this->jacobian[0];
5895  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5896  const VectorizedArrayType inv_det =
5897  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5898  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5899  this->jacobian[0][2][2];
5900 
5901  // J * grad_quad * J^-1 * det(J^-1)
5902  for (unsigned int d = 0; d < dim; ++d)
5903  for (unsigned int comp = 0; comp < n_components; ++comp)
5904  grad_out[comp][d] =
5905  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
5906  inv_t_jac[d][d] * jac[comp][comp] * inv_det;
5907  }
5908  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5909  {
5910  // Affine cell
5911  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5912  this->jacobian[0];
5913  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5914 
5915  // Derivatives are reordered for faces. Need to take this into account
5916  const VectorizedArrayType inv_det =
5917  (is_face && dim == 2 && this->get_face_no() < 2) ?
5918  -determinant(inv_t_jac) :
5919  determinant(inv_t_jac);
5920 
5921  VectorizedArrayType tmp;
5922  // J * grad_quad * J^-1 * det(J^-1)
5923  for (unsigned int comp = 0; comp < n_components; ++comp)
5924  for (unsigned int d = 0; d < dim; ++d)
5925  {
5926  tmp = 0;
5927  for (unsigned int f = 0; f < dim; ++f)
5928  for (unsigned int e = 0; e < dim; ++e)
5929  tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
5930  this->gradients_quad[(f * dim + e) * nqp + q_point];
5931 
5932  grad_out[comp][d] = tmp;
5933  }
5934  }
5935  else
5936  {
5937  // General cell
5938 
5939  // This assert could be removed if we make sure that this is updated
5940  // even though update_hessians or update_jacobian_grads is not passed,
5941  // i.e make the necessary changes in
5942  // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
5943  Assert(this->jacobian_gradients_non_inverse != nullptr,
5945  "update_hessians"));
5946 
5947  const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
5948  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5949  this->jacobian[q_point];
5950  const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
5951 
5952  // Derivatives are reordered for faces. Need to take this into account
5953  const VectorizedArrayType inv_det =
5954  (is_face && dim == 2 && this->get_face_no() < 2) ?
5955  -determinant(inv_t_jac) :
5956  determinant(inv_t_jac);
5957 
5958  VectorizedArrayType tmp;
5959  // J * grad_quad * J^-1 * det(J^-1)
5960  for (unsigned int comp = 0; comp < n_components; ++comp)
5961  for (unsigned int d = 0; d < dim; ++d)
5962  {
5963  tmp = 0;
5964  for (unsigned int f = 0; f < dim; ++f)
5965  for (unsigned int e = 0; e < dim; ++e)
5966  tmp += t_jac[f][comp] * inv_t_jac[d][e] *
5967  this->gradients_quad[(f * dim + e) * nqp + q_point];
5968 
5969  grad_out[comp][d] = tmp * inv_det;
5970  }
5971 
5972  // Contribution from values
5973  {
5974  // Diagonal part of jac_grad
5975 
5976  // Add jac_grad * J^{-1} * values * det(J^{-1})
5977  // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
5978  for (unsigned int i = 0; i < dim; ++i)
5979  for (unsigned int j = 0; j < dim; ++j)
5980  {
5981  tmp = jac_grad[0][i] * inv_t_jac[j][0] *
5982  this->values_quad[q_point];
5983  for (unsigned int f = 1; f < dim; ++f)
5984  tmp += jac_grad[f][i] * inv_t_jac[j][f] *
5985  this->values_quad[f * nqp + q_point];
5986 
5987  grad_out[i][j] += tmp * inv_det;
5988  }
5989 
5990  for (unsigned int i = 0; i < dim; ++i)
5991  for (unsigned int j = 0; j < dim; ++j)
5992  {
5993  tmp = 0;
5994  for (unsigned int f = 0; f < dim; ++f)
5995  for (unsigned int n = 0; n < dim; ++n)
5996  for (unsigned int m = 0; m < dim; ++m)
5997  tmp += inv_t_jac[m][f] * jac_grad[f][m] *
5998  inv_t_jac[j][f] * t_jac[n][i] *
5999  this->values_quad[n * nqp + q_point];
6000  grad_out[i][j] -= tmp * inv_det;
6001  }
6002  }
6003 
6004  {
6005  // Off-diagonal part of jac_grad
6006 
6007  // Add jac_grad * J^{-1} * values * det(J^{-1})
6008  // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
6009  for (unsigned int i = 0; i < dim; ++i)
6010  for (unsigned int j = 0; j < dim; ++j)
6011  {
6012  tmp = 0;
6013  for (unsigned int r = 0, f = dim; r < dim; ++r)
6014  for (unsigned int k = r + 1; k < dim; ++k, ++f)
6015  {
6016  tmp += jac_grad[f][i] *
6017  (inv_t_jac[j][k] *
6018  this->values_quad[r * nqp + q_point] +
6019  inv_t_jac[j][r] *
6020  this->values_quad[k * nqp + q_point]);
6021  for (unsigned int n = 0; n < dim; ++n)
6022  for (unsigned int m = 0; m < dim; ++m)
6023  tmp -= jac_grad[f][m] * t_jac[n][i] *
6024  this->values_quad[n * nqp + q_point] *
6025  (inv_t_jac[m][k] * inv_t_jac[j][r] +
6026  inv_t_jac[m][r] * inv_t_jac[j][k]);
6027  }
6028  grad_out[i][j] += tmp * inv_det;
6029  }
6030  }
6031  }
6032  return grad_out;
6033  }
6034  else
6035  {
6036  return BaseClass::get_gradient(q_point);
6037  }
6038 }
6039 
6040 
6041 
6042 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6043 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6045  get_divergence(const unsigned int q_point) const
6046 {
6047 # ifdef DEBUG
6048  Assert(this->gradients_quad_initialized == true,
6050 # endif
6051  AssertIndexRange(q_point, this->n_quadrature_points);
6052  Assert(this->jacobian != nullptr,
6054  "update_gradients"));
6055 
6056  VectorizedArrayType divergence;
6057  const std::size_t nqp = this->n_quadrature_points;
6058 
6059  if (this->data->element_type ==
6061  {
6062  if (!is_face &&
6063  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6064  {
6065  // Cartesian cell
6066  const VectorizedArrayType inv_det =
6067  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
6068  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
6069  this->jacobian[0][2][2];
6070 
6071  // div * det(J^-1)
6072  divergence = this->gradients_quad[q_point] * inv_det;
6073  for (unsigned int d = 1; d < dim; ++d)
6074  divergence +=
6075  this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
6076  }
6077  else
6078  {
6079  // General cell
6080  // Derivatives are reordered for faces. Need to take this into account
6081  const VectorizedArrayType inv_det =
6082  determinant(
6083  this->jacobian[this->cell_type >
6085  q_point :
6086  0]) *
6087  Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6088 
6089  // div * det(J^-1)
6090  divergence = this->gradients_quad[q_point] * inv_det;
6091  for (unsigned int d = 1; d < dim; ++d)
6092  divergence +=
6093  this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
6094  }
6095  }
6096  else
6097  {
6098  if (!is_face &&
6099  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6100  {
6101  // Cartesian cell
6102  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6103  for (unsigned int d = 1; d < dim; ++d)
6104  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
6105  this->jacobian[0][d][d];
6106  }
6107  else
6108  {
6109  // cell with general/constant Jacobian
6111  this->cell_type == internal::MatrixFreeFunctions::general ?
6112  this->jacobian[q_point] :
6113  this->jacobian[0];
6114  divergence = jac[0][0] * this->gradients_quad[q_point];
6115  for (unsigned int e = 1; e < dim; ++e)
6116  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
6117  for (unsigned int d = 1; d < dim; ++d)
6118  for (unsigned int e = 0; e < dim; ++e)
6119  divergence +=
6120  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
6121  }
6122  }
6123  return divergence;
6124 }
6125 
6126 
6127 
6128 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6131  get_symmetric_gradient(const unsigned int q_point) const
6132 {
6133  // copy from generic function into dim-specialization function
6134  const auto grad = get_gradient(q_point);
6135  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6136  VectorizedArrayType half = Number(0.5);
6137  for (unsigned int d = 0; d < dim; ++d)
6138  symmetrized[d] = grad[d][d];
6139  switch (dim)
6140  {
6141  case 1:
6142  break;
6143  case 2:
6144  symmetrized[2] = grad[0][1] + grad[1][0];
6145  symmetrized[2] *= half;
6146  break;
6147  case 3:
6148  symmetrized[3] = grad[0][1] + grad[1][0];
6149  symmetrized[3] *= half;
6150  symmetrized[4] = grad[0][2] + grad[2][0];
6151  symmetrized[4] *= half;
6152  symmetrized[5] = grad[1][2] + grad[2][1];
6153  symmetrized[5] *= half;
6154  break;
6155  default:
6156  Assert(false, ExcNotImplemented());
6157  }
6159 }
6160 
6161 
6162 
6163 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6164 inline DEAL_II_ALWAYS_INLINE
6165  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6167  const unsigned int q_point) const
6168 {
6169  // copy from generic function into dim-specialization function
6170  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6171  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6172  switch (dim)
6173  {
6174  case 1:
6175  Assert(false,
6176  ExcMessage(
6177  "Computing the curl in 1d is not a useful operation"));
6178  break;
6179  case 2:
6180  curl[0] = grad[1][0] - grad[0][1];
6181  break;
6182  case 3:
6183  curl[0] = grad[2][1] - grad[1][2];
6184  curl[1] = grad[0][2] - grad[2][0];
6185  curl[2] = grad[1][0] - grad[0][1];
6186  break;
6187  default:
6188  Assert(false, ExcNotImplemented());
6189  }
6190  return curl;
6191 }
6192 
6193 
6194 
6195 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6198  get_hessian_diagonal(const unsigned int q_point) const
6199 {
6200  return BaseClass::get_hessian_diagonal(q_point);
6201 }
6202 
6203 
6204 
6205 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6208  const unsigned int q_point) const
6209 {
6210 # ifdef DEBUG
6211  Assert(this->hessians_quad_initialized == true,
6213 # endif
6214  AssertIndexRange(q_point, this->n_quadrature_points);
6215  return BaseClass::get_hessian(q_point);
6216 }
6217 
6218 
6219 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6220 inline DEAL_II_ALWAYS_INLINE void
6223  const unsigned int q_point)
6224 {
6225  if (this->data->element_type ==
6227  {
6228  // Piola transform is required
6229  AssertIndexRange(q_point, this->n_quadrature_points);
6230  Assert(this->J_value != nullptr,
6232  "update_value"));
6233 # ifdef DEBUG
6234  Assert(this->is_reinitialized, ExcNotInitialized());
6235  this->values_quad_submitted = true;
6236 # endif
6237 
6238  const std::size_t nqp = this->n_quadrature_points;
6239  if (!is_face &&
6240  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6241  {
6242  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6243  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6244 
6245  for (unsigned int comp = 0; comp < n_components; ++comp)
6246  this->values_quad[comp * nqp + q_point] =
6247  val_in[comp] * weight * jac[comp][comp];
6248  }
6249  else
6250  {
6251  // Affine or general cell
6252  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6253  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6254  this->jacobian[q_point] :
6255  this->jacobian[0];
6257  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6258  transpose(invert(inv_t_jac)) :
6259  this->jacobian[1];
6260 
6261  // Derivatives are reordered for faces. Need to take this into account
6262  // and 1/inv_det != J_value for faces
6263  const VectorizedArrayType fac =
6264  (!is_face) ?
6265  this->quadrature_weights[q_point] :
6266  (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6267  this->J_value[q_point] :
6268  this->J_value[0] * this->quadrature_weights[q_point]) *
6269  ((dim == 2 && this->get_face_no() < 2) ?
6270  -determinant(inv_t_jac) :
6271  determinant(inv_t_jac)));
6272 
6273  // J^T * u * factor
6274  for (unsigned int comp = 0; comp < n_components; ++comp)
6275  {
6276  this->values_quad[comp * nqp + q_point] =
6277  val_in[0] * jac[0][comp] * fac;
6278  for (unsigned int e = 1; e < dim; ++e)
6279  this->values_quad[comp * nqp + q_point] +=
6280  val_in[e] * jac[e][comp] * fac;
6281  }
6282  }
6283  }
6284  else
6285  {
6286  // No Piola transform
6287  BaseClass::submit_value(val_in, q_point);
6288  }
6289 }
6290 
6291 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6292 inline DEAL_II_ALWAYS_INLINE void
6295  const unsigned int q_point)
6296 {
6297  if (this->data->element_type ==
6299  {
6300  // Piola transform is required
6301 
6302 # ifdef DEBUG
6303  Assert(this->is_reinitialized, ExcNotInitialized());
6304 # endif
6305  AssertIndexRange(q_point, this->n_quadrature_points);
6306  Assert(this->J_value != nullptr,
6308  "update_gradients"));
6309  Assert(this->jacobian != nullptr,
6311  "update_gradients"));
6312 # ifdef DEBUG
6313  this->gradients_quad_submitted = true;
6314 # endif
6315 
6316  const std::size_t nqp = this->n_quadrature_points;
6317  if (!is_face &&
6318  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6319  {
6320  // Cartesian cell
6321  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6322  this->jacobian[0];
6323  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6324  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6325  for (unsigned int d = 0; d < dim; ++d)
6326  for (unsigned int comp = 0; comp < n_components; ++comp)
6327  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6328  grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
6329  }
6330  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6331  {
6332  // Affine cell
6333  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6334  this->jacobian[0];
6335  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6336 
6337  // Derivatives are reordered for faces. Need to take this into account
6338  // and 1/inv_det != J_value for faces
6339  const VectorizedArrayType fac =
6340  (!is_face) ? this->quadrature_weights[q_point] :
6341  this->J_value[0] * this->quadrature_weights[q_point] *
6342  ((dim == 2 && this->get_face_no() < 2) ?
6343  -determinant(inv_t_jac) :
6344  determinant(inv_t_jac));
6345 
6346  // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6347  for (unsigned int comp = 0; comp < n_components; ++comp)
6348  for (unsigned int d = 0; d < dim; ++d)
6349  {
6350  VectorizedArrayType tmp = 0;
6351  for (unsigned int f = 0; f < dim; ++f)
6352  for (unsigned int e = 0; e < dim; ++e)
6353  tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e];
6354 
6355  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6356  tmp * fac;
6357  }
6358  }
6359  else
6360  {
6361  // General cell
6362 
6363  const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
6364  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6365  this->jacobian[q_point];
6366  const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
6367 
6368  // Derivatives are reordered for faces. Need to take this into account
6369  // and 1/inv_det != J_value for faces
6370  const VectorizedArrayType fac =
6371  (!is_face) ?
6372  this->quadrature_weights[q_point] :
6373  this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6374  -determinant(inv_t_jac) :
6375  determinant(inv_t_jac));
6376 
6377  VectorizedArrayType tmp;
6378  // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6379  for (unsigned int comp = 0; comp < n_components; ++comp)
6380  for (unsigned int d = 0; d < dim; ++d)
6381  {
6382  tmp = 0;
6383  for (unsigned int f = 0; f < dim; ++f)
6384  for (unsigned int e = 0; e < dim; ++e)
6385  tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e];
6386 
6387  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6388  tmp * fac;
6389  }
6390 
6391  // Contribution from values
6392  {
6393  // Diagonal part of jac_grad
6394 
6395  // Add jac_grad * J^{-1} * values * factor
6396  // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
6397  for (unsigned int f = 0; f < dim; ++f)
6398  {
6399  tmp = 0;
6400  for (unsigned int i = 0; i < dim; ++i)
6401  for (unsigned int j = 0; j < dim; ++j)
6402  {
6403  tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
6404  for (unsigned int m = 0; m < dim; ++m)
6405  for (unsigned int k = 0; k < dim; ++k)
6406  tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
6407  inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
6408  }
6409  this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
6410  }
6411  }
6412 
6413  {
6414  // Off-diagonal part of jac_grad
6415 
6416  // Add jac_grad * J^{-1} * values * factor
6417  for (unsigned int r = 0, f = dim; r < dim; ++r)
6418  for (unsigned int k = r + 1; k < dim; ++k, ++f)
6419  {
6420  tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
6421  for (unsigned int j = 1; j < dim; ++j)
6422  tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
6423  for (unsigned int i = 1; i < dim; ++i)
6424  for (unsigned int j = 0; j < dim; ++j)
6425  tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
6426  this->values_from_gradients_quad[r * nqp + q_point] +=
6427  tmp * fac;
6428 
6429  tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
6430  for (unsigned int j = 1; j < dim; ++j)
6431  tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
6432  for (unsigned int i = 1; i < dim; ++i)
6433  for (unsigned int j = 0; j < dim; ++j)
6434  tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
6435  this->values_from_gradients_quad[k * nqp + q_point] +=
6436  tmp * fac;
6437  }
6438 
6439  // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
6440  for (unsigned int n = 0; n < dim; ++n)
6441  {
6442  tmp = 0;
6443  for (unsigned int r = 0, f = dim; r < dim; ++r)
6444  for (unsigned int k = r + 1; k < dim; ++k, ++f)
6445  for (unsigned int i = 0; i < dim; ++i)
6446  for (unsigned int j = 0; j < dim; ++j)
6447  for (unsigned int m = 0; m < dim; ++m)
6448  tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
6449  (inv_t_jac[m][k] * inv_t_jac[j][r] +
6450  inv_t_jac[m][r] * inv_t_jac[j][k]);
6451 
6452  this->values_from_gradients_quad[n * nqp + q_point] -=
6453  tmp * fac;
6454  }
6455  }
6456  }
6457  }
6458  else
6459  {
6460  BaseClass::submit_gradient(grad_in, q_point);
6461  }
6462 }
6463 
6464 
6465 
6466 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6467 inline DEAL_II_ALWAYS_INLINE void
6470  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6471  const unsigned int q_point)
6472 {
6473  if (this->data->element_type ==
6475  {
6476  // Piola transform is required
6477  const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
6479  submit_gradient(grad, q_point);
6480  }
6481  else
6482  {
6483  BaseClass::submit_gradient(grad_in, q_point);
6484  }
6485 }
6486 
6487 
6488 
6489 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6490 inline DEAL_II_ALWAYS_INLINE void
6492  submit_divergence(const VectorizedArrayType div_in,
6493  const unsigned int q_point)
6494 {
6495 # ifdef DEBUG
6496  Assert(this->is_reinitialized, ExcNotInitialized());
6497 # endif
6498  AssertIndexRange(q_point, this->n_quadrature_points);
6499  Assert(this->J_value != nullptr,
6501  "update_gradients"));
6502  Assert(this->jacobian != nullptr,
6504  "update_gradients"));
6505 # ifdef DEBUG
6506  this->gradients_quad_submitted = true;
6507 # endif
6508 
6509  const std::size_t nqp = this->n_quadrature_points;
6510  if (this->data->element_type ==
6512  {
6513  // General cell
6514 
6515  // Derivatives are reordered for faces. Need to take this into account
6516  // and 1/inv_det != J_value for faces
6517  const VectorizedArrayType fac =
6518  (!is_face) ?
6519  this->quadrature_weights[q_point] * div_in :
6520  (this->cell_type > internal::MatrixFreeFunctions::affine ?
6521  this->J_value[q_point] :
6522  this->J_value[0] * this->quadrature_weights[q_point]) *
6523  div_in *
6524  determinant(
6525  this->jacobian[this->cell_type >
6527  q_point :
6528  0]) *
6529  Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6530 
6531  for (unsigned int d = 0; d < dim; ++d)
6532  {
6533  this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
6534  for (unsigned int e = d + 1; e < dim; ++e)
6535  {
6536  this->gradients_quad[(dim * d + e) * nqp + q_point] =
6537  VectorizedArrayType();
6538  this->gradients_quad[(dim * e + d) * nqp + q_point] =
6539  VectorizedArrayType();
6540  }
6541  }
6542  this->divergence_is_requested = true;
6543  }
6544  else
6545  {
6546  if (!is_face &&
6547  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6548  {
6549  const VectorizedArrayType fac =
6550  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6551  for (unsigned int d = 0; d < dim; ++d)
6552  {
6553  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6554  (fac * this->jacobian[0][d][d]);
6555  for (unsigned int e = d + 1; e < dim; ++e)
6556  {
6557  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6558  VectorizedArrayType();
6559  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6560  VectorizedArrayType();
6561  }
6562  }
6563  }
6564  else
6565  {
6567  this->cell_type == internal::MatrixFreeFunctions::general ?
6568  this->jacobian[q_point] :
6569  this->jacobian[0];
6570  const VectorizedArrayType fac =
6571  (this->cell_type == internal::MatrixFreeFunctions::general ?
6572  this->J_value[q_point] :
6573  this->J_value[0] * this->quadrature_weights[q_point]) *
6574  div_in;
6575  for (unsigned int d = 0; d < dim; ++d)
6576  {
6577  for (unsigned int e = 0; e < dim; ++e)
6578  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6579  jac[d][e] * fac;
6580  }
6581  }
6582  }
6583 }
6584 
6585 
6586 
6587 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6588 inline DEAL_II_ALWAYS_INLINE void
6592  const unsigned int q_point)
6593 {
6594  AssertThrow(
6595  this->data->element_type !=
6597  ExcNotImplemented());
6598 
6599  // could have used base class operator, but that involves some overhead
6600  // which is inefficient. it is nice to have the symmetric tensor because
6601  // that saves some operations
6602 # ifdef DEBUG
6603  Assert(this->is_reinitialized, ExcNotInitialized());
6604 # endif
6605  AssertIndexRange(q_point, this->n_quadrature_points);
6606  Assert(this->J_value != nullptr,
6608  "update_gradients"));
6609  Assert(this->jacobian != nullptr,
6611  "update_gradients"));
6612 # ifdef DEBUG
6613  this->gradients_quad_submitted = true;
6614 # endif
6615 
6616  const std::size_t nqp = this->n_quadrature_points;
6617  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6618  {
6619  const VectorizedArrayType JxW =
6620  this->J_value[0] * this->quadrature_weights[q_point];
6621  for (unsigned int d = 0; d < dim; ++d)
6622  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6623  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6624  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6625  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6626  {
6627  const VectorizedArrayType value =
6628  sym_grad.access_raw_entry(counter) * JxW;
6629  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6630  value * this->jacobian[0][d][d];
6631  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6632  value * this->jacobian[0][e][e];
6633  }
6634  }
6635  // general/affine cell type
6636  else
6637  {
6638  const VectorizedArrayType JxW =
6639  this->cell_type == internal::MatrixFreeFunctions::general ?
6640  this->J_value[q_point] :
6641  this->J_value[0] * this->quadrature_weights[q_point];
6643  this->cell_type == internal::MatrixFreeFunctions::general ?
6644  this->jacobian[q_point] :
6645  this->jacobian[0];
6646  VectorizedArrayType weighted[dim][dim];
6647  for (unsigned int i = 0; i < dim; ++i)
6648  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6649  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6650  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6651  {
6652  const VectorizedArrayType value =
6653  sym_grad.access_raw_entry(counter) * JxW;
6654  weighted[i][j] = value;
6655  weighted[j][i] = value;
6656  }
6657  for (unsigned int comp = 0; comp < dim; ++comp)
6658  for (unsigned int d = 0; d < dim; ++d)
6659  {
6660  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6661  for (unsigned int e = 1; e < dim; ++e)
6662  new_val += jac[e][d] * weighted[comp][e];
6663  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6664  }
6665  }
6666 }
6667 
6668 
6669 
6670 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6671 inline DEAL_II_ALWAYS_INLINE void
6674  const unsigned int q_point)
6675 {
6677  switch (dim)
6678  {
6679  case 1:
6680  Assert(false,
6681  ExcMessage(
6682  "Testing by the curl in 1d is not a useful operation"));
6683  break;
6684  case 2:
6685  grad[1][0] = curl[0];
6686  grad[0][1] = -curl[0];
6687  break;
6688  case 3:
6689  grad[2][1] = curl[0];
6690  grad[1][2] = -curl[0];
6691  grad[0][2] = curl[1];
6692  grad[2][0] = -curl[1];
6693  grad[1][0] = curl[2];
6694  grad[0][1] = -curl[2];
6695  break;
6696  default:
6697  Assert(false, ExcNotImplemented());
6698  }
6699  submit_gradient(grad, q_point);
6700 }
6701 
6702 
6703 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6704 
6705 
6706 template <typename Number, bool is_face, typename VectorizedArrayType>
6709  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
6710  const unsigned int dof_no,
6711  const unsigned int first_selected_component,
6712  const unsigned int quad_no,
6713  const unsigned int fe_degree,
6714  const unsigned int n_q_points,
6715  const bool is_interior_face,
6716  const unsigned int active_fe_index,
6717  const unsigned int active_quad_index,
6718  const unsigned int face_type)
6719  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6720  matrix_free,
6721  dof_no,
6722  first_selected_component,
6723  quad_no,
6724  fe_degree,
6725  n_q_points,
6726  is_interior_face,
6727  active_fe_index,
6728  active_quad_index,
6729  face_type)
6730 {}
6731 
6732 
6733 
6734 template <typename Number, bool is_face, typename VectorizedArrayType>
6737  const Mapping<1> & mapping,
6738  const FiniteElement<1> &fe,
6739  const Quadrature<1> & quadrature,
6740  const UpdateFlags update_flags,
6741  const unsigned int first_selected_component,
6743  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6744  mapping,
6745  fe,
6746  quadrature,
6747  update_flags,
6748  first_selected_component,
6749  other)
6750 {}
6751 
6752 
6753 
6754 template <typename Number, bool is_face, typename VectorizedArrayType>
6758  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6759 {}
6760 
6761 
6762 
6763 template <typename Number, bool is_face, typename VectorizedArrayType>
6767 {
6769  other);
6770  return *this;
6771 }
6772 
6773 
6774 
6775 template <typename Number, bool is_face, typename VectorizedArrayType>
6776 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6778  const unsigned int dof) const
6779 {
6780  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6781  return this->values_dofs[dof];
6782 }
6783 
6784 
6785 
6786 template <typename Number, bool is_face, typename VectorizedArrayType>
6787 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6789  const unsigned int q_point) const
6790 {
6791 # ifdef DEBUG
6792  Assert(this->values_quad_initialized == true,
6794 # endif
6795  AssertIndexRange(q_point, this->n_quadrature_points);
6796  return this->values_quad[q_point];
6797 }
6798 
6799 
6800 
6801 template <typename Number, bool is_face, typename VectorizedArrayType>
6804  const unsigned int q_point) const
6805 {
6806  // could use the base class gradient, but that involves too many inefficient
6807  // initialization operations on tensors
6808 
6809 # ifdef DEBUG
6810  Assert(this->gradients_quad_initialized == true,
6812 # endif
6813  AssertIndexRange(q_point, this->n_quadrature_points);
6814 
6816  this->cell_type == internal::MatrixFreeFunctions::general ?
6817  this->jacobian[q_point] :
6818  this->jacobian[0];
6819 
6821  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6822 
6823  return grad_out;
6824 }
6825 
6826 
6827 
6828 template <typename Number, bool is_face, typename VectorizedArrayType>
6829 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6831  const unsigned int q_point) const
6832 {
6833  return get_gradient(q_point)[0];
6834 }
6835 
6836 
6837 
6838 template <typename Number, bool is_face, typename VectorizedArrayType>
6839 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6841  get_normal_derivative(const unsigned int q_point) const
6842 {
6843  return BaseClass::get_normal_derivative(q_point)[0];
6844 }
6845 
6846 
6847 
6848 template <typename Number, bool is_face, typename VectorizedArrayType>
6851  const unsigned int q_point) const
6852 {
6853  return BaseClass::get_hessian(q_point)[0];
6854 }
6855 
6856 
6857 
6858 template <typename Number, bool is_face, typename VectorizedArrayType>
6861  get_hessian_diagonal(const unsigned int q_point) const
6862 {
6863  return BaseClass::get_hessian_diagonal(q_point)[0];
6864 }
6865 
6866 
6867 
6868 template <typename Number, bool is_face, typename VectorizedArrayType>
6869 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6871  const unsigned int q_point) const
6872 {
6873  return BaseClass::get_laplacian(q_point)[0];
6874 }
6875 
6876 
6877 
6878 template <typename Number, bool is_face, typename VectorizedArrayType>
6881  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6882 {
6883 # ifdef DEBUG
6884  this->dof_values_initialized = true;
6885  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6886 # endif
6887  this->values_dofs[dof] = val_in;
6888 }
6889 
6890 
6891 
6892 template <typename Number, bool is_face, typename VectorizedArrayType>
6893 inline DEAL_II_ALWAYS_INLINE void
6895  const VectorizedArrayType val_in,
6896  const unsigned int q_point)
6897 {
6898 # ifdef DEBUG
6899  Assert(this->is_reinitialized, ExcNotInitialized());
6900 # endif
6901  AssertIndexRange(q_point, this->n_quadrature_points);
6902 # ifdef DEBUG
6903  this->values_quad_submitted = true;
6904 # endif
6905 
6906  if (this->cell_type == internal::MatrixFreeFunctions::general)
6907  {
6908  const VectorizedArrayType JxW = this->J_value[q_point];
6909  this->values_quad[q_point] = val_in * JxW;
6910  }
6911  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6912  {
6913  const VectorizedArrayType JxW =
6914  this->J_value[0] * this->quadrature_weights[q_point];
6915  this->values_quad[q_point] = val_in * JxW;
6916  }
6917 }
6918 
6919 
6920 
6921 template <typename Number, bool is_face, typename VectorizedArrayType>
6922 inline DEAL_II_ALWAYS_INLINE void
6924  const Tensor<1, 1, VectorizedArrayType> val_in,
6925  const unsigned int q_point)
6926 {
6927  submit_value(val_in[0], q_point);
6928 }
6929 
6930 
6931 
6932 template <typename Number, bool is_face, typename VectorizedArrayType>
6933 inline DEAL_II_ALWAYS_INLINE void
6935  const Tensor<1, 1, VectorizedArrayType> grad_in,
6936  const unsigned int q_point)
6937 {
6938  submit_gradient(grad_in[0], q_point);
6939 }
6940 
6941 
6942 
6943 template <typename Number, bool is_face, typename VectorizedArrayType>
6944 inline DEAL_II_ALWAYS_INLINE void
6946  const VectorizedArrayType grad_in,
6947  const unsigned int q_point)
6948 {
6949 # ifdef DEBUG
6950  Assert(this->is_reinitialized, ExcNotInitialized());
6951 # endif
6952  AssertIndexRange(q_point, this->n_quadrature_points);
6953 # ifdef DEBUG
6954  this->gradients_quad_submitted = true;
6955 # endif
6956 
6958  this->cell_type == internal::MatrixFreeFunctions::general ?
6959  this->jacobian[q_point] :
6960  this->jacobian[0];
6961  const VectorizedArrayType JxW =
6962  this->cell_type == internal::MatrixFreeFunctions::general ?
6963  this->J_value[q_point] :
6964  this->J_value[0] * this->quadrature_weights[q_point];
6965 
6966  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6967 }
6968 
6969 
6970 
6971 template <typename Number, bool is_face, typename VectorizedArrayType>
6972 inline DEAL_II_ALWAYS_INLINE void
6974  const Tensor<2, 1, VectorizedArrayType> grad_in,
6975  const unsigned int q_point)
6976 {
6977  submit_gradient(grad_in[0][0], q_point);
6978 }
6979 
6980 
6981 
6982 template <typename Number, bool is_face, typename VectorizedArrayType>
6983 inline DEAL_II_ALWAYS_INLINE void
6985  submit_normal_derivative(const VectorizedArrayType grad_in,
6986  const unsigned int q_point)
6987 {
6989  grad[0] = grad_in;
6990  BaseClass::submit_normal_derivative(grad, q_point);
6991 }
6992 
6993 
6994 
6995 template <typename Number, bool is_face, typename VectorizedArrayType>
6996 inline DEAL_II_ALWAYS_INLINE void
6999  const unsigned int q_point)
7000 {
7001  BaseClass::submit_normal_derivative(grad_in, q_point);
7002 }
7003 
7004 
7005 template <typename Number, bool is_face, typename VectorizedArrayType>
7006 inline DEAL_II_ALWAYS_INLINE void
7008  const Tensor<2, 1, VectorizedArrayType> hessian_in,
7009  const unsigned int q_point)
7010 {
7012  hessian[0] = hessian_in;
7013  BaseClass::submit_hessian(hessian, q_point);
7014 }
7015 
7016 
7017 template <typename Number, bool is_face, typename VectorizedArrayType>
7018 inline VectorizedArrayType
7020  integrate_value() const
7021 {
7022  return BaseClass::integrate_value()[0];
7023 }
7024 
7025 
7026 
7027 /*-------------------------- FEEvaluation -----------------------------------*/
7028 
7029 
7030 template <int dim,
7031  int fe_degree,
7032  int n_q_points_1d,
7033  int n_components_,
7034  typename Number,
7035  typename VectorizedArrayType>
7036 inline FEEvaluation<dim,
7037  fe_degree,
7038  n_q_points_1d,
7039  n_components_,
7040  Number,
7041  VectorizedArrayType>::
7042  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7043  const unsigned int fe_no,
7044  const unsigned int quad_no,
7045  const unsigned int first_selected_component,
7046  const unsigned int active_fe_index,
7047  const unsigned int active_quad_index)
7048  : BaseClass(matrix_free,
7049  fe_no,
7050  first_selected_component,
7051  quad_no,
7052  fe_degree,
7053  static_n_q_points,
7054  true /*note: this is not a face*/,
7055  active_fe_index,
7056  active_quad_index)
7057  , dofs_per_component(this->data->dofs_per_component_on_cell)
7058  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7059  , n_q_points(this->data->n_q_points)
7060 {
7061  check_template_arguments(fe_no, 0);
7062 }
7063 
7064 
7065 
7066 template <int dim,
7067  int fe_degree,
7068  int n_q_points_1d,
7069  int n_components_,
7070  typename Number,
7071  typename VectorizedArrayType>
7072 inline FEEvaluation<dim,
7073  fe_degree,
7074  n_q_points_1d,
7075  n_components_,
7076  Number,
7077  VectorizedArrayType>::
7078  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7079  const std::pair<unsigned int, unsigned int> & range,
7080  const unsigned int dof_no,
7081  const unsigned int quad_no,
7082  const unsigned int first_selected_component)
7083  : FEEvaluation(matrix_free,
7084  dof_no,
7085  quad_no,
7086  first_selected_component,
7087  matrix_free.get_cell_active_fe_index(range))
7088 {}
7089 
7090 
7091 
7092 template <int dim,
7093  int fe_degree,
7094  int n_q_points_1d,
7095  int n_components_,
7096  typename Number,
7097  typename VectorizedArrayType>
7098 inline FEEvaluation<dim,
7099  fe_degree,
7100  n_q_points_1d,
7101  n_components_,
7102  Number,
7103  VectorizedArrayType>::
7104  FEEvaluation(const Mapping<dim> & mapping,
7105  const FiniteElement<dim> &fe,
7106  const Quadrature<1> & quadrature,
7107  const UpdateFlags update_flags,
7108  const unsigned int first_selected_component)
7109  : BaseClass(mapping,
7110  fe,
7111  quadrature,
7112  update_flags,
7113  first_selected_component,
7114  nullptr)
7115  , dofs_per_component(this->data->dofs_per_component_on_cell)
7116  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7117  , n_q_points(this->data->n_q_points)
7118 {
7119  check_template_arguments(numbers::invalid_unsigned_int, 0);
7120 }
7121 
7122 
7123 
7124 template <int dim,
7125  int fe_degree,
7126  int n_q_points_1d,
7127  int n_components_,
7128  typename Number,
7129  typename VectorizedArrayType>
7130 inline FEEvaluation<dim,
7131  fe_degree,
7132  n_q_points_1d,
7133  n_components_,
7134  Number,
7135  VectorizedArrayType>::
7136  FEEvaluation(const FiniteElement<dim> &fe,
7137  const Quadrature<1> & quadrature,
7138  const UpdateFlags update_flags,
7139  const unsigned int first_selected_component)
7140  : BaseClass(StaticMappingQ1<dim>::mapping,
7141  fe,
7142  quadrature,
7143  update_flags,
7144  first_selected_component,
7145  nullptr)
7146  , dofs_per_component(this->data->dofs_per_component_on_cell)
7147  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7148  , n_q_points(this->data->n_q_points)
7149 {
7150  check_template_arguments(numbers::invalid_unsigned_int, 0);
7151 }
7152 
7153 
7154 
7155 template <int dim,
7156  int fe_degree,
7157  int n_q_points_1d,
7158  int n_components_,
7159  typename Number,
7160  typename VectorizedArrayType>
7161 inline FEEvaluation<dim,
7162  fe_degree,
7163  n_q_points_1d,
7164  n_components_,
7165  Number,
7166  VectorizedArrayType>::
7167  FEEvaluation(const FiniteElement<dim> & fe,
7169  const unsigned int first_selected_component)
7170  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7171  fe,
7172  other.mapped_geometry->get_quadrature(),
7173  other.mapped_geometry->get_fe_values().get_update_flags(),
7174  first_selected_component,
7175  &other)
7176  , dofs_per_component(this->data->dofs_per_component_on_cell)
7177  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7178  , n_q_points(this->data->n_q_points)
7179 {
7180  check_template_arguments(numbers::invalid_unsigned_int, 0);
7181 }
7182 
7183 
7184 
7185 template <int dim,
7186  int fe_degree,
7187  int n_q_points_1d,
7188  int n_components_,
7189  typename Number,
7190  typename VectorizedArrayType>
7191 inline FEEvaluation<dim,
7192  fe_degree,
7193  n_q_points_1d,
7194  n_components_,
7195  Number,
7196  VectorizedArrayType>::FEEvaluation(const FEEvaluation
7197  &other)
7198  : BaseClass(other)
7199  , dofs_per_component(this->data->dofs_per_component_on_cell)
7200  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7201  , n_q_points(this->data->n_q_points)
7202 {
7203  check_template_arguments(numbers::invalid_unsigned_int, 0);
7204 }
7205 
7206 
7207 
7208 template <int dim,
7209  int fe_degree,
7210  int n_q_points_1d,
7211  int n_components_,
7212  typename Number,
7213  typename VectorizedArrayType>
7214 inline FEEvaluation<dim,
7215  fe_degree,
7216  n_q_points_1d,
7217  n_components_,
7218  Number,
7219  VectorizedArrayType> &
7220 FEEvaluation<dim,
7221  fe_degree,
7222  n_q_points_1d,
7223  n_components_,
7224  Number,
7225  VectorizedArrayType>::operator=(const FEEvaluation &other)
7226 {
7227  BaseClass::operator=(other);
7228  check_template_arguments(numbers::invalid_unsigned_int, 0);
7229  return *this;
7230 }
7231 
7232 
7233 
7234 template <int dim,
7235  int fe_degree,
7236  int n_q_points_1d,
7237  int n_components_,
7238  typename Number,
7239  typename VectorizedArrayType>
7240 inline void
7241 FEEvaluation<dim,
7242  fe_degree,
7243  n_q_points_1d,
7244  n_components_,
7245  Number,
7246  VectorizedArrayType>::
7247  check_template_arguments(const unsigned int dof_no,
7248  const unsigned int first_selected_component)
7249 {
7250  (void)dof_no;
7251  (void)first_selected_component;
7252 
7253  Assert(
7254  this->data->dofs_per_component_on_cell > 0,
7255  ExcMessage(
7256  "There is nothing useful you can do with an FEEvaluation object with "
7257  "FE_Nothing, i.e., without DoFs! If you have passed to "
7258  "MatrixFree::reinit() a collection of finite elements also containing "
7259  "FE_Nothing, please check - before creating FEEvaluation - the category "
7260  "of the current range by calling either "
7261  "MatrixFree::get_cell_range_category(range) or "
7262  "MatrixFree::get_face_range_category(range). The returned category "
7263  "is the index of the active FE, which you can use to exclude "
7264  "FE_Nothing."));
7265 
7266 # ifdef DEBUG
7267  // print error message when the dimensions do not match. Propose a possible
7268  // fix
7269  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7270  static_cast<unsigned int>(fe_degree) !=
7271  this->data->data.front().fe_degree) ||
7272  n_q_points != this->n_quadrature_points)
7273  {
7274  std::string message =
7275  "-------------------------------------------------------\n";
7276  message += "Illegal arguments in constructor/wrong template arguments!\n";
7277  message += " Called --> FEEvaluation<dim,";
7278  message += Utilities::int_to_string(fe_degree) + ",";
7279  message += Utilities::int_to_string(n_q_points_1d);
7280  message += "," + Utilities::int_to_string(n_components);
7281  message += ",Number>(data";
7282  if (first_selected_component != numbers::invalid_unsigned_int)
7283  {
7284  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7285  message += Utilities::int_to_string(this->quad_no) + ", ";
7286  message += Utilities::int_to_string(first_selected_component);
7287  }
7288  message += ")\n";
7289 
7290  // check whether some other vector component has the correct number of
7291  // points
7292  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7293  proposed_fe_comp = numbers::invalid_unsigned_int,
7294  proposed_quad_comp = numbers::invalid_unsigned_int;
7295  if (dof_no != numbers::invalid_unsigned_int)
7296  {
7297  if (static_cast<unsigned int>(fe_degree) ==
7298  this->data->data.front().fe_degree)
7299  {
7300  proposed_dof_comp = dof_no;
7301  proposed_fe_comp = first_selected_component;
7302  }
7303  else
7304  for (unsigned int no = 0; no < this->matrix_free->n_components();
7305  ++no)
7306  for (unsigned int nf = 0;
7307  nf < this->matrix_free->n_base_elements(no);
7308  ++nf)
7309  if (this->matrix_free
7310  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7311  .data.front()
7312  .fe_degree == static_cast<unsigned int>(fe_degree))
7313  {
7314  proposed_dof_comp = no;
7315  proposed_fe_comp = nf;
7316  break;
7317  }
7318  if (n_q_points ==
7319  this->mapping_data->descriptor[this->active_quad_index]
7320  .n_q_points)
7321  proposed_quad_comp = this->quad_no;
7322  else
7323  for (unsigned int no = 0;
7324  no < this->matrix_free->get_mapping_info().cell_data.size();
7325  ++no)
7326  if (this->matrix_free->get_mapping_info()
7327  .cell_data[no]
7328  .descriptor[this->active_quad_index]
7329  .n_q_points == n_q_points)
7330  {
7331  proposed_quad_comp = no;
7332  break;
7333  }
7334  }
7335  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7336  proposed_quad_comp != numbers::invalid_unsigned_int)
7337  {
7338  if (proposed_dof_comp != first_selected_component)
7339  message += "Wrong vector component selection:\n";
7340  else
7341  message += "Wrong quadrature formula selection:\n";
7342  message += " Did you mean FEEvaluation<dim,";
7343  message += Utilities::int_to_string(fe_degree) + ",";
7344  message += Utilities::int_to_string(n_q_points_1d);
7345  message += "," + Utilities::int_to_string(n_components);
7346  message += ",Number>(data";
7347  if (dof_no != numbers::invalid_unsigned_int)
7348  {
7349  message +=
7350  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7351  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7352  message += Utilities::int_to_string(proposed_fe_comp);
7353  }
7354  message += ")?\n";
7355  std::string correct_pos;
7356  if (proposed_dof_comp != dof_no)
7357  correct_pos = " ^ ";
7358  else
7359  correct_pos = " ";
7360  if (proposed_quad_comp != this->quad_no)
7361  correct_pos += " ^ ";
7362  else
7363  correct_pos += " ";
7364  if (proposed_fe_comp != first_selected_component)
7365  correct_pos += " ^\n";
7366  else
7367  correct_pos += " \n";
7368  message += " " +
7369  correct_pos;
7370  }
7371  // ok, did not find the numbers specified by the template arguments in
7372  // the given list. Suggest correct template arguments
7373  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7374  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7375  message += "Wrong template arguments:\n";
7376  message += " Did you mean FEEvaluation<dim,";
7377  message +=
7378  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7379  message += Utilities::int_to_string(proposed_n_q_points_1d);
7380  message += "," + Utilities::int_to_string(n_components);
7381  message += ",Number>(data";
7382  if (dof_no != numbers::invalid_unsigned_int)
7383  {
7384  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7385  message += Utilities::int_to_string(this->quad_no);
7386  message += ", " + Utilities::int_to_string(first_selected_component);
7387  }
7388  message += ")?\n";
7389  std::string correct_pos;
7390  if (this->data->data.front().fe_degree !=
7391  static_cast<unsigned int>(fe_degree))
7392  correct_pos = " ^";
7393  else
7394  correct_pos = " ";
7395  if (proposed_n_q_points_1d != n_q_points_1d)
7396  correct_pos += " ^\n";
7397  else
7398  correct_pos += " \n";
7399  message += " " + correct_pos;
7400 
7401  Assert(static_cast<unsigned int>(fe_degree) ==
7402  this->data->data.front().fe_degree &&
7403  n_q_points == this->n_quadrature_points,
7404  ExcMessage(message));
7405  }
7406  if (dof_no != numbers::invalid_unsigned_int)
7408  n_q_points,
7409  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7410 # endif
7411 }
7412 
7413 
7414 
7415 template <int dim,
7416  int fe_degree,
7417  int n_q_points_1d,
7418  int n_components_,
7419  typename Number,
7420  typename VectorizedArrayType>
7421 inline void
7422 FEEvaluation<dim,
7423  fe_degree,
7424  n_q_points_1d,
7425  n_components_,
7426  Number,
7427  VectorizedArrayType>::reinit(const unsigned int cell_index)
7428 {
7429  Assert(this->matrix_free != nullptr,
7430  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7431  " Integer indexing is not possible."));
7432 
7433  Assert(this->dof_info != nullptr, ExcNotInitialized());
7434  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7435  this->cell = cell_index;
7436  this->cell_type =
7437  this->matrix_free->get_mapping_info().get_cell_type(cell_index);
7438 
7439  const unsigned int offsets =
7440  this->mapping_data->data_index_offsets[cell_index];
7441  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7442  this->J_value = &this->mapping_data->JxW_values[offsets];
7443  if (!this->mapping_data->jacobian_gradients[0].empty())
7444  {
7445  this->jacobian_gradients =
7446  this->mapping_data->jacobian_gradients[0].data() + offsets;
7447  this->jacobian_gradients_non_inverse =
7448  this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
7449  }
7450 
7451  if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
7452  VectorizedArrayType::size())
7453  {
7455  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
7456  this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7457  }
7458  else
7459  {
7460  unsigned int i = 0;
7461  for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7462  ++i)
7463  this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7464  for (; i < VectorizedArrayType::size(); ++i)
7465  this->cell_ids[i] = numbers::invalid_unsigned_int;
7466  }
7467 
7468  if (this->mapping_data->quadrature_points.empty() == false)
7469  this->quadrature_points =
7470  &this->mapping_data->quadrature_points
7471  [this->mapping_data->quadrature_point_offsets[this->cell]];
7472 
7473 # ifdef DEBUG
7474  this->is_reinitialized = true;
7475  this->dof_values_initialized = false;
7476  this->values_quad_initialized = false;
7477  this->gradients_quad_initialized = false;
7478  this->hessians_quad_initialized = false;
7479 # endif
7480 }
7481 
7482 
7483 
7484 template <int dim,
7485  int fe_degree,
7486  int n_q_points_1d,
7487  int n_components_,
7488  typename Number,
7489  typename VectorizedArrayType>
7490 inline void
7491 FEEvaluation<dim,
7492  fe_degree,
7493  n_q_points_1d,
7494  n_components_,
7495  Number,
7496  VectorizedArrayType>::
7497  reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
7498 {
7499  Assert(this->dof_info != nullptr, ExcNotInitialized());
7500  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7501 
7502  this->cell = numbers::invalid_unsigned_int;
7503  this->cell_ids = cell_ids;
7504 
7505  // determine type of cell batch
7507 
7508  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7509  {
7510  const unsigned int cell_index = cell_ids[v];
7511 
7513  continue;
7514 
7515  this->cell_type =
7516  std::max(this->cell_type,
7517  this->matrix_free->get_mapping_info().get_cell_type(
7518  cell_index / VectorizedArrayType::size()));
7519  }
7520 
7521  // allocate memory for internal data storage
7522  if (this->mapped_geometry == nullptr)
7523  this->mapped_geometry =
7524  std::make_shared<internal::MatrixFreeFunctions::
7525  MappingDataOnTheFly<dim, VectorizedArrayType>>();
7526 
7527  auto &mapping_storage = this->mapped_geometry->get_data_storage();
7528 
7529  auto &this_jacobian_data = mapping_storage.jacobians[0];
7530  auto &this_J_value_data = mapping_storage.JxW_values;
7531  auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7532  auto &this_jacobian_gradients_non_inverse_data =
7533  mapping_storage.jacobian_gradients_non_inverse[0];
7534  auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7535 
7537  {
7538  if (this->mapping_data->jacobians[0].size() > 0)
7539  this_jacobian_data.resize_fast(2);
7540 
7541  if (this->mapping_data->JxW_values.size() > 0)
7542  this_J_value_data.resize_fast(1);
7543 
7544  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7545  this_jacobian_gradients_data.resize_fast(1);
7546 
7547  if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7548  this_jacobian_gradients_non_inverse_data.resize_fast(1);
7549 
7550  if (this->mapping_data->quadrature_points.size() > 0)
7551  this_quadrature_points_data.resize_fast(1);
7552  }
7553  else
7554  {
7555  if (this->mapping_data->jacobians[0].size() > 0)
7556  this_jacobian_data.resize_fast(this->n_quadrature_points);
7557 
7558  if (this->mapping_data->JxW_values.size() > 0)
7559  this_J_value_data.resize_fast(this->n_quadrature_points);
7560 
7561  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7562  this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7563 
7564  if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7565  this_jacobian_gradients_non_inverse_data.resize_fast(
7566  this->n_quadrature_points);
7567 
7568  if (this->mapping_data->quadrature_points.size() > 0)
7569  this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7570  }
7571 
7572  // set pointers to internal data storage
7573  this->jacobian = this_jacobian_data.data();
7574  this->J_value = this_J_value_data.data();
7575  this->jacobian_gradients = this_jacobian_gradients_data.data();
7576  this->jacobian_gradients_non_inverse =
7577  this_jacobian_gradients_non_inverse_data.data();
7578  this->quadrature_points = this_quadrature_points_data.data();
7579 
7580  // fill internal data storage lane by lane
7581  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7582  {
7583&#