Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
fe_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 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 
281 
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 
582 
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 
1956  static constexpr unsigned int static_n_q_points =
1957  Utilities::pow(n_q_points_1d, dim);
1958 
1966  static constexpr unsigned int static_dofs_per_component =
1967  Utilities::pow(fe_degree + 1, dim);
1968 
1976  static constexpr unsigned int tensor_dofs_per_cell =
1978 
1986  static constexpr unsigned int static_dofs_per_cell =
1988 
2025  const unsigned int dof_no = 0,
2026  const unsigned int quad_no = 0,
2027  const unsigned int first_selected_component = 0,
2028  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2029  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2030 
2039  const std::pair<unsigned int, unsigned int> & range,
2040  const unsigned int dof_no = 0,
2041  const unsigned int quad_no = 0,
2042  const unsigned int first_selected_component = 0);
2043 
2072  FEEvaluation(const Mapping<dim> & mapping,
2073  const FiniteElement<dim> &fe,
2074  const Quadrature<1> & quadrature,
2075  const UpdateFlags update_flags,
2076  const unsigned int first_selected_component = 0);
2077 
2084  const Quadrature<1> & quadrature,
2085  const UpdateFlags update_flags,
2086  const unsigned int first_selected_component = 0);
2087 
2100  const unsigned int first_selected_component = 0);
2101 
2109 
2116  FEEvaluation &
2117  operator=(const FEEvaluation &other);
2118 
2127  void
2128  reinit(const unsigned int cell_batch_index);
2129 
2136  void
2137  reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids);
2138 
2151  template <bool level_dof_access>
2152  void
2154 
2165  void
2167 
2171  static bool
2172  fast_evaluation_supported(const unsigned int given_degree,
2173  const unsigned int give_n_q_points_1d);
2174 
2184  void
2186 
2191  DEAL_II_DEPRECATED void
2192  evaluate(const bool evaluate_values,
2193  const bool evaluate_gradients,
2194  const bool evaluate_hessians = false);
2195 
2208  void
2209  evaluate(const VectorizedArrayType * values_array,
2210  const EvaluationFlags::EvaluationFlags evaluation_flag);
2211 
2216  DEAL_II_DEPRECATED void
2217  evaluate(const VectorizedArrayType *values_array,
2218  const bool evaluate_values,
2219  const bool evaluate_gradients,
2220  const bool evaluate_hessians = false);
2221 
2235  template <typename VectorType>
2236  void
2237  gather_evaluate(const VectorType & input_vector,
2238  const EvaluationFlags::EvaluationFlags evaluation_flag);
2239 
2243  template <typename VectorType>
2244  DEAL_II_DEPRECATED void
2245  gather_evaluate(const VectorType &input_vector,
2246  const bool evaluate_values,
2247  const bool evaluate_gradients,
2248  const bool evaluate_hessians = false);
2249 
2260  void
2262 
2266  DEAL_II_DEPRECATED void
2267  integrate(const bool integrate_values, const bool integrate_gradients);
2268 
2280  void
2282  VectorizedArrayType * values_array,
2283  const bool sum_into_values = false);
2284 
2288  DEAL_II_DEPRECATED void
2289  integrate(const bool integrate_values,
2290  const bool integrate_gradients,
2291  VectorizedArrayType *values_array);
2292 
2306  template <typename VectorType>
2307  void
2309  VectorType & output_vector);
2310 
2314  template <typename VectorType>
2315  DEAL_II_DEPRECATED void
2316  integrate_scatter(const bool integrate_values,
2317  const bool integrate_gradients,
2318  VectorType &output_vector);
2319 
2326  dof_indices() const;
2327 
2334  const unsigned int dofs_per_component;
2335 
2342  const unsigned int dofs_per_cell;
2343 
2351  const unsigned int n_q_points;
2352 
2353 private:
2358  void
2359  check_template_arguments(const unsigned int fe_no,
2360  const unsigned int first_selected_component);
2361 };
2362 
2363 
2364 
2400 template <int dim,
2401  int fe_degree,
2402  int n_q_points_1d = fe_degree + 1,
2403  int n_components_ = 1,
2404  typename Number = double,
2405  typename VectorizedArrayType = VectorizedArray<Number>>
2407  n_components_,
2408  Number,
2409  true,
2410  VectorizedArrayType>
2411 {
2412  static_assert(
2413  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2414  "Type of Number and of VectorizedArrayType do not match.");
2415 
2416 public:
2420  using BaseClass =
2422 
2426  using number_type = Number;
2427 
2434 
2441 
2445  static constexpr unsigned int dimension = dim;
2446 
2451  static constexpr unsigned int n_components = n_components_;
2452 
2460  static constexpr unsigned int static_n_q_points =
2461  Utilities::pow(n_q_points_1d, dim - 1);
2462 
2469  static constexpr unsigned int static_n_q_points_cell =
2470  Utilities::pow(n_q_points_1d, dim);
2471 
2478  static constexpr unsigned int static_dofs_per_component =
2479  Utilities::pow(fe_degree + 1, dim);
2480 
2487  static constexpr unsigned int tensor_dofs_per_cell =
2489 
2496  static constexpr unsigned int static_dofs_per_cell =
2498 
2542  const bool is_interior_face = true,
2543  const unsigned int dof_no = 0,
2544  const unsigned int quad_no = 0,
2545  const unsigned int first_selected_component = 0,
2546  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2547  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
2548  const unsigned int face_type = numbers::invalid_unsigned_int);
2549 
2559  const std::pair<unsigned int, unsigned int> & range,
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 
2575  void
2576  reinit(const unsigned int face_batch_number);
2577 
2585  void
2586  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2587 
2591  static bool
2592  fast_evaluation_supported(const unsigned int given_degree,
2593  const unsigned int give_n_q_points_1d);
2594 
2605  void
2607 
2611  DEAL_II_DEPRECATED void
2612  evaluate(const bool evaluate_values, const bool evaluate_gradients);
2613 
2626  void
2627  evaluate(const VectorizedArrayType * values_array,
2628  const EvaluationFlags::EvaluationFlags evaluation_flag);
2629 
2633  DEAL_II_DEPRECATED void
2634  evaluate(const VectorizedArrayType *values_array,
2635  const bool evaluate_values,
2636  const bool evaluate_gradients);
2637 
2649  template <typename VectorType>
2650  void
2651  gather_evaluate(const VectorType & input_vector,
2652  const EvaluationFlags::EvaluationFlags evaluation_flag);
2653 
2657  template <typename VectorType>
2658  DEAL_II_DEPRECATED void
2659  gather_evaluate(const VectorType &input_vector,
2660  const bool evaluate_values,
2661  const bool evaluate_gradients);
2662 
2672  void
2674 
2678  DEAL_II_DEPRECATED void
2679  integrate(const bool integrate_values, const bool integrate_gradients);
2680 
2689  void
2691  VectorizedArrayType * values_array);
2692 
2696  DEAL_II_DEPRECATED void
2697  integrate(const bool integrate_values,
2698  const bool integrate_gradients,
2699  VectorizedArrayType *values_array);
2700 
2712  template <typename VectorType>
2713  void
2715  VectorType & output_vector);
2716 
2720  template <typename VectorType>
2721  void
2722  integrate_scatter(const bool integrate_values,
2723  const bool integrate_gradients,
2724  VectorType &output_vector);
2725 
2732  dof_indices() const;
2733 
2740  const unsigned int dofs_per_component;
2741 
2748  const unsigned int dofs_per_cell;
2749 
2757  const unsigned int n_q_points;
2758 };
2759 
2760 
2761 
2762 namespace internal
2763 {
2764  namespace MatrixFreeFunctions
2765  {
2766  // a helper function to compute the number of DoFs of a DGP element at
2767  // compile time, depending on the degree
2768  template <int dim, int degree>
2770  {
2771  // this division is always without remainder
2772  static constexpr unsigned int value =
2773  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2774  };
2775 
2776  // base specialization: 1d elements have 'degree+1' degrees of freedom
2777  template <int degree>
2778  struct DGP_dofs_per_component<1, degree>
2779  {
2780  static constexpr unsigned int value = degree + 1;
2781  };
2782  } // namespace MatrixFreeFunctions
2783 } // namespace internal
2784 
2785 
2786 /*----------------------- Inline functions ----------------------------------*/
2787 
2788 #ifndef DOXYGEN
2789 
2790 
2791 namespace internal
2792 {
2793  // Extract all internal data pointers and indices in a single function that
2794  // get passed on to the constructor of FEEvaluationData, avoiding to look
2795  // things up multiple times
2796  template <bool is_face,
2797  int dim,
2798  typename Number,
2799  typename VectorizedArrayType>
2801  InitializationData
2802  extract_initialization_data(
2804  const unsigned int dof_no,
2805  const unsigned int first_selected_component,
2806  const unsigned int quad_no,
2807  const unsigned int fe_degree,
2808  const unsigned int n_q_points,
2809  const unsigned int active_fe_index_given,
2810  const unsigned int active_quad_index_given,
2811  const unsigned int face_type)
2812  {
2814  InitializationData init_data;
2815 
2816  init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2817  init_data.mapping_data =
2818  &internal::MatrixFreeFunctions::
2819  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2820  matrix_free.get_mapping_info(), quad_no);
2821 
2822  init_data.active_fe_index =
2823  fe_degree != numbers::invalid_unsigned_int ?
2824  init_data.dof_info->fe_index_from_degree(first_selected_component,
2825  fe_degree) :
2826  (active_fe_index_given != numbers::invalid_unsigned_int ?
2827  active_fe_index_given :
2828  0);
2829  init_data.active_quad_index =
2830  fe_degree == numbers::invalid_unsigned_int ?
2831  (active_quad_index_given != numbers::invalid_unsigned_int ?
2832  active_quad_index_given :
2833  std::min<unsigned int>(init_data.active_fe_index,
2834  init_data.mapping_data->descriptor.size() -
2835  1)) :
2836  init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2837 
2838  init_data.shape_info = &matrix_free.get_shape_info(
2839  dof_no,
2840  quad_no,
2841  init_data.dof_info->component_to_base_index[first_selected_component],
2842  init_data.active_fe_index,
2843  init_data.active_quad_index);
2844  init_data.descriptor =
2845  &init_data.mapping_data->descriptor
2846  [is_face ?
2847  (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2848  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2849  init_data.active_quad_index];
2850 
2851  return init_data;
2852  }
2853 } // namespace internal
2854 
2855 
2856 
2857 /*----------------------- FEEvaluationBase ----------------------------------*/
2858 
2859 template <int dim,
2860  int n_components_,
2861  typename Number,
2862  bool is_face,
2863  typename VectorizedArrayType>
2864 inline FEEvaluationBase<dim,
2865  n_components_,
2866  Number,
2867  is_face,
2868  VectorizedArrayType>::
2869  FEEvaluationBase(
2871  const unsigned int dof_no,
2872  const unsigned int first_selected_component,
2873  const unsigned int quad_no,
2874  const unsigned int fe_degree,
2875  const unsigned int n_q_points,
2876  const bool is_interior_face,
2877  const unsigned int active_fe_index,
2878  const unsigned int active_quad_index,
2879  const unsigned int face_type)
2880  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2881  internal::extract_initialization_data<is_face>(matrix_free,
2882  dof_no,
2883  first_selected_component,
2884  quad_no,
2885  fe_degree,
2886  n_q_points,
2887  active_fe_index,
2888  active_quad_index,
2889  face_type),
2890  is_interior_face,
2891  quad_no,
2892  first_selected_component)
2893  , scratch_data_array(matrix_free.acquire_scratch_data())
2894  , matrix_free(&matrix_free)
2895 {
2896  this->set_data_pointers(scratch_data_array, n_components_);
2897  Assert(
2898  this->dof_info->start_components.back() == 1 ||
2899  static_cast<int>(n_components_) <=
2900  static_cast<int>(
2901  this->dof_info->start_components
2902  [this->dof_info->component_to_base_index[first_selected_component] +
2903  1]) -
2904  first_selected_component,
2905  ExcMessage(
2906  "You tried to construct a vector-valued evaluator with " +
2907  std::to_string(n_components) +
2908  " components. However, "
2909  "the current base element has only " +
2911  this->dof_info->start_components
2912  [this->dof_info->component_to_base_index[first_selected_component] +
2913  1] -
2914  first_selected_component) +
2915  " components left when starting from local element index " +
2917  first_selected_component -
2918  this->dof_info->start_components
2919  [this->dof_info->component_to_base_index[first_selected_component]]) +
2920  " (global index " + std::to_string(first_selected_component) + ")"));
2921 
2922  // do not check for correct dimensions of data fields here, should be done
2923  // in derived classes
2924 }
2925 
2926 
2927 
2928 template <int dim,
2929  int n_components_,
2930  typename Number,
2931  bool is_face,
2932  typename VectorizedArrayType>
2933 inline FEEvaluationBase<dim,
2934  n_components_,
2935  Number,
2936  is_face,
2937  VectorizedArrayType>::
2938  FEEvaluationBase(
2939  const Mapping<dim> & mapping,
2940  const FiniteElement<dim> &fe,
2941  const Quadrature<1> & quadrature,
2942  const UpdateFlags update_flags,
2943  const unsigned int first_selected_component,
2945  : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2946  other != nullptr &&
2947  other->mapped_geometry->get_quadrature() == quadrature ?
2948  other->mapped_geometry :
2949  std::make_shared<internal::MatrixFreeFunctions::
2950  MappingDataOnTheFly<dim, VectorizedArrayType>>(
2951  mapping,
2952  quadrature,
2953  update_flags),
2954  n_components_,
2955  first_selected_component)
2956  , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2957  , matrix_free(nullptr)
2958 {
2959  const unsigned int base_element_number =
2960  fe.component_to_base_index(first_selected_component).first;
2961  Assert(fe.element_multiplicity(base_element_number) == 1 ||
2962  fe.element_multiplicity(base_element_number) -
2963  first_selected_component >=
2964  n_components_,
2965  ExcMessage("The underlying element must at least contain as many "
2966  "components as requested by this class"));
2967  (void)base_element_number;
2968 
2969  Assert(this->data == nullptr, ExcInternalError());
2970  this->data =
2972  Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2973  fe,
2974  fe.component_to_base_index(first_selected_component).first);
2975 
2976  this->set_data_pointers(scratch_data_array, n_components_);
2977 }
2978 
2979 
2980 
2981 template <int dim,
2982  int n_components_,
2983  typename Number,
2984  bool is_face,
2985  typename VectorizedArrayType>
2986 inline FEEvaluationBase<dim,
2987  n_components_,
2988  Number,
2989  is_face,
2990  VectorizedArrayType>::
2991  FEEvaluationBase(const FEEvaluationBase<dim,
2992  n_components_,
2993  Number,
2994  is_face,
2995  VectorizedArrayType> &other)
2996  : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
2997  , scratch_data_array(other.matrix_free == nullptr ?
2998  new AlignedVector<VectorizedArrayType>() :
2999  other.matrix_free->acquire_scratch_data())
3000  , matrix_free(other.matrix_free)
3001 {
3002  if (other.matrix_free == nullptr)
3003  {
3004  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3005  this->data =
3007  *other.data);
3008 
3009  // Create deep copy of mapped geometry for use in parallel
3010  this->mapped_geometry =
3011  std::make_shared<internal::MatrixFreeFunctions::
3012  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3013  other.mapped_geometry->get_fe_values().get_mapping(),
3014  other.mapped_geometry->get_quadrature(),
3015  other.mapped_geometry->get_fe_values().get_update_flags());
3016  this->mapping_data = &this->mapped_geometry->get_data_storage();
3017  this->cell = 0;
3018 
3019  this->jacobian =
3020  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3021  this->J_value =
3022  this->mapped_geometry->get_data_storage().JxW_values.begin();
3023  this->jacobian_gradients =
3024  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3025  this->quadrature_points =
3026  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3027  }
3028 
3029  this->set_data_pointers(scratch_data_array, n_components_);
3030 }
3031 
3032 
3033 
3034 template <int dim,
3035  int n_components_,
3036  typename Number,
3037  bool is_face,
3038  typename VectorizedArrayType>
3039 inline FEEvaluationBase<dim,
3040  n_components_,
3041  Number,
3042  is_face,
3043  VectorizedArrayType> &
3045 operator=(const FEEvaluationBase<dim,
3046  n_components_,
3047  Number,
3048  is_face,
3049  VectorizedArrayType> &other)
3050 {
3051  // release old memory
3052  if (matrix_free == nullptr)
3053  {
3054  delete this->data;
3055  delete scratch_data_array;
3056  }
3057  else
3058  {
3059  matrix_free->release_scratch_data(scratch_data_array);
3060  }
3061 
3063 
3064  matrix_free = other.matrix_free;
3065 
3066  if (other.matrix_free == nullptr)
3067  {
3068  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3069  this->data =
3071  *other.data);
3072  scratch_data_array = new AlignedVector<VectorizedArrayType>();
3073 
3074  // Create deep copy of mapped geometry for use in parallel
3075  this->mapped_geometry =
3076  std::make_shared<internal::MatrixFreeFunctions::
3077  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3078  other.mapped_geometry->get_fe_values().get_mapping(),
3079  other.mapped_geometry->get_quadrature(),
3080  other.mapped_geometry->get_fe_values().get_update_flags());
3081  this->cell = 0;
3082  this->mapping_data = &this->mapped_geometry->get_data_storage();
3083  this->jacobian =
3084  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3085  this->J_value =
3086  this->mapped_geometry->get_data_storage().JxW_values.begin();
3087  this->jacobian_gradients =
3088  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3089  this->quadrature_points =
3090  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3091  }
3092  else
3093  {
3094  scratch_data_array = matrix_free->acquire_scratch_data();
3095  }
3096 
3097  this->set_data_pointers(scratch_data_array, n_components_);
3098 
3099  return *this;
3100 }
3101 
3102 
3103 
3104 template <int dim,
3105  int n_components_,
3106  typename Number,
3107  bool is_face,
3108  typename VectorizedArrayType>
3109 inline FEEvaluationBase<dim,
3110  n_components_,
3111  Number,
3112  is_face,
3113  VectorizedArrayType>::~FEEvaluationBase()
3114 {
3115  if (matrix_free != nullptr)
3116  {
3117  try
3118  {
3119  matrix_free->release_scratch_data(scratch_data_array);
3120  }
3121  catch (...)
3122  {}
3123  }
3124  else
3125  {
3126  delete scratch_data_array;
3127  delete this->data;
3128  }
3129 }
3130 
3131 
3132 
3133 template <int dim,
3134  int n_components_,
3135  typename Number,
3136  bool is_face,
3137  typename VectorizedArrayType>
3140  get_matrix_free() const
3141 {
3142  Assert(matrix_free != nullptr,
3143  ExcMessage(
3144  "FEEvaluation was not initialized with a MatrixFree object!"));
3145  return *matrix_free;
3146 }
3147 
3148 
3149 
3150 namespace internal
3151 {
3152  // given a block vector return the underlying vector type
3153  // including constness (specified by bool)
3154  template <typename VectorType, bool>
3155  struct ConstBlockVectorSelector;
3156 
3157  template <typename VectorType>
3158  struct ConstBlockVectorSelector<VectorType, true>
3159  {
3160  using BaseVectorType = const typename VectorType::BlockType;
3161  };
3162 
3163  template <typename VectorType>
3164  struct ConstBlockVectorSelector<VectorType, false>
3165  {
3166  using BaseVectorType = typename VectorType::BlockType;
3167  };
3168 
3169  // allows to select between block vectors and non-block vectors, which
3170  // allows to use a unified interface for extracting blocks on block vectors
3171  // and doing nothing on usual vectors
3172  template <typename VectorType, bool>
3173  struct BlockVectorSelector;
3174 
3175  template <typename VectorType>
3176  struct BlockVectorSelector<VectorType, true>
3177  {
3178  using BaseVectorType = typename ConstBlockVectorSelector<
3179  VectorType,
3180  std::is_const<VectorType>::value>::BaseVectorType;
3181 
3182  static BaseVectorType *
3183  get_vector_component(VectorType &vec, const unsigned int component)
3184  {
3185  AssertIndexRange(component, vec.n_blocks());
3186  return &vec.block(component);
3187  }
3188  };
3189 
3190  template <typename VectorType>
3191  struct BlockVectorSelector<VectorType, false>
3192  {
3193  using BaseVectorType = VectorType;
3194 
3195  static BaseVectorType *
3196  get_vector_component(VectorType &vec, const unsigned int component)
3197  {
3198  // FEEvaluation allows to combine several vectors from a scalar
3199  // FiniteElement into a "vector-valued" FEEvaluation object with
3200  // multiple components. These components can be extracted with the other
3201  // get_vector_component functions. If we do not get a vector of vectors
3202  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3203  // must make sure that we do not duplicate the components in input
3204  // and/or duplicate the resulting integrals. In such a case, we should
3205  // only get the zeroth component in the vector contained set nullptr for
3206  // the others which allows us to catch unintended use in
3207  // read_write_operation.
3208  if (component == 0)
3209  return &vec;
3210  else
3211  return nullptr;
3212  }
3213  };
3214 
3215  template <typename VectorType>
3216  struct BlockVectorSelector<std::vector<VectorType>, false>
3217  {
3218  using BaseVectorType = VectorType;
3219 
3220  static BaseVectorType *
3221  get_vector_component(std::vector<VectorType> &vec,
3222  const unsigned int component)
3223  {
3224  AssertIndexRange(component, vec.size());
3225  return &vec[component];
3226  }
3227  };
3228 
3229  template <typename VectorType>
3230  struct BlockVectorSelector<const std::vector<VectorType>, false>
3231  {
3232  using BaseVectorType = const VectorType;
3233 
3234  static const BaseVectorType *
3235  get_vector_component(const std::vector<VectorType> &vec,
3236  const unsigned int component)
3237  {
3238  AssertIndexRange(component, vec.size());
3239  return &vec[component];
3240  }
3241  };
3242 
3243  template <typename VectorType>
3244  struct BlockVectorSelector<std::vector<VectorType *>, false>
3245  {
3246  using BaseVectorType = VectorType;
3247 
3248  static BaseVectorType *
3249  get_vector_component(std::vector<VectorType *> &vec,
3250  const unsigned int component)
3251  {
3252  AssertIndexRange(component, vec.size());
3253  return vec[component];
3254  }
3255  };
3256 
3257  template <typename VectorType>
3258  struct BlockVectorSelector<const std::vector<VectorType *>, false>
3259  {
3260  using BaseVectorType = const VectorType;
3261 
3262  static const BaseVectorType *
3263  get_vector_component(const std::vector<VectorType *> &vec,
3264  const unsigned int component)
3265  {
3266  AssertIndexRange(component, vec.size());
3267  return vec[component];
3268  }
3269  };
3270 } // namespace internal
3271 
3272 
3273 
3274 template <int dim,
3275  int n_components_,
3276  typename Number,
3277  bool is_face,
3278  typename VectorizedArrayType>
3279 template <typename VectorType, typename VectorOperation>
3280 inline void
3283  const VectorOperation & operation,
3284  const std::array<VectorType *, n_components_> &src,
3285  const std::array<
3287  n_components_> & src_sm,
3288  const std::bitset<VectorizedArrayType::size()> &mask,
3289  const bool apply_constraints) const
3290 {
3291  // Case 1: No MatrixFree object given, simple case because we do not need to
3292  // process constraints and need not care about vectorization -> go to
3293  // separate function
3294  if (this->matrix_free == nullptr)
3295  {
3296  read_write_operation_global(operation, src);
3297  return;
3298  }
3299 
3300  Assert(this->dof_info != nullptr, ExcNotInitialized());
3301  Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3302  if (this->n_fe_components == 1)
3303  for (unsigned int comp = 0; comp < n_components; ++comp)
3304  {
3305  Assert(src[comp] != nullptr,
3306  ExcMessage("The finite element underlying this FEEvaluation "
3307  "object is scalar, but you requested " +
3308  std::to_string(n_components) +
3309  " components via the template argument in "
3310  "FEEvaluation. In that case, you must pass an "
3311  "std::vector<VectorType> or a BlockVector to " +
3312  "read_dof_values and distribute_local_to_global."));
3314  *this->matrix_free,
3315  *this->dof_info);
3316  }
3317  else
3318  {
3320  *this->matrix_free,
3321  *this->dof_info);
3322  }
3323 
3324  // Case 2: contiguous indices which use reduced storage of indices and can
3325  // use vectorized load/store operations -> go to separate function
3326  if (this->cell != numbers::invalid_unsigned_int)
3327  {
3329  this->cell,
3330  this->dof_info->index_storage_variants[this->dof_access_index].size());
3331  if (this->dof_info->index_storage_variants
3332  [is_face ? this->dof_access_index :
3334  [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::
3335  IndexStorageVariants::contiguous)
3336  {
3337  read_write_operation_contiguous(operation, src, src_sm, mask);
3338  return;
3339  }
3340  }
3341 
3342  // Case 3: standard operation with one index per degree of freedom -> go on
3343  // here
3344  constexpr unsigned int n_lanes = VectorizedArrayType::size();
3345 
3346  std::array<unsigned int, VectorizedArrayType::size()> cells =
3347  this->get_cell_ids();
3348 
3349  for (unsigned int v = 0; v < n_lanes; ++v)
3350  if (mask[v] == false)
3351  cells[v] = numbers::invalid_unsigned_int;
3352 
3353  bool has_hn_constraints = false;
3354 
3355  if (is_face == false)
3356  {
3357  for (unsigned int v = 0; v < n_lanes; ++v)
3358  if (cells[v] != numbers::invalid_unsigned_int &&
3359  this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3360  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3361  this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3364  this->dof_info->hanging_node_constraint_masks_comp
3365  [this->active_fe_index][this->first_selected_component])
3366  has_hn_constraints = true;
3367  }
3368 
3369  std::integral_constant<bool,
3370  internal::is_vectorizable<VectorType, Number>::value>
3371  vector_selector;
3372 
3373  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3374  std::array<VectorizedArrayType *, n_components> values_dofs;
3375  for (unsigned int c = 0; c < n_components; ++c)
3376  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3377  c * dofs_per_component;
3378 
3379  if (this->cell != numbers::invalid_unsigned_int &&
3380  this->dof_info->index_storage_variants
3381  [is_face ? this->dof_access_index :
3383  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3384  IndexStorageVariants::interleaved &&
3385  (has_hn_constraints == false))
3386  {
3387  const unsigned int *dof_indices =
3388  this->dof_info->dof_indices_interleaved.data() +
3389  this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes]
3390  .first +
3391  this->dof_info
3392  ->component_dof_indices_offset[this->active_fe_index]
3393  [this->first_selected_component] *
3394  n_lanes;
3395  if (n_components == 1 || this->n_fe_components == 1)
3396  for (unsigned int i = 0; i < dofs_per_component;
3397  ++i, dof_indices += n_lanes)
3398  for (unsigned int comp = 0; comp < n_components; ++comp)
3399  operation.process_dof_gather(dof_indices,
3400  *src[comp],
3401  0,
3402  values_dofs[comp][i],
3403  vector_selector);
3404  else
3405  for (unsigned int comp = 0; comp < n_components; ++comp)
3406  for (unsigned int i = 0; i < dofs_per_component;
3407  ++i, dof_indices += n_lanes)
3408  operation.process_dof_gather(
3409  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3410  return;
3411  }
3412 
3413  // Allocate pointers, then initialize all of them to nullptrs and
3414  // below overwrite the ones we actually use:
3415  std::array<const unsigned int *, n_lanes> dof_indices;
3416  dof_indices.fill(nullptr);
3417 
3418  // Assign the appropriate cell ids for face/cell case and get the pointers
3419  // to the dof indices of the cells on all lanes
3420 
3421  bool has_constraints = false;
3422  const unsigned int n_components_read =
3423  this->n_fe_components > 1 ? n_components : 1;
3424 
3425  if (is_face)
3426  {
3427  for (unsigned int v = 0; v < n_lanes; ++v)
3428  {
3429  if (cells[v] == numbers::invalid_unsigned_int)
3430  continue;
3431 
3432  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
3433  ExcInternalError());
3434  const std::pair<unsigned int, unsigned int> *my_index_start =
3435  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3436  this->first_selected_component];
3437 
3438  // check whether any of the SIMD lanes has constraints, i.e., the
3439  // constraint indicator which is the second entry of row_starts
3440  // increments on this cell
3441  if (my_index_start[n_components_read].second !=
3442  my_index_start[0].second)
3443  has_constraints = true;
3444 
3445  dof_indices[v] =
3446  this->dof_info->dof_indices.data() + my_index_start[0].first;
3447  }
3448  }
3449  else
3450  {
3451  for (unsigned int v = 0; v < n_lanes; ++v)
3452  {
3453  if (cells[v] == numbers::invalid_unsigned_int)
3454  continue;
3455 
3456  const std::pair<unsigned int, unsigned int> *my_index_start =
3457  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3458  this->first_selected_component];
3459  if (my_index_start[n_components_read].second !=
3460  my_index_start[0].second)
3461  has_constraints = true;
3462 
3463  if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3464  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3465  this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3468  this->dof_info->hanging_node_constraint_masks_comp
3469  [this->active_fe_index][this->first_selected_component])
3470  has_hn_constraints = true;
3471 
3472  Assert(my_index_start[n_components_read].first ==
3473  my_index_start[0].first ||
3474  my_index_start[0].first < this->dof_info->dof_indices.size(),
3475  ExcIndexRange(0,
3476  my_index_start[0].first,
3477  this->dof_info->dof_indices.size()));
3478  dof_indices[v] =
3479  this->dof_info->dof_indices.data() + my_index_start[0].first;
3480  }
3481  }
3482 
3483  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3484  return i != numbers::invalid_unsigned_int;
3485  }) < n_lanes)
3486  for (unsigned int comp = 0; comp < n_components; ++comp)
3487  for (unsigned int i = 0; i < dofs_per_component; ++i)
3488  operation.process_empty(values_dofs[comp][i]);
3489 
3490  // Case where we have no constraints throughout the whole cell: Can go
3491  // through the list of DoFs directly
3492  if (!has_constraints && apply_constraints)
3493  {
3494  if (n_components == 1 || this->n_fe_components == 1)
3495  {
3496  for (unsigned int v = 0; v < n_lanes; ++v)
3497  {
3498  if (cells[v] == numbers::invalid_unsigned_int)
3499  continue;
3500 
3501  for (unsigned int i = 0; i < dofs_per_component; ++i)
3502  for (unsigned int comp = 0; comp < n_components; ++comp)
3503  operation.process_dof(dof_indices[v][i],
3504  *src[comp],
3505  values_dofs[comp][i][v]);
3506  }
3507  }
3508  else
3509  {
3510  for (unsigned int comp = 0; comp < n_components; ++comp)
3511  for (unsigned int v = 0; v < n_lanes; ++v)
3512  {
3513  if (cells[v] == numbers::invalid_unsigned_int)
3514  continue;
3515 
3516  for (unsigned int i = 0; i < dofs_per_component; ++i)
3517  operation.process_dof(
3518  dof_indices[v][comp * dofs_per_component + i],
3519  *src[0],
3520  values_dofs[comp][i][v]);
3521  }
3522  }
3523  return;
3524  }
3525 
3526  // In the case where there are some constraints to be resolved, loop over
3527  // all vector components that are filled and then over local dofs. ind_local
3528  // holds local number on cell, index iterates over the elements of
3529  // index_local_to_global and dof_indices points to the global indices stored
3530  // in index_local_to_global
3531 
3532  for (unsigned int v = 0; v < n_lanes; ++v)
3533  {
3534  if (cells[v] == numbers::invalid_unsigned_int)
3535  continue;
3536 
3537  const unsigned int cell_index = cells[v];
3538  const unsigned int cell_dof_index =
3539  cell_index * this->n_fe_components + this->first_selected_component;
3540  const unsigned int n_components_read =
3541  this->n_fe_components > 1 ? n_components : 1;
3542  unsigned int index_indicators =
3543  this->dof_info->row_starts[cell_dof_index].second;
3544  unsigned int next_index_indicators =
3545  this->dof_info->row_starts[cell_dof_index + 1].second;
3546 
3547  // For read_dof_values_plain, redirect the dof_indices field to the
3548  // unconstrained indices
3549  if (apply_constraints == false &&
3550  (this->dof_info->row_starts[cell_dof_index].second !=
3551  this->dof_info->row_starts[cell_dof_index + n_components_read]
3552  .second ||
3553  ((this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3554  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3555  this->dof_info->hanging_node_constraint_masks[cell_index] !=
3558  this->dof_info->hanging_node_constraint_masks_comp
3559  [this->active_fe_index][this->first_selected_component])))
3560  {
3561  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
3563  ExcNotInitialized());
3564  dof_indices[v] =
3565  this->dof_info->plain_dof_indices.data() +
3566  this->dof_info
3567  ->component_dof_indices_offset[this->active_fe_index]
3568  [this->first_selected_component] +
3569  this->dof_info->row_starts_plain_indices[cell_index];
3570  next_index_indicators = index_indicators;
3571  }
3572 
3573  if (n_components == 1 || this->n_fe_components == 1)
3574  {
3575  unsigned int ind_local = 0;
3576  for (; index_indicators != next_index_indicators; ++index_indicators)
3577  {
3578  const std::pair<unsigned short, unsigned short> indicator =
3579  this->dof_info->constraint_indicator[index_indicators];
3580  // run through values up to next constraint
3581  for (unsigned int j = 0; j < indicator.first; ++j)
3582  for (unsigned int comp = 0; comp < n_components; ++comp)
3583  operation.process_dof(dof_indices[v][j],
3584  *src[comp],
3585  values_dofs[comp][ind_local + j][v]);
3586 
3587  ind_local += indicator.first;
3588  dof_indices[v] += indicator.first;
3589 
3590  // constrained case: build the local value as a linear
3591  // combination of the global value according to constraints
3592  Number value[n_components];
3593  for (unsigned int comp = 0; comp < n_components; ++comp)
3594  operation.pre_constraints(values_dofs[comp][ind_local][v],
3595  value[comp]);
3596 
3597  const Number *data_val =
3598  this->matrix_free->constraint_pool_begin(indicator.second);
3599  const Number *end_pool =
3600  this->matrix_free->constraint_pool_end(indicator.second);
3601  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3602  for (unsigned int comp = 0; comp < n_components; ++comp)
3603  operation.process_constraint(*dof_indices[v],
3604  *data_val,
3605  *src[comp],
3606  value[comp]);
3607 
3608  for (unsigned int comp = 0; comp < n_components; ++comp)
3609  operation.post_constraints(value[comp],
3610  values_dofs[comp][ind_local][v]);
3611  ind_local++;
3612  }
3613 
3614  AssertIndexRange(ind_local, dofs_per_component + 1);
3615 
3616  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3617  for (unsigned int comp = 0; comp < n_components; ++comp)
3618  operation.process_dof(*dof_indices[v],
3619  *src[comp],
3620  values_dofs[comp][ind_local][v]);
3621  }
3622  else
3623  {
3624  // case with vector-valued finite elements where all components are
3625  // included in one single vector. Assumption: first come all entries
3626  // to the first component, then all entries to the second one, and
3627  // so on. This is ensured by the way MatrixFree reads out the
3628  // indices.
3629  for (unsigned int comp = 0; comp < n_components; ++comp)
3630  {
3631  unsigned int ind_local = 0;
3632 
3633  // check whether there is any constraint on the current cell
3634  for (; index_indicators != next_index_indicators;
3635  ++index_indicators)
3636  {
3637  const std::pair<unsigned short, unsigned short> indicator =
3638  this->dof_info->constraint_indicator[index_indicators];
3639 
3640  // run through values up to next constraint
3641  for (unsigned int j = 0; j < indicator.first; ++j)
3642  operation.process_dof(dof_indices[v][j],
3643  *src[0],
3644  values_dofs[comp][ind_local + j][v]);
3645  ind_local += indicator.first;
3646  dof_indices[v] += indicator.first;
3647 
3648  // constrained case: build the local value as a linear
3649  // combination of the global value according to constraints
3650  Number value;
3651  operation.pre_constraints(values_dofs[comp][ind_local][v],
3652  value);
3653 
3654  const Number *data_val =
3655  this->matrix_free->constraint_pool_begin(indicator.second);
3656  const Number *end_pool =
3657  this->matrix_free->constraint_pool_end(indicator.second);
3658 
3659  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3660  operation.process_constraint(*dof_indices[v],
3661  *data_val,
3662  *src[0],
3663  value);
3664 
3665  operation.post_constraints(value,
3666  values_dofs[comp][ind_local][v]);
3667  ind_local++;
3668  }
3669 
3670  AssertIndexRange(ind_local, dofs_per_component + 1);
3671 
3672  // get the dof values past the last constraint
3673  for (; ind_local < dofs_per_component;
3674  ++dof_indices[v], ++ind_local)
3675  {
3676  AssertIndexRange(*dof_indices[v], src[0]->size());
3677  operation.process_dof(*dof_indices[v],
3678  *src[0],
3679  values_dofs[comp][ind_local][v]);
3680  }
3681 
3682  if (apply_constraints == true && comp + 1 < n_components)
3683  next_index_indicators =
3684  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
3685  }
3686  }
3687  }
3688 }
3689 
3690 
3691 
3692 template <int dim,
3693  int n_components_,
3694  typename Number,
3695  bool is_face,
3696  typename VectorizedArrayType>
3697 template <typename VectorType, typename VectorOperation>
3698 inline void
3701  const VectorOperation & operation,
3702  const std::array<VectorType *, n_components_> &src) const
3703 {
3704  Assert(!local_dof_indices.empty(), ExcNotInitialized());
3705 
3706  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3707  unsigned int index = this->first_selected_component * dofs_per_component;
3708  for (unsigned int comp = 0; comp < n_components; ++comp)
3709  {
3710  for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3711  {
3712  operation.process_empty(
3713  this->values_dofs[comp * dofs_per_component + i]);
3714  operation.process_dof_global(
3715  local_dof_indices[this->data->lexicographic_numbering[index]],
3716  *src[0],
3717  this->values_dofs[comp * dofs_per_component + i][0]);
3718  }
3719  }
3720 }
3721 
3722 
3723 
3724 template <int dim,
3725  int n_components_,
3726  typename Number,
3727  bool is_face,
3728  typename VectorizedArrayType>
3729 template <typename VectorType, typename VectorOperation>
3730 inline void
3733  const VectorOperation & operation,
3734  const std::array<VectorType *, n_components_> &src,
3735  const std::array<
3737  n_components_> & vectors_sm,
3738  const std::bitset<VectorizedArrayType::size()> &mask) const
3739 {
3740  // This functions processes the functions read_dof_values,
3741  // distribute_local_to_global, and set_dof_values with the same code for
3742  // contiguous cell indices (DG case). The distinction between these three
3743  // cases is made by the input VectorOperation that either reads values from
3744  // a vector and puts the data into the local data field or write local data
3745  // into the vector. Certain operations are no-ops for the given use case.
3746 
3747  std::integral_constant<bool,
3748  internal::is_vectorizable<VectorType, Number>::value>
3749  vector_selector;
3751  is_face ? this->dof_access_index :
3753  const unsigned int n_lanes = mask.count();
3754 
3755  const std::vector<unsigned int> &dof_indices_cont =
3756  this->dof_info->dof_indices_contiguous[ind];
3757 
3758  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3759  std::array<VectorizedArrayType *, n_components> values_dofs;
3760  for (unsigned int c = 0; c < n_components; ++c)
3761  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3762  c * dofs_per_component;
3763 
3765 
3766  // Simple case: We have contiguous storage, so we can simply copy out the
3767  // data
3768  if ((this->dof_info->index_storage_variants[ind][this->cell] ==
3770  interleaved_contiguous &&
3771  n_lanes == VectorizedArrayType::size()) &&
3772  !(is_face &&
3773  this->dof_access_index ==
3775  this->is_interior_face() == false) &&
3776  !(!is_face && !this->is_interior_face()))
3777  {
3778  const unsigned int dof_index =
3779  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3780  this->dof_info
3781  ->component_dof_indices_offset[this->active_fe_index]
3782  [this->first_selected_component] *
3783  VectorizedArrayType::size();
3784  if (n_components == 1 || this->n_fe_components == 1)
3785  for (unsigned int comp = 0; comp < n_components; ++comp)
3786  operation.process_dofs_vectorized(dofs_per_component,
3787  dof_index,
3788  *src[comp],
3789  values_dofs[comp],
3790  vector_selector);
3791  else
3792  operation.process_dofs_vectorized(dofs_per_component * n_components,
3793  dof_index,
3794  *src[0],
3795  values_dofs[0],
3796  vector_selector);
3797  return;
3798  }
3799 
3800  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3801  this->get_cell_or_face_ids();
3802 
3803  // More general case: Must go through the components one by one and apply
3804  // some transformations
3805  const unsigned int n_filled_lanes =
3806  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
3807 
3808  const bool is_ecl =
3809  (this->dof_access_index ==
3811  this->is_interior_face() == false) ||
3812  (!is_face && !this->is_interior_face());
3813 
3814  if (vectors_sm[0] != nullptr)
3815  {
3816  const auto compute_vector_ptrs = [&](const unsigned int comp) {
3817  std::array<typename VectorType::value_type *,
3818  VectorizedArrayType::size()>
3819  vector_ptrs = {};
3820 
3821  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3822  {
3823  if (mask[v] == false)
3824  {
3825  vector_ptrs[v] = nullptr;
3826  continue;
3827  }
3828 
3830  ExcNotImplemented());
3831  Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
3832  ExcIndexRange(
3833  ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
3834  Assert(cells[v] <
3835  this->dof_info->dof_indices_contiguous_sm[ind].size(),
3836  ExcIndexRange(
3837  cells[v],
3838  0,
3839  this->dof_info->dof_indices_contiguous_sm[ind].size()));
3840 
3841  const auto &temp =
3842  this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
3843 
3844  if (temp.first != numbers::invalid_unsigned_int)
3845  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3846  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3847  this->dof_info->component_dof_indices_offset
3848  [this->active_fe_index][this->first_selected_component]);
3849  else
3850  vector_ptrs[v] = nullptr;
3851  }
3852  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3853  ++v)
3854  vector_ptrs[v] = nullptr;
3855 
3856  return vector_ptrs;
3857  };
3858 
3859  if (n_filled_lanes == VectorizedArrayType::size() &&
3860  n_lanes == VectorizedArrayType::size() && !is_ecl)
3861  {
3862  if (n_components == 1 || this->n_fe_components == 1)
3863  {
3864  for (unsigned int comp = 0; comp < n_components; ++comp)
3865  {
3866  auto vector_ptrs = compute_vector_ptrs(comp);
3867  operation.process_dofs_vectorized_transpose(
3868  dofs_per_component,
3869  vector_ptrs,
3870  values_dofs[comp],
3871  vector_selector);
3872  }
3873  }
3874  else
3875  {
3876  auto vector_ptrs = compute_vector_ptrs(0);
3877  operation.process_dofs_vectorized_transpose(dofs_per_component *
3878  n_components,
3879  vector_ptrs,
3880  &values_dofs[0][0],
3881  vector_selector);
3882  }
3883  }
3884  else
3885  for (unsigned int comp = 0; comp < n_components; ++comp)
3886  {
3887  auto vector_ptrs = compute_vector_ptrs(
3888  (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3889 
3890  for (unsigned int i = 0; i < dofs_per_component; ++i)
3891  operation.process_empty(values_dofs[comp][i]);
3892 
3893  if (n_components == 1 || this->n_fe_components == 1)
3894  {
3895  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3896  if (mask[v] == true)
3897  for (unsigned int i = 0; i < dofs_per_component; ++i)
3898  operation.process_dof(vector_ptrs[v][i],
3899  values_dofs[comp][i][v]);
3900  }
3901  else
3902  {
3903  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3904  if (mask[v] == true)
3905  for (unsigned int i = 0; i < dofs_per_component; ++i)
3906  operation.process_dof(
3907  vector_ptrs[v][i + comp * dofs_per_component],
3908  values_dofs[comp][i][v]);
3909  }
3910  }
3911  return;
3912  }
3913 
3914  unsigned int dof_indices[VectorizedArrayType::size()];
3915 
3916  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3917  {
3919  dof_indices[v] =
3920  dof_indices_cont[cells[v]] +
3921  this->dof_info
3922  ->component_dof_indices_offset[this->active_fe_index]
3923  [this->first_selected_component] *
3924  this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
3925  }
3926 
3927  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3928  dof_indices[v] = numbers::invalid_unsigned_int;
3929 
3930  // In the case with contiguous cell indices, we know that there are no
3931  // constraints and that the indices within each element are contiguous
3932  if (n_filled_lanes == VectorizedArrayType::size() &&
3933  n_lanes == VectorizedArrayType::size() && !is_ecl)
3934  {
3935  if (this->dof_info->index_storage_variants[ind][this->cell] ==
3937  contiguous)
3938  {
3939  if (n_components == 1 || this->n_fe_components == 1)
3940  for (unsigned int comp = 0; comp < n_components; ++comp)
3941  operation.process_dofs_vectorized_transpose(dofs_per_component,
3942  dof_indices,
3943  *src[comp],
3944  values_dofs[comp],
3945  vector_selector);
3946  else
3947  operation.process_dofs_vectorized_transpose(dofs_per_component *
3948  n_components,
3949  dof_indices,
3950  *src[0],
3951  &values_dofs[0][0],
3952  vector_selector);
3953  }
3954  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
3956  interleaved_contiguous_strided)
3957  {
3958  if (n_components == 1 || this->n_fe_components == 1)
3959  for (unsigned int i = 0; i < dofs_per_component; ++i)
3960  {
3961  for (unsigned int comp = 0; comp < n_components; ++comp)
3962  operation.process_dof_gather(dof_indices,
3963  *src[comp],
3964  i * VectorizedArrayType::size(),
3965  values_dofs[comp][i],
3966  vector_selector);
3967  }
3968  else
3969  for (unsigned int comp = 0; comp < n_components; ++comp)
3970  for (unsigned int i = 0; i < dofs_per_component; ++i)
3971  {
3972  operation.process_dof_gather(dof_indices,
3973  *src[0],
3974  (comp * dofs_per_component + i) *
3975  VectorizedArrayType::size(),
3976  values_dofs[comp][i],
3977  vector_selector);
3978  }
3979  }
3980  else
3981  {
3982  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
3984  IndexStorageVariants::interleaved_contiguous_mixed_strides,
3985  ExcNotImplemented());
3986  const unsigned int *offsets =
3987  &this->dof_info->dof_indices_interleave_strides
3988  [ind][VectorizedArrayType::size() * this->cell];
3989  if (n_components == 1 || this->n_fe_components == 1)
3990  for (unsigned int i = 0; i < dofs_per_component; ++i)
3991  {
3992  for (unsigned int comp = 0; comp < n_components; ++comp)
3993  operation.process_dof_gather(dof_indices,
3994  *src[comp],
3995  0,
3996  values_dofs[comp][i],
3997  vector_selector);
3999  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4000  dof_indices[v] += offsets[v];
4001  }
4002  else
4003  for (unsigned int comp = 0; comp < n_components; ++comp)
4004  for (unsigned int i = 0; i < dofs_per_component; ++i)
4005  {
4006  operation.process_dof_gather(dof_indices,
4007  *src[0],
4008  0,
4009  values_dofs[comp][i],
4010  vector_selector);
4012  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4013  dof_indices[v] += offsets[v];
4014  }
4015  }
4016  }
4017  else
4018  for (unsigned int comp = 0; comp < n_components; ++comp)
4019  {
4020  for (unsigned int i = 0; i < dofs_per_component; ++i)
4021  operation.process_empty(values_dofs[comp][i]);
4022  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4024  contiguous)
4025  {
4026  if (n_components == 1 || this->n_fe_components == 1)
4027  {
4028  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4029  if (mask[v] == true)
4030  for (unsigned int i = 0; i < dofs_per_component; ++i)
4031  operation.process_dof(dof_indices[v] + i,
4032  *src[comp],
4033  values_dofs[comp][i][v]);
4034  }
4035  else
4036  {
4037  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4038  if (mask[v] == true)
4039  for (unsigned int i = 0; i < dofs_per_component; ++i)
4040  operation.process_dof(dof_indices[v] + i +
4041  comp * dofs_per_component,
4042  *src[0],
4043  values_dofs[comp][i][v]);
4044  }
4045  }
4046  else
4047  {
4048  const unsigned int *offsets =
4049  &this->dof_info->dof_indices_interleave_strides
4050  [ind][VectorizedArrayType::size() * this->cell];
4051  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4052  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4053  if (n_components == 1 || this->n_fe_components == 1)
4054  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4055  {
4056  if (mask[v] == true)
4057  for (unsigned int i = 0; i < dofs_per_component; ++i)
4058  operation.process_dof(dof_indices[v] + i * offsets[v],
4059  *src[comp],
4060  values_dofs[comp][i][v]);
4061  }
4062  else
4063  {
4064  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4065  if (mask[v] == true)
4066  for (unsigned int i = 0; i < dofs_per_component; ++i)
4067  operation.process_dof(dof_indices[v] +
4068  (i + comp * dofs_per_component) *
4069  offsets[v],
4070  *src[0],
4071  values_dofs[comp][i][v]);
4072  }
4073  }
4074  }
4075 }
4076 
4077 namespace internal
4078 {
4079  template <typename Number,
4080  typename VectorType,
4081  typename std::enable_if<!IsBlockVector<VectorType>::value,
4082  VectorType>::type * = nullptr>
4083  decltype(std::declval<VectorType>().begin())
4084  get_beginning(VectorType &vec)
4085  {
4086  return vec.begin();
4087  }
4088 
4089  template <typename Number,
4090  typename VectorType,
4091  typename std::enable_if<IsBlockVector<VectorType>::value,
4092  VectorType>::type * = nullptr>
4093  typename VectorType::value_type *
4094  get_beginning(VectorType &)
4095  {
4096  return nullptr;
4097  }
4098 
4099  template <typename VectorType,
4100  typename std::enable_if<has_shared_vector_data<VectorType>,
4101  VectorType>::type * = nullptr>
4102  const std::vector<ArrayView<const typename VectorType::value_type>> *
4103  get_shared_vector_data(VectorType & vec,
4104  const bool is_valid_mode_for_sm,
4105  const unsigned int active_fe_index,
4107  {
4108  // note: no hp is supported
4109  if (is_valid_mode_for_sm &&
4110  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4111  .size() > 0 &&
4112  active_fe_index == 0)
4113  return &vec.shared_vector_data();
4114  else
4115  return nullptr;
4116  }
4117 
4118  template <typename VectorType,
4119  typename std::enable_if<!has_shared_vector_data<VectorType>,
4120  VectorType>::type * = nullptr>
4121  const std::vector<ArrayView<const typename VectorType::value_type>> *
4122  get_shared_vector_data(VectorType &,
4123  const bool,
4124  const unsigned int,
4126  {
4127  return nullptr;
4128  }
4129 
4130  template <int n_components, typename VectorType>
4131  std::pair<
4132  std::array<typename internal::BlockVectorSelector<
4133  VectorType,
4134  IsBlockVector<VectorType>::value>::BaseVectorType *,
4135  n_components>,
4136  std::array<
4137  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4138  VectorType,
4139  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4140  n_components>>
4141  get_vector_data(VectorType & src,
4142  const unsigned int first_index,
4143  const bool is_valid_mode_for_sm,
4144  const unsigned int active_fe_index,
4146  {
4147  // select between block vectors and non-block vectors. Note that the number
4148  // of components is checked in the internal data
4149  std::pair<
4150  std::array<typename internal::BlockVectorSelector<
4151  VectorType,
4152  IsBlockVector<VectorType>::value>::BaseVectorType *,
4153  n_components>,
4154  std::array<
4155  const std::vector<
4156  ArrayView<const typename internal::BlockVectorSelector<
4157  VectorType,
4158  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4159  n_components>>
4160  src_data;
4161 
4162  for (unsigned int d = 0; d < n_components; ++d)
4163  src_data.first[d] = internal::BlockVectorSelector<
4164  VectorType,
4165  IsBlockVector<VectorType>::value>::get_vector_component(src,
4166  d +
4167  first_index);
4168 
4169  for (unsigned int d = 0; d < n_components; ++d)
4170  src_data.second[d] = get_shared_vector_data(
4171  const_cast<typename internal::BlockVectorSelector<
4172  typename std::remove_const<VectorType>::type,
4174  BaseVectorType &>(*src_data.first[d]),
4175  is_valid_mode_for_sm,
4176  active_fe_index,
4177  dof_info);
4178 
4179  return src_data;
4180  }
4181 } // namespace internal
4182 
4183 
4184 
4185 template <int dim,
4186  int n_components_,
4187  typename Number,
4188  bool is_face,
4189  typename VectorizedArrayType>
4190 inline void
4193 {
4194  if (this->dof_info == nullptr ||
4195  this->dof_info->hanging_node_constraint_masks.size() == 0 ||
4196  this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4197  this->dof_info->hanging_node_constraint_masks_comp
4198  [this->active_fe_index][this->first_selected_component] == false)
4199  return; // nothing to do with faces
4200 
4201  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4202  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4203  constraint_mask;
4204 
4205  bool hn_available = false;
4206 
4207  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
4208  this->get_cell_ids();
4209 
4210  for (unsigned int v = 0; v < n_lanes; ++v)
4211  {
4212  if (cells[v] == numbers::invalid_unsigned_int)
4213  {
4214  constraint_mask[v] = internal::MatrixFreeFunctions::
4216  continue;
4217  }
4218 
4219  const unsigned int cell_index = cells[v];
4220  const auto mask =
4221  this->dof_info->hanging_node_constraint_masks[cell_index];
4222  constraint_mask[v] = mask;
4223 
4224  hn_available |= (mask != internal::MatrixFreeFunctions::
4226  }
4227 
4228  if (hn_available == false)
4229  return; // no hanging node on cell batch -> nothing to do
4230 
4232  apply(n_components,
4233  this->data->data.front().fe_degree,
4234  this->get_shape_info(),
4235  transpose,
4236  constraint_mask,
4237  this->values_dofs);
4238 }
4239 
4240 
4241 
4242 template <int dim,
4243  int n_components_,
4244  typename Number,
4245  bool is_face,
4246  typename VectorizedArrayType>
4247 template <typename VectorType>
4248 inline void
4250  read_dof_values(const VectorType & src,
4251  const unsigned int first_index,
4252  const std::bitset<VectorizedArrayType::size()> &mask)
4253 {
4254  const auto src_data = internal::get_vector_data<n_components_>(
4255  src,
4256  first_index,
4257  this->dof_access_index ==
4259  this->active_fe_index,
4260  this->dof_info);
4261 
4263  read_write_operation(reader, src_data.first, src_data.second, mask, true);
4264 
4265  apply_hanging_node_constraints(false);
4266 
4267 # ifdef DEBUG
4268  this->dof_values_initialized = true;
4269 # endif
4270 }
4271 
4272 
4273 
4274 template <int dim,
4275  int n_components_,
4276  typename Number,
4277  bool is_face,
4278  typename VectorizedArrayType>
4279 template <typename VectorType>
4280 inline void
4282  read_dof_values_plain(const VectorType & src,
4283  const unsigned int first_index,
4284  const std::bitset<VectorizedArrayType::size()> &mask)
4285 {
4286  const auto src_data = internal::get_vector_data<n_components_>(
4287  src,
4288  first_index,
4289  this->dof_access_index ==
4291  this->active_fe_index,
4292  this->dof_info);
4293 
4295  read_write_operation(reader, src_data.first, src_data.second, mask, false);
4296 
4297 # ifdef DEBUG
4298  this->dof_values_initialized = true;
4299 # endif
4300 }
4301 
4302 
4303 
4304 template <int dim,
4305  int n_components_,
4306  typename Number,
4307  bool is_face,
4308  typename VectorizedArrayType>
4309 template <typename VectorType>
4310 inline void
4313  VectorType & dst,
4314  const unsigned int first_index,
4315  const std::bitset<VectorizedArrayType::size()> &mask) const
4316 {
4317 # ifdef DEBUG
4318  Assert(this->dof_values_initialized == true,
4320 # endif
4321 
4322  apply_hanging_node_constraints(true);
4323 
4324  const auto dst_data = internal::get_vector_data<n_components_>(
4325  dst,
4326  first_index,
4327  this->dof_access_index ==
4329  this->active_fe_index,
4330  this->dof_info);
4331 
4333  distributor;
4334  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4335 }
4336 
4337 
4338 
4339 template <int dim,
4340  int n_components_,
4341  typename Number,
4342  bool is_face,
4343  typename VectorizedArrayType>
4344 template <typename VectorType>
4345 inline void
4347  set_dof_values(VectorType & dst,
4348  const unsigned int first_index,
4349  const std::bitset<VectorizedArrayType::size()> &mask) const
4350 {
4351 # ifdef DEBUG
4352  Assert(this->dof_values_initialized == true,
4354 # endif
4355 
4356  const auto dst_data = internal::get_vector_data<n_components_>(
4357  dst,
4358  first_index,
4359  this->dof_access_index ==
4361  this->active_fe_index,
4362  this->dof_info);
4363 
4365  read_write_operation(setter, dst_data.first, dst_data.second, mask);
4366 }
4367 
4368 
4369 
4370 template <int dim,
4371  int n_components_,
4372  typename Number,
4373  bool is_face,
4374  typename VectorizedArrayType>
4375 template <typename VectorType>
4376 inline void
4379  VectorType & dst,
4380  const unsigned int first_index,
4381  const std::bitset<VectorizedArrayType::size()> &mask) const
4382 {
4383 # ifdef DEBUG
4384  Assert(this->dof_values_initialized == true,
4386 # endif
4387 
4388  const auto dst_data = internal::get_vector_data<n_components_>(
4389  dst,
4390  first_index,
4391  this->dof_access_index ==
4393  this->active_fe_index,
4394  this->dof_info);
4395 
4397  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4398 }
4399 
4400 
4401 
4402 /*------------------------------ access to data fields ----------------------*/
4403 
4404 
4405 
4406 template <int dim,
4407  int n_components_,
4408  typename Number,
4409  bool is_face,
4410  typename VectorizedArrayType>
4413  get_dof_value(const unsigned int dof) const
4414 {
4415  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4416  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4418  for (unsigned int comp = 0; comp < n_components; ++comp)
4419  return_value[comp] = this->values_dofs[comp * dofs + dof];
4420  return return_value;
4421 }
4422 
4423 
4424 
4425 template <int dim,
4426  int n_components_,
4427  typename Number,
4428  bool is_face,
4429  typename VectorizedArrayType>
4432  get_value(const unsigned int q_point) const
4433 {
4434 # ifdef DEBUG
4435  Assert(this->values_quad_initialized == true,
4437 # endif
4438 
4439  AssertIndexRange(q_point, this->n_quadrature_points);
4440  const std::size_t nqp = this->n_quadrature_points;
4442  for (unsigned int comp = 0; comp < n_components; ++comp)
4443  return_value[comp] = this->values_quad[comp * nqp + q_point];
4444  return return_value;
4445 }
4446 
4447 
4448 
4449 template <int dim,
4450  int n_components_,
4451  typename Number,
4452  bool is_face,
4453  typename VectorizedArrayType>
4454 inline DEAL_II_ALWAYS_INLINE
4457  get_gradient(const unsigned int q_point) const
4458 {
4459 # ifdef DEBUG
4460  Assert(this->gradients_quad_initialized == true,
4462 # endif
4463 
4464  AssertIndexRange(q_point, this->n_quadrature_points);
4465  Assert(this->jacobian != nullptr,
4467  "update_gradients"));
4468  const std::size_t nqp = this->n_quadrature_points;
4470 
4471  // Cartesian cell
4472  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4473  {
4474  for (unsigned int d = 0; d < dim; ++d)
4475  for (unsigned int comp = 0; comp < n_components; ++comp)
4476  grad_out[comp][d] =
4477  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4478  this->jacobian[0][d][d];
4479  }
4480  // cell with general/affine Jacobian
4481  else
4482  {
4484  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4485  q_point :
4486  0];
4487  for (unsigned int comp = 0; comp < n_components; ++comp)
4488  for (unsigned int d = 0; d < dim; ++d)
4489  {
4490  grad_out[comp][d] =
4491  jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4492  for (unsigned int e = 1; e < dim; ++e)
4493  grad_out[comp][d] +=
4494  jac[d][e] *
4495  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4496  }
4497  }
4498  return grad_out;
4499 }
4500 
4501 
4502 
4503 template <int dim,
4504  int n_components_,
4505  typename Number,
4506  bool is_face,
4507  typename VectorizedArrayType>
4510  get_normal_derivative(const unsigned int q_point) const
4511 {
4512  AssertIndexRange(q_point, this->n_quadrature_points);
4513 # ifdef DEBUG
4514  Assert(this->gradients_quad_initialized == true,
4516 # endif
4517 
4518  Assert(this->normal_x_jacobian != nullptr,
4520  "update_gradients"));
4521 
4522  const std::size_t nqp = this->n_quadrature_points;
4524 
4525  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4526  for (unsigned int comp = 0; comp < n_components; ++comp)
4527  grad_out[comp] =
4528  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4529  (this->normal_x_jacobian[0][dim - 1]);
4530  else
4531  {
4532  const std::size_t index =
4533  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4534  for (unsigned int comp = 0; comp < n_components; ++comp)
4535  {
4536  grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4537  this->normal_x_jacobian[index][0];
4538  for (unsigned int d = 1; d < dim; ++d)
4539  grad_out[comp] +=
4540  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4541  this->normal_x_jacobian[index][d];
4542  }
4543  }
4544  return grad_out;
4545 }
4546 
4547 
4548 
4549 namespace internal
4550 {
4551  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4552  // store the lower diagonal because of symmetry
4553  template <typename VectorizedArrayType>
4554  inline void
4555  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4556  const VectorizedArrayType *const hessians,
4557  const unsigned int,
4558  VectorizedArrayType (&tmp)[1][1])
4559  {
4560  tmp[0][0] = jac[0][0] * hessians[0];
4561  }
4562 
4563  template <typename VectorizedArrayType>
4564  inline void
4565  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4566  const VectorizedArrayType *const hessians,
4567  const unsigned int nqp,
4568  VectorizedArrayType (&tmp)[2][2])
4569  {
4570  for (unsigned int d = 0; d < 2; ++d)
4571  {
4572  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4573  tmp[1][d] =
4574  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4575  }
4576  }
4577 
4578  template <typename VectorizedArrayType>
4579  inline void
4580  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4581  const VectorizedArrayType *const hessians,
4582  const unsigned int nqp,
4583  VectorizedArrayType (&tmp)[3][3])
4584  {
4585  for (unsigned int d = 0; d < 3; ++d)
4586  {
4587  tmp[0][d] =
4588  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4589  jac[d][2] * hessians[4 * nqp]);
4590  tmp[1][d] =
4591  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4592  jac[d][2] * hessians[5 * nqp]);
4593  tmp[2][d] =
4594  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4595  jac[d][2] * hessians[2 * nqp]);
4596  }
4597  }
4598 } // namespace internal
4599 
4600 
4601 
4602 template <int dim,
4603  int n_components_,
4604  typename Number,
4605  bool is_face,
4606  typename VectorizedArrayType>
4609  get_hessian(const unsigned int q_point) const
4610 {
4611 # ifdef DEBUG
4612  Assert(this->hessians_quad_initialized == true,
4614 # endif
4615  AssertIndexRange(q_point, this->n_quadrature_points);
4616 
4617  Assert(this->jacobian != nullptr,
4619  "update_hessian"));
4621  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4622  0 :
4623  q_point];
4624 
4626 
4627  const std::size_t nqp = this->n_quadrature_points;
4628  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4629 
4630  // Cartesian cell
4631  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4632  {
4633  for (unsigned int comp = 0; comp < n_components; ++comp)
4634  {
4635  for (unsigned int d = 0; d < dim; ++d)
4636  hessian_out[comp][d][d] =
4637  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4638  (jac[d][d] * jac[d][d]);
4639  switch (dim)
4640  {
4641  case 1:
4642  break;
4643  case 2:
4644  hessian_out[comp][0][1] =
4645  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4646  (jac[0][0] * jac[1][1]);
4647  break;
4648  case 3:
4649  hessian_out[comp][0][1] =
4650  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4651  (jac[0][0] * jac[1][1]);
4652  hessian_out[comp][0][2] =
4653  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4654  (jac[0][0] * jac[2][2]);
4655  hessian_out[comp][1][2] =
4656  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4657  (jac[1][1] * jac[2][2]);
4658  break;
4659  default:
4660  Assert(false, ExcNotImplemented());
4661  }
4662  for (unsigned int d = 0; d < dim; ++d)
4663  for (unsigned int e = d + 1; e < dim; ++e)
4664  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4665  }
4666  }
4667  // cell with general Jacobian, but constant within the cell
4668  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4669  {
4670  for (unsigned int comp = 0; comp < n_components; ++comp)
4671  {
4672  VectorizedArrayType tmp[dim][dim];
4673  internal::hessian_unit_times_jac(
4674  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4675 
4676  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4677  for (unsigned int d = 0; d < dim; ++d)
4678  for (unsigned int e = d; e < dim; ++e)
4679  {
4680  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4681  for (unsigned int f = 1; f < dim; ++f)
4682  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4683  }
4684 
4685  // no J' * grad(u) part here because the Jacobian is constant
4686  // throughout the cell and hence, its derivative is zero
4687 
4688  // take symmetric part
4689  for (unsigned int d = 0; d < dim; ++d)
4690  for (unsigned int e = d + 1; e < dim; ++e)
4691  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4692  }
4693  }
4694  // cell with general Jacobian
4695  else
4696  {
4697  const auto &jac_grad = this->jacobian_gradients[q_point];
4698  for (unsigned int comp = 0; comp < n_components; ++comp)
4699  {
4700  VectorizedArrayType tmp[dim][dim];
4701  internal::hessian_unit_times_jac(
4702  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4703 
4704  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4705  for (unsigned int d = 0; d < dim; ++d)
4706  for (unsigned int e = d; e < dim; ++e)
4707  {
4708  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4709  for (unsigned int f = 1; f < dim; ++f)
4710  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4711  }
4712 
4713  // add diagonal part of J' * grad(u)
4714  for (unsigned int d = 0; d < dim; ++d)
4715  for (unsigned int e = 0; e < dim; ++e)
4716  hessian_out[comp][d][d] +=
4717  jac_grad[d][e] *
4718  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4719 
4720  // add off-diagonal part of J' * grad(u)
4721  for (unsigned int d = 0, count = dim; d < dim; ++d)
4722  for (unsigned int e = d + 1; e < dim; ++e, ++count)
4723  for (unsigned int f = 0; f < dim; ++f)
4724  hessian_out[comp][d][e] +=
4725  jac_grad[count][f] *
4726  this->gradients_quad[(comp * dim + f) * nqp + q_point];
4727 
4728  // take symmetric part
4729  for (unsigned int d = 0; d < dim; ++d)
4730  for (unsigned int e = d + 1; e < dim; ++e)
4731  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4732  }
4733  }
4734  return hessian_out;
4735 }
4736 
4737 
4738 
4739 template <int dim,
4740  int n_components_,
4741  typename Number,
4742  bool is_face,
4743  typename VectorizedArrayType>
4746  get_hessian_diagonal(const unsigned int q_point) const
4747 {
4748  Assert(!is_face, ExcNotImplemented());
4749 # ifdef DEBUG
4750  Assert(this->hessians_quad_initialized == true,
4752 # endif
4753  AssertIndexRange(q_point, this->n_quadrature_points);
4754 
4755  Assert(this->jacobian != nullptr, ExcNotImplemented());
4757  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4758  0 :
4759  q_point];
4760 
4761  const std::size_t nqp = this->n_quadrature_points;
4762  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4764 
4765  // Cartesian cell
4766  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4767  {
4768  for (unsigned int comp = 0; comp < n_components; ++comp)
4769  for (unsigned int d = 0; d < dim; ++d)
4770  hessian_out[comp][d] =
4771  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4772  (jac[d][d] * jac[d][d]);
4773  }
4774  // cell with general Jacobian, but constant within the cell
4775  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4776  {
4777  for (unsigned int comp = 0; comp < n_components; ++comp)
4778  {
4779  // compute laplacian before the gradient because it needs to access
4780  // unscaled gradient data
4781  VectorizedArrayType tmp[dim][dim];
4782  internal::hessian_unit_times_jac(
4783  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4784 
4785  // compute only the trace part of hessian, J * tmp = J *
4786  // hess_unit(u) * J^T
4787  for (unsigned int d = 0; d < dim; ++d)
4788  {
4789  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4790  for (unsigned int f = 1; f < dim; ++f)
4791  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4792  }
4793  }
4794  }
4795  // cell with general Jacobian
4796  else
4797  {
4798  const auto &jac_grad = this->jacobian_gradients[q_point];
4799  for (unsigned int comp = 0; comp < n_components; ++comp)
4800  {
4801  // compute laplacian before the gradient because it needs to access
4802  // unscaled gradient data
4803  VectorizedArrayType tmp[dim][dim];
4804  internal::hessian_unit_times_jac(
4805  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4806 
4807  // compute only the trace part of hessian, J * tmp = J *
4808  // hess_unit(u) * J^T
4809  for (unsigned int d = 0; d < dim; ++d)
4810  {
4811  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4812  for (unsigned int f = 1; f < dim; ++f)
4813  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4814  }
4815 
4816  for (unsigned int d = 0; d < dim; ++d)
4817  for (unsigned int e = 0; e < dim; ++e)
4818  hessian_out[comp][d] +=
4819  jac_grad[d][e] *
4820  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4821  }
4822  }
4823  return hessian_out;
4824 }
4825 
4826 
4827 
4828 template <int dim,
4829  int n_components_,
4830  typename Number,
4831  bool is_face,
4832  typename VectorizedArrayType>
4835  get_laplacian(const unsigned int q_point) const
4836 {
4837  Assert(is_face == false, ExcNotImplemented());
4838 # ifdef DEBUG
4839  Assert(this->hessians_quad_initialized == true,
4841 # endif
4842  AssertIndexRange(q_point, this->n_quadrature_points);
4843 
4845  const auto hess_diag = get_hessian_diagonal(q_point);
4846  for (unsigned int comp = 0; comp < n_components; ++comp)
4847  {
4848  laplacian_out[comp] = hess_diag[comp][0];
4849  for (unsigned int d = 1; d < dim; ++d)
4850  laplacian_out[comp] += hess_diag[comp][d];
4851  }
4852  return laplacian_out;
4853 }
4854 
4855 
4856 
4857 template <int dim,
4858  int n_components_,
4859  typename Number,
4860  bool is_face,
4861  typename VectorizedArrayType>
4862 inline DEAL_II_ALWAYS_INLINE void
4865  const unsigned int dof)
4866 {
4867 # ifdef DEBUG
4868  this->dof_values_initialized = true;
4869 # endif
4870  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4871  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4872  for (unsigned int comp = 0; comp < n_components; ++comp)
4873  this->values_dofs[comp * dofs + dof] = val_in[comp];
4874 }
4875 
4876 
4877 
4878 template <int dim,
4879  int n_components_,
4880  typename Number,
4881  bool is_face,
4882  typename VectorizedArrayType>
4883 inline DEAL_II_ALWAYS_INLINE void
4886  const unsigned int q_point)
4887 {
4888 # ifdef DEBUG
4889  Assert(this->is_reinitialized, ExcNotInitialized());
4890 # endif
4891  AssertIndexRange(q_point, this->n_quadrature_points);
4892  Assert(this->J_value != nullptr,
4894  "update_values"));
4895 # ifdef DEBUG
4896  this->values_quad_submitted = true;
4897 # endif
4898 
4899  const std::size_t nqp = this->n_quadrature_points;
4900  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4901  {
4902  const VectorizedArrayType JxW =
4903  this->J_value[0] * this->quadrature_weights[q_point];
4904  for (unsigned int comp = 0; comp < n_components; ++comp)
4905  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4906  }
4907  else
4908  {
4909  const VectorizedArrayType JxW = this->J_value[q_point];
4910  for (unsigned int comp = 0; comp < n_components; ++comp)
4911  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4912  }
4913 }
4914 
4915 
4916 
4917 template <int dim,
4918  int n_components_,
4919  typename Number,
4920  bool is_face,
4921  typename VectorizedArrayType>
4922 inline DEAL_II_ALWAYS_INLINE void
4925  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
4926  const unsigned int q_point)
4927 {
4928 # ifdef DEBUG
4929  Assert(this->is_reinitialized, ExcNotInitialized());
4930 # endif
4931  AssertIndexRange(q_point, this->n_quadrature_points);
4932  Assert(this->J_value != nullptr,
4934  "update_gradients"));
4935  Assert(this->jacobian != nullptr,
4937  "update_gradients"));
4938 # ifdef DEBUG
4939  this->gradients_quad_submitted = true;
4940 # endif
4941 
4942  const std::size_t nqp = this->n_quadrature_points;
4943  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4944  {
4945  const VectorizedArrayType JxW =
4946  this->J_value[0] * this->quadrature_weights[q_point];
4947  for (unsigned int d = 0; d < dim; ++d)
4948  {
4949  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
4950  for (unsigned int comp = 0; comp < n_components; ++comp)
4951  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
4952  grad_in[comp][d] * factor;
4953  }
4954  }
4955  else
4956  {
4958  this->cell_type > internal::MatrixFreeFunctions::affine ?
4959  this->jacobian[q_point] :
4960  this->jacobian[0];
4961  const VectorizedArrayType JxW =
4962  this->cell_type > internal::MatrixFreeFunctions::affine ?
4963  this->J_value[q_point] :
4964  this->J_value[0] * this->quadrature_weights[q_point];
4965  for (unsigned int comp = 0; comp < n_components; ++comp)
4966  for (unsigned int d = 0; d < dim; ++d)
4967  {
4968  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
4969  for (unsigned int e = 1; e < dim; ++e)
4970  new_val += (jac[e][d] * grad_in[comp][e]);
4971  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
4972  new_val * JxW;
4973  }
4974  }
4975 }
4976 
4977 
4978 
4979 template <int dim,
4980  int n_components_,
4981  typename Number,
4982  bool is_face,
4983  typename VectorizedArrayType>
4984 inline DEAL_II_ALWAYS_INLINE void
4988  const unsigned int q_point)
4989 {
4990  AssertIndexRange(q_point, this->n_quadrature_points);
4991  Assert(this->normal_x_jacobian != nullptr,
4993  "update_gradients"));
4994 # ifdef DEBUG
4995  this->gradients_quad_submitted = true;
4996 # endif
4997 
4998  const std::size_t nqp = this->n_quadrature_points;
4999  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5000  for (unsigned int comp = 0; comp < n_components; ++comp)
5001  {
5002  for (unsigned int d = 0; d < dim - 1; ++d)
5003  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5004  VectorizedArrayType();
5005  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5006  grad_in[comp] *
5007  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5008  this->quadrature_weights[q_point]);
5009  }
5010  else
5011  {
5012  const unsigned int index =
5013  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5015  this->normal_x_jacobian[index];
5016  for (unsigned int comp = 0; comp < n_components; ++comp)
5017  {
5018  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5019  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5020  factor = factor * this->quadrature_weights[q_point];
5021  for (unsigned int d = 0; d < dim; ++d)
5022  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5023  factor * jac[d];
5024  }
5025  }
5026 }
5027 
5028 
5029 
5030 template <int dim,
5031  int n_components_,
5032  typename Number,
5033  bool is_face,
5034  typename VectorizedArrayType>
5035 inline DEAL_II_ALWAYS_INLINE void
5038  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5039  hessian_in,
5040  const unsigned int q_point)
5041 {
5042 # ifdef DEBUG
5043  Assert(this->is_reinitialized, ExcNotInitialized());
5044 # endif
5045  AssertIndexRange(q_point, this->n_quadrature_points);
5046  Assert(this->J_value != nullptr,
5048  "update_hessians"));
5049  Assert(this->jacobian != nullptr,
5051  "update_hessians"));
5052 # ifdef DEBUG
5053  this->hessians_quad_submitted = true;
5054 # endif
5055 
5056  // compute hessian_unit = J^T * hessian_in(u) * J
5057  const std::size_t nqp = this->n_quadrature_points;
5058  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5059  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5060  {
5061  const VectorizedArrayType JxW =
5062  this->J_value[0] * this->quadrature_weights[q_point];
5063 
5064  // diagonal part
5065  for (unsigned int d = 0; d < dim; ++d)
5066  {
5067  const auto jac_d = this->jacobian[0][d][d];
5068  const VectorizedArrayType factor = jac_d * jac_d * JxW;
5069  for (unsigned int comp = 0; comp < n_components; ++comp)
5070  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5071  hessian_in[comp][d][d] * factor;
5072  }
5073 
5074  // off diagonal part
5075  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5076  for (unsigned int e = 0; e < d; ++e, ++off_dia)
5077  {
5078  const auto jac_d = this->jacobian[0][d][d];
5079  const auto jac_e = this->jacobian[0][e][e];
5080  const VectorizedArrayType factor = jac_d * jac_e * JxW;
5081  for (unsigned int comp = 0; comp < n_components; ++comp)
5082  this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5083  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5084  }
5085  }
5086  // cell with general Jacobian, but constant within the cell
5087  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5088  {
5089  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5090  const VectorizedArrayType JxW =
5091  this->J_value[0] * this->quadrature_weights[q_point];
5092  for (unsigned int comp = 0; comp < n_components; ++comp)
5093  {
5094  // 1. tmp = hessian_in(u) * J
5095  VectorizedArrayType tmp[dim][dim];
5096  for (unsigned int i = 0; i < dim; ++i)
5097  for (unsigned int j = 0; j < dim; ++j)
5098  {
5099  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5100  for (unsigned int k = 1; k < dim; ++k)
5101  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5102  }
5103 
5104  // 2. hessian_unit = J^T * tmp
5105  VectorizedArrayType tmp2[dim][dim];
5106  for (unsigned int i = 0; i < dim; ++i)
5107  for (unsigned int j = 0; j < dim; ++j)
5108  {
5109  tmp2[i][j] = jac[0][i] * tmp[0][j];
5110  for (unsigned int k = 1; k < dim; ++k)
5111  tmp2[i][j] += jac[k][i] * tmp[k][j];
5112  }
5113 
5114  // diagonal part
5115  for (unsigned int d = 0; d < dim; ++d)
5116  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5117  tmp2[d][d] * JxW;
5118 
5119  // off diagonal part
5120  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5121  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5122  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5123  (tmp2[d][e] + tmp2[e][d]) * JxW;
5124  }
5125  }
5126  else
5127  {
5128  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5129  const VectorizedArrayType JxW = this->J_value[q_point];
5130  const auto &jac_grad = this->jacobian_gradients[q_point];
5131  for (unsigned int comp = 0; comp < n_components; ++comp)
5132  {
5133  // 1. tmp = hessian_in(u) * J
5134  VectorizedArrayType tmp[dim][dim];
5135  for (unsigned int i = 0; i < dim; ++i)
5136  for (unsigned int j = 0; j < dim; ++j)
5137  {
5138  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5139  for (unsigned int k = 1; k < dim; ++k)
5140  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5141  }
5142 
5143  // 2. hessian_unit = J^T * tmp
5144  VectorizedArrayType tmp2[dim][dim];
5145  for (unsigned int i = 0; i < dim; ++i)
5146  for (unsigned int j = 0; j < dim; ++j)
5147  {
5148  tmp2[i][j] = jac[0][i] * tmp[0][j];
5149  for (unsigned int k = 1; k < dim; ++k)
5150  tmp2[i][j] += jac[k][i] * tmp[k][j];
5151  }
5152 
5153  // diagonal part
5154  for (unsigned int d = 0; d < dim; ++d)
5155  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5156  tmp2[d][d] * JxW;
5157 
5158  // off diagonal part
5159  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5160  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5161  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5162  (tmp2[d][e] + tmp2[e][d]) * JxW;
5163 
5164  // 3. gradient_unit = J' ** hessian_in
5165  for (unsigned int d = 0; d < dim; ++d)
5166  {
5167  VectorizedArrayType sum = 0;
5168  for (unsigned int e = 0; e < dim; ++e)
5169  sum += hessian_in[comp][e][e] * jac_grad[e][d];
5170  for (unsigned int e = 0, count = dim; e < dim; ++e)
5171  for (unsigned int f = e + 1; f < dim; ++f, ++count)
5172  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5173  jac_grad[count][d];
5174  this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
5175  q_point] = sum * JxW;
5176  }
5177  }
5178  }
5179 }
5180 
5181 
5182 
5183 template <int dim,
5184  int n_components_,
5185  typename Number,
5186  bool is_face,
5187  typename VectorizedArrayType>
5190  integrate_value() const
5191 {
5192 # ifdef DEBUG
5193  Assert(this->is_reinitialized, ExcNotInitialized());
5194  Assert(this->values_quad_submitted == true,
5196 # endif
5197 
5199  const std::size_t nqp = this->n_quadrature_points;
5200  for (unsigned int q = 0; q < nqp; ++q)
5201  for (unsigned int comp = 0; comp < n_components; ++comp)
5202  return_value[comp] += this->values_quad[comp * nqp + q];
5203  return (return_value);
5204 }
5205 
5206 
5207 
5208 /*----------------------- FEEvaluationAccess --------------------------------*/
5209 
5210 
5211 template <int dim,
5212  int n_components_,
5213  typename Number,
5214  bool is_face,
5215  typename VectorizedArrayType>
5216 inline FEEvaluationAccess<dim,
5217  n_components_,
5218  Number,
5219  is_face,
5220  VectorizedArrayType>::
5221  FEEvaluationAccess(
5223  const unsigned int dof_no,
5224  const unsigned int first_selected_component,
5225  const unsigned int quad_no,
5226  const unsigned int fe_degree,
5227  const unsigned int n_q_points,
5228  const bool is_interior_face,
5229  const unsigned int active_fe_index,
5230  const unsigned int active_quad_index,
5231  const unsigned int face_type)
5232  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5233  matrix_free,
5234  dof_no,
5235  first_selected_component,
5236  quad_no,
5237  fe_degree,
5238  n_q_points,
5239  is_interior_face,
5240  active_fe_index,
5241  active_quad_index,
5242  face_type)
5243 {}
5244 
5245 
5246 
5247 template <int dim,
5248  int n_components_,
5249  typename Number,
5250  bool is_face,
5251  typename VectorizedArrayType>
5252 inline FEEvaluationAccess<dim,
5253  n_components_,
5254  Number,
5255  is_face,
5256  VectorizedArrayType>::
5257  FEEvaluationAccess(
5258  const Mapping<dim> & mapping,
5259  const FiniteElement<dim> &fe,
5260  const Quadrature<1> & quadrature,
5261  const UpdateFlags update_flags,
5262  const unsigned int first_selected_component,
5264  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5265  mapping,
5266  fe,
5267  quadrature,
5268  update_flags,
5269  first_selected_component,
5270  other)
5271 {}
5272 
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(const FEEvaluationAccess<dim,
5286  n_components_,
5287  Number,
5288  is_face,
5289  VectorizedArrayType> &other)
5290  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5291  other)
5292 {}
5293 
5294 
5295 
5296 template <int dim,
5297  int n_components_,
5298  typename Number,
5299  bool is_face,
5300  typename VectorizedArrayType>
5301 inline FEEvaluationAccess<dim,
5302  n_components_,
5303  Number,
5304  is_face,
5305  VectorizedArrayType> &
5307 operator=(const FEEvaluationAccess<dim,
5308  n_components_,
5309  Number,
5310  is_face,
5311  VectorizedArrayType> &other)
5312 {
5313  this->FEEvaluationBase<dim,
5314  n_components_,
5315  Number,
5316  is_face,
5317  VectorizedArrayType>::operator=(other);
5318  return *this;
5319 }
5320 
5321 
5322 
5323 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5324 
5325 
5326 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5330  const unsigned int dof_no,
5331  const unsigned int first_selected_component,
5332  const unsigned int quad_no,
5333  const unsigned int fe_degree,
5334  const unsigned int n_q_points,
5335  const bool is_interior_face,
5336  const unsigned int active_fe_index,
5337  const unsigned int active_quad_index,
5338  const unsigned int face_type)
5339  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5340  matrix_free,
5341  dof_no,
5342  first_selected_component,
5343  quad_no,
5344  fe_degree,
5345  n_q_points,
5346  is_interior_face,
5347  active_fe_index,
5348  active_quad_index,
5349  face_type)
5350 {}
5351 
5352 
5353 
5354 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5357  const Mapping<dim> & mapping,
5358  const FiniteElement<dim> &fe,
5359  const Quadrature<1> & quadrature,
5360  const UpdateFlags update_flags,
5361  const unsigned int first_selected_component,
5363  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5364  mapping,
5365  fe,
5366  quadrature,
5367  update_flags,
5368  first_selected_component,
5369  other)
5370 {}
5371 
5372 
5373 
5374 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5378  &other)
5379  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5380 {}
5381 
5382 
5383 
5384 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5388 {
5389  this
5390  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5391  other);
5392  return *this;
5393 }
5394 
5395 
5396 
5397 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5398 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5400  const unsigned int dof) const
5401 {
5402  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5403  return this->values_dofs[dof];
5404 }
5405 
5406 
5407 
5408 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5409 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5411  const unsigned int q_point) const
5412 {
5413 # ifdef DEBUG
5414  Assert(this->values_quad_initialized == true,
5416 # endif
5417  AssertIndexRange(q_point, this->n_quadrature_points);
5418  return this->values_quad[q_point];
5419 }
5420 
5421 
5422 
5423 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5424 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5426  get_normal_derivative(const unsigned int q_point) const
5427 {
5428  return BaseClass::get_normal_derivative(q_point)[0];
5429 }
5430 
5431 
5432 
5433 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5436  const unsigned int q_point) const
5437 {
5438  // could use the base class gradient, but that involves too many expensive
5439  // initialization operations on tensors
5440 
5441 # ifdef DEBUG
5442  Assert(this->gradients_quad_initialized == true,
5444 # endif
5445  AssertIndexRange(q_point, this->n_quadrature_points);
5446 
5447  Assert(this->jacobian != nullptr,
5449  "update_gradients"));
5450 
5452 
5453  const std::size_t nqp = this->n_quadrature_points;
5454  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5455  {
5456  for (unsigned int d = 0; d < dim; ++d)
5457  grad_out[d] =
5458  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5459  }
5460  // cell with general/affine Jacobian
5461  else
5462  {
5464  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5465  q_point :
5466  0];
5467  for (unsigned int d = 0; d < dim; ++d)
5468  {
5469  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5470  for (unsigned int e = 1; e < dim; ++e)
5471  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5472  }
5473  }
5474  return grad_out;
5475 }
5476 
5477 
5478 
5479 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5482  const unsigned int q_point) const
5483 {
5484  return BaseClass::get_hessian(q_point)[0];
5485 }
5486 
5487 
5488 
5489 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5492  get_hessian_diagonal(const unsigned int q_point) const
5493 {
5494  return BaseClass::get_hessian_diagonal(q_point)[0];
5495 }
5496 
5497 
5498 
5499 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5500 inline VectorizedArrayType
5502  const unsigned int q_point) const
5503 {
5504  return BaseClass::get_laplacian(q_point)[0];
5505 }
5506 
5507 
5508 
5509 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5510 inline void DEAL_II_ALWAYS_INLINE
5512  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5513 {
5514 # ifdef DEBUG
5515  this->dof_values_initialized = true;
5516  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5517 # endif
5518  this->values_dofs[dof] = val_in;
5519 }
5520 
5521 
5522 
5523 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5524 inline void DEAL_II_ALWAYS_INLINE
5526  const VectorizedArrayType val_in,
5527  const unsigned int q_point)
5528 {
5529 # ifdef DEBUG
5530  Assert(this->is_reinitialized, ExcNotInitialized());
5531 # endif
5532  AssertIndexRange(q_point, this->n_quadrature_points);
5533  Assert(this->J_value != nullptr,
5535  "update_value"));
5536 # ifdef DEBUG
5537  this->values_quad_submitted = true;
5538 # endif
5539 
5540  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5541  {
5542  const VectorizedArrayType JxW =
5543  this->J_value[0] * this->quadrature_weights[q_point];
5544  this->values_quad[q_point] = val_in * JxW;
5545  }
5546  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5547  {
5548  this->values_quad[q_point] = val_in * this->J_value[q_point];
5549  }
5550 }
5551 
5552 
5553 
5554 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5555 inline DEAL_II_ALWAYS_INLINE void
5557  const Tensor<1, 1, VectorizedArrayType> val_in,
5558  const unsigned int q_point)
5559 {
5560  submit_value(val_in[0], q_point);
5561 }
5562 
5563 
5564 
5565 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5566 inline DEAL_II_ALWAYS_INLINE void
5568  submit_normal_derivative(const VectorizedArrayType grad_in,
5569  const unsigned int q_point)
5570 {
5572  grad[0] = grad_in;
5573  BaseClass::submit_normal_derivative(grad, q_point);
5574 }
5575 
5576 
5577 
5578 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5579 inline DEAL_II_ALWAYS_INLINE void
5582  const unsigned int q_point)
5583 {
5584 # ifdef DEBUG
5585  Assert(this->is_reinitialized, ExcNotInitialized());
5586 # endif
5587  AssertIndexRange(q_point, this->n_quadrature_points);
5588  Assert(this->J_value != nullptr,
5590  "update_gradients"));
5591  Assert(this->jacobian != nullptr,
5593  "update_gradients"));
5594 # ifdef DEBUG
5595  this->gradients_quad_submitted = true;
5596 # endif
5597 
5598  const std::size_t nqp = this->n_quadrature_points;
5599  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5600  {
5601  const VectorizedArrayType JxW =
5602  this->J_value[0] * this->quadrature_weights[q_point];
5603  for (unsigned int d = 0; d < dim; ++d)
5604  this->gradients_quad[d * nqp + q_point] =
5605  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5606  }
5607  // general/affine cell type
5608  else
5609  {
5611  this->cell_type > internal::MatrixFreeFunctions::affine ?
5612  this->jacobian[q_point] :
5613  this->jacobian[0];
5614  const VectorizedArrayType JxW =
5615  this->cell_type > internal::MatrixFreeFunctions::affine ?
5616  this->J_value[q_point] :
5617  this->J_value[0] * this->quadrature_weights[q_point];
5618  for (unsigned int d = 0; d < dim; ++d)
5619  {
5620  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5621  for (unsigned int e = 1; e < dim; ++e)
5622  new_val += jac[e][d] * grad_in[e];
5623  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
5624  }
5625  }
5626 }
5627 
5628 
5629 
5630 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5631 inline DEAL_II_ALWAYS_INLINE void
5634  const unsigned int q_point)
5635 {
5637  hessian[0] = hessian_in;
5638  BaseClass::submit_hessian(hessian, q_point);
5639 }
5640 
5641 
5642 
5643 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5644 inline VectorizedArrayType
5646  integrate_value() const
5647 {
5648  return BaseClass::integrate_value()[0];
5649 }
5650 
5651 
5652 
5653 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5654 
5655 
5656 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5660  const unsigned int dof_no,
5661  const unsigned int first_selected_component,
5662  const unsigned int quad_no,
5663  const unsigned int fe_degree,
5664  const unsigned int n_q_points,
5665  const bool is_interior_face,
5666  const unsigned int active_fe_index,
5667  const unsigned int active_quad_index,
5668  const unsigned int face_type)
5669  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5670  matrix_free,
5671  dof_no,
5672  first_selected_component,
5673  quad_no,
5674  fe_degree,
5675  n_q_points,
5676  is_interior_face,
5677  active_fe_index,
5678  active_quad_index,
5679  face_type)
5680 {}
5681 
5682 
5683 
5684 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5687  const Mapping<dim> & mapping,
5688  const FiniteElement<dim> &fe,
5689  const Quadrature<1> & quadrature,
5690  const UpdateFlags update_flags,
5691  const unsigned int first_selected_component,
5693  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5694  mapping,
5695  fe,
5696  quadrature,
5697  update_flags,
5698  first_selected_component,
5699  other)
5700 {}
5701 
5702 
5703 
5704 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5708  &other)
5709  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5710 {}
5711 
5712 
5713 
5714 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5718  &other)
5719 {
5721  operator=(other);
5722  return *this;
5723 }
5724 
5725 
5726 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5729  const unsigned int q_point) const
5730 {
5731  if (this->data->element_type ==
5733  {
5734  // Piola transform is required
5735 # ifdef DEBUG
5736  Assert(this->values_quad_initialized == true,
5738 # endif
5739 
5740  AssertIndexRange(q_point, this->n_quadrature_points);
5741  Assert(this->J_value != nullptr,
5743  "update_values"));
5744  const std::size_t nqp = this->n_quadrature_points;
5746 
5747  if (!is_face &&
5748  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5749  {
5750  // Cartesian cell
5751  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5752  const VectorizedArrayType inv_det =
5753  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5754  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5755  this->jacobian[0][2][2];
5756 
5757  // J * u * det(J^-1)
5758  for (unsigned int comp = 0; comp < n_components; ++comp)
5759  value_out[comp] = this->values_quad[comp * nqp + q_point] *
5760  jac[comp][comp] * inv_det;
5761  }
5762  else
5763  {
5764  // Affine or general cell
5765  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5766  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5767  this->jacobian[q_point] :
5768  this->jacobian[0];
5770  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5771  transpose(invert(inv_t_jac)) :
5772  this->jacobian[1];
5773 
5774  // Derivatives are reordered for faces. Need to take this into account
5775  const VectorizedArrayType inv_det =
5776  (is_face && dim == 2 && this->get_face_no() < 2) ?
5777  -determinant(inv_t_jac) :
5778  determinant(inv_t_jac);
5779  // J * u * det(J^-1)
5780  for (unsigned int comp = 0; comp < n_components; ++comp)
5781  {
5782  value_out[comp] =
5783  this->values_quad[q_point] * jac[comp][0] * inv_det;
5784  for (unsigned int e = 1; e < dim; ++e)
5785  value_out[comp] +=
5786  this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
5787  }
5788  }
5789  return value_out;
5790  }
5791  else
5792  {
5793  // No Piola needed
5794  return BaseClass::get_value(q_point);
5795  }
5796 }
5797 
5798 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5801  get_gradient(const unsigned int q_point) const
5802 {
5803  if (this->data->element_type ==
5805  {
5806  // Piola transform is required
5807 # ifdef DEBUG
5808  Assert(this->gradients_quad_initialized == true,
5810 # endif
5811 
5812  AssertIndexRange(q_point, this->n_quadrature_points);
5813  Assert(this->jacobian != nullptr,
5815  "update_gradients"));
5816  const std::size_t nqp = this->n_quadrature_points;
5818 
5819  if (!is_face &&
5820  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5821  {
5822  // Cartesian cell
5823  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5824  this->jacobian[0];
5825  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5826  const VectorizedArrayType inv_det =
5827  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5828  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5829  this->jacobian[0][2][2];
5830 
5831  // J * grad_quad * J^-1 * det(J^-1)
5832  for (unsigned int d = 0; d < dim; ++d)
5833  for (unsigned int comp = 0; comp < n_components; ++comp)
5834  grad_out[comp][d] =
5835  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
5836  inv_t_jac[d][d] * jac[comp][comp] * inv_det;
5837  }
5838  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5839  {
5840  // Affine cell
5841  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5842  this->jacobian[0];
5843  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5844 
5845  // Derivatives are reordered for faces. Need to take this into account
5846  const VectorizedArrayType inv_det =
5847  (is_face && dim == 2 && this->get_face_no() < 2) ?
5848  -determinant(inv_t_jac) :
5849  determinant(inv_t_jac);
5850 
5851  VectorizedArrayType tmp;
5852  // J * grad_quad * J^-1 * det(J^-1)
5853  for (unsigned int comp = 0; comp < n_components; ++comp)
5854  for (unsigned int d = 0; d < dim; ++d)
5855  {
5856  tmp = 0;
5857  for (unsigned int f = 0; f < dim; ++f)
5858  for (unsigned int e = 0; e < dim; ++e)
5859  tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
5860  this->gradients_quad[(f * dim + e) * nqp + q_point];
5861 
5862  grad_out[comp][d] = tmp;
5863  }
5864  }
5865  else
5866  {
5867  // General cell
5868  // Here we need the jacobian gradient and not the inverse which is
5869  // stored in this->jacobian_gradients
5870  AssertThrow(false, ExcNotImplemented());
5871  }
5872  return grad_out;
5873  }
5874  else
5875  {
5876  return BaseClass::get_gradient(q_point);
5877  }
5878 }
5879 
5880 
5881 
5882 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5883 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5885  get_divergence(const unsigned int q_point) const
5886 {
5887 # ifdef DEBUG
5888  Assert(this->gradients_quad_initialized == true,
5890 # endif
5891  AssertIndexRange(q_point, this->n_quadrature_points);
5892  Assert(this->jacobian != nullptr,
5894  "update_gradients"));
5895 
5896  VectorizedArrayType divergence;
5897  const std::size_t nqp = this->n_quadrature_points;
5898 
5899  if (this->data->element_type ==
5901  {
5902  if (!is_face &&
5903  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5904  {
5905  // Cartesian cell
5906  const VectorizedArrayType inv_det =
5907  (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5908  this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5909  this->jacobian[0][2][2];
5910 
5911  // div * det(J^-1)
5912  divergence = this->gradients_quad[q_point] * inv_det;
5913  for (unsigned int d = 1; d < dim; ++d)
5914  divergence +=
5915  this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5916  }
5917  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5918  {
5919  // Affine cell
5920  // Derivatives are reordered for faces. Need to take this into account
5921  const VectorizedArrayType inv_det =
5922  (is_face && dim == 2 && this->get_face_no() < 2) ?
5923  -determinant(this->jacobian[0]) :
5924  determinant(this->jacobian[0]);
5925 
5926  // div * det(J^-1)
5927  divergence = this->gradients_quad[q_point] * inv_det;
5928  for (unsigned int d = 1; d < dim; ++d)
5929  divergence +=
5930  this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5931  }
5932  else
5933  {
5934  // General cell
5935  Assert(false, ExcNotImplemented());
5936  }
5937  }
5938  else
5939  {
5940  if (!is_face &&
5941  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5942  {
5943  // Cartesian cell
5944  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
5945  for (unsigned int d = 1; d < dim; ++d)
5946  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
5947  this->jacobian[0][d][d];
5948  }
5949  else
5950  {
5951  // cell with general/constant Jacobian
5953  this->cell_type == internal::MatrixFreeFunctions::general ?
5954  this->jacobian[q_point] :
5955  this->jacobian[0];
5956  divergence = jac[0][0] * this->gradients_quad[q_point];
5957  for (unsigned int e = 1; e < dim; ++e)
5958  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
5959  for (unsigned int d = 1; d < dim; ++d)
5960  for (unsigned int e = 0; e < dim; ++e)
5961  divergence +=
5962  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
5963  }
5964  }
5965  return divergence;
5966 }
5967 
5968 
5969 
5970 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5973  get_symmetric_gradient(const unsigned int q_point) const
5974 {
5975  // copy from generic function into dim-specialization function
5976  const auto grad = get_gradient(q_point);
5977  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5978  VectorizedArrayType half = Number(0.5);
5979  for (unsigned int d = 0; d < dim; ++d)
5980  symmetrized[d] = grad[d][d];
5981  switch (dim)
5982  {
5983  case 1:
5984  break;
5985  case 2:
5986  symmetrized[2] = grad[0][1] + grad[1][0];
5987  symmetrized[2] *= half;
5988  break;
5989  case 3:
5990  symmetrized[3] = grad[0][1] + grad[1][0];
5991  symmetrized[3] *= half;
5992  symmetrized[4] = grad[0][2] + grad[2][0];
5993  symmetrized[4] *= half;
5994  symmetrized[5] = grad[1][2] + grad[2][1];
5995  symmetrized[5] *= half;
5996  break;
5997  default:
5998  Assert(false, ExcNotImplemented());
5999  }
6001 }
6002 
6003 
6004 
6005 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6006 inline DEAL_II_ALWAYS_INLINE
6007  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6009  const unsigned int q_point) const
6010 {
6011  // copy from generic function into dim-specialization function
6012  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6013  Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6014  switch (dim)
6015  {
6016  case 1:
6017  Assert(false,
6018  ExcMessage(
6019  "Computing the curl in 1d is not a useful operation"));
6020  break;
6021  case 2:
6022  curl[0] = grad[1][0] - grad[0][1];
6023  break;
6024  case 3:
6025  curl[0] = grad[2][1] - grad[1][2];
6026  curl[1] = grad[0][2] - grad[2][0];
6027  curl[2] = grad[1][0] - grad[0][1];
6028  break;
6029  default:
6030  Assert(false, ExcNotImplemented());
6031  }
6032  return curl;
6033 }
6034 
6035 
6036 
6037 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6040  get_hessian_diagonal(const unsigned int q_point) const
6041 {
6042  return BaseClass::get_hessian_diagonal(q_point);
6043 }
6044 
6045 
6046 
6047 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6050  const unsigned int q_point) const
6051 {
6052 # ifdef DEBUG
6053  Assert(this->hessians_quad_initialized == true,
6055 # endif
6056  AssertIndexRange(q_point, this->n_quadrature_points);
6057  return BaseClass::get_hessian(q_point);
6058 }
6059 
6060 
6061 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6062 inline DEAL_II_ALWAYS_INLINE void
6065  const unsigned int q_point)
6066 {
6067  if (this->data->element_type ==
6069  {
6070  // Piola transform is required
6071  AssertIndexRange(q_point, this->n_quadrature_points);
6072  Assert(this->J_value != nullptr,
6074  "update_value"));
6075 # ifdef DEBUG
6076  Assert(this->is_reinitialized, ExcNotInitialized());
6077  this->values_quad_submitted = true;
6078 # endif
6079 
6080  const std::size_t nqp = this->n_quadrature_points;
6081  if (!is_face &&
6082  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6083  {
6084  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6085  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6086 
6087  for (unsigned int comp = 0; comp < n_components; ++comp)
6088  this->values_quad[comp * nqp + q_point] =
6089  val_in[comp] * weight * jac[comp][comp];
6090  }
6091  else
6092  {
6093  // Affine or general cell
6094  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6095  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6096  this->jacobian[q_point] :
6097  this->jacobian[0];
6099  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6100  transpose(invert(inv_t_jac)) :
6101  this->jacobian[1];
6102 
6103  // Derivatives are reordered for faces. Need to take this into account
6104  // and 1/inv_det != J_value for faces
6105  const VectorizedArrayType fac =
6106  (!is_face) ?
6107  this->quadrature_weights[q_point] :
6108  (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6109  this->J_value[q_point] :
6110  this->J_value[0] * this->quadrature_weights[q_point]) *
6111  ((dim == 2 && this->get_face_no() < 2) ?
6112  -determinant(inv_t_jac) :
6113  determinant(inv_t_jac)));
6114 
6115  // J^T * u * factor
6116  for (unsigned int comp = 0; comp < n_components; ++comp)
6117  {
6118  this->values_quad[comp * nqp + q_point] =
6119  val_in[0] * jac[0][comp] * fac;
6120  for (unsigned int e = 1; e < dim; ++e)
6121  this->values_quad[comp * nqp + q_point] +=
6122  val_in[e] * jac[e][comp] * fac;
6123  }
6124  }
6125  }
6126  else
6127  {
6128  // No Piola transform
6129  BaseClass::submit_value(val_in, q_point);
6130  }
6131 }
6132 
6133 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6134 inline DEAL_II_ALWAYS_INLINE void
6137  const unsigned int q_point)
6138 {
6139  if (this->data->element_type ==
6141  {
6142  // Piola transform is required
6143 
6144 # ifdef DEBUG
6145  Assert(this->is_reinitialized, ExcNotInitialized());
6146 # endif
6147  AssertIndexRange(q_point, this->n_quadrature_points);
6148  Assert(this->J_value != nullptr,
6150  "update_gradients"));
6151  Assert(this->jacobian != nullptr,
6153  "update_gradients"));
6154 # ifdef DEBUG
6155  this->gradients_quad_submitted = true;
6156 # endif
6157 
6158  const std::size_t nqp = this->n_quadrature_points;
6159  if (!is_face &&
6160  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6161  {
6162  // Cartesian cell
6163  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6164  this->jacobian[0];
6165  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6166  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6167  for (unsigned int d = 0; d < dim; ++d)
6168  for (unsigned int comp = 0; comp < n_components; ++comp)
6169  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6170  grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
6171  }
6172  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6173  {
6174  // Affine cell
6175  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6176  this->jacobian[0];
6177  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6178 
6179  // Derivatives are reordered for faces. Need to take this into account
6180  // and 1/inv_det != J_value for faces
6181  const VectorizedArrayType fac =
6182  (!is_face) ? this->quadrature_weights[q_point] :
6183  this->J_value[0] * this->quadrature_weights[q_point] *
6184  ((dim == 2 && this->get_face_no() < 2) ?
6185  -determinant(inv_t_jac) :
6186  determinant(inv_t_jac));
6187 
6188  // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6189  for (unsigned int comp = 0; comp < n_components; ++comp)
6190  for (unsigned int d = 0; d < dim; ++d)
6191  {
6192  VectorizedArrayType tmp = 0;
6193  for (unsigned int f = 0; f < dim; ++f)
6194  for (unsigned int e = 0; e < dim; ++e)
6195  tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
6196 
6197  this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
6198  }
6199  }
6200  else
6201  {
6202  // General cell
6203  AssertThrow(false, ExcNotImplemented());
6204  }
6205  }
6206  else
6207  {
6208  BaseClass::submit_gradient(grad_in, q_point);
6209  }
6210 }
6211 
6212 
6213 
6214 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6215 inline DEAL_II_ALWAYS_INLINE void
6218  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6219  const unsigned int q_point)
6220 {
6221  if (this->data->element_type ==
6223  {
6224  // Piola transform is required
6225  const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
6227  submit_gradient(grad, q_point);
6228  }
6229  else
6230  {
6231  BaseClass::submit_gradient(grad_in, q_point);
6232  }
6233 }
6234 
6235 
6236 
6237 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6238 inline DEAL_II_ALWAYS_INLINE void
6240  submit_divergence(const VectorizedArrayType div_in,
6241  const unsigned int q_point)
6242 {
6243 # ifdef DEBUG
6244  Assert(this->is_reinitialized, ExcNotInitialized());
6245 # endif
6246  AssertIndexRange(q_point, this->n_quadrature_points);
6247  Assert(this->J_value != nullptr,
6249  "update_gradients"));
6250  Assert(this->jacobian != nullptr,
6252  "update_gradients"));
6253 # ifdef DEBUG
6254  this->gradients_quad_submitted = true;
6255 # endif
6256 
6257  const std::size_t nqp = this->n_quadrature_points;
6258  if (this->data->element_type ==
6260  {
6261  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6262  {
6263  // Affine cell
6264 
6265  // Derivatives are reordered for faces. Need to take this into account
6266  // and 1/inv_det != J_value for faces
6267  const VectorizedArrayType fac =
6268  ((!is_face) ?
6269  1 :
6270  this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
6271  -determinant(this->jacobian[0]) :
6272  determinant(this->jacobian[0]))) *
6273  this->quadrature_weights[q_point] * div_in;
6274 
6275  for (unsigned int d = 0; d < dim; ++d)
6276  {
6277  this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
6278  for (unsigned int e = d + 1; e < dim; ++e)
6279  {
6280  this->gradients_quad[(dim * d + e) * nqp + q_point] =
6281  VectorizedArrayType();
6282  this->gradients_quad[(dim * e + d) * nqp + q_point] =
6283  VectorizedArrayType();
6284  }
6285  }
6286  }
6287  else
6288  {
6289  // General cell
6290  AssertThrow(false, ExcNotImplemented());
6291  }
6292  }
6293  else
6294  {
6295  if (!is_face &&
6296  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6297  {
6298  const VectorizedArrayType fac =
6299  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6300  for (unsigned int d = 0; d < dim; ++d)
6301  {
6302  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6303  (fac * this->jacobian[0][d][d]);
6304  for (unsigned int e = d + 1; e < dim; ++e)
6305  {
6306  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6307  VectorizedArrayType();
6308  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6309  VectorizedArrayType();
6310  }
6311  }
6312  }
6313  else
6314  {
6316  this->cell_type == internal::MatrixFreeFunctions::general ?
6317  this->jacobian[q_point] :
6318  this->jacobian[0];
6319  const VectorizedArrayType fac =
6320  (this->cell_type == internal::MatrixFreeFunctions::general ?
6321  this->J_value[q_point] :
6322  this->J_value[0] * this->quadrature_weights[q_point]) *
6323  div_in;
6324  for (unsigned int d = 0; d < dim; ++d)
6325  {
6326  for (unsigned int e = 0; e < dim; ++e)
6327  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6328  jac[d][e] * fac;
6329  }
6330  }
6331  }
6332 }
6333 
6334 
6335 
6336 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6337 inline DEAL_II_ALWAYS_INLINE void
6341  const unsigned int q_point)
6342 {
6343  AssertThrow(
6344  this->data->element_type !=
6346  ExcNotImplemented());
6347 
6348  // could have used base class operator, but that involves some overhead
6349  // which is inefficient. it is nice to have the symmetric tensor because
6350  // that saves some operations
6351 # ifdef DEBUG
6352  Assert(this->is_reinitialized, ExcNotInitialized());
6353 # endif
6354  AssertIndexRange(q_point, this->n_quadrature_points);
6355  Assert(this->J_value != nullptr,
6357  "update_gradients"));
6358  Assert(this->jacobian != nullptr,
6360  "update_gradients"));
6361 # ifdef DEBUG
6362  this->gradients_quad_submitted = true;
6363 # endif
6364 
6365  const std::size_t nqp = this->n_quadrature_points;
6366  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6367  {
6368  const VectorizedArrayType JxW =
6369  this->J_value[0] * this->quadrature_weights[q_point];
6370  for (unsigned int d = 0; d < dim; ++d)
6371  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6372  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6373  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6374  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6375  {
6376  const VectorizedArrayType value =
6377  sym_grad.access_raw_entry(counter) * JxW;
6378  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6379  value * this->jacobian[0][d][d];
6380  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6381  value * this->jacobian[0][e][e];
6382  }
6383  }
6384  // general/affine cell type
6385  else
6386  {
6387  const VectorizedArrayType JxW =
6388  this->cell_type == internal::MatrixFreeFunctions::general ?
6389  this->J_value[q_point] :
6390  this->J_value[0] * this->quadrature_weights[q_point];
6392  this->cell_type == internal::MatrixFreeFunctions::general ?
6393  this->jacobian[q_point] :
6394  this->jacobian[0];
6395  VectorizedArrayType weighted[dim][dim];
6396  for (unsigned int i = 0; i < dim; ++i)
6397  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6398  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6399  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6400  {
6401  const VectorizedArrayType value =
6402  sym_grad.access_raw_entry(counter) * JxW;
6403  weighted[i][j] = value;
6404  weighted[j][i] = value;
6405  }
6406  for (unsigned int comp = 0; comp < dim; ++comp)
6407  for (unsigned int d = 0; d < dim; ++d)
6408  {
6409  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6410  for (unsigned int e = 1; e < dim; ++e)
6411  new_val += jac[e][d] * weighted[comp][e];
6412  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6413  }
6414  }
6415 }
6416 
6417 
6418 
6419 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6420 inline DEAL_II_ALWAYS_INLINE void
6423  const unsigned int q_point)
6424 {
6426  switch (dim)
6427  {
6428  case 1:
6429  Assert(false,
6430  ExcMessage(
6431  "Testing by the curl in 1d is not a useful operation"));
6432  break;
6433  case 2:
6434  grad[1][0] = curl[0];
6435  grad[0][1] = -curl[0];
6436  break;
6437  case 3:
6438  grad[2][1] = curl[0];
6439  grad[1][2] = -curl[0];
6440  grad[0][2] = curl[1];
6441  grad[2][0] = -curl[1];
6442  grad[1][0] = curl[2];
6443  grad[0][1] = -curl[2];
6444  break;
6445  default:
6446  Assert(false, ExcNotImplemented());
6447  }
6448  submit_gradient(grad, q_point);
6449 }
6450 
6451 
6452 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6453 
6454 
6455 template <typename Number, bool is_face, typename VectorizedArrayType>
6458  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
6459  const unsigned int dof_no,
6460  const unsigned int first_selected_component,
6461  const unsigned int quad_no,
6462  const unsigned int fe_degree,
6463  const unsigned int n_q_points,
6464  const bool is_interior_face,
6465  const unsigned int active_fe_index,
6466  const unsigned int active_quad_index,
6467  const unsigned int face_type)
6468  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6469  matrix_free,
6470  dof_no,
6471  first_selected_component,
6472  quad_no,
6473  fe_degree,
6474  n_q_points,
6475  is_interior_face,
6476  active_fe_index,
6477  active_quad_index,
6478  face_type)
6479 {}
6480 
6481 
6482 
6483 template <typename Number, bool is_face, typename VectorizedArrayType>
6486  const Mapping<1> & mapping,
6487  const FiniteElement<1> &fe,
6488  const Quadrature<1> & quadrature,
6489  const UpdateFlags update_flags,
6490  const unsigned int first_selected_component,
6492  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6493  mapping,
6494  fe,
6495  quadrature,
6496  update_flags,
6497  first_selected_component,
6498  other)
6499 {}
6500 
6501 
6502 
6503 template <typename Number, bool is_face, typename VectorizedArrayType>
6507  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6508 {}
6509 
6510 
6511 
6512 template <typename Number, bool is_face, typename VectorizedArrayType>
6516 {
6518  other);
6519  return *this;
6520 }
6521 
6522 
6523 
6524 template <typename Number, bool is_face, typename VectorizedArrayType>
6525 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6527  const unsigned int dof) const
6528 {
6529  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6530  return this->values_dofs[dof];
6531 }
6532 
6533 
6534 
6535 template <typename Number, bool is_face, typename VectorizedArrayType>
6536 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6538  const unsigned int q_point) const
6539 {
6540 # ifdef DEBUG
6541  Assert(this->values_quad_initialized == true,
6543 # endif
6544  AssertIndexRange(q_point, this->n_quadrature_points);
6545  return this->values_quad[q_point];
6546 }
6547 
6548 
6549 
6550 template <typename Number, bool is_face, typename VectorizedArrayType>
6553  const unsigned int q_point) const
6554 {
6555  // could use the base class gradient, but that involves too many inefficient
6556  // initialization operations on tensors
6557 
6558 # ifdef DEBUG
6559  Assert(this->gradients_quad_initialized == true,
6561 # endif
6562  AssertIndexRange(q_point, this->n_quadrature_points);
6563 
6565  this->cell_type == internal::MatrixFreeFunctions::general ?
6566  this->jacobian[q_point] :
6567  this->jacobian[0];
6568 
6570  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6571 
6572  return grad_out;
6573 }
6574 
6575 
6576 
6577 template <typename Number, bool is_face, typename VectorizedArrayType>
6578 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6580  const unsigned int q_point) const
6581 {
6582  return get_gradient(q_point)[0];
6583 }
6584 
6585 
6586 
6587 template <typename Number, bool is_face, typename VectorizedArrayType>
6588 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6590  get_normal_derivative(const unsigned int q_point) const
6591 {
6592  return BaseClass::get_normal_derivative(q_point)[0];
6593 }
6594 
6595 
6596 
6597 template <typename Number, bool is_face, typename VectorizedArrayType>
6600  const unsigned int q_point) const
6601 {
6602  return BaseClass::get_hessian(q_point)[0];
6603 }
6604 
6605 
6606 
6607 template <typename Number, bool is_face, typename VectorizedArrayType>
6610  get_hessian_diagonal(const unsigned int q_point) const
6611 {
6612  return BaseClass::get_hessian_diagonal(q_point)[0];
6613 }
6614 
6615 
6616 
6617 template <typename Number, bool is_face, typename VectorizedArrayType>
6618 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6620  const unsigned int q_point) const
6621 {
6622  return BaseClass::get_laplacian(q_point)[0];
6623 }
6624 
6625 
6626 
6627 template <typename Number, bool is_face, typename VectorizedArrayType>
6630  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6631 {
6632 # ifdef DEBUG
6633  this->dof_values_initialized = true;
6634  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6635 # endif
6636  this->values_dofs[dof] = val_in;
6637 }
6638 
6639 
6640 
6641 template <typename Number, bool is_face, typename VectorizedArrayType>
6642 inline DEAL_II_ALWAYS_INLINE void
6644  const VectorizedArrayType val_in,
6645  const unsigned int q_point)
6646 {
6647 # ifdef DEBUG
6648  Assert(this->is_reinitialized, ExcNotInitialized());
6649 # endif
6650  AssertIndexRange(q_point, this->n_quadrature_points);
6651 # ifdef DEBUG
6652  this->values_quad_submitted = true;
6653 # endif
6654 
6655  if (this->cell_type == internal::MatrixFreeFunctions::general)
6656  {
6657  const VectorizedArrayType JxW = this->J_value[q_point];
6658  this->values_quad[q_point] = val_in * JxW;
6659  }
6660  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6661  {
6662  const VectorizedArrayType JxW =
6663  this->J_value[0] * this->quadrature_weights[q_point];
6664  this->values_quad[q_point] = val_in * JxW;
6665  }
6666 }
6667 
6668 
6669 
6670 template <typename Number, bool is_face, typename VectorizedArrayType>
6671 inline DEAL_II_ALWAYS_INLINE void
6673  const Tensor<1, 1, VectorizedArrayType> val_in,
6674  const unsigned int q_point)
6675 {
6676  submit_value(val_in[0], q_point);
6677 }
6678 
6679 
6680 
6681 template <typename Number, bool is_face, typename VectorizedArrayType>
6682 inline DEAL_II_ALWAYS_INLINE void
6684  const Tensor<1, 1, VectorizedArrayType> grad_in,
6685  const unsigned int q_point)
6686 {
6687  submit_gradient(grad_in[0], q_point);
6688 }
6689 
6690 
6691 
6692 template <typename Number, bool is_face, typename VectorizedArrayType>
6693 inline DEAL_II_ALWAYS_INLINE void
6695  const VectorizedArrayType grad_in,
6696  const unsigned int q_point)
6697 {
6698 # ifdef DEBUG
6699  Assert(this->is_reinitialized, ExcNotInitialized());
6700 # endif
6701  AssertIndexRange(q_point, this->n_quadrature_points);
6702 # ifdef DEBUG
6703  this->gradients_quad_submitted = true;
6704 # endif
6705 
6707  this->cell_type == internal::MatrixFreeFunctions::general ?
6708  this->jacobian[q_point] :
6709  this->jacobian[0];
6710  const VectorizedArrayType JxW =
6711  this->cell_type == internal::MatrixFreeFunctions::general ?
6712  this->J_value[q_point] :
6713  this->J_value[0] * this->quadrature_weights[q_point];
6714 
6715  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6716 }
6717 
6718 
6719 
6720 template <typename Number, bool is_face, typename VectorizedArrayType>
6721 inline DEAL_II_ALWAYS_INLINE void
6723  const Tensor<2, 1, VectorizedArrayType> grad_in,
6724  const unsigned int q_point)
6725 {
6726  submit_gradient(grad_in[0][0], q_point);
6727 }
6728 
6729 
6730 
6731 template <typename Number, bool is_face, typename VectorizedArrayType>
6732 inline DEAL_II_ALWAYS_INLINE void
6734  submit_normal_derivative(const VectorizedArrayType grad_in,
6735  const unsigned int q_point)
6736 {
6738  grad[0] = grad_in;
6739  BaseClass::submit_normal_derivative(grad, q_point);
6740 }
6741 
6742 
6743 
6744 template <typename Number, bool is_face, typename VectorizedArrayType>
6745 inline DEAL_II_ALWAYS_INLINE void
6748  const unsigned int q_point)
6749 {
6750  BaseClass::submit_normal_derivative(grad_in, q_point);
6751 }
6752 
6753 
6754 template <typename Number, bool is_face, typename VectorizedArrayType>
6755 inline DEAL_II_ALWAYS_INLINE void
6757  const Tensor<2, 1, VectorizedArrayType> hessian_in,
6758  const unsigned int q_point)
6759 {
6761  hessian[0] = hessian_in;
6762  BaseClass::submit_hessian(hessian, q_point);
6763 }
6764 
6765 
6766 template <typename Number, bool is_face, typename VectorizedArrayType>
6767 inline VectorizedArrayType
6769  integrate_value() const
6770 {
6771  return BaseClass::integrate_value()[0];
6772 }
6773 
6774 
6775 
6776 /*-------------------------- FEEvaluation -----------------------------------*/
6777 
6778 
6779 template <int dim,
6780  int fe_degree,
6781  int n_q_points_1d,
6782  int n_components_,
6783  typename Number,
6784  typename VectorizedArrayType>
6785 inline FEEvaluation<dim,
6786  fe_degree,
6787  n_q_points_1d,
6788  n_components_,
6789  Number,
6790  VectorizedArrayType>::
6791  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6792  const unsigned int fe_no,
6793  const unsigned int quad_no,
6794  const unsigned int first_selected_component,
6795  const unsigned int active_fe_index,
6796  const unsigned int active_quad_index)
6797  : BaseClass(matrix_free,
6798  fe_no,
6799  first_selected_component,
6800  quad_no,
6801  fe_degree,
6802  static_n_q_points,
6803  true /*note: this is not a face*/,
6804  active_fe_index,
6805  active_quad_index)
6806  , dofs_per_component(this->data->dofs_per_component_on_cell)
6807  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6808  , n_q_points(this->data->n_q_points)
6809 {
6810  check_template_arguments(fe_no, 0);
6811 }
6812 
6813 
6814 
6815 template <int dim,
6816  int fe_degree,
6817  int n_q_points_1d,
6818  int n_components_,
6819  typename Number,
6820  typename VectorizedArrayType>
6821 inline FEEvaluation<dim,
6822  fe_degree,
6823  n_q_points_1d,
6824  n_components_,
6825  Number,
6826  VectorizedArrayType>::
6827  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6828  const std::pair<unsigned int, unsigned int> & range,
6829  const unsigned int dof_no,
6830  const unsigned int quad_no,
6831  const unsigned int first_selected_component)
6832  : FEEvaluation(matrix_free,
6833  dof_no,
6834  quad_no,
6835  first_selected_component,
6836  matrix_free.get_cell_active_fe_index(range))
6837 {}
6838 
6839 
6840 
6841 template <int dim,
6842  int fe_degree,
6843  int n_q_points_1d,
6844  int n_components_,
6845  typename Number,
6846  typename VectorizedArrayType>
6847 inline FEEvaluation<dim,
6848  fe_degree,
6849  n_q_points_1d,
6850  n_components_,
6851  Number,
6852  VectorizedArrayType>::
6853  FEEvaluation(const Mapping<dim> & mapping,
6854  const FiniteElement<dim> &fe,
6855  const Quadrature<1> & quadrature,
6856  const UpdateFlags update_flags,
6857  const unsigned int first_selected_component)
6858  : BaseClass(mapping,
6859  fe,
6860  quadrature,
6861  update_flags,
6862  first_selected_component,
6863  nullptr)
6864  , dofs_per_component(this->data->dofs_per_component_on_cell)
6865  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6866  , n_q_points(this->data->n_q_points)
6867 {
6868  check_template_arguments(numbers::invalid_unsigned_int, 0);
6869 }
6870 
6871 
6872 
6873 template <int dim,
6874  int fe_degree,
6875  int n_q_points_1d,
6876  int n_components_,
6877  typename Number,
6878  typename VectorizedArrayType>
6879 inline FEEvaluation<dim,
6880  fe_degree,
6881  n_q_points_1d,
6882  n_components_,
6883  Number,
6884  VectorizedArrayType>::
6885  FEEvaluation(const FiniteElement<dim> &fe,
6886  const Quadrature<1> & quadrature,
6887  const UpdateFlags update_flags,
6888  const unsigned int first_selected_component)
6889  : BaseClass(StaticMappingQ1<dim>::mapping,
6890  fe,
6891  quadrature,
6892  update_flags,
6893  first_selected_component,
6894  nullptr)
6895  , dofs_per_component(this->data->dofs_per_component_on_cell)
6896  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6897  , n_q_points(this->data->n_q_points)
6898 {
6899  check_template_arguments(numbers::invalid_unsigned_int, 0);
6900 }
6901 
6902 
6903 
6904 template <int dim,
6905  int fe_degree,
6906  int n_q_points_1d,
6907  int n_components_,
6908  typename Number,
6909  typename VectorizedArrayType>
6910 inline FEEvaluation<dim,
6911  fe_degree,
6912  n_q_points_1d,
6913  n_components_,
6914  Number,
6915  VectorizedArrayType>::
6916  FEEvaluation(const FiniteElement<dim> & fe,
6918  const unsigned int first_selected_component)
6919  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6920  fe,
6921  other.mapped_geometry->get_quadrature(),
6922  other.mapped_geometry->get_fe_values().get_update_flags(),
6923  first_selected_component,
6924  &other)
6925  , dofs_per_component(this->data->dofs_per_component_on_cell)
6926  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6927  , n_q_points(this->data->n_q_points)
6928 {
6929  check_template_arguments(numbers::invalid_unsigned_int, 0);
6930 }
6931 
6932 
6933 
6934 template <int dim,
6935  int fe_degree,
6936  int n_q_points_1d,
6937  int n_components_,
6938  typename Number,
6939  typename VectorizedArrayType>
6940 inline FEEvaluation<dim,
6941  fe_degree,
6942  n_q_points_1d,
6943  n_components_,
6944  Number,
6945  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6946  &other)
6947  : BaseClass(other)
6948  , dofs_per_component(this->data->dofs_per_component_on_cell)
6949  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6950  , n_q_points(this->data->n_q_points)
6951 {
6952  check_template_arguments(numbers::invalid_unsigned_int, 0);
6953 }
6954 
6955 
6956 
6957 template <int dim,
6958  int fe_degree,
6959  int n_q_points_1d,
6960  int n_components_,
6961  typename Number,
6962  typename VectorizedArrayType>
6963 inline FEEvaluation<dim,
6964  fe_degree,
6965  n_q_points_1d,
6966  n_components_,
6967  Number,
6968  VectorizedArrayType> &
6969 FEEvaluation<dim,
6970  fe_degree,
6971  n_q_points_1d,
6972  n_components_,
6973  Number,
6974  VectorizedArrayType>::operator=(const FEEvaluation &other)
6975 {
6976  BaseClass::operator=(other);
6977  check_template_arguments(numbers::invalid_unsigned_int, 0);
6978  return *this;
6979 }
6980 
6981 
6982 
6983 template <int dim,
6984  int fe_degree,
6985  int n_q_points_1d,
6986  int n_components_,
6987  typename Number,
6988  typename VectorizedArrayType>
6989 inline void
6990 FEEvaluation<dim,
6991  fe_degree,
6992  n_q_points_1d,
6993  n_components_,
6994  Number,
6995  VectorizedArrayType>::
6996  check_template_arguments(const unsigned int dof_no,
6997  const unsigned int first_selected_component)
6998 {
6999  (void)dof_no;
7000  (void)first_selected_component;
7001 
7002  Assert(
7003  this->data->dofs_per_component_on_cell > 0,
7004  ExcMessage(
7005  "There is nothing useful you can do with an FEEvaluation object with "
7006  "FE_Nothing, i.e., without DoFs! If you have passed to "
7007  "MatrixFree::reinit() a collection of finite elements also containing "
7008  "FE_Nothing, please check - before creating FEEvaluation - the category "
7009  "of the current range by calling either "
7010  "MatrixFree::get_cell_range_category(range) or "
7011  "MatrixFree::get_face_range_category(range). The returned category "
7012  "is the index of the active FE, which you can use to exclude "
7013  "FE_Nothing."));
7014 
7015 # ifdef DEBUG
7016  // print error message when the dimensions do not match. Propose a possible
7017  // fix
7018  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7019  static_cast<unsigned int>(fe_degree) !=
7020  this->data->data.front().fe_degree) ||
7021  n_q_points != this->n_quadrature_points)
7022  {
7023  std::string message =
7024  "-------------------------------------------------------\n";
7025  message += "Illegal arguments in constructor/wrong template arguments!\n";
7026  message += " Called --> FEEvaluation<dim,";
7027  message += Utilities::int_to_string(fe_degree) + ",";
7028  message += Utilities::int_to_string(n_q_points_1d);
7029  message += "," + Utilities::int_to_string(n_components);
7030  message += ",Number>(data";
7031  if (first_selected_component != numbers::invalid_unsigned_int)
7032  {
7033  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7034  message += Utilities::int_to_string(this->quad_no) + ", ";
7035  message += Utilities::int_to_string(first_selected_component);
7036  }
7037  message += ")\n";
7038 
7039  // check whether some other vector component has the correct number of
7040  // points
7041  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7042  proposed_fe_comp = numbers::invalid_unsigned_int,
7043  proposed_quad_comp = numbers::invalid_unsigned_int;
7044  if (dof_no != numbers::invalid_unsigned_int)
7045  {
7046  if (static_cast<unsigned int>(fe_degree) ==
7047  this->data->data.front().fe_degree)
7048  {
7049  proposed_dof_comp = dof_no;
7050  proposed_fe_comp = first_selected_component;
7051  }
7052  else
7053  for (unsigned int no = 0; no < this->matrix_free->n_components();
7054  ++no)
7055  for (unsigned int nf = 0;
7056  nf < this->matrix_free->n_base_elements(no);
7057  ++nf)
7058  if (this->matrix_free
7059  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7060  .data.front()
7061  .fe_degree == static_cast<unsigned int>(fe_degree))
7062  {
7063  proposed_dof_comp = no;
7064  proposed_fe_comp = nf;
7065  break;
7066  }
7067  if (n_q_points ==
7068  this->mapping_data->descriptor[this->active_quad_index]
7069  .n_q_points)
7070  proposed_quad_comp = this->quad_no;
7071  else
7072  for (unsigned int no = 0;
7073  no < this->matrix_free->get_mapping_info().cell_data.size();
7074  ++no)
7075  if (this->matrix_free->get_mapping_info()
7076  .cell_data[no]
7077  .descriptor[this->active_quad_index]
7078  .n_q_points == n_q_points)
7079  {
7080  proposed_quad_comp = no;
7081  break;
7082  }
7083  }
7084  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7085  proposed_quad_comp != numbers::invalid_unsigned_int)
7086  {
7087  if (proposed_dof_comp != first_selected_component)
7088  message += "Wrong vector component selection:\n";
7089  else
7090  message += "Wrong quadrature formula selection:\n";
7091  message += " Did you mean FEEvaluation<dim,";
7092  message += Utilities::int_to_string(fe_degree) + ",";
7093  message += Utilities::int_to_string(n_q_points_1d);
7094  message += "," + Utilities::int_to_string(n_components);
7095  message += ",Number>(data";
7096  if (dof_no != numbers::invalid_unsigned_int)
7097  {
7098  message +=
7099  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7100  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7101  message += Utilities::int_to_string(proposed_fe_comp);
7102  }
7103  message += ")?\n";
7104  std::string correct_pos;
7105  if (proposed_dof_comp != dof_no)
7106  correct_pos = " ^ ";
7107  else
7108  correct_pos = " ";
7109  if (proposed_quad_comp != this->quad_no)
7110  correct_pos += " ^ ";
7111  else
7112  correct_pos += " ";
7113  if (proposed_fe_comp != first_selected_component)
7114  correct_pos += " ^\n";
7115  else
7116  correct_pos += " \n";
7117  message += " " +
7118  correct_pos;
7119  }
7120  // ok, did not find the numbers specified by the template arguments in
7121  // the given list. Suggest correct template arguments
7122  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7123  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7124  message += "Wrong template arguments:\n";
7125  message += " Did you mean FEEvaluation<dim,";
7126  message +=
7127  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7128  message += Utilities::int_to_string(proposed_n_q_points_1d);
7129  message += "," + Utilities::int_to_string(n_components);
7130  message += ",Number>(data";
7131  if (dof_no != numbers::invalid_unsigned_int)
7132  {
7133  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7134  message += Utilities::int_to_string(this->quad_no);
7135  message += ", " + Utilities::int_to_string(first_selected_component);
7136  }
7137  message += ")?\n";
7138  std::string correct_pos;
7139  if (this->data->data.front().fe_degree !=
7140  static_cast<unsigned int>(fe_degree))
7141  correct_pos = " ^";
7142  else
7143  correct_pos = " ";
7144  if (proposed_n_q_points_1d != n_q_points_1d)
7145  correct_pos += " ^\n";
7146  else
7147  correct_pos += " \n";
7148  message += " " + correct_pos;
7149 
7150  Assert(static_cast<unsigned int>(fe_degree) ==
7151  this->data->data.front().fe_degree &&
7152  n_q_points == this->n_quadrature_points,
7153  ExcMessage(message));
7154  }
7155  if (dof_no != numbers::invalid_unsigned_int)
7157  n_q_points,
7158  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7159 # endif
7160 }
7161 
7162 
7163 
7164 template <int dim,
7165  int fe_degree,
7166  int n_q_points_1d,
7167  int n_components_,
7168  typename Number,
7169  typename VectorizedArrayType>
7170 inline void
7171 FEEvaluation<dim,
7172  fe_degree,
7173  n_q_points_1d,
7174  n_components_,
7175  Number,
7176  VectorizedArrayType>::reinit(const unsigned int cell_index)
7177 {
7178  Assert(this->mapped_geometry == nullptr,
7179  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7180  " Integer indexing is not possible"));
7181  if (this->mapped_geometry != nullptr)
7182  return;
7183 
7184  Assert(this->dof_info != nullptr, ExcNotInitialized());
7185  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7186  this->cell = cell_index;
7187  this->cell_type =
7188  this->matrix_free->get_mapping_info().get_cell_type(cell_index);
7189 
7190  const unsigned int offsets =
7191  this->mapping_data->data_index_offsets[cell_index];
7192  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7193  this->J_value = &this->mapping_data->JxW_values[offsets];
7194  this->jacobian_gradients =
7195  this->mapping_data->jacobian_gradients[0].data() + offsets;
7196 
7197  unsigned int i = 0;
7198  for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7199  ++i)
7200  this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7201  for (; i < VectorizedArrayType::size(); ++i)
7202  this->cell_ids[i] = numbers::invalid_unsigned_int;
7203 
7204  if (this->mapping_data->quadrature_points.empty() == false)
7205  this->quadrature_points =
7206  &this->mapping_data->quadrature_points
7207  [this->mapping_data->quadrature_point_offsets[this->cell]];
7208 
7209 # ifdef DEBUG
7210  this->is_reinitialized = true;
7211  this->dof_values_initialized = false;
7212  this->values_quad_initialized = false;
7213  this->gradients_quad_initialized = false;
7214  this->hessians_quad_initialized = false;
7215 # endif
7216 }
7217 
7218 
7219 
7220 template <int dim,
7221  int fe_degree,
7222  int n_q_points_1d,
7223  int n_components_,
7224  typename Number,
7225  typename VectorizedArrayType>
7226 inline void
7227 FEEvaluation<dim,
7228  fe_degree,
7229  n_q_points_1d,
7230  n_components_,
7231  Number,
7232  VectorizedArrayType>::
7233  reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
7234 {
7235  Assert(this->dof_info != nullptr, ExcNotInitialized());
7236  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7237 
7238  this->cell = numbers::invalid_unsigned_int;
7239  this->cell_ids = cell_ids;
7240 
7241  // determine type of cell batch
7243 
7244  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7245  {
7246  const unsigned int cell_index = cell_ids[v];
7247 
7249  continue;
7250 
7251  this->cell_type =
7252  std::max(this->cell_type,
7253  this->matrix_free->get_mapping_info().get_cell_type(
7254  cell_index / VectorizedArrayType::size()));
7255  }
7256 
7257  // allocate memory for internal data storage
7258  if (this->mapped_geometry == nullptr)
7259  this->mapped_geometry =
7260  std::make_shared<internal::MatrixFreeFunctions::
7261  MappingDataOnTheFly<dim, VectorizedArrayType>>();
7262 
7263  auto &mapping_storage = this->mapped_geometry->get_data_storage();
7264 
7265  auto &this_jacobian_data = mapping_storage.jacobians[0];
7266  auto &this_J_value_data = mapping_storage.JxW_values;
7267  auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7268  auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7269 
7271  {
7272  if (this->mapping_data->jacobians[0].size() > 0)
7273  this_jacobian_data.resize_fast(2);
7274 
7275  if (this->mapping_data->JxW_values.size() > 0)
7276  this_J_value_data.resize_fast(1);
7277 
7278  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7279  this_jacobian_gradients_data.resize_fast(1);
7280 
7281  if (this->mapping_data->quadrature_points.size() > 0)
7282  this_quadrature_points_data.resize_fast(1);
7283  }
7284  else
7285  {
7286  if (this->mapping_data->jacobians[0].size() > 0)
7287  this_jacobian_data.resize_fast(this->n_quadrature_points);
7288 
7289  if (this->mapping_data->JxW_values.size() > 0)
7290  this_J_value_data.resize_fast(this->n_quadrature_points);
7291 
7292  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7293  this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7294 
7295  if (this->mapping_data->quadrature_points.size() > 0)
7296  this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7297  }
7298 
7299  // set pointers to internal data storage
7300  this->jacobian = this_jacobian_data.data();
7301  this->J_value = this_J_value_data.data();
7302  this->jacobian_gradients = this_jacobian_gradients_data.data();
7303  this->quadrature_points = this_quadrature_points_data.data();
7304 
7305  // fill internal data storage lane by lane
7306  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7307  {
7308  const unsigned int cell_index = cell_ids[v];
7309 
7311  continue;
7312 
7313  const unsigned int cell_batch_index =
7314  cell_index / VectorizedArrayType::size();
7315  const unsigned int offsets =
7316  this->mapping_data->data_index_offsets[cell_batch_index];
7317  const unsigned int lane = cell_index % VectorizedArrayType::size();
7318 
7319  if (this->cell_type <=
7321  {
7322  // case that all cells are Cartesian or affine
7323  const unsigned int q = 0;
7324 
7325  if (this->mapping_data->JxW_values.size() > 0)
7326  this_J_value_data[q][v] =
7327  this->mapping_data->JxW_values[offsets + q][lane];
7328 
7329  if (this->mapping_data->jacobians[0].size() > 0)
7330  for (unsigned int q = 0; q < 2; ++q)
7331  for (unsigned int i = 0; i < dim; ++i)
7332  for (unsigned int j = 0; j < dim; ++j)
7333  this_jacobian_data[q][i][j][v] =
7334  this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7335 
7336  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7337  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7338  for (unsigned int j = 0; j < dim; ++j)
7339  this_jacobian_gradients_data[q][i][j][v] =
7340  this->mapping_data
7341  ->jacobian_gradients[0][offsets + q][i][j][lane];
7342 
7343  if (this->mapping_data->quadrature_points.size() > 0)
7344  for (unsigned int i = 0; i < dim; ++i)
7345  this_quadrature_points_data[q][i][v] =
7346  this->mapping_data->quadrature_points
7347  [this->mapping_data
7348  ->quadrature_point_offsets[cell_batch_index] +
7349  q][i][lane];
7350  }
7351  else
7352  {
7353  // general case that at least one cell is not Cartesian or affine
7354  const auto cell_type =
7355  this->matrix_free->get_mapping_info().get_cell_type(
7356  cell_batch_index);
7357 
7358  for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
7359  {
7360  const unsigned int q_src =
7361  (cell_type <=
7363  0 :
7364  q;
7365 
7366  if (this->mapping_data->JxW_values.size() > 0)
7367  this_J_value_data[q][v] =
7368  this->mapping_data->JxW_values[offsets + q_src][lane];
7369 
7370  if (this->mapping_data->jacobians[0].size() > 0)
7371  for (unsigned int i = 0; i < dim; ++i)
7372  for (unsigned int j = 0; j < dim; ++j)
7373  this_jacobian_data[q][i][j][v] =
7374  this->mapping_data
7375  ->jacobians[0][offsets + q_src][i][j][lane];
7376 
7377  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7378  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7379  for (unsigned int j = 0; j < dim; ++j)
7380  this_jacobian_gradients_data[q][i][j][v] =
7381  this->mapping_data
7382  ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7383 
7384  if (this->mapping_data->quadrature_points.size() > 0)
7385  {
7386  if (cell_type <=
7388  {
7389  // affine case: quadrature points are not available but
7390  // have to be computed from the corner point and the
7391  // Jacobian
7393  this->mapping_data->quadrature_points
7394  [this->mapping_data
7395  ->quadrature_point_offsets[cell_batch_index] +
7396  0];
7397 
7399  this->mapping_data->jacobians[0][offsets + 1];
7401  for (unsigned int d = 0; d < dim; ++d)
7402  point[d] +=
7403  jac[d][d] *
7404  static_cast<Number>(
7405  this->descriptor->quadrature.point(q)[d]);
7406  else
7407  for (unsigned int d = 0; d < dim; ++d)
7408  for (unsigned int e = 0; e < dim; ++e)
7409  point[d] +=
7410  jac[d][e] *
7411  static_cast<Number>(
7412  this->descriptor->quadrature.point(q)[e]);
7413 
7414  for (unsigned int i = 0; i < dim; ++i)
7415  this_quadrature_points_data[q][i][v] = point[i][lane];
7416  }
7417  else
7418  {
7419  // general case: quadrature points are available
7420  for (unsigned int i = 0; i < dim; ++i)
7421  this_quadrature_points_data[q][i][v] =
7422  this->mapping_data->quadrature_points
7423  [this->mapping_data
7424  ->quadrature_point_offsets[cell_batch_index] +
7425  q][i][lane];
7426  }
7427  }
7428  }
7429  }
7430  }
7431 
7432 # ifdef DEBUG
7433  this->is_reinitialized = true;
7434  this->dof_values_initialized = false;
7435  this->values_quad_initialized = false;
7436  this->gradients_quad_initialized = false;
7437  this->hessians_quad_initialized = false;
7438 # endif
7439 }
7440 
7441 
7442 
7443 template <int dim,
7444  int fe_degree,
7445  int n_q_points_1d,
7446  int n_components_,
7447  typename Number,
7448  typename VectorizedArrayType>
7449 template <bool level_dof_access>
7450 inline void
7451 FEEvaluation<dim,
7452  fe_degree,
7453  n_q_points_1d,
7454  n_components_,
7455  Number,
7456  VectorizedArrayType>::
7458 {
7459  Assert(this->matrix_free == nullptr,
7460  ExcMessage("Cannot use initialization from cell iterator if "
7461  "initialized from MatrixFree object. Use variant for "
7462  "on the fly computation with arguments as for FEValues "
7463  "instead"));
7464  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7465  this->mapped_geometry->reinit(
7466  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7467  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7468  if (level_dof_access)
7469  cell->get_mg_dof_indices(this->local_dof_indices);
7470  else
7471  cell->get_dof_indices(this->local_dof_indices);
7472 
7473 # ifdef DEBUG
7474  this->is_reinitialized = true;
7475 # endif
7476 }
7477 
7478 
7479 
7480 template <int dim,
7481  int fe_degree,
7482  int n_q_points_1d,
7483  int n_components_,
7484  typename Number,
7485  typename VectorizedArrayType>
7486 inline void
7487 FEEvaluation<dim,
7488  fe_degree,
7489  n_q_points_1d,
7490  n_components_,
7491  Number,
7492  VectorizedArrayType>::
7493  reinit(const typename Triangulation<dim>::cell_iterator &cell)
7494 {
7495  Assert(this->matrix_free == 0,
7496  ExcMessage("Cannot use initialization from cell iterator if "
7497  "initialized from MatrixFree object. Use variant for "
7498  "on the fly computation with arguments as for FEValues "
7499  "instead"));
7500  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7501  this->mapped_geometry->reinit(cell);
7502 
7503 # ifdef DEBUG
7504  this->is_reinitialized = true;
7505 # endif
7506 }
7507 
7508 
7509 
7510 template <int dim,
7511  int fe_degree,
7512  int n_q_points_1d,
7513  int n_components_,
7514  typename Number,
7515  typename VectorizedArrayType>
7516 inline void
7517 FEEvaluation<dim,
7518  fe_degree,
7519  n_q_points_1d,
7520  n_components_,
7521  Number,
7522  VectorizedArrayType>::evaluate(const bool evaluate_values,
7523  const bool evaluate_gradients,
7524  const bool evaluate_hessians)
7525 {
7526 # ifdef DEBUG
7527  Assert(this->dof_values_initialized == true,
7529 # endif
7530  evaluate(this->values_dofs,
7531  evaluate_values,
7532  evaluate_gradients,
7533  evaluate_hessians);
7534 }
7535 
7536 
7537 template <int dim,
7538  int fe_degree,
7539  int n_q_points_1d,
7540  int n_components_,
7541  typename Number,
7542  typename VectorizedArrayType>
7543 inline void
7544 FEEvaluation<dim,
7545  fe_degree,
7546  n_q_points_1d,
7547  n_components_,
7548  Number,
7549  VectorizedArrayType>::
7550  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7551 {
7552 # ifdef DEBUG
7553  Assert(this->dof_values_initialized == true,
7555 # endif
7556  evaluate(this->values_dofs, evaluation_flags);
7557 }
7558 
7559 
7560 
7561 template <int dim,
7562  int fe_degree,
7563  int n_q_points_1d,
7564  int n_components_,
7565  typename Number,
7566  typename VectorizedArrayType>
7567 inline void
7568 FEEvaluation<dim,
7569  fe_degree,
7570  n_q_points_1d,
7571  n_components_,
7572  Number,
7573  VectorizedArrayType>::evaluate(const VectorizedArrayType
7574  * values_array,
7575  const bool evaluate_values,
7576  const bool evaluate_gradients,
7577  const bool evaluate_hessians)
7578 {
7580  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7581  ((evaluate_gradients) ? EvaluationFlags::gradients :
7583  ((evaluate_hessians) ? EvaluationFlags::hessians :
7585 
7586  evaluate(values_array, flag);
7587 }
7588 
7589 
7590