Reference documentation for deal.II version 9.4.0
\(\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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY 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_EARLY void
2679  integrate(const bool integrate_values, const bool integrate_gradients);
2680 
2689  void
2691  VectorizedArrayType * values_array);
2692 
2696  DEAL_II_DEPRECATED_EARLY 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  Assert(mask.count() == n_lanes,
3346  ExcNotImplemented("Masking currently not implemented for "
3347  "non-contiguous DoF storage"));
3348 
3349  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3350  this->get_cell_ids();
3351 
3352  bool has_hn_constraints = false;
3353 
3354  if (is_face == false)
3355  {
3356  for (unsigned int v = 0; v < n_lanes; ++v)
3357  if (cells[v] != numbers::invalid_unsigned_int &&
3358  this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3359  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3360  this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3363  this->dof_info->hanging_node_constraint_masks_comp
3364  [this->active_fe_index][this->first_selected_component])
3365  has_hn_constraints = true;
3366  }
3367 
3368  std::integral_constant<bool,
3369  internal::is_vectorizable<VectorType, Number>::value>
3370  vector_selector;
3371 
3372  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3373  std::array<VectorizedArrayType *, n_components> values_dofs;
3374  for (unsigned int c = 0; c < n_components; ++c)
3375  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3376  c * dofs_per_component;
3377 
3378  if (this->cell != numbers::invalid_unsigned_int &&
3379  this->dof_info->index_storage_variants
3380  [is_face ? this->dof_access_index :
3382  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3383  IndexStorageVariants::interleaved &&
3384  (has_hn_constraints == false))
3385  {
3386  const unsigned int *dof_indices =
3387  this->dof_info->dof_indices_interleaved.data() +
3388  this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes]
3389  .first +
3390  this->dof_info
3391  ->component_dof_indices_offset[this->active_fe_index]
3392  [this->first_selected_component] *
3393  n_lanes;
3394  if (n_components == 1 || this->n_fe_components == 1)
3395  for (unsigned int i = 0; i < dofs_per_component;
3396  ++i, dof_indices += n_lanes)
3397  for (unsigned int comp = 0; comp < n_components; ++comp)
3398  operation.process_dof_gather(dof_indices,
3399  *src[comp],
3400  0,
3401  values_dofs[comp][i],
3402  vector_selector);
3403  else
3404  for (unsigned int comp = 0; comp < n_components; ++comp)
3405  for (unsigned int i = 0; i < dofs_per_component;
3406  ++i, dof_indices += n_lanes)
3407  operation.process_dof_gather(
3408  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3409  return;
3410  }
3411 
3412  // Allocate pointers, then initialize all of them to nullptrs and
3413  // below overwrite the ones we actually use:
3414  std::array<const unsigned int *, n_lanes> dof_indices;
3415  dof_indices.fill(nullptr);
3416 
3417  // Assign the appropriate cell ids for face/cell case and get the pointers
3418  // to the dof indices of the cells on all lanes
3419 
3420  bool has_constraints = false;
3421  const unsigned int n_components_read =
3422  this->n_fe_components > 1 ? n_components : 1;
3423 
3424  if (is_face)
3425  {
3426  for (unsigned int v = 0; v < n_lanes; ++v)
3427  {
3428  if (cells[v] == numbers::invalid_unsigned_int)
3429  continue;
3430 
3431  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
3432  ExcInternalError());
3433  const std::pair<unsigned int, unsigned int> *my_index_start =
3434  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3435  this->first_selected_component];
3436 
3437  // check whether any of the SIMD lanes has constraints, i.e., the
3438  // constraint indicator which is the second entry of row_starts
3439  // increments on this cell
3440  if (my_index_start[n_components_read].second !=
3441  my_index_start[0].second)
3442  has_constraints = true;
3443 
3444  dof_indices[v] =
3445  this->dof_info->dof_indices.data() + my_index_start[0].first;
3446  }
3447  }
3448  else
3449  {
3450  for (unsigned int v = 0; v < n_lanes; ++v)
3451  {
3452  if (cells[v] == numbers::invalid_unsigned_int)
3453  continue;
3454 
3455  const std::pair<unsigned int, unsigned int> *my_index_start =
3456  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3457  this->first_selected_component];
3458  if (my_index_start[n_components_read].second !=
3459  my_index_start[0].second)
3460  has_constraints = true;
3461 
3462  if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3463  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3464  this->dof_info->hanging_node_constraint_masks[cells[v]] !=
3467  this->dof_info->hanging_node_constraint_masks_comp
3468  [this->active_fe_index][this->first_selected_component])
3469  has_hn_constraints = true;
3470 
3471  Assert(my_index_start[n_components_read].first ==
3472  my_index_start[0].first ||
3473  my_index_start[0].first < this->dof_info->dof_indices.size(),
3474  ExcIndexRange(0,
3475  my_index_start[0].first,
3476  this->dof_info->dof_indices.size()));
3477  dof_indices[v] =
3478  this->dof_info->dof_indices.data() + my_index_start[0].first;
3479  }
3480  }
3481 
3482  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3483  return i != numbers::invalid_unsigned_int;
3484  }) < n_lanes)
3485  for (unsigned int comp = 0; comp < n_components; ++comp)
3486  for (unsigned int i = 0; i < dofs_per_component; ++i)
3487  operation.process_empty(values_dofs[comp][i]);
3488 
3489  // Case where we have no constraints throughout the whole cell: Can go
3490  // through the list of DoFs directly
3491  if (!has_constraints && apply_constraints)
3492  {
3493  if (n_components == 1 || this->n_fe_components == 1)
3494  {
3495  for (unsigned int v = 0; v < n_lanes; ++v)
3496  {
3497  if (cells[v] == numbers::invalid_unsigned_int)
3498  continue;
3499 
3500  for (unsigned int i = 0; i < dofs_per_component; ++i)
3501  for (unsigned int comp = 0; comp < n_components; ++comp)
3502  operation.process_dof(dof_indices[v][i],
3503  *src[comp],
3504  values_dofs[comp][i][v]);
3505  }
3506  }
3507  else
3508  {
3509  for (unsigned int comp = 0; comp < n_components; ++comp)
3510  for (unsigned int v = 0; v < n_lanes; ++v)
3511  {
3512  if (cells[v] == numbers::invalid_unsigned_int)
3513  continue;
3514 
3515  for (unsigned int i = 0; i < dofs_per_component; ++i)
3516  operation.process_dof(
3517  dof_indices[v][comp * dofs_per_component + i],
3518  *src[0],
3519  values_dofs[comp][i][v]);
3520  }
3521  }
3522  return;
3523  }
3524 
3525  // In the case where there are some constraints to be resolved, loop over
3526  // all vector components that are filled and then over local dofs. ind_local
3527  // holds local number on cell, index iterates over the elements of
3528  // index_local_to_global and dof_indices points to the global indices stored
3529  // in index_local_to_global
3530 
3531  for (unsigned int v = 0; v < n_lanes; ++v)
3532  {
3533  if (cells[v] == numbers::invalid_unsigned_int)
3534  continue;
3535 
3536  const unsigned int cell_index = cells[v];
3537  const unsigned int cell_dof_index =
3538  cell_index * this->n_fe_components + this->first_selected_component;
3539  const unsigned int n_components_read =
3540  this->n_fe_components > 1 ? n_components : 1;
3541  unsigned int index_indicators =
3542  this->dof_info->row_starts[cell_dof_index].second;
3543  unsigned int next_index_indicators =
3544  this->dof_info->row_starts[cell_dof_index + 1].second;
3545 
3546  // For read_dof_values_plain, redirect the dof_indices field to the
3547  // unconstrained indices
3548  if (apply_constraints == false &&
3549  (this->dof_info->row_starts[cell_dof_index].second !=
3550  this->dof_info->row_starts[cell_dof_index + n_components_read]
3551  .second ||
3552  ((this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3553  this->dof_info->hanging_node_constraint_masks_comp.size() > 0 &&
3554  this->dof_info->hanging_node_constraint_masks[cell_index] !=
3557  this->dof_info->hanging_node_constraint_masks_comp
3558  [this->active_fe_index][this->first_selected_component])))
3559  {
3560  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
3562  ExcNotInitialized());
3563  dof_indices[v] =
3564  this->dof_info->plain_dof_indices.data() +
3565  this->dof_info
3566  ->component_dof_indices_offset[this->active_fe_index]
3567  [this->first_selected_component] +
3568  this->dof_info->row_starts_plain_indices[cell_index];
3569  next_index_indicators = index_indicators;
3570  }
3571 
3572  if (n_components == 1 || this->n_fe_components == 1)
3573  {
3574  unsigned int ind_local = 0;
3575  for (; index_indicators != next_index_indicators; ++index_indicators)
3576  {
3577  const std::pair<unsigned short, unsigned short> indicator =
3578  this->dof_info->constraint_indicator[index_indicators];
3579  // run through values up to next constraint
3580  for (unsigned int j = 0; j < indicator.first; ++j)
3581  for (unsigned int comp = 0; comp < n_components; ++comp)
3582  operation.process_dof(dof_indices[v][j],
3583  *src[comp],
3584  values_dofs[comp][ind_local + j][v]);
3585 
3586  ind_local += indicator.first;
3587  dof_indices[v] += indicator.first;
3588 
3589  // constrained case: build the local value as a linear
3590  // combination of the global value according to constraints
3591  Number value[n_components];
3592  for (unsigned int comp = 0; comp < n_components; ++comp)
3593  operation.pre_constraints(values_dofs[comp][ind_local][v],
3594  value[comp]);
3595 
3596  const Number *data_val =
3597  this->matrix_free->constraint_pool_begin(indicator.second);
3598  const Number *end_pool =
3599  this->matrix_free->constraint_pool_end(indicator.second);
3600  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3601  for (unsigned int comp = 0; comp < n_components; ++comp)
3602  operation.process_constraint(*dof_indices[v],
3603  *data_val,
3604  *src[comp],
3605  value[comp]);
3606 
3607  for (unsigned int comp = 0; comp < n_components; ++comp)
3608  operation.post_constraints(value[comp],
3609  values_dofs[comp][ind_local][v]);
3610  ind_local++;
3611  }
3612 
3613  AssertIndexRange(ind_local, dofs_per_component + 1);
3614 
3615  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3616  for (unsigned int comp = 0; comp < n_components; ++comp)
3617  operation.process_dof(*dof_indices[v],
3618  *src[comp],
3619  values_dofs[comp][ind_local][v]);
3620  }
3621  else
3622  {
3623  // case with vector-valued finite elements where all components are
3624  // included in one single vector. Assumption: first come all entries
3625  // to the first component, then all entries to the second one, and
3626  // so on. This is ensured by the way MatrixFree reads out the
3627  // indices.
3628  for (unsigned int comp = 0; comp < n_components; ++comp)
3629  {
3630  unsigned int ind_local = 0;
3631 
3632  // check whether there is any constraint on the current cell
3633  for (; index_indicators != next_index_indicators;
3634  ++index_indicators)
3635  {
3636  const std::pair<unsigned short, unsigned short> indicator =
3637  this->dof_info->constraint_indicator[index_indicators];
3638 
3639  // run through values up to next constraint
3640  for (unsigned int j = 0; j < indicator.first; ++j)
3641  operation.process_dof(dof_indices[v][j],
3642  *src[0],
3643  values_dofs[comp][ind_local + j][v]);
3644  ind_local += indicator.first;
3645  dof_indices[v] += indicator.first;
3646 
3647  // constrained case: build the local value as a linear
3648  // combination of the global value according to constraints
3649  Number value;
3650  operation.pre_constraints(values_dofs[comp][ind_local][v],
3651  value);
3652 
3653  const Number *data_val =
3654  this->matrix_free->constraint_pool_begin(indicator.second);
3655  const Number *end_pool =
3656  this->matrix_free->constraint_pool_end(indicator.second);
3657 
3658  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3659  operation.process_constraint(*dof_indices[v],
3660  *data_val,
3661  *src[0],
3662  value);
3663 
3664  operation.post_constraints(value,
3665  values_dofs[comp][ind_local][v]);
3666  ind_local++;
3667  }
3668 
3669  AssertIndexRange(ind_local, dofs_per_component + 1);
3670 
3671  // get the dof values past the last constraint
3672  for (; ind_local < dofs_per_component;
3673  ++dof_indices[v], ++ind_local)
3674  {
3675  AssertIndexRange(*dof_indices[v], src[0]->size());
3676  operation.process_dof(*dof_indices[v],
3677  *src[0],
3678  values_dofs[comp][ind_local][v]);
3679  }
3680 
3681  if (apply_constraints == true && comp + 1 < n_components)
3682  next_index_indicators =
3683  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
3684  }
3685  }
3686  }
3687 }
3688 
3689 
3690 
3691 template <int dim,
3692  int n_components_,
3693  typename Number,
3694  bool is_face,
3695  typename VectorizedArrayType>
3696 template <typename VectorType, typename VectorOperation>
3697 inline void
3700  const VectorOperation & operation,
3701  const std::array<VectorType *, n_components_> &src) const
3702 {
3703  Assert(!local_dof_indices.empty(), ExcNotInitialized());
3704 
3705  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3706  unsigned int index = this->first_selected_component * dofs_per_component;
3707  for (unsigned int comp = 0; comp < n_components; ++comp)
3708  {
3709  for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3710  {
3711  operation.process_empty(
3712  this->values_dofs[comp * dofs_per_component + i]);
3713  operation.process_dof_global(
3714  local_dof_indices[this->data->lexicographic_numbering[index]],
3715  *src[0],
3716  this->values_dofs[comp * dofs_per_component + i][0]);
3717  }
3718  }
3719 }
3720 
3721 
3722 
3723 template <int dim,
3724  int n_components_,
3725  typename Number,
3726  bool is_face,
3727  typename VectorizedArrayType>
3728 template <typename VectorType, typename VectorOperation>
3729 inline void
3732  const VectorOperation & operation,
3733  const std::array<VectorType *, n_components_> &src,
3734  const std::array<
3736  n_components_> & vectors_sm,
3737  const std::bitset<VectorizedArrayType::size()> &mask) const
3738 {
3739  // This functions processes the functions read_dof_values,
3740  // distribute_local_to_global, and set_dof_values with the same code for
3741  // contiguous cell indices (DG case). The distinction between these three
3742  // cases is made by the input VectorOperation that either reads values from
3743  // a vector and puts the data into the local data field or write local data
3744  // into the vector. Certain operations are no-ops for the given use case.
3745 
3746  std::integral_constant<bool,
3747  internal::is_vectorizable<VectorType, Number>::value>
3748  vector_selector;
3750  is_face ? this->dof_access_index :
3752  const unsigned int n_lanes = mask.count();
3753 
3754  const std::vector<unsigned int> &dof_indices_cont =
3755  this->dof_info->dof_indices_contiguous[ind];
3756 
3757  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3758  std::array<VectorizedArrayType *, n_components> values_dofs;
3759  for (unsigned int c = 0; c < n_components; ++c)
3760  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3761  c * dofs_per_component;
3762 
3764 
3765  // Simple case: We have contiguous storage, so we can simply copy out the
3766  // data
3767  if ((this->dof_info->index_storage_variants[ind][this->cell] ==
3769  interleaved_contiguous &&
3770  n_lanes == VectorizedArrayType::size()) &&
3771  !(is_face &&
3772  this->dof_access_index ==
3774  this->is_interior_face() == false) &&
3775  !(!is_face && !this->is_interior_face()))
3776  {
3777  const unsigned int dof_index =
3778  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3779  this->dof_info
3780  ->component_dof_indices_offset[this->active_fe_index]
3781  [this->first_selected_component] *
3782  VectorizedArrayType::size();
3783  if (n_components == 1 || this->n_fe_components == 1)
3784  for (unsigned int comp = 0; comp < n_components; ++comp)
3785  operation.process_dofs_vectorized(dofs_per_component,
3786  dof_index,
3787  *src[comp],
3788  values_dofs[comp],
3789  vector_selector);
3790  else
3791  operation.process_dofs_vectorized(dofs_per_component * n_components,
3792  dof_index,
3793  *src[0],
3794  values_dofs[0],
3795  vector_selector);
3796  return;
3797  }
3798 
3799  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3800  this->get_cell_or_face_ids();
3801 
3802  // More general case: Must go through the components one by one and apply
3803  // some transformations
3804  const unsigned int n_filled_lanes =
3805  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
3806 
3807  const bool is_ecl =
3808  (this->dof_access_index ==
3810  this->is_interior_face() == false) ||
3811  (!is_face && !this->is_interior_face());
3812 
3813  if (vectors_sm[0] != nullptr)
3814  {
3815  const auto compute_vector_ptrs = [&](const unsigned int comp) {
3816  std::array<typename VectorType::value_type *,
3817  VectorizedArrayType::size()>
3818  vector_ptrs = {};
3819 
3820  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3821  {
3822  if (mask[v] == false)
3823  {
3824  vector_ptrs[v] = nullptr;
3825  continue;
3826  }
3827 
3829  ExcNotImplemented());
3830  Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
3831  ExcIndexRange(
3832  ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
3833  Assert(cells[v] <
3834  this->dof_info->dof_indices_contiguous_sm[ind].size(),
3835  ExcIndexRange(
3836  cells[v],
3837  0,
3838  this->dof_info->dof_indices_contiguous_sm[ind].size()));
3839 
3840  const auto &temp =
3841  this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
3842 
3843  if (temp.first != numbers::invalid_unsigned_int)
3844  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3845  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3846  this->dof_info->component_dof_indices_offset
3847  [this->active_fe_index][this->first_selected_component]);
3848  else
3849  vector_ptrs[v] = nullptr;
3850  }
3851  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3852  ++v)
3853  vector_ptrs[v] = nullptr;
3854 
3855  return vector_ptrs;
3856  };
3857 
3858  if (n_filled_lanes == VectorizedArrayType::size() &&
3859  n_lanes == VectorizedArrayType::size() && !is_ecl)
3860  {
3861  if (n_components == 1 || this->n_fe_components == 1)
3862  {
3863  for (unsigned int comp = 0; comp < n_components; ++comp)
3864  {
3865  auto vector_ptrs = compute_vector_ptrs(comp);
3866  operation.process_dofs_vectorized_transpose(
3867  dofs_per_component,
3868  vector_ptrs,
3869  values_dofs[comp],
3870  vector_selector);
3871  }
3872  }
3873  else
3874  {
3875  auto vector_ptrs = compute_vector_ptrs(0);
3876  operation.process_dofs_vectorized_transpose(dofs_per_component *
3877  n_components,
3878  vector_ptrs,
3879  &values_dofs[0][0],
3880  vector_selector);
3881  }
3882  }
3883  else
3884  for (unsigned int comp = 0; comp < n_components; ++comp)
3885  {
3886  auto vector_ptrs = compute_vector_ptrs(
3887  (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3888 
3889  for (unsigned int i = 0; i < dofs_per_component; ++i)
3890  operation.process_empty(values_dofs[comp][i]);
3891 
3892  if (n_components == 1 || this->n_fe_components == 1)
3893  {
3894  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3895  if (mask[v] == true)
3896  for (unsigned int i = 0; i < dofs_per_component; ++i)
3897  operation.process_dof(vector_ptrs[v][i],
3898  values_dofs[comp][i][v]);
3899  }
3900  else
3901  {
3902  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3903  if (mask[v] == true)
3904  for (unsigned int i = 0; i < dofs_per_component; ++i)
3905  operation.process_dof(
3906  vector_ptrs[v][i + comp * dofs_per_component],
3907  values_dofs[comp][i][v]);
3908  }
3909  }
3910  return;
3911  }
3912 
3913  unsigned int dof_indices[VectorizedArrayType::size()];
3914 
3915  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3916  {
3918  dof_indices[v] =
3919  dof_indices_cont[cells[v]] +
3920  this->dof_info
3921  ->component_dof_indices_offset[this->active_fe_index]
3922  [this->first_selected_component] *
3923  this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
3924  }
3925 
3926  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
3927  dof_indices[v] = numbers::invalid_unsigned_int;
3928 
3929  // In the case with contiguous cell indices, we know that there are no
3930  // constraints and that the indices within each element are contiguous
3931  if (n_filled_lanes == VectorizedArrayType::size() &&
3932  n_lanes == VectorizedArrayType::size() && !is_ecl)
3933  {
3934  if (this->dof_info->index_storage_variants[ind][this->cell] ==
3936  contiguous)
3937  {
3938  if (n_components == 1 || this->n_fe_components == 1)
3939  for (unsigned int comp = 0; comp < n_components; ++comp)
3940  operation.process_dofs_vectorized_transpose(dofs_per_component,
3941  dof_indices,
3942  *src[comp],
3943  values_dofs[comp],
3944  vector_selector);
3945  else
3946  operation.process_dofs_vectorized_transpose(dofs_per_component *
3947  n_components,
3948  dof_indices,
3949  *src[0],
3950  &values_dofs[0][0],
3951  vector_selector);
3952  }
3953  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
3955  interleaved_contiguous_strided)
3956  {
3957  if (n_components == 1 || this->n_fe_components == 1)
3958  for (unsigned int i = 0; i < dofs_per_component; ++i)
3959  {
3960  for (unsigned int comp = 0; comp < n_components; ++comp)
3961  operation.process_dof_gather(dof_indices,
3962  *src[comp],
3963  i * VectorizedArrayType::size(),
3964  values_dofs[comp][i],
3965  vector_selector);
3966  }
3967  else
3968  for (unsigned int comp = 0; comp < n_components; ++comp)
3969  for (unsigned int i = 0; i < dofs_per_component; ++i)
3970  {
3971  operation.process_dof_gather(dof_indices,
3972  *src[0],
3973  (comp * dofs_per_component + i) *
3974  VectorizedArrayType::size(),
3975  values_dofs[comp][i],
3976  vector_selector);
3977  }
3978  }
3979  else
3980  {
3981  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
3983  IndexStorageVariants::interleaved_contiguous_mixed_strides,
3984  ExcNotImplemented());
3985  const unsigned int *offsets =
3986  &this->dof_info->dof_indices_interleave_strides
3987  [ind][VectorizedArrayType::size() * this->cell];
3988  if (n_components == 1 || this->n_fe_components == 1)
3989  for (unsigned int i = 0; i < dofs_per_component; ++i)
3990  {
3991  for (unsigned int comp = 0; comp < n_components; ++comp)
3992  operation.process_dof_gather(dof_indices,
3993  *src[comp],
3994  0,
3995  values_dofs[comp][i],
3996  vector_selector);
3998  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
3999  dof_indices[v] += offsets[v];
4000  }
4001  else
4002  for (unsigned int comp = 0; comp < n_components; ++comp)
4003  for (unsigned int i = 0; i < dofs_per_component; ++i)
4004  {
4005  operation.process_dof_gather(dof_indices,
4006  *src[0],
4007  0,
4008  values_dofs[comp][i],
4009  vector_selector);
4011  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4012  dof_indices[v] += offsets[v];
4013  }
4014  }
4015  }
4016  else
4017  for (unsigned int comp = 0; comp < n_components; ++comp)
4018  {
4019  for (unsigned int i = 0; i < dofs_per_component; ++i)
4020  operation.process_empty(values_dofs[comp][i]);
4021  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4023  contiguous)
4024  {
4025  if (n_components == 1 || this->n_fe_components == 1)
4026  {
4027  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4028  if (mask[v] == true)
4029  for (unsigned int i = 0; i < dofs_per_component; ++i)
4030  operation.process_dof(dof_indices[v] + i,
4031  *src[comp],
4032  values_dofs[comp][i][v]);
4033  }
4034  else
4035  {
4036  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4037  if (mask[v] == true)
4038  for (unsigned int i = 0; i < dofs_per_component; ++i)
4039  operation.process_dof(dof_indices[v] + i +
4040  comp * dofs_per_component,
4041  *src[0],
4042  values_dofs[comp][i][v]);
4043  }
4044  }
4045  else
4046  {
4047  const unsigned int *offsets =
4048  &this->dof_info->dof_indices_interleave_strides
4049  [ind][VectorizedArrayType::size() * this->cell];
4050  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4051  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4052  if (n_components == 1 || this->n_fe_components == 1)
4053  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4054  {
4055  if (mask[v] == true)
4056  for (unsigned int i = 0; i < dofs_per_component; ++i)
4057  operation.process_dof(dof_indices[v] + i * offsets[v],
4058  *src[comp],
4059  values_dofs[comp][i][v]);
4060  }
4061  else
4062  {
4063  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4064  if (mask[v] == true)
4065  for (unsigned int i = 0; i < dofs_per_component; ++i)
4066  operation.process_dof(dof_indices[v] +
4067  (i + comp * dofs_per_component) *
4068  offsets[v],
4069  *src[0],
4070  values_dofs[comp][i][v]);
4071  }
4072  }
4073  }
4074 }
4075 
4076 namespace internal
4077 {
4078  template <typename Number,
4079  typename VectorType,
4080  typename std::enable_if<!IsBlockVector<VectorType>::value,
4081  VectorType>::type * = nullptr>
4082  decltype(std::declval<VectorType>().begin())
4083  get_beginning(VectorType &vec)
4084  {
4085  return vec.begin();
4086  }
4087 
4088  template <typename Number,
4089  typename VectorType,
4090  typename std::enable_if<IsBlockVector<VectorType>::value,
4091  VectorType>::type * = nullptr>
4092  typename VectorType::value_type *
4093  get_beginning(VectorType &)
4094  {
4095  return nullptr;
4096  }
4097 
4098  template <typename VectorType,
4099  typename std::enable_if<has_shared_vector_data<VectorType>,
4100  VectorType>::type * = nullptr>
4101  const std::vector<ArrayView<const typename VectorType::value_type>> *
4102  get_shared_vector_data(VectorType & vec,
4103  const bool is_valid_mode_for_sm,
4104  const unsigned int active_fe_index,
4106  {
4107  // note: no hp is supported
4108  if (is_valid_mode_for_sm &&
4109  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4110  .size() > 0 &&
4111  active_fe_index == 0)
4112  return &vec.shared_vector_data();
4113  else
4114  return nullptr;
4115  }
4116 
4117  template <typename VectorType,
4118  typename std::enable_if<!has_shared_vector_data<VectorType>,
4119  VectorType>::type * = nullptr>
4120  const std::vector<ArrayView<const typename VectorType::value_type>> *
4121  get_shared_vector_data(VectorType &,
4122  const bool,
4123  const unsigned int,
4125  {
4126  return nullptr;
4127  }
4128 
4129  template <int n_components, typename VectorType>
4130  std::pair<
4131  std::array<typename internal::BlockVectorSelector<
4132  VectorType,
4133  IsBlockVector<VectorType>::value>::BaseVectorType *,
4134  n_components>,
4135  std::array<
4136  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4137  VectorType,
4138  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4139  n_components>>
4140  get_vector_data(VectorType & src,
4141  const unsigned int first_index,
4142  const bool is_valid_mode_for_sm,
4143  const unsigned int active_fe_index,
4145  {
4146  // select between block vectors and non-block vectors. Note that the number
4147  // of components is checked in the internal data
4148  std::pair<
4149  std::array<typename internal::BlockVectorSelector<
4150  VectorType,
4151  IsBlockVector<VectorType>::value>::BaseVectorType *,
4152  n_components>,
4153  std::array<
4154  const std::vector<
4155  ArrayView<const typename internal::BlockVectorSelector<
4156  VectorType,
4157  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4158  n_components>>
4159  src_data;
4160 
4161  for (unsigned int d = 0; d < n_components; ++d)
4162  src_data.first[d] = internal::BlockVectorSelector<
4163  VectorType,
4164  IsBlockVector<VectorType>::value>::get_vector_component(src,
4165  d +
4166  first_index);
4167 
4168  for (unsigned int d = 0; d < n_components; ++d)
4169  src_data.second[d] = get_shared_vector_data(
4170  const_cast<typename internal::BlockVectorSelector<
4171  typename std::remove_const<VectorType>::type,
4173  BaseVectorType &>(*src_data.first[d]),
4174  is_valid_mode_for_sm,
4175  active_fe_index,
4176  dof_info);
4177 
4178  return src_data;
4179  }
4180 } // namespace internal
4181 
4182 
4183 
4184 template <int dim,
4185  int n_components_,
4186  typename Number,
4187  bool is_face,
4188  typename VectorizedArrayType>
4189 inline void
4192 {
4193  if (this->dof_info == nullptr ||
4194  this->dof_info->hanging_node_constraint_masks.size() == 0 ||
4195  this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4196  this->dof_info->hanging_node_constraint_masks_comp
4197  [this->active_fe_index][this->first_selected_component] == false)
4198  return; // nothing to do with faces
4199 
4200  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4201  std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4202  constraint_mask;
4203 
4204  bool hn_available = false;
4205 
4206  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
4207  this->get_cell_ids();
4208 
4209  for (unsigned int v = 0; v < n_lanes; ++v)
4210  {
4211  if (cells[v] == numbers::invalid_unsigned_int)
4212  {
4213  constraint_mask[v] = internal::MatrixFreeFunctions::
4215  continue;
4216  }
4217 
4218  const unsigned int cell_index = cells[v];
4219  const auto mask =
4220  this->dof_info->hanging_node_constraint_masks[cell_index];
4221  constraint_mask[v] = mask;
4222 
4223  hn_available |= (mask != internal::MatrixFreeFunctions::
4225  }
4226 
4227  if (hn_available == false)
4228  return; // no hanging node on cell batch -> nothing to do
4229 
4231  apply(n_components,
4232  this->data->data.front().fe_degree,
4233  this->get_shape_info(),
4234  transpose,
4235  constraint_mask,
4236  this->values_dofs);
4237 }
4238 
4239 
4240 
4241 template <int dim,
4242  int n_components_,
4243  typename Number,
4244  bool is_face,
4245  typename VectorizedArrayType>
4246 template <typename VectorType>
4247 inline void
4249  read_dof_values(const VectorType & src,
4250  const unsigned int first_index,
4251  const std::bitset<VectorizedArrayType::size()> &mask)
4252 {
4253  const auto src_data = internal::get_vector_data<n_components_>(
4254  src,
4255  first_index,
4256  this->dof_access_index ==
4258  this->active_fe_index,
4259  this->dof_info);
4260 
4262  read_write_operation(reader, src_data.first, src_data.second, mask, true);
4263 
4264  apply_hanging_node_constraints(false);
4265 
4266 # ifdef DEBUG
4267  this->dof_values_initialized = true;
4268 # endif
4269 }
4270 
4271 
4272 
4273 template <int dim,
4274  int n_components_,
4275  typename Number,
4276  bool is_face,
4277  typename VectorizedArrayType>
4278 template <typename VectorType>
4279 inline void
4281  read_dof_values_plain(const VectorType & src,
4282  const unsigned int first_index,
4283  const std::bitset<VectorizedArrayType::size()> &mask)
4284 {
4285  const auto src_data = internal::get_vector_data<n_components_>(
4286  src,
4287  first_index,
4288  this->dof_access_index ==
4290  this->active_fe_index,
4291  this->dof_info);
4292 
4294  read_write_operation(reader, src_data.first, src_data.second, mask, false);
4295 
4296 # ifdef DEBUG
4297  this->dof_values_initialized = true;
4298 # endif
4299 }
4300 
4301 
4302 
4303 template <int dim,
4304  int n_components_,
4305  typename Number,
4306  bool is_face,
4307  typename VectorizedArrayType>
4308 template <typename VectorType>
4309 inline void
4312  VectorType & dst,
4313  const unsigned int first_index,
4314  const std::bitset<VectorizedArrayType::size()> &mask) const
4315 {
4316 # ifdef DEBUG
4317  Assert(this->dof_values_initialized == true,
4319 # endif
4320 
4321  apply_hanging_node_constraints(true);
4322 
4323  const auto dst_data = internal::get_vector_data<n_components_>(
4324  dst,
4325  first_index,
4326  this->dof_access_index ==
4328  this->active_fe_index,
4329  this->dof_info);
4330 
4332  distributor;
4333  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4334 }
4335 
4336 
4337 
4338 template <int dim,
4339  int n_components_,
4340  typename Number,
4341  bool is_face,
4342  typename VectorizedArrayType>
4343 template <typename VectorType>
4344 inline void
4346  set_dof_values(VectorType & dst,
4347  const unsigned int first_index,
4348  const std::bitset<VectorizedArrayType::size()> &mask) const
4349 {
4350 # ifdef DEBUG
4351  Assert(this->dof_values_initialized == true,
4353 # endif
4354 
4355  const auto dst_data = internal::get_vector_data<n_components_>(
4356  dst,
4357  first_index,
4358  this->dof_access_index ==
4360  this->active_fe_index,
4361  this->dof_info);
4362 
4364  read_write_operation(setter, dst_data.first, dst_data.second, mask);
4365 }
4366 
4367 
4368 
4369 template <int dim,
4370  int n_components_,
4371  typename Number,
4372  bool is_face,
4373  typename VectorizedArrayType>
4374 template <typename VectorType>
4375 inline void
4378  VectorType & dst,
4379  const unsigned int first_index,
4380  const std::bitset<VectorizedArrayType::size()> &mask) const
4381 {
4382 # ifdef DEBUG
4383  Assert(this->dof_values_initialized == true,
4385 # endif
4386 
4387  const auto dst_data = internal::get_vector_data<n_components_>(
4388  dst,
4389  first_index,
4390  this->dof_access_index ==
4392  this->active_fe_index,
4393  this->dof_info);
4394 
4396  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4397 }
4398 
4399 
4400 
4401 /*------------------------------ access to data fields ----------------------*/
4402 
4403 
4404 
4405 template <int dim,
4406  int n_components_,
4407  typename Number,
4408  bool is_face,
4409  typename VectorizedArrayType>
4412  get_dof_value(const unsigned int dof) const
4413 {
4414  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4415  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4417  for (unsigned int comp = 0; comp < n_components; ++comp)
4418  return_value[comp] = this->values_dofs[comp * dofs + dof];
4419  return return_value;
4420 }
4421 
4422 
4423 
4424 template <int dim,
4425  int n_components_,
4426  typename Number,
4427  bool is_face,
4428  typename VectorizedArrayType>
4431  get_value(const unsigned int q_point) const
4432 {
4433 # ifdef DEBUG
4434  Assert(this->values_quad_initialized == true,
4436 # endif
4437 
4438  AssertIndexRange(q_point, this->n_quadrature_points);
4439  const std::size_t nqp = this->n_quadrature_points;
4441  for (unsigned int comp = 0; comp < n_components; ++comp)
4442  return_value[comp] = this->values_quad[comp * nqp + q_point];
4443  return return_value;
4444 }
4445 
4446 
4447 
4448 template <int dim,
4449  int n_components_,
4450  typename Number,
4451  bool is_face,
4452  typename VectorizedArrayType>
4453 inline DEAL_II_ALWAYS_INLINE
4456  get_gradient(const unsigned int q_point) const
4457 {
4458 # ifdef DEBUG
4459  Assert(this->gradients_quad_initialized == true,
4461 # endif
4462 
4463  AssertIndexRange(q_point, this->n_quadrature_points);
4464  Assert(this->jacobian != nullptr,
4466  "update_gradients"));
4467  const std::size_t nqp = this->n_quadrature_points;
4469 
4470  // Cartesian cell
4471  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4472  {
4473  for (unsigned int d = 0; d < dim; ++d)
4474  for (unsigned int comp = 0; comp < n_components; ++comp)
4475  grad_out[comp][d] =
4476  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4477  this->jacobian[0][d][d];
4478  }
4479  // cell with general/affine Jacobian
4480  else
4481  {
4483  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4484  q_point :
4485  0];
4486  for (unsigned int comp = 0; comp < n_components; ++comp)
4487  for (unsigned int d = 0; d < dim; ++d)
4488  {
4489  grad_out[comp][d] =
4490  jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4491  for (unsigned int e = 1; e < dim; ++e)
4492  grad_out[comp][d] +=
4493  jac[d][e] *
4494  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4495  }
4496  }
4497  return grad_out;
4498 }
4499 
4500 
4501 
4502 template <int dim,
4503  int n_components_,
4504  typename Number,
4505  bool is_face,
4506  typename VectorizedArrayType>
4509  get_normal_derivative(const unsigned int q_point) const
4510 {
4511  AssertIndexRange(q_point, this->n_quadrature_points);
4512 # ifdef DEBUG
4513  Assert(this->gradients_quad_initialized == true,
4515 # endif
4516 
4517  Assert(this->normal_x_jacobian != nullptr,
4519  "update_gradients"));
4520 
4521  const std::size_t nqp = this->n_quadrature_points;
4523 
4524  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4525  for (unsigned int comp = 0; comp < n_components; ++comp)
4526  grad_out[comp] =
4527  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4528  (this->normal_x_jacobian[0][dim - 1]);
4529  else
4530  {
4531  const std::size_t index =
4532  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4533  for (unsigned int comp = 0; comp < n_components; ++comp)
4534  {
4535  grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4536  this->normal_x_jacobian[index][0];
4537  for (unsigned int d = 1; d < dim; ++d)
4538  grad_out[comp] +=
4539  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4540  this->normal_x_jacobian[index][d];
4541  }
4542  }
4543  return grad_out;
4544 }
4545 
4546 
4547 
4548 namespace internal
4549 {
4550  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4551  // store the lower diagonal because of symmetry
4552  template <typename VectorizedArrayType>
4553  inline void
4554  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4555  const VectorizedArrayType *const hessians,
4556  const unsigned int,
4557  VectorizedArrayType (&tmp)[1][1])
4558  {
4559  tmp[0][0] = jac[0][0] * hessians[0];
4560  }
4561 
4562  template <typename VectorizedArrayType>
4563  inline void
4564  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4565  const VectorizedArrayType *const hessians,
4566  const unsigned int nqp,
4567  VectorizedArrayType (&tmp)[2][2])
4568  {
4569  for (unsigned int d = 0; d < 2; ++d)
4570  {
4571  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4572  tmp[1][d] =
4573  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4574  }
4575  }
4576 
4577  template <typename VectorizedArrayType>
4578  inline void
4579  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4580  const VectorizedArrayType *const hessians,
4581  const unsigned int nqp,
4582  VectorizedArrayType (&tmp)[3][3])
4583  {
4584  for (unsigned int d = 0; d < 3; ++d)
4585  {
4586  tmp[0][d] =
4587  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4588  jac[d][2] * hessians[4 * nqp]);
4589  tmp[1][d] =
4590  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4591  jac[d][2] * hessians[5 * nqp]);
4592  tmp[2][d] =
4593  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4594  jac[d][2] * hessians[2 * nqp]);
4595  }
4596  }
4597 } // namespace internal
4598 
4599 
4600 
4601 template <int dim,
4602  int n_components_,
4603  typename Number,
4604  bool is_face,
4605  typename VectorizedArrayType>
4608  get_hessian(const unsigned int q_point) const
4609 {
4610 # ifdef DEBUG
4611  Assert(this->hessians_quad_initialized == true,
4613 # endif
4614  AssertIndexRange(q_point, this->n_quadrature_points);
4615 
4616  Assert(this->jacobian != nullptr,
4618  "update_hessian"));
4620  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4621  0 :
4622  q_point];
4623 
4625 
4626  const std::size_t nqp = this->n_quadrature_points;
4627  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4628 
4629  // Cartesian cell
4630  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4631  {
4632  for (unsigned int comp = 0; comp < n_components; ++comp)
4633  {
4634  for (unsigned int d = 0; d < dim; ++d)
4635  hessian_out[comp][d][d] =
4636  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4637  (jac[d][d] * jac[d][d]);
4638  switch (dim)
4639  {
4640  case 1:
4641  break;
4642  case 2:
4643  hessian_out[comp][0][1] =
4644  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4645  (jac[0][0] * jac[1][1]);
4646  break;
4647  case 3:
4648  hessian_out[comp][0][1] =
4649  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4650  (jac[0][0] * jac[1][1]);
4651  hessian_out[comp][0][2] =
4652  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4653  (jac[0][0] * jac[2][2]);
4654  hessian_out[comp][1][2] =
4655  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4656  (jac[1][1] * jac[2][2]);
4657  break;
4658  default:
4659  Assert(false, ExcNotImplemented());
4660  }
4661  for (unsigned int d = 0; d < dim; ++d)
4662  for (unsigned int e = d + 1; e < dim; ++e)
4663  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4664  }
4665  }
4666  // cell with general Jacobian, but constant within the cell
4667  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4668  {
4669  for (unsigned int comp = 0; comp < n_components; ++comp)
4670  {
4671  VectorizedArrayType tmp[dim][dim];
4672  internal::hessian_unit_times_jac(
4673  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4674 
4675  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4676  for (unsigned int d = 0; d < dim; ++d)
4677  for (unsigned int e = d; e < dim; ++e)
4678  {
4679  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4680  for (unsigned int f = 1; f < dim; ++f)
4681  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4682  }
4683 
4684  // no J' * grad(u) part here because the Jacobian is constant
4685  // throughout the cell and hence, its derivative is zero
4686 
4687  // take symmetric part
4688  for (unsigned int d = 0; d < dim; ++d)
4689  for (unsigned int e = d + 1; e < dim; ++e)
4690  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4691  }
4692  }
4693  // cell with general Jacobian
4694  else
4695  {
4696  const auto &jac_grad = this->jacobian_gradients[q_point];
4697  for (unsigned int comp = 0; comp < n_components; ++comp)
4698  {
4699  VectorizedArrayType tmp[dim][dim];
4700  internal::hessian_unit_times_jac(
4701  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4702 
4703  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4704  for (unsigned int d = 0; d < dim; ++d)
4705  for (unsigned int e = d; e < dim; ++e)
4706  {
4707  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4708  for (unsigned int f = 1; f < dim; ++f)
4709  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4710  }
4711 
4712  // add diagonal part of J' * grad(u)
4713  for (unsigned int d = 0; d < dim; ++d)
4714  for (unsigned int e = 0; e < dim; ++e)
4715  hessian_out[comp][d][d] +=
4716  jac_grad[d][e] *
4717  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4718 
4719  // add off-diagonal part of J' * grad(u)
4720  for (unsigned int d = 0, count = dim; d < dim; ++d)
4721  for (unsigned int e = d + 1; e < dim; ++e, ++count)
4722  for (unsigned int f = 0; f < dim; ++f)
4723  hessian_out[comp][d][e] +=
4724  jac_grad[count][f] *
4725  this->gradients_quad[(comp * dim + f) * nqp + q_point];
4726 
4727  // take symmetric part
4728  for (unsigned int d = 0; d < dim; ++d)
4729  for (unsigned int e = d + 1; e < dim; ++e)
4730  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4731  }
4732  }
4733  return hessian_out;
4734 }
4735 
4736 
4737 
4738 template <int dim,
4739  int n_components_,
4740  typename Number,
4741  bool is_face,
4742  typename VectorizedArrayType>
4745  get_hessian_diagonal(const unsigned int q_point) const
4746 {
4747  Assert(!is_face, ExcNotImplemented());
4748 # ifdef DEBUG
4749  Assert(this->hessians_quad_initialized == true,
4751 # endif
4752  AssertIndexRange(q_point, this->n_quadrature_points);
4753 
4754  Assert(this->jacobian != nullptr, ExcNotImplemented());
4756  this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4757  0 :
4758  q_point];
4759 
4760  const std::size_t nqp = this->n_quadrature_points;
4761  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4763 
4764  // Cartesian cell
4765  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4766  {
4767  for (unsigned int comp = 0; comp < n_components; ++comp)
4768  for (unsigned int d = 0; d < dim; ++d)
4769  hessian_out[comp][d] =
4770  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4771  (jac[d][d] * jac[d][d]);
4772  }
4773  // cell with general Jacobian, but constant within the cell
4774  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4775  {
4776  for (unsigned int comp = 0; comp < n_components; ++comp)
4777  {
4778  // compute laplacian before the gradient because it needs to access
4779  // unscaled gradient data
4780  VectorizedArrayType tmp[dim][dim];
4781  internal::hessian_unit_times_jac(
4782  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4783 
4784  // compute only the trace part of hessian, J * tmp = J *
4785  // hess_unit(u) * J^T
4786  for (unsigned int d = 0; d < dim; ++d)
4787  {
4788  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4789  for (unsigned int f = 1; f < dim; ++f)
4790  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4791  }
4792  }
4793  }
4794  // cell with general Jacobian
4795  else
4796  {
4797  const auto &jac_grad = this->jacobian_gradients[q_point];
4798  for (unsigned int comp = 0; comp < n_components; ++comp)
4799  {
4800  // compute laplacian before the gradient because it needs to access
4801  // unscaled gradient data
4802  VectorizedArrayType tmp[dim][dim];
4803  internal::hessian_unit_times_jac(
4804  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4805 
4806  // compute only the trace part of hessian, J * tmp = J *
4807  // hess_unit(u) * J^T
4808  for (unsigned int d = 0; d < dim; ++d)
4809  {
4810  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4811  for (unsigned int f = 1; f < dim; ++f)
4812  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4813  }
4814 
4815  for (unsigned int d = 0; d < dim; ++d)
4816  for (unsigned int e = 0; e < dim; ++e)
4817  hessian_out[comp][d] +=
4818  jac_grad[d][e] *
4819  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4820  }
4821  }
4822  return hessian_out;
4823 }
4824 
4825 
4826 
4827 template <int dim,
4828  int n_components_,
4829  typename Number,
4830  bool is_face,
4831  typename VectorizedArrayType>
4834  get_laplacian(const unsigned int q_point) const
4835 {
4836  Assert(is_face == false, ExcNotImplemented());
4837 # ifdef DEBUG
4838  Assert(this->hessians_quad_initialized == true,
4840 # endif
4841  AssertIndexRange(q_point, this->n_quadrature_points);
4842 
4844  const auto hess_diag = get_hessian_diagonal(q_point);
4845  for (unsigned int comp = 0; comp < n_components; ++comp)
4846  {
4847  laplacian_out[comp] = hess_diag[comp][0];
4848  for (unsigned int d = 1; d < dim; ++d)
4849  laplacian_out[comp] += hess_diag[comp][d];
4850  }
4851  return laplacian_out;
4852 }
4853 
4854 
4855 
4856 template <int dim,
4857  int n_components_,
4858  typename Number,
4859  bool is_face,
4860  typename VectorizedArrayType>
4861 inline DEAL_II_ALWAYS_INLINE void
4864  const unsigned int dof)
4865 {
4866 # ifdef DEBUG
4867  this->dof_values_initialized = true;
4868 # endif
4869  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4870  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4871  for (unsigned int comp = 0; comp < n_components; ++comp)
4872  this->values_dofs[comp * dofs + dof] = val_in[comp];
4873 }
4874 
4875 
4876 
4877 template <int dim,
4878  int n_components_,
4879  typename Number,
4880  bool is_face,
4881  typename VectorizedArrayType>
4882 inline DEAL_II_ALWAYS_INLINE void
4885  const unsigned int q_point)
4886 {
4887 # ifdef DEBUG
4888  Assert(this->is_reinitialized, ExcNotInitialized());
4889 # endif
4890  AssertIndexRange(q_point, this->n_quadrature_points);
4891  Assert(this->J_value != nullptr,
4893  "update_values"));
4894 # ifdef DEBUG
4895  this->values_quad_submitted = true;
4896 # endif
4897 
4898  const std::size_t nqp = this->n_quadrature_points;
4899  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4900  {
4901  const VectorizedArrayType JxW =
4902  this->J_value[0] * this->quadrature_weights[q_point];
4903  for (unsigned int comp = 0; comp < n_components; ++comp)
4904  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4905  }
4906  else
4907  {
4908  const VectorizedArrayType JxW = this->J_value[q_point];
4909  for (unsigned int comp = 0; comp < n_components; ++comp)
4910  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
4911  }
4912 }
4913 
4914 
4915 
4916 template <int dim,
4917  int n_components_,
4918  typename Number,
4919  bool is_face,
4920  typename VectorizedArrayType>
4921 inline DEAL_II_ALWAYS_INLINE void
4924  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
4925  const unsigned int q_point)
4926 {
4927 # ifdef DEBUG
4928  Assert(this->is_reinitialized, ExcNotInitialized());
4929 # endif
4930  AssertIndexRange(q_point, this->n_quadrature_points);
4931  Assert(this->J_value != nullptr,
4933  "update_gradients"));
4934  Assert(this->jacobian != nullptr,
4936  "update_gradients"));
4937 # ifdef DEBUG
4938  this->gradients_quad_submitted = true;
4939 # endif
4940 
4941  const std::size_t nqp = this->n_quadrature_points;
4942  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4943  {
4944  const VectorizedArrayType JxW =
4945  this->J_value[0] * this->quadrature_weights[q_point];
4946  for (unsigned int d = 0; d < dim; ++d)
4947  {
4948  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
4949  for (unsigned int comp = 0; comp < n_components; ++comp)
4950  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
4951  grad_in[comp][d] * factor;
4952  }
4953  }
4954  else
4955  {
4957  this->cell_type > internal::MatrixFreeFunctions::affine ?
4958  this->jacobian[q_point] :
4959  this->jacobian[0];
4960  const VectorizedArrayType JxW =
4961  this->cell_type > internal::MatrixFreeFunctions::affine ?
4962  this->J_value[q_point] :
4963  this->J_value[0] * this->quadrature_weights[q_point];
4964  for (unsigned int comp = 0; comp < n_components; ++comp)
4965  for (unsigned int d = 0; d < dim; ++d)
4966  {
4967  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
4968  for (unsigned int e = 1; e < dim; ++e)
4969  new_val += (jac[e][d] * grad_in[comp][e]);
4970  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
4971  new_val * JxW;
4972  }
4973  }
4974 }
4975 
4976 
4977 
4978 template <int dim,
4979  int n_components_,
4980  typename Number,
4981  bool is_face,
4982  typename VectorizedArrayType>
4983 inline DEAL_II_ALWAYS_INLINE void
4987  const unsigned int q_point)
4988 {
4989  AssertIndexRange(q_point, this->n_quadrature_points);
4990  Assert(this->normal_x_jacobian != nullptr,
4992  "update_gradients"));
4993 # ifdef DEBUG
4994  this->gradients_quad_submitted = true;
4995 # endif
4996 
4997  const std::size_t nqp = this->n_quadrature_points;
4998  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4999  for (unsigned int comp = 0; comp < n_components; ++comp)
5000  {
5001  for (unsigned int d = 0; d < dim - 1; ++d)
5002  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5003  VectorizedArrayType();
5004  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5005  grad_in[comp] *
5006  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5007  this->quadrature_weights[q_point]);
5008  }
5009  else
5010  {
5011  const unsigned int index =
5012  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5014  this->normal_x_jacobian[index];
5015  for (unsigned int comp = 0; comp < n_components; ++comp)
5016  {
5017  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5018  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5019  factor = factor * this->quadrature_weights[q_point];
5020  for (unsigned int d = 0; d < dim; ++d)
5021  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5022  factor * jac[d];
5023  }
5024  }
5025 }
5026 
5027 
5028 
5029 template <int dim,
5030  int n_components_,
5031  typename Number,
5032  bool is_face,
5033  typename VectorizedArrayType>
5034 inline DEAL_II_ALWAYS_INLINE void
5037  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5038  hessian_in,
5039  const unsigned int q_point)
5040 {
5041 # ifdef DEBUG
5042  Assert(this->is_reinitialized, ExcNotInitialized());
5043 # endif
5044  AssertIndexRange(q_point, this->n_quadrature_points);
5045  Assert(this->J_value != nullptr,
5047  "update_hessians"));
5048  Assert(this->jacobian != nullptr,
5050  "update_hessians"));
5051 # ifdef DEBUG
5052  this->hessians_quad_submitted = true;
5053 # endif
5054 
5055  // compute hessian_unit = J^T * hessian_in(u) * J
5056  const std::size_t nqp = this->n_quadrature_points;
5057  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5058  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5059  {
5060  const VectorizedArrayType JxW =
5061  this->J_value[0] * this->quadrature_weights[q_point];
5062 
5063  // diagonal part
5064  for (unsigned int d = 0; d < dim; ++d)
5065  {
5066  const auto jac_d = this->jacobian[0][d][d];
5067  const VectorizedArrayType factor = jac_d * jac_d * JxW;
5068  for (unsigned int comp = 0; comp < n_components; ++comp)
5069  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5070  hessian_in[comp][d][d] * factor;
5071  }
5072 
5073  // off diagonal part
5074  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5075  for (unsigned int e = 0; e < d; ++e, ++off_dia)
5076  {
5077  const auto jac_d = this->jacobian[0][d][d];
5078  const auto jac_e = this->jacobian[0][e][e];
5079  const VectorizedArrayType factor = jac_d * jac_e * JxW;
5080  for (unsigned int comp = 0; comp < n_components; ++comp)
5081  this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5082  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5083  }
5084  }
5085  // cell with general Jacobian, but constant within the cell
5086  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5087  {
5088  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5089  const VectorizedArrayType JxW =
5090  this->J_value[0] * this->quadrature_weights[q_point];
5091  for (unsigned int comp = 0; comp < n_components; ++comp)
5092  {
5093  // 1. tmp = hessian_in(u) * J
5094  VectorizedArrayType tmp[dim][dim];
5095  for (unsigned int i = 0; i < dim; ++i)
5096  for (unsigned int j = 0; j < dim; ++j)
5097  {
5098  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5099  for (unsigned int k = 1; k < dim; ++k)
5100  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5101  }
5102 
5103  // 2. hessian_unit = J^T * tmp
5104  VectorizedArrayType tmp2[dim][dim];
5105  for (unsigned int i = 0; i < dim; ++i)
5106  for (unsigned int j = 0; j < dim; ++j)
5107  {
5108  tmp2[i][j] = jac[0][i] * tmp[0][j];
5109  for (unsigned int k = 1; k < dim; ++k)
5110  tmp2[i][j] += jac[k][i] * tmp[k][j];
5111  }
5112 
5113  // diagonal part
5114  for (unsigned int d = 0; d < dim; ++d)
5115  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5116  tmp2[d][d] * JxW;
5117 
5118  // off diagonal part
5119  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5120  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5121  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5122  (tmp2[d][e] + tmp2[e][d]) * JxW;
5123  }
5124  }
5125  else
5126  {
5127  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5128  const VectorizedArrayType JxW = this->J_value[q_point];
5129  const auto &jac_grad = this->jacobian_gradients[q_point];
5130  for (unsigned int comp = 0; comp < n_components; ++comp)
5131  {
5132  // 1. tmp = hessian_in(u) * J
5133  VectorizedArrayType tmp[dim][dim];
5134  for (unsigned int i = 0; i < dim; ++i)
5135  for (unsigned int j = 0; j < dim; ++j)
5136  {
5137  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5138  for (unsigned int k = 1; k < dim; ++k)
5139  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5140  }
5141 
5142  // 2. hessian_unit = J^T * tmp
5143  VectorizedArrayType tmp2[dim][dim];
5144  for (unsigned int i = 0; i < dim; ++i)
5145  for (unsigned int j = 0; j < dim; ++j)
5146  {
5147  tmp2[i][j] = jac[0][i] * tmp[0][j];
5148  for (unsigned int k = 1; k < dim; ++k)
5149  tmp2[i][j] += jac[k][i] * tmp[k][j];
5150  }
5151 
5152  // diagonal part
5153  for (unsigned int d = 0; d < dim; ++d)
5154  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5155  tmp2[d][d] * JxW;
5156 
5157  // off diagonal part
5158  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5159  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5160  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5161  (tmp2[d][e] + tmp2[e][d]) * JxW;
5162 
5163  // 3. gradient_unit = J' ** hessian_in
5164  for (unsigned int d = 0; d < dim; ++d)
5165  {
5166  VectorizedArrayType sum = 0;
5167  for (unsigned int e = 0; e < dim; ++e)
5168  sum += hessian_in[comp][e][e] * jac_grad[e][d];
5169  for (unsigned int e = 0, count = dim; e < dim; ++e)
5170  for (unsigned int f = e + 1; f < dim; ++f, ++count)
5171  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5172  jac_grad[count][d];
5173  this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
5174  q_point] = sum * JxW;
5175  }
5176  }
5177  }
5178 }
5179 
5180 
5181 
5182 template <int dim,
5183  int n_components_,
5184  typename Number,
5185  bool is_face,
5186  typename VectorizedArrayType>
5189  integrate_value() const
5190 {
5191 # ifdef DEBUG
5192  Assert(this->is_reinitialized, ExcNotInitialized());
5193  Assert(this->values_quad_submitted == true,
5195 # endif
5196 
5198  const std::size_t nqp = this->n_quadrature_points;
5199  for (unsigned int q = 0; q < nqp; ++q)
5200  for (unsigned int comp = 0; comp < n_components; ++comp)
5201  return_value[comp] += this->values_quad[comp * nqp + q];
5202  return (return_value);
5203 }
5204 
5205 
5206 
5207 /*----------------------- FEEvaluationAccess --------------------------------*/
5208 
5209 
5210 template <int dim,
5211  int n_components_,
5212  typename Number,
5213  bool is_face,
5214  typename VectorizedArrayType>
5215 inline FEEvaluationAccess<dim,
5216  n_components_,
5217  Number,
5218  is_face,
5219  VectorizedArrayType>::
5220  FEEvaluationAccess(
5222  const unsigned int dof_no,
5223  const unsigned int first_selected_component,
5224  const unsigned int quad_no,
5225  const unsigned int fe_degree,
5226  const unsigned int n_q_points,
5227  const bool is_interior_face,
5228  const unsigned int active_fe_index,
5229  const unsigned int active_quad_index,
5230  const unsigned int face_type)
5231  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5232  matrix_free,
5233  dof_no,
5234  first_selected_component,
5235  quad_no,
5236  fe_degree,
5237  n_q_points,
5238  is_interior_face,
5239  active_fe_index,
5240  active_quad_index,
5241  face_type)
5242 {}
5243 
5244 
5245 
5246 template <int dim,
5247  int n_components_,
5248  typename Number,
5249  bool is_face,
5250  typename VectorizedArrayType>
5251 inline FEEvaluationAccess<dim,
5252  n_components_,
5253  Number,
5254  is_face,
5255  VectorizedArrayType>::
5256  FEEvaluationAccess(
5257  const Mapping<dim> & mapping,
5258  const FiniteElement<dim> &fe,
5259  const Quadrature<1> & quadrature,
5260  const UpdateFlags update_flags,
5261  const unsigned int first_selected_component,
5263  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5264  mapping,
5265  fe,
5266  quadrature,
5267  update_flags,
5268  first_selected_component,
5269  other)
5270 {}
5271 
5272 
5273 
5274 template <int dim,
5275  int n_components_,
5276  typename Number,
5277  bool is_face,
5278  typename VectorizedArrayType>
5279 inline FEEvaluationAccess<dim,
5280  n_components_,
5281  Number,
5282  is_face,
5283  VectorizedArrayType>::
5284  FEEvaluationAccess(const FEEvaluationAccess<dim,
5285  n_components_,
5286  Number,
5287  is_face,
5288  VectorizedArrayType> &other)
5289  : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5290  other)
5291 {}
5292 
5293 
5294 
5295 template <int dim,
5296  int n_components_,
5297  typename Number,
5298  bool is_face,
5299  typename VectorizedArrayType>
5300 inline FEEvaluationAccess<dim,
5301  n_components_,
5302  Number,
5303  is_face,
5304  VectorizedArrayType> &
5306 operator=(const FEEvaluationAccess<dim,
5307  n_components_,
5308  Number,
5309  is_face,
5310  VectorizedArrayType> &other)
5311 {
5312  this->FEEvaluationBase<dim,
5313  n_components_,
5314  Number,
5315  is_face,
5316  VectorizedArrayType>::operator=(other);
5317  return *this;
5318 }
5319 
5320 
5321 
5322 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5323 
5324 
5325 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5329  const unsigned int dof_no,
5330  const unsigned int first_selected_component,
5331  const unsigned int quad_no,
5332  const unsigned int fe_degree,
5333  const unsigned int n_q_points,
5334  const bool is_interior_face,
5335  const unsigned int active_fe_index,
5336  const unsigned int active_quad_index,
5337  const unsigned int face_type)
5338  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5339  matrix_free,
5340  dof_no,
5341  first_selected_component,
5342  quad_no,
5343  fe_degree,
5344  n_q_points,
5345  is_interior_face,
5346  active_fe_index,
5347  active_quad_index,
5348  face_type)
5349 {}
5350 
5351 
5352 
5353 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5356  const Mapping<dim> & mapping,
5357  const FiniteElement<dim> &fe,
5358  const Quadrature<1> & quadrature,
5359  const UpdateFlags update_flags,
5360  const unsigned int first_selected_component,
5362  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5363  mapping,
5364  fe,
5365  quadrature,
5366  update_flags,
5367  first_selected_component,
5368  other)
5369 {}
5370 
5371 
5372 
5373 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5377  &other)
5378  : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5379 {}
5380 
5381 
5382 
5383 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5387 {
5388  this
5389  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5390  other);
5391  return *this;
5392 }
5393 
5394 
5395 
5396 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5397 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5399  const unsigned int dof) const
5400 {
5401  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5402  return this->values_dofs[dof];
5403 }
5404 
5405 
5406 
5407 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5408 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5410  const unsigned int q_point) const
5411 {
5412 # ifdef DEBUG
5413  Assert(this->values_quad_initialized == true,
5415 # endif
5416  AssertIndexRange(q_point, this->n_quadrature_points);
5417  return this->values_quad[q_point];
5418 }
5419 
5420 
5421 
5422 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5423 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5425  get_normal_derivative(const unsigned int q_point) const
5426 {
5427  return BaseClass::get_normal_derivative(q_point)[0];
5428 }
5429 
5430 
5431 
5432 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5435  const unsigned int q_point) const
5436 {
5437  // could use the base class gradient, but that involves too many expensive
5438  // initialization operations on tensors
5439 
5440 # ifdef DEBUG
5441  Assert(this->gradients_quad_initialized == true,
5443 # endif
5444  AssertIndexRange(q_point, this->n_quadrature_points);
5445 
5446  Assert(this->jacobian != nullptr,
5448  "update_gradients"));
5449 
5451 
5452  const std::size_t nqp = this->n_quadrature_points;
5453  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5454  {
5455  for (unsigned int d = 0; d < dim; ++d)
5456  grad_out[d] =
5457  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5458  }
5459  // cell with general/affine Jacobian
5460  else
5461  {
5463  this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5464  q_point :
5465  0];
5466  for (unsigned int d = 0; d < dim; ++d)
5467  {
5468  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5469  for (unsigned int e = 1; e < dim; ++e)
5470  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5471  }
5472  }
5473  return grad_out;
5474 }
5475 
5476 
5477 
5478 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5481  const unsigned int q_point) const
5482 {
5483  return BaseClass::get_hessian(q_point)[0];
5484 }
5485 
5486 
5487 
5488 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5491  get_hessian_diagonal(const unsigned int q_point) const
5492 {
5493  return BaseClass::get_hessian_diagonal(q_point)[0];
5494 }
5495 
5496 
5497 
5498 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5499 inline VectorizedArrayType
5501  const unsigned int q_point) const
5502 {
5503  return BaseClass::get_laplacian(q_point)[0];
5504 }
5505 
5506 
5507 
5508 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5509 inline void DEAL_II_ALWAYS_INLINE
5511  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5512 {
5513 # ifdef DEBUG
5514  this->dof_values_initialized = true;
5515  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5516 # endif
5517  this->values_dofs[dof] = val_in;
5518 }
5519 
5520 
5521 
5522 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5523 inline void DEAL_II_ALWAYS_INLINE
5525  const VectorizedArrayType val_in,
5526  const unsigned int q_point)
5527 {
5528 # ifdef DEBUG
5529  Assert(this->is_reinitialized, ExcNotInitialized());
5530 # endif
5531  AssertIndexRange(q_point, this->n_quadrature_points);
5532  Assert(this->J_value != nullptr,
5534  "update_value"));
5535 # ifdef DEBUG
5536  this->values_quad_submitted = true;
5537 # endif
5538 
5539  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5540  {
5541  const VectorizedArrayType JxW =
5542  this->J_value[0] * this->quadrature_weights[q_point];
5543  this->values_quad[q_point] = val_in * JxW;
5544  }
5545  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5546  {
5547  this->values_quad[q_point] = val_in * this->J_value[q_point];
5548  }
5549 }
5550 
5551 
5552 
5553 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5554 inline DEAL_II_ALWAYS_INLINE void
5556  const Tensor<1, 1, VectorizedArrayType> val_in,
5557  const unsigned int q_point)
5558 {
5559  submit_value(val_in[0], q_point);
5560 }
5561 
5562 
5563 
5564 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5565 inline DEAL_II_ALWAYS_INLINE void
5567  submit_normal_derivative(const VectorizedArrayType grad_in,
5568  const unsigned int q_point)
5569 {
5571  grad[0] = grad_in;
5572  BaseClass::submit_normal_derivative(grad, q_point);
5573 }
5574 
5575 
5576 
5577 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5578 inline DEAL_II_ALWAYS_INLINE void
5581  const unsigned int q_point)
5582 {
5583 # ifdef DEBUG
5584  Assert(this->is_reinitialized, ExcNotInitialized());
5585 # endif
5586  AssertIndexRange(q_point, this->n_quadrature_points);
5587  Assert(this->J_value != nullptr,
5589  "update_gradients"));
5590  Assert(this->jacobian != nullptr,
5592  "update_gradients"));
5593 # ifdef DEBUG
5594  this->gradients_quad_submitted = true;
5595 # endif
5596 
5597  const std::size_t nqp = this->n_quadrature_points;
5598  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5599  {
5600  const VectorizedArrayType JxW =
5601  this->J_value[0] * this->quadrature_weights[q_point];
5602  for (unsigned int d = 0; d < dim; ++d)
5603  this->gradients_quad[d * nqp + q_point] =
5604  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5605  }
5606  // general/affine cell type
5607  else
5608  {
5610  this->cell_type > internal::MatrixFreeFunctions::affine ?
5611  this->jacobian[q_point] :
5612  this->jacobian[0];
5613  const VectorizedArrayType JxW =
5614  this->cell_type > internal::MatrixFreeFunctions::affine ?
5615  this->J_value[q_point] :
5616  this->J_value[0] * this->quadrature_weights[q_point];
5617  for (unsigned int d = 0; d < dim; ++d)
5618  {
5619  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5620  for (unsigned int e = 1; e < dim; ++e)
5621  new_val += jac[e][d] * grad_in[e];
5622  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
5623  }
5624  }
5625 }
5626 
5627 
5628 
5629 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5630 inline DEAL_II_ALWAYS_INLINE void
5633  const unsigned int q_point)
5634 {
5636  hessian[0] = hessian_in;
5637  BaseClass::submit_hessian(hessian, q_point);
5638 }
5639 
5640 
5641 
5642 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5643 inline VectorizedArrayType
5645  integrate_value() const
5646 {
5647  return BaseClass::integrate_value()[0];
5648 }
5649 
5650 
5651 
5652 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5653 
5654 
5655 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5659  const unsigned int dof_no,
5660  const unsigned int first_selected_component,
5661  const unsigned int quad_no,
5662  const unsigned int fe_degree,
5663  const unsigned int n_q_points,
5664  const bool is_interior_face,
5665  const unsigned int active_fe_index,
5666  const unsigned int active_quad_index,
5667  const unsigned int face_type)
5668  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5669  matrix_free,
5670  dof_no,
5671  first_selected_component,
5672  quad_no,
5673  fe_degree,
5674  n_q_points,
5675  is_interior_face,
5676  active_fe_index,
5677  active_quad_index,
5678  face_type)
5679 {}
5680 
5681 
5682 
5683 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5686  const Mapping<dim> & mapping,
5687  const FiniteElement<dim> &fe,
5688  const Quadrature<1> & quadrature,
5689  const UpdateFlags update_flags,
5690  const unsigned int first_selected_component,
5692  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5693  mapping,
5694  fe,
5695  quadrature,
5696  update_flags,
5697  first_selected_component,
5698  other)
5699 {}
5700 
5701 
5702 
5703 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5707  &other)
5708  : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5709 {}
5710 
5711 
5712 
5713 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5717  &other)
5718 {
5720  operator=(other);
5721  return *this;
5722 }
5723 
5724 
5725 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5728  const unsigned int q_point) const
5729 {
5730  if (this->data->element_type ==
5732  {
5733  // Piola transform is required
5734 # ifdef DEBUG
5735  Assert(this->values_quad_initialized == true,
5737 # endif
5738 
5739  AssertIndexRange(q_point, this->n_quadrature_points);
5740  Assert(this->J_value != nullptr,
5742  "update_values"));
5743  const std::size_t nqp = this->n_quadrature_points;
5745 
5746  if (!is_face &&
5747  this->cell_type == internal::MatrixFreeFunctions::cartesian)
5748  {
5749  // Cartesian cell
5751  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, ::VectorizedArray<Number>> &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, ::VectorizedArray<Number>> &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  {
6085  this->jacobian[1];
6086  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6087 
6088  for (unsigned int comp = 0; comp < n_components; ++comp)
6089  this->values_quad[comp * nqp + q_point] =
6090  val_in[comp] * weight * jac[comp][comp];
6091  }
6092  else
6093  {
6094  // Affine or general cell
6095  const Tensor<2, dim, ::VectorizedArray<Number>> &inv_t_jac =
6096  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6097  this->jacobian[q_point] :
6098  this->jacobian[0];
6100  (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6101  transpose(invert(inv_t_jac)) :
6102  this->jacobian[1];
6103 
6104  // Derivatives are reordered for faces. Need to take this into account
6105  // and 1/inv_det != J_value for faces
6106  const VectorizedArrayType fac =
6107  (!is_face) ?
6108  this->quadrature_weights[q_point] :
6109  (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6110  this->J_value[q_point] :
6111  this->J_value[0] * this->quadrature_weights[q_point]) *
6112  ((dim == 2 && this->get_face_no() < 2) ?
6113  -determinant(inv_t_jac) :
6114  determinant(inv_t_jac)));
6115 
6116  // J^T * u * factor
6117  for (unsigned int comp = 0; comp < n_components; ++comp)
6118  {
6119  this->values_quad[comp * nqp + q_point] =
6120  val_in[0] * jac[0][comp] * fac;
6121  for (unsigned int e = 1; e < dim; ++e)
6122  this->values_quad[comp * nqp + q_point] +=
6123  val_in[e] * jac[e][comp] * fac;
6124  }
6125  }
6126  }
6127  else
6128  {
6129  // No Piola transform
6130  BaseClass::submit_value(val_in, q_point);
6131  }
6132 }
6133 
6134 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6135 inline DEAL_II_ALWAYS_INLINE void
6138  const unsigned int q_point)
6139 {
6140  if (this->data->element_type ==
6142  {
6143  // Piola transform is required
6144 
6145 # ifdef DEBUG
6146  Assert(this->is_reinitialized, ExcNotInitialized());
6147 # endif
6148  AssertIndexRange(q_point, this->n_quadrature_points);
6149  Assert(this->J_value != nullptr,
6151  "update_gradients"));
6152  Assert(this->jacobian != nullptr,
6154  "update_gradients"));
6155 # ifdef DEBUG
6156  this->gradients_quad_submitted = true;
6157 # endif
6158 
6159  const std::size_t nqp = this->n_quadrature_points;
6160  if (!is_face &&
6161  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6162  {
6163  // Cartesian cell
6164  const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6165  this->jacobian[0];
6166  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6167  const VectorizedArrayType weight = this->quadrature_weights[q_point];
6168  for (unsigned int d = 0; d < dim; ++d)
6169  for (unsigned int comp = 0; comp < n_components; ++comp)
6170  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6171  grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
6172  }
6173  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6174  {
6175  // Affine cell
6176  const Tensor<2, dim, ::VectorizedArray<Number>> &inv_t_jac =
6177  this->jacobian[0];
6178  const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6179 
6180  // Derivatives are reordered for faces. Need to take this into account
6181  // and 1/inv_det != J_value for faces
6182  const VectorizedArrayType fac =
6183  (!is_face) ? this->quadrature_weights[q_point] :
6184  this->J_value[0] * this->quadrature_weights[q_point] *
6185  ((dim == 2 && this->get_face_no() < 2) ?
6186  -determinant(inv_t_jac) :
6187  determinant(inv_t_jac));
6188 
6189  // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6190  for (unsigned int comp = 0; comp < n_components; ++comp)
6191  for (unsigned int d = 0; d < dim; ++d)
6192  {
6193  VectorizedArrayType tmp = 0;
6194  for (unsigned int f = 0; f < dim; ++f)
6195  for (unsigned int e = 0; e < dim; ++e)
6196  tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
6197 
6198  this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
6199  }
6200  }
6201  else
6202  {
6203  // General cell
6204  AssertThrow(false, ExcNotImplemented());
6205  }
6206  }
6207  else
6208  {
6209  BaseClass::submit_gradient(grad_in, q_point);
6210  }
6211 }
6212 
6213 
6214 
6215 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6216 inline DEAL_II_ALWAYS_INLINE void
6219  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6220  const unsigned int q_point)
6221 {
6222  if (this->data->element_type ==
6224  {
6225  // Piola transform is required
6226  const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
6228  submit_gradient(grad, q_point);
6229  }
6230  else
6231  {
6232  BaseClass::submit_gradient(grad_in, q_point);
6233  }
6234 }
6235 
6236 
6237 
6238 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6239 inline DEAL_II_ALWAYS_INLINE void
6241  submit_divergence(const VectorizedArrayType div_in,
6242  const unsigned int q_point)
6243 {
6244 # ifdef DEBUG
6245  Assert(this->is_reinitialized, ExcNotInitialized());
6246 # endif
6247  AssertIndexRange(q_point, this->n_quadrature_points);
6248  Assert(this->J_value != nullptr,
6250  "update_gradients"));
6251  Assert(this->jacobian != nullptr,
6253  "update_gradients"));
6254 # ifdef DEBUG
6255  this->gradients_quad_submitted = true;
6256 # endif
6257 
6258  const std::size_t nqp = this->n_quadrature_points;
6259  if (this->data->element_type ==
6261  {
6262  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6263  {
6264  // Affine cell
6265 
6266  // Derivatives are reordered for faces. Need to take this into account
6267  // and 1/inv_det != J_value for faces
6268  const VectorizedArrayType fac =
6269  ((!is_face) ?
6270  1 :
6271  this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
6272  -determinant(this->jacobian[0]) :
6273  determinant(this->jacobian[0]))) *
6274  this->quadrature_weights[q_point] * div_in;
6275 
6276  for (unsigned int d = 0; d < dim; ++d)
6277  {
6278  this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
6279  for (unsigned int e = d + 1; e < dim; ++e)
6280  {
6281  this->gradients_quad[(dim * d + e) * nqp + q_point] =
6282  VectorizedArrayType();
6283  this->gradients_quad[(dim * e + d) * nqp + q_point] =
6284  VectorizedArrayType();
6285  }
6286  }
6287  }
6288  else
6289  {
6290  // General cell
6291  AssertThrow(false, ExcNotImplemented());
6292  }
6293  }
6294  else
6295  {
6296  if (!is_face &&
6297  this->cell_type == internal::MatrixFreeFunctions::cartesian)
6298  {
6299  const VectorizedArrayType fac =
6300  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6301  for (unsigned int d = 0; d < dim; ++d)
6302  {
6303  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6304  (fac * this->jacobian[0][d][d]);
6305  for (unsigned int e = d + 1; e < dim; ++e)
6306  {
6307  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6308  VectorizedArrayType();
6309  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6310  VectorizedArrayType();
6311  }
6312  }
6313  }
6314  else
6315  {
6317  this->cell_type == internal::MatrixFreeFunctions::general ?
6318  this->jacobian[q_point] :
6319  this->jacobian[0];
6320  const VectorizedArrayType fac =
6321  (this->cell_type == internal::MatrixFreeFunctions::general ?
6322  this->J_value[q_point] :
6323  this->J_value[0] * this->quadrature_weights[q_point]) *
6324  div_in;
6325  for (unsigned int d = 0; d < dim; ++d)
6326  {
6327  for (unsigned int e = 0; e < dim; ++e)
6328  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6329  jac[d][e] * fac;
6330  }
6331  }
6332  }
6333 }
6334 
6335 
6336 
6337 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6338 inline DEAL_II_ALWAYS_INLINE void
6342  const unsigned int q_point)
6343 {
6344  AssertThrow(
6345  this->data->element_type !=
6347  ExcNotImplemented());
6348 
6349  // could have used base class operator, but that involves some overhead
6350  // which is inefficient. it is nice to have the symmetric tensor because
6351  // that saves some operations
6352 # ifdef DEBUG
6353  Assert(this->is_reinitialized, ExcNotInitialized());
6354 # endif
6355  AssertIndexRange(q_point, this->n_quadrature_points);
6356  Assert(this->J_value != nullptr,
6358  "update_gradients"));
6359  Assert(this->jacobian != nullptr,
6361  "update_gradients"));
6362 # ifdef DEBUG
6363  this->gradients_quad_submitted = true;
6364 # endif
6365 
6366  const std::size_t nqp = this->n_quadrature_points;
6367  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6368  {
6369  const VectorizedArrayType JxW =
6370  this->J_value[0] * this->quadrature_weights[q_point];
6371  for (unsigned int d = 0; d < dim; ++d)
6372  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6373  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6374  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6375  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6376  {
6377  const VectorizedArrayType value =
6378  sym_grad.access_raw_entry(counter) * JxW;
6379  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6380  value * this->jacobian[0][d][d];
6381  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6382  value * this->jacobian[0][e][e];
6383  }
6384  }
6385  // general/affine cell type
6386  else
6387  {
6388  const VectorizedArrayType JxW =
6389  this->cell_type == internal::MatrixFreeFunctions::general ?
6390  this->J_value[q_point] :
6391  this->J_value[0] * this->quadrature_weights[q_point];
6393  this->cell_type == internal::MatrixFreeFunctions::general ?
6394  this->jacobian[q_point] :
6395  this->jacobian[0];
6396  VectorizedArrayType weighted[dim][dim];
6397  for (unsigned int i = 0; i < dim; ++i)
6398  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6399  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6400  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6401  {
6402  const VectorizedArrayType value =
6403  sym_grad.access_raw_entry(counter) * JxW;
6404  weighted[i][j] = value;
6405  weighted[j][i] = value;
6406  }
6407  for (unsigned int comp = 0; comp < dim; ++comp)
6408  for (unsigned int d = 0; d < dim; ++d)
6409  {
6410  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6411  for (unsigned int e = 1; e < dim; ++e)
6412  new_val += jac[e][d] * weighted[comp][e];
6413  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6414  }
6415  }
6416 }
6417 
6418 
6419 
6420 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6421 inline DEAL_II_ALWAYS_INLINE void
6424  const unsigned int q_point)
6425 {
6427  switch (dim)
6428  {
6429  case 1:
6430  Assert(false,
6431  ExcMessage(
6432  "Testing by the curl in 1d is not a useful operation"));
6433  break;
6434  case 2:
6435  grad[1][0] = curl[0];
6436  grad[0][1] = -curl[0];
6437  break;
6438  case 3:
6439  grad[2][1] = curl[0];
6440  grad[1][2] = -curl[0];
6441  grad[0][2] = curl[1];
6442  grad[2][0] = -curl[1];
6443  grad[1][0] = curl[2];
6444  grad[0][1] = -curl[2];
6445  break;
6446  default:
6447  Assert(false, ExcNotImplemented());
6448  }
6449  submit_gradient(grad, q_point);
6450 }
6451 
6452 
6453 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6454 
6455 
6456 template <typename Number, bool is_face, typename VectorizedArrayType>
6459  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
6460  const unsigned int dof_no,
6461  const unsigned int first_selected_component,
6462  const unsigned int quad_no,
6463  const unsigned int fe_degree,
6464  const unsigned int n_q_points,
6465  const bool is_interior_face,
6466  const unsigned int active_fe_index,
6467  const unsigned int active_quad_index,
6468  const unsigned int face_type)
6469  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6470  matrix_free,
6471  dof_no,
6472  first_selected_component,
6473  quad_no,
6474  fe_degree,
6475  n_q_points,
6476  is_interior_face,
6477  active_fe_index,
6478  active_quad_index,
6479  face_type)
6480 {}
6481 
6482 
6483 
6484 template <typename Number, bool is_face, typename VectorizedArrayType>
6487  const Mapping<1> & mapping,
6488  const FiniteElement<1> &fe,
6489  const Quadrature<1> & quadrature,
6490  const UpdateFlags update_flags,
6491  const unsigned int first_selected_component,
6493  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6494  mapping,
6495  fe,
6496  quadrature,
6497  update_flags,
6498  first_selected_component,
6499  other)
6500 {}
6501 
6502 
6503 
6504 template <typename Number, bool is_face, typename VectorizedArrayType>
6508  : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6509 {}
6510 
6511 
6512 
6513 template <typename Number, bool is_face, typename VectorizedArrayType>
6517 {
6519  other);
6520  return *this;
6521 }
6522 
6523 
6524 
6525 template <typename Number, bool is_face, typename VectorizedArrayType>
6526 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6528  const unsigned int dof) const
6529 {
6530  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6531  return this->values_dofs[dof];
6532 }
6533 
6534 
6535 
6536 template <typename Number, bool is_face, typename VectorizedArrayType>
6537 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6539  const unsigned int q_point) const
6540 {
6541 # ifdef DEBUG
6542  Assert(this->values_quad_initialized == true,
6544 # endif
6545  AssertIndexRange(q_point, this->n_quadrature_points);
6546  return this->values_quad[q_point];
6547 }
6548 
6549 
6550 
6551 template <typename Number, bool is_face, typename VectorizedArrayType>
6554  const unsigned int q_point) const
6555 {
6556  // could use the base class gradient, but that involves too many inefficient
6557  // initialization operations on tensors
6558 
6559 # ifdef DEBUG
6560  Assert(this->gradients_quad_initialized == true,
6562 # endif
6563  AssertIndexRange(q_point, this->n_quadrature_points);
6564 
6566  this->cell_type == internal::MatrixFreeFunctions::general ?
6567  this->jacobian[q_point] :
6568  this->jacobian[0];
6569 
6571  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6572 
6573  return grad_out;
6574 }
6575 
6576 
6577 
6578 template <typename Number, bool is_face, typename VectorizedArrayType>
6579 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6581  const unsigned int q_point) const
6582 {
6583  return get_gradient(q_point)[0];
6584 }
6585 
6586 
6587 
6588 template <typename Number, bool is_face, typename VectorizedArrayType>
6589 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6591  get_normal_derivative(const unsigned int q_point) const
6592 {
6593  return BaseClass::get_normal_derivative(q_point)[0];
6594 }
6595 
6596 
6597 
6598 template <typename Number, bool is_face, typename VectorizedArrayType>
6601  const unsigned int q_point) const
6602 {
6603  return BaseClass::get_hessian(q_point)[0];
6604 }
6605 
6606 
6607 
6608 template <typename Number, bool is_face, typename VectorizedArrayType>
6611  get_hessian_diagonal(const unsigned int q_point) const
6612 {
6613  return BaseClass::get_hessian_diagonal(q_point)[0];
6614 }
6615 
6616 
6617 
6618 template <typename Number, bool is_face, typename VectorizedArrayType>
6619 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6621  const unsigned int q_point) const
6622 {
6623  return BaseClass::get_laplacian(q_point)[0];
6624 }
6625 
6626 
6627 
6628 template <typename Number, bool is_face, typename VectorizedArrayType>
6631  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6632 {
6633 # ifdef DEBUG
6634  this->dof_values_initialized = true;
6635  AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6636 # endif
6637  this->values_dofs[dof] = val_in;
6638 }
6639 
6640 
6641 
6642 template <typename Number, bool is_face, typename VectorizedArrayType>
6643 inline DEAL_II_ALWAYS_INLINE void
6645  const VectorizedArrayType val_in,
6646  const unsigned int q_point)
6647 {
6648 # ifdef DEBUG
6649  Assert(this->is_reinitialized, ExcNotInitialized());
6650 # endif
6651  AssertIndexRange(q_point, this->n_quadrature_points);
6652 # ifdef DEBUG
6653  this->values_quad_submitted = true;
6654 # endif
6655 
6656  if (this->cell_type == internal::MatrixFreeFunctions::general)
6657  {
6658  const VectorizedArrayType JxW = this->J_value[q_point];
6659  this->values_quad[q_point] = val_in * JxW;
6660  }
6661  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6662  {
6663  const VectorizedArrayType JxW =
6664  this->J_value[0] * this->quadrature_weights[q_point];
6665  this->values_quad[q_point] = val_in * JxW;
6666  }
6667 }
6668 
6669 
6670 
6671 template <typename Number, bool is_face, typename VectorizedArrayType>
6672 inline DEAL_II_ALWAYS_INLINE void
6674  const Tensor<1, 1, VectorizedArrayType> val_in,
6675  const unsigned int q_point)
6676 {
6677  submit_value(val_in[0], q_point);
6678 }
6679 
6680 
6681 
6682 template <typename Number, bool is_face, typename VectorizedArrayType>
6683 inline DEAL_II_ALWAYS_INLINE void
6685  const Tensor<1, 1, VectorizedArrayType> grad_in,
6686  const unsigned int q_point)
6687 {
6688  submit_gradient(grad_in[0], q_point);
6689 }
6690 
6691 
6692 
6693 template <typename Number, bool is_face, typename VectorizedArrayType>
6694 inline DEAL_II_ALWAYS_INLINE void
6696  const VectorizedArrayType grad_in,
6697  const unsigned int q_point)
6698 {
6699 # ifdef DEBUG
6700  Assert(this->is_reinitialized, ExcNotInitialized());
6701 # endif
6702  AssertIndexRange(q_point, this->n_quadrature_points);
6703 # ifdef DEBUG
6704  this->gradients_quad_submitted = true;
6705 # endif
6706 
6708  this->cell_type == internal::MatrixFreeFunctions::general ?
6709  this->jacobian[q_point] :
6710  this->jacobian[0];
6711  const VectorizedArrayType JxW =
6712  this->cell_type == internal::MatrixFreeFunctions::general ?
6713  this->J_value[q_point] :
6714  this->J_value[0] * this->quadrature_weights[q_point];
6715 
6716  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6717 }
6718 
6719 
6720 
6721 template <typename Number, bool is_face, typename VectorizedArrayType>
6722 inline DEAL_II_ALWAYS_INLINE void
6724  const Tensor<2, 1, VectorizedArrayType> grad_in,
6725  const unsigned int q_point)
6726 {
6727  submit_gradient(grad_in[0][0], q_point);
6728 }
6729 
6730 
6731 
6732 template <typename Number, bool is_face, typename VectorizedArrayType>
6733 inline DEAL_II_ALWAYS_INLINE void
6735  submit_normal_derivative(const VectorizedArrayType grad_in,
6736  const unsigned int q_point)
6737 {
6739  grad[0] = grad_in;
6740  BaseClass::submit_normal_derivative(grad, q_point);
6741 }
6742 
6743 
6744 
6745 template <typename Number, bool is_face, typename VectorizedArrayType>
6746 inline DEAL_II_ALWAYS_INLINE void
6749  const unsigned int q_point)
6750 {
6751  BaseClass::submit_normal_derivative(grad_in, q_point);
6752 }
6753 
6754 
6755 template <typename Number, bool is_face, typename VectorizedArrayType>
6756 inline DEAL_II_ALWAYS_INLINE void
6758  const Tensor<2, 1, VectorizedArrayType> hessian_in,
6759  const unsigned int q_point)
6760 {
6762  hessian[0] = hessian_in;
6763  BaseClass::submit_hessian(hessian, q_point);
6764 }
6765 
6766 
6767 template <typename Number, bool is_face, typename VectorizedArrayType>
6768 inline VectorizedArrayType
6770  integrate_value() const
6771 {
6772  return BaseClass::integrate_value()[0];
6773 }
6774 
6775 
6776 
6777 /*-------------------------- FEEvaluation -----------------------------------*/
6778 
6779 
6780 template <int dim,
6781  int fe_degree,
6782  int n_q_points_1d,
6783  int n_components_,
6784  typename Number,
6785  typename VectorizedArrayType>
6786 inline FEEvaluation<dim,
6787  fe_degree,
6788  n_q_points_1d,
6789  n_components_,
6790  Number,
6791  VectorizedArrayType>::
6792  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6793  const unsigned int fe_no,
6794  const unsigned int quad_no,
6795  const unsigned int first_selected_component,
6796  const unsigned int active_fe_index,
6797  const unsigned int active_quad_index)
6798  : BaseClass(matrix_free,
6799  fe_no,
6800  first_selected_component,
6801  quad_no,
6802  fe_degree,
6803  static_n_q_points,
6804  true /*note: this is not a face*/,
6805  active_fe_index,
6806  active_quad_index)
6807  , dofs_per_component(this->data->dofs_per_component_on_cell)
6808  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6809  , n_q_points(this->data->n_q_points)
6810 {
6811  check_template_arguments(fe_no, 0);
6812 }
6813 
6814 
6815 
6816 template <int dim,
6817  int fe_degree,
6818  int n_q_points_1d,
6819  int n_components_,
6820  typename Number,
6821  typename VectorizedArrayType>
6822 inline FEEvaluation<dim,
6823  fe_degree,
6824  n_q_points_1d,
6825  n_components_,
6826  Number,
6827  VectorizedArrayType>::
6828  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6829  const std::pair<unsigned int, unsigned int> & range,
6830  const unsigned int dof_no,
6831  const unsigned int quad_no,
6832  const unsigned int first_selected_component)
6833  : FEEvaluation(matrix_free,
6834  dof_no,
6835  quad_no,
6836  first_selected_component,
6837  matrix_free.get_cell_active_fe_index(range))
6838 {}
6839 
6840 
6841 
6842 template <int dim,
6843  int fe_degree,
6844  int n_q_points_1d,
6845  int n_components_,
6846  typename Number,
6847  typename VectorizedArrayType>
6848 inline FEEvaluation<dim,
6849  fe_degree,
6850  n_q_points_1d,
6851  n_components_,
6852  Number,
6853  VectorizedArrayType>::
6854  FEEvaluation(const Mapping<dim> & mapping,
6855  const FiniteElement<dim> &fe,
6856  const Quadrature<1> & quadrature,
6857  const UpdateFlags update_flags,
6858  const unsigned int first_selected_component)
6859  : BaseClass(mapping,
6860  fe,
6861  quadrature,
6862  update_flags,
6863  first_selected_component,
6864  nullptr)
6865  , dofs_per_component(this->data->dofs_per_component_on_cell)
6866  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6867  , n_q_points(this->data->n_q_points)
6868 {
6869  check_template_arguments(numbers::invalid_unsigned_int, 0);
6870 }
6871 
6872 
6873 
6874 template <int dim,
6875  int fe_degree,
6876  int n_q_points_1d,
6877  int n_components_,
6878  typename Number,
6879  typename VectorizedArrayType>
6880 inline FEEvaluation<dim,
6881  fe_degree,
6882  n_q_points_1d,
6883  n_components_,
6884  Number,
6885  VectorizedArrayType>::
6886  FEEvaluation(const FiniteElement<dim> &fe,
6887  const Quadrature<1> & quadrature,
6888  const UpdateFlags update_flags,
6889  const unsigned int first_selected_component)
6890  : BaseClass(StaticMappingQ1<dim>::mapping,
6891  fe,
6892  quadrature,
6893  update_flags,
6894  first_selected_component,
6895  nullptr)
6896  , dofs_per_component(this->data->dofs_per_component_on_cell)
6897  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6898  , n_q_points(this->data->n_q_points)
6899 {
6900  check_template_arguments(numbers::invalid_unsigned_int, 0);
6901 }
6902 
6903 
6904 
6905 template <int dim,
6906  int fe_degree,
6907  int n_q_points_1d,
6908  int n_components_,
6909  typename Number,
6910  typename VectorizedArrayType>
6911 inline FEEvaluation<dim,
6912  fe_degree,
6913  n_q_points_1d,
6914  n_components_,
6915  Number,
6916  VectorizedArrayType>::
6917  FEEvaluation(const FiniteElement<dim> & fe,
6919  const unsigned int first_selected_component)
6920  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6921  fe,
6922  other.mapped_geometry->get_quadrature(),
6923  other.mapped_geometry->get_fe_values().get_update_flags(),
6924  first_selected_component,
6925  &other)
6926  , dofs_per_component(this->data->dofs_per_component_on_cell)
6927  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6928  , n_q_points(this->data->n_q_points)
6929 {
6930  check_template_arguments(numbers::invalid_unsigned_int, 0);
6931 }
6932 
6933 
6934 
6935 template <int dim,
6936  int fe_degree,
6937  int n_q_points_1d,
6938  int n_components_,
6939  typename Number,
6940  typename VectorizedArrayType>
6941 inline FEEvaluation<dim,
6942  fe_degree,
6943  n_q_points_1d,
6944  n_components_,
6945  Number,
6946  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6947  &other)
6948  : BaseClass(other)
6949  , dofs_per_component(this->data->dofs_per_component_on_cell)
6950  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6951  , n_q_points(this->data->n_q_points)
6952 {
6953  check_template_arguments(numbers::invalid_unsigned_int, 0);
6954 }
6955 
6956 
6957 
6958 template <int dim,
6959  int fe_degree,
6960  int n_q_points_1d,
6961  int n_components_,
6962  typename Number,
6963  typename VectorizedArrayType>
6964 inline FEEvaluation<dim,
6965  fe_degree,
6966  n_q_points_1d,
6967  n_components_,
6968  Number,
6969  VectorizedArrayType> &
6970 FEEvaluation<dim,
6971  fe_degree,
6972  n_q_points_1d,
6973  n_components_,
6974  Number,
6975  VectorizedArrayType>::operator=(const FEEvaluation &other)
6976 {
6977  BaseClass::operator=(other);
6978  check_template_arguments(numbers::invalid_unsigned_int, 0);
6979  return *this;
6980 }
6981 
6982 
6983 
6984 template <int dim,
6985  int fe_degree,
6986  int n_q_points_1d,
6987  int n_components_,
6988  typename Number,
6989  typename VectorizedArrayType>
6990 inline void
6991 FEEvaluation<dim,
6992  fe_degree,
6993  n_q_points_1d,
6994  n_components_,
6995  Number,
6996  VectorizedArrayType>::
6997  check_template_arguments(const unsigned int dof_no,
6998  const unsigned int first_selected_component)
6999 {
7000  (void)dof_no;
7001  (void)first_selected_component;
7002 
7003  Assert(
7004  this->data->dofs_per_component_on_cell > 0,
7005  ExcMessage(
7006  "There is nothing useful you can do with an FEEvaluation object with "
7007  "FE_Nothing, i.e., without DoFs! If you have passed to "
7008  "MatrixFree::reinit() a collection of finite elements also containing "
7009  "FE_Nothing, please check - before creating FEEvaluation - the category "
7010  "of the current range by calling either "
7011  "MatrixFree::get_cell_range_category(range) or "
7012  "MatrixFree::get_face_range_category(range). The returned category "
7013  "is the index of the active FE, which you can use to exclude "
7014  "FE_Nothing."));
7015 
7016 # ifdef DEBUG
7017  // print error message when the dimensions do not match. Propose a possible
7018  // fix
7019  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7020  static_cast<unsigned int>(fe_degree) !=
7021  this->data->data.front().fe_degree) ||
7022  n_q_points != this->n_quadrature_points)
7023  {
7024  std::string message =
7025  "-------------------------------------------------------\n";
7026  message += "Illegal arguments in constructor/wrong template arguments!\n";
7027  message += " Called --> FEEvaluation<dim,";
7028  message += Utilities::int_to_string(fe_degree) + ",";
7029  message += Utilities::int_to_string(n_q_points_1d);
7030  message += "," + Utilities::int_to_string(n_components);
7031  message += ",Number>(data";
7032  if (first_selected_component != numbers::invalid_unsigned_int)
7033  {
7034  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7035  message += Utilities::int_to_string(this->quad_no) + ", ";
7036  message += Utilities::int_to_string(first_selected_component);
7037  }
7038  message += ")\n";
7039 
7040  // check whether some other vector component has the correct number of
7041  // points
7042  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7043  proposed_fe_comp = numbers::invalid_unsigned_int,
7044  proposed_quad_comp = numbers::invalid_unsigned_int;
7045  if (dof_no != numbers::invalid_unsigned_int)
7046  {
7047  if (static_cast<unsigned int>(fe_degree) ==
7048  this->data->data.front().fe_degree)
7049  {
7050  proposed_dof_comp = dof_no;
7051  proposed_fe_comp = first_selected_component;
7052  }
7053  else
7054  for (unsigned int no = 0; no < this->matrix_free->n_components();
7055  ++no)
7056  for (unsigned int nf = 0;
7057  nf < this->matrix_free->n_base_elements(no);
7058  ++nf)
7059  if (this->matrix_free
7060  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7061  .data.front()
7062  .fe_degree == static_cast<unsigned int>(fe_degree))
7063  {
7064  proposed_dof_comp = no;
7065  proposed_fe_comp = nf;
7066  break;
7067  }
7068  if (n_q_points ==
7069  this->mapping_data->descriptor[this->active_quad_index]
7070  .n_q_points)
7071  proposed_quad_comp = this->quad_no;
7072  else
7073  for (unsigned int no = 0;
7074  no < this->matrix_free->get_mapping_info().cell_data.size();
7075  ++no)
7076  if (this->matrix_free->get_mapping_info()
7077  .cell_data[no]
7078  .descriptor[this->active_quad_index]
7079  .n_q_points == n_q_points)
7080  {
7081  proposed_quad_comp = no;
7082  break;
7083  }
7084  }
7085  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7086  proposed_quad_comp != numbers::invalid_unsigned_int)
7087  {
7088  if (proposed_dof_comp != first_selected_component)
7089  message += "Wrong vector component selection:\n";
7090  else
7091  message += "Wrong quadrature formula selection:\n";
7092  message += " Did you mean FEEvaluation<dim,";
7093  message += Utilities::int_to_string(fe_degree) + ",";
7094  message += Utilities::int_to_string(n_q_points_1d);
7095  message += "," + Utilities::int_to_string(n_components);
7096  message += ",Number>(data";
7097  if (dof_no != numbers::invalid_unsigned_int)
7098  {
7099  message +=
7100  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7101  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7102  message += Utilities::int_to_string(proposed_fe_comp);
7103  }
7104  message += ")?\n";
7105  std::string correct_pos;
7106  if (proposed_dof_comp != dof_no)
7107  correct_pos = " ^ ";
7108  else
7109  correct_pos = " ";
7110  if (proposed_quad_comp != this->quad_no)
7111  correct_pos += " ^ ";
7112  else
7113  correct_pos += " ";
7114  if (proposed_fe_comp != first_selected_component)
7115  correct_pos += " ^\n";
7116  else
7117  correct_pos += " \n";
7118  message += " " +
7119  correct_pos;
7120  }
7121  // ok, did not find the numbers specified by the template arguments in
7122  // the given list. Suggest correct template arguments
7123  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7124  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7125  message += "Wrong template arguments:\n";
7126  message += " Did you mean FEEvaluation<dim,";
7127  message +=
7128  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7129  message += Utilities::int_to_string(proposed_n_q_points_1d);
7130  message += "," + Utilities::int_to_string(n_components);
7131  message += ",Number>(data";
7132  if (dof_no != numbers::invalid_unsigned_int)
7133  {
7134  message += ", " + Utilities::int_to_string(dof_no) + ", ";
7135  message += Utilities::int_to_string(this->quad_no);
7136  message += ", " + Utilities::int_to_string(first_selected_component);
7137  }
7138  message += ")?\n";
7139  std::string correct_pos;
7140  if (this->data->data.front().fe_degree !=
7141  static_cast<unsigned int>(fe_degree))
7142  correct_pos = " ^";
7143  else
7144  correct_pos = " ";
7145  if (proposed_n_q_points_1d != n_q_points_1d)
7146  correct_pos += " ^\n";
7147  else
7148  correct_pos += " \n";
7149  message += " " + correct_pos;
7150 
7151  Assert(static_cast<unsigned int>(fe_degree) ==
7152  this->data->data.front().fe_degree &&
7153  n_q_points == this->n_quadrature_points,
7154  ExcMessage(message));
7155  }
7156  if (dof_no != numbers::invalid_unsigned_int)
7158  n_q_points,
7159  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7160 # endif
7161 }
7162 
7163 
7164 
7165 template <int dim,
7166  int fe_degree,
7167  int n_q_points_1d,
7168  int n_components_,
7169  typename Number,
7170  typename VectorizedArrayType>
7171 inline void
7172 FEEvaluation<dim,
7173  fe_degree,
7174  n_q_points_1d,
7175  n_components_,
7176  Number,
7177  VectorizedArrayType>::reinit(const unsigned int cell_index)
7178 {
7179  Assert(this->mapped_geometry == nullptr,
7180  ExcMessage("FEEvaluation was initialized without a matrix-free object."
7181  " Integer indexing is not possible"));
7182  if (this->mapped_geometry != nullptr)
7183  return;
7184 
7185  Assert(this->dof_info != nullptr, ExcNotInitialized());
7186  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7187  this->cell = cell_index;
7188  this->cell_type =
7189  this->matrix_free->get_mapping_info().get_cell_type(cell_index);
7190 
7191  const unsigned int offsets =
7192  this->mapping_data->data_index_offsets[cell_index];
7193  this->jacobian = &this->mapping_data->jacobians[0][offsets];
7194  this->J_value = &this->mapping_data->JxW_values[offsets];
7195  this->jacobian_gradients =
7196  this->mapping_data->jacobian_gradients[0].data() + offsets;
7197 
7198  unsigned int i = 0;
7199  for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7200  ++i)
7201  this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7202  for (; i < VectorizedArrayType::size(); ++i)
7203  this->cell_ids[i] = numbers::invalid_unsigned_int;
7204 
7205  if (this->mapping_data->quadrature_points.empty() == false)
7206  this->quadrature_points =
7207  &this->mapping_data->quadrature_points
7208  [this->mapping_data->quadrature_point_offsets[this->cell]];
7209 
7210 # ifdef DEBUG
7211  this->is_reinitialized = true;
7212  this->dof_values_initialized = false;
7213  this->values_quad_initialized = false;
7214  this->gradients_quad_initialized = false;
7215  this->hessians_quad_initialized = false;
7216 # endif
7217 }
7218 
7219 
7220 
7221 template <int dim,
7222  int fe_degree,
7223  int n_q_points_1d,
7224  int n_components_,
7225  typename Number,
7226  typename VectorizedArrayType>
7227 inline void
7228 FEEvaluation<dim,
7229  fe_degree,
7230  n_q_points_1d,
7231  n_components_,
7232  Number,
7233  VectorizedArrayType>::
7234  reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
7235 {
7236  Assert(this->dof_info != nullptr, ExcNotInitialized());
7237  Assert(this->mapping_data != nullptr, ExcNotInitialized());
7238 
7239  this->cell = numbers::invalid_unsigned_int;
7240  this->cell_ids = cell_ids;
7241 
7242  // determine type of cell batch
7244 
7245  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7246  {
7247  const unsigned int cell_index = cell_ids[v];
7248 
7250  continue;
7251 
7252  this->cell_type =
7253  std::max(this->cell_type,
7254  this->matrix_free->get_mapping_info().get_cell_type(
7255  cell_index / VectorizedArrayType::size()));
7256  }
7257 
7258  // allocate memory for internal data storage
7259  if (this->mapped_geometry == nullptr)
7260  this->mapped_geometry =
7261  std::make_shared<internal::MatrixFreeFunctions::
7262  MappingDataOnTheFly<dim, VectorizedArrayType>>();
7263 
7264  auto &mapping_storage = this->mapped_geometry->get_data_storage();
7265 
7266  auto &this_jacobian_data = mapping_storage.jacobians[0];
7267  auto &this_J_value_data = mapping_storage.JxW_values;
7268  auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7269  auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7270 
7272  {
7273  if (this->mapping_data->jacobians[0].size() > 0)
7274  this_jacobian_data.resize_fast(2);
7275 
7276  if (this->mapping_data->JxW_values.size() > 0)
7277  this_J_value_data.resize_fast(1);
7278 
7279  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7280  this_jacobian_gradients_data.resize_fast(1);
7281 
7282  if (this->mapping_data->quadrature_points.size() > 0)
7283  this_quadrature_points_data.resize_fast(1);
7284  }
7285  else
7286  {
7287  if (this->mapping_data->jacobians[0].size() > 0)
7288  this_jacobian_data.resize_fast(this->n_quadrature_points);
7289 
7290  if (this->mapping_data->JxW_values.size() > 0)
7291  this_J_value_data.resize_fast(this->n_quadrature_points);
7292 
7293  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7294  this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7295 
7296  if (this->mapping_data->quadrature_points.size() > 0)
7297  this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7298  }
7299 
7300  // set pointers to internal data storage
7301  this->jacobian = this_jacobian_data.data();
7302  this->J_value = this_J_value_data.data();
7303  this->jacobian_gradients = this_jacobian_gradients_data.data();
7304  this->quadrature_points = this_quadrature_points_data.data();
7305 
7306  // fill internal data storage lane by lane
7307  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7308  {
7309  const unsigned int cell_index = cell_ids[v];
7310 
7312  continue;
7313 
7314  const unsigned int cell_batch_index =
7315  cell_index / VectorizedArrayType::size();
7316  const unsigned int offsets =
7317  this->mapping_data->data_index_offsets[cell_batch_index];
7318  const unsigned int lane = cell_index % VectorizedArrayType::size();
7319 
7320  if (this->cell_type <=
7322  {
7323  // case that all cells are Cartesian or affine
7324  const unsigned int q = 0;
7325 
7326  if (this->mapping_data->JxW_values.size() > 0)
7327  this_J_value_data[q][v] =
7328  this->mapping_data->JxW_values[offsets + q][lane];
7329 
7330  if (this->mapping_data->jacobians[0].size() > 0)
7331  for (unsigned int q = 0; q < 2; ++q)
7332  for (unsigned int i = 0; i < dim; ++i)
7333  for (unsigned int j = 0; j < dim; ++j)
7334  this_jacobian_data[q][i][j][v] =
7335  this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7336 
7337  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7338  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7339  for (unsigned int j = 0; j < dim; ++j)
7340  this_jacobian_gradients_data[q][i][j][v] =
7341  this->mapping_data
7342  ->jacobian_gradients[0][offsets + q][i][j][lane];
7343 
7344  if (this->mapping_data->quadrature_points.size() > 0)
7345  for (unsigned int i = 0; i < dim; ++i)
7346  this_quadrature_points_data[q][i][v] =
7347  this->mapping_data->quadrature_points
7348  [this->mapping_data
7349  ->quadrature_point_offsets[cell_batch_index] +
7350  q][i][lane];
7351  }
7352  else
7353  {
7354  // general case that at least one cell is not Cartesian or affine
7355  const auto cell_type =
7356  this->matrix_free->get_mapping_info().get_cell_type(
7357  cell_batch_index);
7358 
7359  for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
7360  {
7361  const unsigned int q_src =
7362  (cell_type <=
7364  0 :
7365  q;
7366 
7367  if (this->mapping_data->JxW_values.size() > 0)
7368  this_J_value_data[q][v] =
7369  this->mapping_data->JxW_values[offsets + q_src][lane];
7370 
7371  if (this->mapping_data->jacobians[0].size() > 0)
7372  for (unsigned int i = 0; i < dim; ++i)
7373  for (unsigned int j = 0; j < dim; ++j)
7374  this_jacobian_data[q][i][j][v] =
7375  this->mapping_data
7376  ->jacobians[0][offsets + q_src][i][j][lane];
7377 
7378  if (this->mapping_data->jacobian_gradients[0].size() > 0)
7379  for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7380  for (unsigned int j = 0; j < dim; ++j)
7381  this_jacobian_gradients_data[q][i][j][v] =
7382  this->mapping_data
7383  ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7384 
7385  if (this->mapping_data->quadrature_points.size() > 0)
7386  {
7387  if (cell_type <=
7389  {
7390  // affine case: quadrature points are not available but
7391  // have to be computed from the corner point and the
7392  // Jacobian
7394  this->mapping_data->quadrature_points
7395  [this->mapping_data
7396  ->quadrature_point_offsets[cell_batch_index] +
7397  0];
7398 
7400  this->mapping_data->jacobians[0][offsets + 1];
7402  for (unsigned int d = 0; d < dim; ++d)
7403  point[d] +=
7404  jac[d][d] *
7405  static_cast<Number>(
7406  this->descriptor->quadrature.point(q)[d]);
7407  else
7408  for (unsigned int d = 0; d < dim; ++d)
7409  for (unsigned int e = 0; e < dim; ++e)
7410  point[d] +=
7411  jac[d][e] *
7412  static_cast<Number>(
7413  this->descriptor->quadrature.point(q)[e]);
7414 
7415  for (unsigned int i = 0; i < dim; ++i)
7416  this_quadrature_points_data[q][i][v] = point[i][lane];
7417  }
7418  else
7419  {
7420  // general case: quadrature points are available
7421  for (unsigned int i = 0; i < dim; ++i)
7422  this_quadrature_points_data[q][i][v] =
7423  this->mapping_data->quadrature_points
7424  [this->mapping_data
7425  ->quadrature_point_offsets[cell_batch_index] +
7426  q][i][lane];
7427  }
7428  }
7429  }
7430  }
7431  }
7432 
7433 # ifdef DEBUG
7434  this->is_reinitialized = true;
7435  this->dof_values_initialized = false;
7436  this->values_quad_initialized = false;
7437  this->gradients_quad_initialized = false;
7438  this->hessians_quad_initialized = false;
7439 # endif
7440 }
7441 
7442 
7443 
7444 template <int dim,
7445  int fe_degree,
7446  int n_q_points_1d,
7447  int n_components_,
7448  typename Number,
7449  typename VectorizedArrayType>
7450 template <bool level_dof_access>
7451 inline void
7452 FEEvaluation<dim,
7453  fe_degree,
7454  n_q_points_1d,
7455  n_components_,
7456  Number,
7457  VectorizedArrayType>::
7459 {
7460  Assert(this->matrix_free == nullptr,
7461  ExcMessage("Cannot use initialization from cell iterator if "
7462  "initialized from MatrixFree object. Use variant for "
7463  "on the fly computation with arguments as for FEValues "
7464  "instead"));
7465  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7466  this->mapped_geometry->reinit(
7467  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7468  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7469  if (level_dof_access)
7470  cell->get_mg_dof_indices(this->local_dof_indices);
7471  else
7472  cell->get_dof_indices(this->local_dof_indices);
7473 
7474 # ifdef DEBUG
7475  this->is_reinitialized = true;
7476 # endif
7477 }
7478 
7479 
7480 
7481 template <int dim,
7482  int fe_degree,
7483  int n_q_points_1d,
7484  int n_components_,
7485  typename Number,
7486  typename VectorizedArrayType>
7487 inline void
7488 FEEvaluation<dim,
7489  fe_degree,
7490  n_q_points_1d,
7491  n_components_,
7492  Number,
7493  VectorizedArrayType>::
7494  reinit(const typename Triangulation<dim>::cell_iterator &cell)
7495 {
7496  Assert(this->matrix_free == 0,
7497  ExcMessage("Cannot use initialization from cell iterator if "
7498  "initialized from MatrixFree object. Use variant for "
7499  "on the fly computation with arguments as for FEValues "
7500  "instead"));
7501  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7502  this->mapped_geometry->reinit(cell);
7503 
7504 # ifdef DEBUG
7505  this->is_reinitialized = true;
7506 # endif
7507 }
7508 
7509 
7510 
7511 template <int dim,
7512  int fe_degree,
7513  int n_q_points_1d,
7514  int n_components_,
7515  typename Number,
7516  typename VectorizedArrayType>
7517 inline void
7518 FEEvaluation<dim,
7519  fe_degree,
7520  n_q_points_1d,
7521  n_components_,
7522  Number,
7523  VectorizedArrayType>::evaluate(const bool evaluate_values,
7524  const bool evaluate_gradients,
7525  const bool evaluate_hessians)
7526 {
7527 # ifdef DEBUG
7528  Assert(this->dof_values_initialized == true,
7530 # endif
7531  evaluate(this->values_dofs,
7532  evaluate_values,
7533  evaluate_gradients,
7534  evaluate_hessians);
7535 }
7536 
7537 
7538 template <int dim,
7539  int fe_degree,
7540  int n_q_points_1d,
7541  int n_components_,
7542  typename Number,
7543  typename VectorizedArrayType>
7544 inline void
7545 FEEvaluation<dim,
7546  fe_degree,
7547  n_q_points_1d,
7548  n_components_,
7549  Number,
7550  VectorizedArrayType>::
7551  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7552 {
7553 # ifdef DEBUG
7554  Assert(this->dof_values_initialized == true,
7556 # endif
7557  evaluate(this->values_dofs, evaluation_flags);
7558 }
7559 
7560 
7561 
7562 template <int dim,
7563  int fe_degree,
7564  int n_q_points_1d,
7565  int n_components_,
7566  typename Number,
7567  typename VectorizedArrayType>
7568 inline void
7569 FEEvaluation<dim,
7570  fe_degree,
7571  n_q_points_1d,
7572  n_components_,
7573  Number,
7574  VectorizedArrayType>::evaluate(const VectorizedArrayType
7575  * values_array,
7576  const bool evaluate_values,
7577  const bool evaluate_gradients,
7578  const bool evaluate_hessians)
7579 {
7581  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7582  ((evaluate_gradients) ? EvaluationFlags::gradients :
7584  ((evaluate_hessians) ? EvaluationFlags::hessians :
7586 
7587  evaluate(values_array, flag);
7588 }
7589 
7590 
7591 
7592 template <int dim,
7593  int fe_degree,
7594  int n_q_points_1d,
7595  int n_components_,
7596  typename Number,
7597  typename VectorizedArrayType>
7598 inline void
7599 FEEvaluation<dim,
7600  fe_degree,
7601  n_q_points_1d,
7602  n_components_,
7603  Number,
7604  VectorizedArrayType>::
7605  evaluate(const VectorizedArrayType * values_array,
7606  const EvaluationFlags::EvaluationFlags evaluation_flag)
7607 {
7608  const bool hessians_on_general_cells =
7609  evaluation_flag & EvaluationFlags::hessians &&
7610  (this->cell_type > internal::MatrixFreeFunctions::affine);
7611  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7612  if (hessians_on_general_cells)
7613  evaluation_flag_actual |= EvaluationFlags::gradients;
7614 
7615  if (fe_degree > -1)
7616  {
7618  evaluate(n_components, evaluation_flag_actual, values_array, *this);
7619  }
7620  else
7621  {
7623  n_components,
7624  evaluation_flag_actual,
7625  const_cast<VectorizedArrayType *>(values_array),
7626  *this);
7627  }
7628 
7629 # ifdef DEBUG
7630  if (evaluation_flag_actual & EvaluationFlags::values)
7631  this->values_quad_initialized = true;
7632  if (evaluation_flag_actual & EvaluationFlags::gradients)
7633  this->gradients_quad_initialized = true;
7634  if (evaluation_flag_actual & EvaluationFlags::hessians)
7635  this->hessians_quad_initialized = true;
7636 # endif
7637 }
7638 
7639 
7640 
7641 template <int dim,
7642  int fe_degree,
7643  int n_q_points_1d,
7644  int n_components_,
7645  typename Number,
7646  typename VectorizedArrayType>
7647 template <typename VectorType>
7648 inline void
7649 FEEvaluation<
7650  dim,
7651  fe_degree,
7652  n_q_points_1d,
7653  n_components_,
7654  Number,
7655  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7656  const bool evaluate_values,
7657  const bool evaluate_gradients,
7658  const bool evaluate_hessians)
7659 {
7661  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7662  ((evaluate_gradients) ? EvaluationFlags::gradients :
7664  ((evaluate_hessians) ? EvaluationFlags::hessians :
7666 
7667  gather_evaluate(input_vector, flag);
7668 }
7669 
7670 
7671 namespace internal
7672 {
7676  template <typename Number,
7677  typename VectorizedArrayType,
7678  typename VectorType,
7679  typename EvaluatorType,
7680  typename std::enable_if<internal::has_begin<VectorType> &&
7682  VectorType>::type * = nullptr>
7683  VectorizedArrayType *
7684  check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
7685  {
7686  // for user-defined cell batches this functionality is not supported
7687  if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
7688  return nullptr;
7689 
7690  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7691  const auto & dof_info = fe_eval.get_dof_info();
7692 
7693  // If the index storage is interleaved and contiguous and the vector
7694  // storage has the correct alignment, we can directly pass the pointer
7695  // into the vector to the evaluate() and integrate() calls, without
7696  // reading the vector entries into a separate data field. This saves some
7697  // operations.
7698  if (std::is_same<typename VectorType::value_type, Number>::value &&
7699  dof_info.index_storage_variants
7702  interleaved_contiguous &&
7703  reinterpret_cast<std::size_t>(
7704  vector.begin() +
7705  dof_info.dof_indices_contiguous
7707  [cell * VectorizedArrayType::size()]) %
7708  sizeof(VectorizedArrayType) ==
7709  0)
7710  {
7711  return reinterpret_cast<VectorizedArrayType *>(
7712  vector.begin() +
7713  dof_info.dof_indices_contiguous
7715  [cell * VectorizedArrayType::size()] +
7717  [fe_eval.get_active_fe_index()]
7718  [fe_eval.get_first_selected_component()] *
7719  VectorizedArrayType::size());
7720  }
7721  else
7722  return nullptr;
7723  }
7724 
7728  template <typename Number,
7729  typename VectorizedArrayType,
7730  typename VectorType,
7731  typename EvaluatorType,
7732  typename std::enable_if<!internal::has_begin<VectorType> ||
7734  VectorType>::type * = nullptr>
7735  VectorizedArrayType *
7736  check_vector_access_inplace(const EvaluatorType &, VectorType &)
7737  {
7738  return nullptr;
7739  }
7740 } // namespace internal
7741 
7742 
7743 
7744 template <int dim,
7745  int fe_degree,
7746  int n_q_points_1d,
7747  int n_components_,
7748  typename Number,
7749  typename VectorizedArrayType>
7750 template <typename VectorType>
7751 inline void
7752 FEEvaluation<dim,
7753  fe_degree,
7754  n_q_points_1d,
7755  n_components_,
7756  Number,
7757  VectorizedArrayType>::
7758  gather_evaluate(const VectorType & input_vector,
7759  const EvaluationFlags::EvaluationFlags evaluation_flag)
7760 {
7761  const VectorizedArrayType *src_ptr =
7762  internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7763  *this, input_vector);
7764  if (src_ptr != nullptr)
7765  evaluate(src_ptr, evaluation_flag);
7766  else
7767  {
7768  this->read_dof_values(input_vector);
7769  evaluate(this->begin_dof_values(), evaluation_flag);
7770  }
7771 }
7772 
7773 
7774 
7775 template <int dim,
7776  int fe_degree,
7777  int n_q_points_1d,
7778  int n_components_,
7779  typename Number,
7780  typename VectorizedArrayType>
7781 inline void
7782 FEEvaluation<dim,
7783  fe_degree,
7784  n_q_points_1d,
7785  n_components_,
7786  Number,
7787  VectorizedArrayType>::integrate(const bool integrate_values,
7788  const bool integrate_gradients)
7789 {
7790  integrate(integrate_values, integrate_gradients, this->values_dofs);
7791 
7792 # ifdef DEBUG
7793  this->dof_values_initialized = true;
7794 # endif
7795 }
7796 
7797 
7798 
7799 template <int dim,
7800  int fe_degree,
7801  int n_q_points_1d,
7802  int n_components_,
7803  typename Number,
7804  typename VectorizedArrayType>
7805 inline void
7806 FEEvaluation<dim,
7807  fe_degree,
7808  n_q_points_1d,
7809  n_components_,
7810  Number,
7811  VectorizedArrayType>::
7812  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7813 {
7814  integrate(integration_flag, this->values_dofs);
7815 
7816 # ifdef DEBUG
7817  this->dof_values_initialized = true;
7818 # endif
7819 }
7820 
7821 
7822 
7823 template <int dim,
7824  int fe_degree,
7825  int n_q_points_1d,
7826  int n_components_,
7827  typename Number,
7828  typename VectorizedArrayType>
7829 inline void
7830 FEEvaluation<dim,
7831  fe_degree,
7832  n_q_points_1d,
7833  n_components_,
7834  Number,
7835  VectorizedArrayType>::integrate(const bool integrate_values,
7836  const bool integrate_gradients,
7837  VectorizedArrayType *values_array)
7838 {
7840  (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
7841  (integrate_gradients ? EvaluationFlags::gradients :
7843  integrate(flag, values_array);
7844 }
7845 
7846 
7847 
7848 template <int dim,
7849  int fe_degree,
7850  int n_q_points_1d,
7851  int n_components_,
7852  typename Number,
7853  typename VectorizedArrayType>
7854 inline void
7855 FEEvaluation<dim,
7856  fe_degree,
7857  n_q_points_1d,
7858  n_components_,
7859  Number,
7860  VectorizedArrayType>::
7861  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7862  VectorizedArrayType * values_array,
7863  const bool sum_into_values_array)
7864 {
7865 # ifdef DEBUG
7866  if (integration_flag & EvaluationFlags::values)
7867  Assert(this->values_quad_submitted == true,
7869  if (integration_flag & EvaluationFlags::gradients)
7870  Assert(this->gradients_quad_submitted == true,
7872  if ((integration_flag & EvaluationFlags::hessians) != 0u)
7873  Assert(this->hessians_quad_submitted == true,
7875 # endif
7876  Assert(this->matrix_free != nullptr ||
7877  this->mapped_geometry->is_initialized(),
7878  ExcNotInitialized());
7879 
7880  Assert(
7881  (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
7883  ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7884  "EvaluationFlags::hessians are supported."));
7885 
7886  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7887  if (integration_flag & EvaluationFlags::hessians &&
7888  (this->cell_type > internal::MatrixFreeFunctions::affine))
7889  {
7890  unsigned int size = n_components * dim * n_q_points;
7891  if ((integration_flag & EvaluationFlags::gradients) != 0u)
7892  {
7893  for (unsigned int i = 0; i < size; ++i)
7894  this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7895  }
7896  else
7897  {
7898  for (unsigned int i = 0; i < size; ++i)
7899  this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7900  integration_flag_actual |= EvaluationFlags::gradients;
7901  }
7902  }
7903 
7904  if (fe_degree > -1)
7905  {
7907  integrate(n_components,
7908  integration_flag_actual,
7909  values_array,
7910  *this,
7911  sum_into_values_array);
7912  }
7913  else
7914  {
7916  n_components,
7917  integration_flag_actual,
7918  values_array,
7919  *this,
7920  sum_into_values_array);
7921  }
7922 
7923 # ifdef DEBUG
7924  this->dof_values_initialized = true;
7925 # endif
7926 }
7927 
7928 
7929 
7930 template <int dim,
7931  int fe_degree,
7932  int n_q_points_1d,
7933  int n_components_,
7934  typename Number,
7935  typename VectorizedArrayType>
7936 template <typename VectorType>
7937 inline void
7938 FEEvaluation<
7939  dim,
7940  fe_degree,
7941  n_q_points_1d,
7942  n_components_,
7943  Number,
7944  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7945  const bool integrate_gradients,
7946  VectorType &destination)
7947 {
7949  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7950  ((integrate_gradients) ? EvaluationFlags::gradients :
7952 
7953  integrate_scatter(flag, destination);
7954 }
7955 
7956 
7957 
7958 template <int dim,
7959  int fe_degree,
7960  int n_q_points_1d,
7961  int n_components_,
7962  typename Number,
7963  typename VectorizedArrayType>
7964 template <typename VectorType>
7965 inline void
7966 FEEvaluation<dim,
7967  fe_degree,
7968  n_q_points_1d,
7969  n_components_,
7970  Number,
7971  VectorizedArrayType>::
7972  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
7973  VectorType & destination)
7974 {
7975  VectorizedArrayType *dst_ptr =
7976  internal::check_vector_access_inplace<Number, VectorizedArrayType>(
7977  *this, destination);
7978  if (dst_ptr != nullptr)
7979  integrate(integration_flag, dst_ptr, true);
7980  else
7981  {
7982  integrate(integration_flag, this->begin_dof_values());
7983  this->distribute_local_to_global(destination);
7984  }
7985 }
7986 
7987 
7988 
7989 template <int dim,
7990  int fe_degree,
7991  int n_q_points_1d,
7992  int n_components_,
7993  typename Number,
7994  typename VectorizedArrayType>
7996 FEEvaluation<dim,
7997  fe_degree,
7998  n_q_points_1d,
7999  n_components_,
8000  Number,
8001  VectorizedArrayType>::dof_indices() const
8002 {
8003  return {0U, dofs_per_cell};
8004 }
8005 
8006 
8007 
8008 /*-------------------------- FEFaceEvaluation ---------------------------*/
8009 
8010 
8011 
8012 template <int dim,
8013  int fe_degree,
8014  int n_q_points_1d,
8015  int n_components_,
8016  typename Number,
8017  typename VectorizedArrayType>
8018 inline FEFaceEvaluation<dim,
8019  fe_degree,
8020  n_q_points_1d,
8021  n_components_,
8022  Number,
8023  VectorizedArrayType>::
8024  FEFaceEvaluation(
8026  const bool is_interior_face,
8027  const unsigned int dof_no,
8028  const unsigned int quad_no,
8029  const unsigned int first_selected_component,
8030  const unsigned int active_fe_index,
8031  const unsigned int active_quad_index,
8032  const unsigned int face_type)
8033  : BaseClass(matrix_free,
8034  dof_no,
8035  first_selected_component,
8036  quad_no,
8037  fe_degree,
8038  static_n_q_points,
8039  is_interior_face,
8040  active_fe_index,
8041  active_quad_index,
8042  face_type)
8043  , dofs_per_component(this->data->dofs_per_component_on_cell)
8044  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
8045  , n_q_points(this->n_quadrature_points)
8046 {}
8047 
8048 
8049 
8050 template <int dim,
8051  int fe_degree,
8052  int n_q_points_1d,
8053  int n_components_,
8054  typename Number,
8055  typename VectorizedArrayType>
8056 inline FEFaceEvaluation<dim,
8057  fe_degree,
8058  n_q_points_1d,
8059  n_components_,
8060  Number,
8061  VectorizedArrayType>::
8062  FEFaceEvaluation(
8064  const std::pair<unsigned int, unsigned int> & range,
8065  const bool is_interior_face,
8066  const unsigned int dof_no,
8067  const unsigned int quad_no,
8068  const unsigned int first_selected_component)
8069  : FEFaceEvaluation(matrix_free,
8070  is_interior_face,
8071  dof_no,
8072  quad_no,
8073  first_selected_component,
8074  matrix_free.get_face_active_fe_index(range,
8075  is_interior_face),
8077  matrix_free.get_face_info(range.first).face_type)
8078 {}
8079 
8080 
8081 
8082 template <int dim,
8083  int fe_degree,
8084  int n_q_points_1d,
8085  int n_components_,
8086  typename Number,
8087  typename VectorizedArrayType>
8088 inline void
8089 FEFaceEvaluation<dim,
8090  fe_degree,
8091  n_q_points_1d,
8092  n_components_,
8093  Number,
8094  VectorizedArrayType>::reinit(const unsigned int face_index)
8095 {
8096  Assert(this->mapped_geometry == nullptr,
8097  ExcMessage("FEEvaluation was initialized without a matrix-free object."
8098  " Integer indexing is not possible"));
8099  if (this->mapped_geometry != nullptr)
8100  return;
8101 
8102  this->cell = face_index;
8103  this->dof_access_index =
8104  this->is_interior_face() ?
8107  Assert(this->mapping_data != nullptr, ExcNotInitialized());
8108 
8109  if (face_index >=
8110  this->matrix_free->get_task_info().face_partition_data.back() &&
8111  face_index <
8112  this->matrix_free->get_task_info().boundary_partition_data.back())
8113  Assert(this->is_interior_face(),
8114  ExcMessage(
8115  "Boundary faces do not have a neighbor. When looping over "
8116  "boundary faces use FEFaceEvaluation with the parameter "
8117  "is_interior_face set to true. "));
8118 
8119  this->reinit_face(this->matrix_free->get_face_info(face_index));
8120 
8121  unsigned int i = 0;
8122  for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
8123  ++i)
8124  this->face_ids[i] = face_index * VectorizedArrayType::size() + i;
8125  for (; i < VectorizedArrayType::size(); ++i)
8126  this->face_ids[i] = numbers::invalid_unsigned_int;
8127 
8128  this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
8129  const unsigned int offsets =
8130  this->mapping_data->data_index_offsets[face_index];
8131  this->J_value = &this->mapping_data->JxW_values[offsets];
8132  this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8133  this->jacobian =
8134  &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
8135  this->normal_x_jacobian =
8136  &this->mapping_data
8137  ->normals_times_jacobians[!this->is_interior_face()][offsets];
8138  this->jacobian_gradients =
8139  this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8140  offsets;
8141 
8142  if (this->mapping_data->quadrature_point_offsets.empty() == false)
8143  {
8144  AssertIndexRange(this->cell,
8145  this->mapping_data->quadrature_point_offsets.size());
8146  this->quadrature_points =
8147  this->mapping_data->quadrature_points.data() +
8148  this->mapping_data->quadrature_point_offsets[this->cell];
8149  }
8150 
8151 # ifdef DEBUG
8152  this->is_reinitialized = true;
8153  this->dof_values_initialized = false;
8154  this->values_quad_initialized = false;
8155  this->gradients_quad_initialized = false;
8156  this->hessians_quad_initialized = false;
8157 # endif
8158 }
8159 
8160 
8161 
8162 template <int dim,
8163  int fe_degree,
8164  int n_q_points_1d,
8165  int n_components_,
8166  typename Number,
8167  typename VectorizedArrayType>
8168 inline void
8169 FEFaceEvaluation<dim,
8170  fe_degree,
8171  n_q_points_1d,
8172  n_components_,
8173  Number,
8174  VectorizedArrayType>::reinit(const unsigned int cell_index,
8175  const unsigned int face_number)
8176 {
8177  Assert(
8178  this->quad_no <
8179  this->matrix_free->get_mapping_info().face_data_by_cells.size(),
8180  ExcMessage(
8181  "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8184  this->matrix_free->get_mapping_info().cell_type.size());
8185  Assert(this->mapped_geometry == nullptr,
8186  ExcMessage("FEEvaluation was initialized without a matrix-free object."
8187  " Integer indexing is not possible"));
8188  if (this->mapped_geometry != nullptr)
8189  return;
8190  Assert(this->matrix_free != nullptr, ExcNotInitialized());
8191 
8192  this->cell_type = this->matrix_free->get_mapping_info()
8193  .faces_by_cells_type[cell_index][face_number];
8194  this->cell = cell_index;
8195  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
8196  this->dof_access_index =
8198 
8199  constexpr unsigned int n_lanes = VectorizedArrayType::size();
8200 
8201  if (this->is_interior_face() == false)
8202  {
8203  // for this case, we need to look into the FaceInfo field that collects
8204  // information from both sides of a face once for the global mesh, and
8205  // pick the face id that is not the local one (cell_this).
8206  for (unsigned int i = 0; i < n_lanes; ++i)
8207  {
8208  // compute actual (non vectorized) cell ID
8209  const unsigned int cell_this = cell_index * n_lanes + i;
8210  // compute face ID
8211  unsigned int face_index =
8212  this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
8213  face_number,
8214  i);
8215 
8216  this->face_ids[i] = face_index;
8217 
8218  if (face_index == numbers::invalid_unsigned_int)
8219  {
8220  this->cell_ids[i] = numbers::invalid_unsigned_int;
8221  this->face_numbers[i] = static_cast<std::uint8_t>(-1);
8222  this->face_orientations[i] = static_cast<std::uint8_t>(-1);
8223  continue; // invalid face ID: no neighbor on boundary
8224  }
8225 
8226  const auto &faces =
8227  this->matrix_free->get_face_info(face_index / n_lanes);
8228  // get cell ID on both sides of face
8229  auto cell_m = faces.cells_interior[face_index % n_lanes];
8230  auto cell_p = faces.cells_exterior[face_index % n_lanes];
8231 
8232  const bool face_identifies_as_interior = cell_m != cell_this;
8233 
8234  Assert(cell_m == cell_this || cell_p == cell_this,
8235  ExcInternalError());
8236 
8237  // compare the IDs with the given cell ID
8238  if (face_identifies_as_interior)
8239  {
8240  this->cell_ids[i] = cell_m; // neighbor has the other ID
8241  this->face_numbers[i] = faces.interior_face_no;
8242  }
8243  else
8244  {
8245  this->cell_ids[i] = cell_p;
8246  this->face_numbers[i] = faces.exterior_face_no;
8247  }
8248 
8249  const bool orientation_interior_face = faces.face_orientation >= 8;
8250  unsigned int face_orientation = faces.face_orientation % 8;
8251  if (face_identifies_as_interior != orientation_interior_face)
8252  {
8253  constexpr std::array<std::uint8_t, 8> table{
8254  {0, 1, 2, 3, 6, 5, 4, 7}};
8255  face_orientation = table[face_orientation];
8256  }
8257  this->face_orientations[i] = face_orientation;
8258  }
8259  }
8260  else
8261  {
8262  this->face_orientations[0] = 0;
8263  this->face_numbers[0] = face_number;
8264  for (unsigned int i = 0; i < n_lanes; ++i)
8265  this->cell_ids[i] = cell_index * n_lanes + i;
8266  for (unsigned int i = 0; i < n_lanes; ++i)
8267  this->face_ids[i] =
8268  this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
8269  face_number,
8270  i);
8271  }
8272 
8273  const unsigned int offsets =
8274  this->matrix_free->get_mapping_info()
8275  .face_data_by_cells[this->quad_no]
8276  .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
8277  face_number];
8278  AssertIndexRange(offsets,
8279  this->matrix_free->get_mapping_info()
8280  .face_data_by_cells[this->quad_no]
8281  .JxW_values.size());
8282  this->J_value = &this->matrix_free->get_mapping_info()
8283  .face_data_by_cells[this->quad_no]
8284  .JxW_values[offsets];
8285  this->normal_vectors = &this->matrix_free->get_mapping_info()
8286  .face_data_by_cells[this->quad_no]
8287  .normal_vectors[offsets];
8288  this->jacobian = &this->matrix_free->get_mapping_info()
8289  .face_data_by_cells[this->quad_no]
8290  .jacobians[!this->is_interior_face()][offsets];
8291  this->normal_x_jacobian =
8292  &this->matrix_free->get_mapping_info()
8293  .face_data_by_cells[this->quad_no]
8294  .normals_times_jacobians[!this->is_interior_face()][offsets];
8295  this->jacobian_gradients =
8296  this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8297  offsets;
8298 
8299  if (this->matrix_free->get_mapping_info()
8300  .face_data_by_cells[this->quad_no]
8301  .quadrature_point_offsets.empty() == false)
8302  {
8303  const unsigned int index =
8304  this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
8305  AssertIndexRange(index,
8306  this->matrix_free->get_mapping_info()
8307  .face_data_by_cells[this->quad_no]
8308  .quadrature_point_offsets.size());
8309  this->quadrature_points = this->matrix_free->get_mapping_info()
8310  .face_data_by_cells[this->quad_no]
8311  .quadrature_points.data() +
8312  this->matrix_free->get_mapping_info()
8313  .face_data_by_cells[this->quad_no]
8314  .quadrature_point_offsets[index];
8315  }
8316 
8317 # ifdef DEBUG
8318  this->is_reinitialized = true;
8319  this->dof_values_initialized = false;
8320  this->values_quad_initialized = false;
8321  this->gradients_quad_initialized = false;
8322  this->hessians_quad_initialized = false;
8323 # endif
8324 }
8325 
8326 
8327 
8328 template <int dim,
8329  int fe_degree,
8330  int n_q_points_1d,
8331  int n_components,
8332  typename Number,
8333  typename VectorizedArrayType>
8334 inline void
8335 FEFaceEvaluation<dim,
8336  fe_degree,
8337  n_q_points_1d,
8338  n_components,
8339  Number,
8340  VectorizedArrayType>::evaluate(const bool evaluate_values,
8341  const bool evaluate_gradients)
8342 {
8343 # ifdef DEBUG
8344  Assert(this->dof_values_initialized, ExcNotInitialized());
8345 # endif
8346 
8347  evaluate(this->values_dofs, evaluate_values, evaluate_gradients);
8348 }
8349 
8350 
8351 
8352 template <int dim,
8353  int fe_degree,
8354  int n_q_points_1d,
8355  int n_components,
8356  typename Number,
8357  typename VectorizedArrayType>
8358 inline void
8359 FEFaceEvaluation<dim,
8360  fe_degree,
8361  n_q_points_1d,
8362  n_components,
8363  Number,
8364  VectorizedArrayType>::
8365  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
8366 {
8367 # ifdef DEBUG
8368  Assert(this->dof_values_initialized, ExcNotInitialized());
8369 # endif
8370 
8371  evaluate(this->values_dofs, evaluation_flag);
8372 }
8373 
8374 
8375 
8376 template <int dim,
8377  int fe_degree,
8378  int n_q_points_1d,
8379  int n_components,
8380  typename Number,
8381  typename VectorizedArrayType>
8382 inline void
8383 FEFaceEvaluation<dim,
8384  fe_degree,
8385  n_q_points_1d,
8386  n_components,
8387  Number,
8388  VectorizedArrayType>::evaluate(const VectorizedArrayType
8389  * values_array,
8390  const bool evaluate_values,
8391  const bool evaluate_gradients)
8392 {
8394  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8395  ((evaluate_gradients) ? EvaluationFlags::gradients :
8397 
8398  evaluate(values_array, flag);
8399 }
8400 
8401 
8402 
8403 template <int dim,
8404  int fe_degree,
8405  int n_q_points_1d,
8406  int n_components,
8407  typename Number,
8408  typename VectorizedArrayType>
8409 inline void
8410 FEFaceEvaluation<dim,
8411  fe_degree,
8412  n_q_points_1d,
8413  n_components,
8414  Number,
8415  VectorizedArrayType>::
8416  evaluate(const VectorizedArrayType * values_array,
8417  const EvaluationFlags::EvaluationFlags evaluation_flag)
8418 {
8419  Assert((evaluation_flag &
8422  ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8423  "and EvaluationFlags::hessians are supported."));
8424 
8425  const bool hessians_on_general_cells =
8426  evaluation_flag & EvaluationFlags::hessians &&
8427  (this->cell_type > internal::MatrixFreeFunctions::affine);
8428  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
8429  if (hessians_on_general_cells)
8430  evaluation_flag_actual |= EvaluationFlags::gradients;
8431 
8432  if (fe_degree > -1)
8434  template run<fe_degree, n_q_points_1d>(n_components,
8435  evaluation_flag_actual,
8436  values_array,
8437  *this);
8438  else
8440  n_components, evaluation_flag_actual, values_array, *this);
8441 
8442 # ifdef DEBUG
8443  if (evaluation_flag_actual & EvaluationFlags::values)
8444  this->values_quad_initialized = true;
8445  if (evaluation_flag_actual & EvaluationFlags::gradients)
8446  this->gradients_quad_initialized = true;
8447  if ((evaluation_flag_actual & EvaluationFlags::hessians) != 0u)
8448  this->hessians_quad_initialized = true;
8449 # endif
8450 }
8451 
8452 
8453 
8454 template <int dim,
8455  int fe_degree,
8456  int n_q_points_1d,
8457  int n_components,
8458  typename Number,
8459  typename VectorizedArrayType>
8460 inline void
8461 FEFaceEvaluation<dim,
8462  fe_degree,
8463  n_q_points_1d,
8464  n_components,
8465  Number,
8466  VectorizedArrayType>::
8467  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
8468 {
8469  integrate(integration_flag, this->values_dofs);
8470 
8471 # ifdef DEBUG
8472  this->dof_values_initialized = true;
8473 # endif
8474 }
8475 
8476 
8477 
8478 template <int dim,
8479  int fe_degree,
8480  int n_q_points_1d,
8481  int n_components,
8482  typename Number,
8483  typename VectorizedArrayType>
8484 inline void
8485 FEFaceEvaluation<dim,
8486  fe_degree,
8487  n_q_points_1d,
8488  n_components,
8489  Number,
8490  VectorizedArrayType>::integrate(const bool integrate_values,
8491  const bool integrate_gradients)
8492 {
8493  integrate(integrate_values, integrate_gradients, this->values_dofs);
8494 
8495 # ifdef DEBUG
8496  this->dof_values_initialized = true;
8497 # endif
8498 }
8499 
8500 
8501 
8502 template <int dim,
8503  int fe_degree,
8504  int n_q_points_1d,
8505  int n_components,
8506  typename Number,
8507  typename VectorizedArrayType>
8508 inline void
8509 FEFaceEvaluation<dim,
8510  fe_degree,
8511  n_q_points_1d,
8512  n_components,
8513  Number,
8514  VectorizedArrayType>::integrate(const bool integrate_values,
8515  const bool integrate_gradients,
8516  VectorizedArrayType
8517  *values_array)
8518 {
8520  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8521  ((integrate_gradients) ? EvaluationFlags::gradients :
8523 
8524  integrate(flag, values_array);
8525 }
8526 
8527 
8528 
8529 template <int dim,
8530  int fe_degree,
8531  int n_q_points_1d,
8532  int n_components,
8533  typename Number,
8534  typename VectorizedArrayType>
8535 inline void
8536 FEFaceEvaluation<dim,
8537  fe_degree,
8538  n_q_points_1d,
8539  n_components,
8540  Number,
8541  VectorizedArrayType>::
8542  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
8543  VectorizedArrayType * values_array)
8544 {
8545  Assert((integration_flag &
8548  ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8549  "and EvaluationFlags::hessians are supported."));
8550 
8551  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
8552  if (integration_flag & EvaluationFlags::hessians &&
8553  (this->cell_type > internal::MatrixFreeFunctions::affine))
8554  {
8555  unsigned int size = n_components * dim * n_q_points;
8556  if ((integration_flag & EvaluationFlags::gradients) != 0u)
8557  {
8558  for (unsigned int i = 0; i < size; ++i)
8559  this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
8560  }
8561  else
8562  {
8563  for (unsigned int i = 0; i < size; ++i)
8564  this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
8565  integration_flag_actual |= EvaluationFlags::gradients;
8566  }
8567  }
8568 
8569  if (fe_degree > -1)
8571  template run<fe_degree, n_q_points_1d>(n_components,
8572  integration_flag_actual,
8573  values_array,
8574  *this);
8575  else
8577  n_components, integration_flag_actual, values_array, *this);
8578 }
8579 
8580 
8581 
8582 template <int dim,
8583  int fe_degree,
8584  int n_q_points_1d,
8585  int n_components_,
8586  typename Number,
8587  typename VectorizedArrayType>
8588 template <typename VectorType>
8589 inline void
8591  dim,
8592  fe_degree,
8593  n_q_points_1d,
8594  n_components_,
8595  Number,
8596  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
8597  const bool evaluate_values,
8598  const bool evaluate_gradients)
8599 {
8601  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8602  ((evaluate_gradients) ? EvaluationFlags::gradients :
8604 
8605  gather_evaluate(input_vector, flag);
8606 }
8607 
8608 
8609 
8610 template <int dim,
8611  int fe_degree,
8612  int n_q_points_1d,
8613  int n_components_,
8614  typename Number,
8615  typename VectorizedArrayType>
8616 template <typename VectorType>
8617 inline void
8618 FEFaceEvaluation<dim,
8619  fe_degree,
8620  n_q_points_1d,
8621  n_components_,
8622  Number,
8623  VectorizedArrayType>::
8624  gather_evaluate(const VectorType & input_vector,
8625  const EvaluationFlags::EvaluationFlags evaluation_flag)
8626 {
8627  Assert((evaluation_flag &
8630  ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8631  "and EvaluationFlags::hessians are supported."));
8632 
8633  const auto shared_vector_data = internal::get_shared_vector_data(
8634  input_vector,
8635  this->dof_access_index ==
8637  this->active_fe_index,
8638  this->dof_info);
8639 
8640  if (this->data->data.front().fe_degree > 0 &&
8641  fast_evaluation_supported(this->data->data.front().fe_degree,
8642  this->data->data.front().n_q_points_1d) &&
8644  dim,
8645  typename VectorType::value_type,
8646  VectorizedArrayType>::
8647  supports(evaluation_flag,
8648  *this->data,
8649  internal::get_beginning<typename VectorType::value_type>(
8650  input_vector),
8651  this->dof_info->index_storage_variants[this->dof_access_index]
8652  [this->cell]))
8653  {
8654  if (fe_degree > -1)
8655  {
8657  dim,
8658  typename VectorType::value_type,
8659  VectorizedArrayType>::template run<fe_degree,
8660  n_q_points_1d>(
8661  n_components,
8662  evaluation_flag,
8663  internal::get_beginning<typename VectorType::value_type>(
8664  input_vector),
8665  shared_vector_data,
8666  *this);
8667  }
8668  else
8669  {
8671  dim,
8672  typename VectorType::value_type,
8673  VectorizedArrayType>::evaluate(n_components,
8674  evaluation_flag,
8675  internal::get_beginning<
8676  typename VectorType::value_type>(
8677  input_vector),
8678  shared_vector_data,
8679  *this);
8680  }
8681  }
8682  else
8683  {
8684  this->read_dof_values(input_vector);
8685  this->evaluate(evaluation_flag);
8686  }
8687 
8688 # ifdef DEBUG
8689  if (evaluation_flag & EvaluationFlags::values)
8690  this->values_quad_initialized = true;
8691  if (evaluation_flag & EvaluationFlags::gradients)
8692  this->gradients_quad_initialized = true;
8693  if (evaluation_flag & EvaluationFlags::hessians)
8694  this->hessians_quad_initialized = true;
8695 # endif
8696 }
8697 
8698 
8699 
8700 template <int dim,
8701  int fe_degree,
8702  int n_q_points_1d,
8703  int n_components_,
8704  typename Number,
8705  typename VectorizedArrayType>
8706 template <typename VectorType>
8707 inline void
8709  dim,
8710  fe_degree,
8711  n_q_points_1d,
8712  n_components_,
8713  Number,
8714  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
8715  const bool integrate_gradients,
8716  VectorType &destination)
8717 {
8719  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8720  ((integrate_gradients) ? EvaluationFlags::gradients :
8722 
8723  integrate_scatter(flag, destination);
8724 }
8725 
8726 
8727 
8728 template <int dim,
8729  int fe_degree,
8730  int n_q_points_1d,
8731  int n_components_,
8732  typename Number,
8733  typename VectorizedArrayType>
8734 template <typename VectorType>
8735 inline void
8736 FEFaceEvaluation<dim,
8737  fe_degree,
8738  n_q_points_1d,
8739  n_components_,
8740  Number,
8741  VectorizedArrayType>::
8742  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
8743  VectorType & destination)
8744 {
8745  Assert((this->dof_access_index ==
8747  this->is_interior_face() == false) == false,
8748  ExcNotImplemented());
8749 
8750  const auto shared_vector_data = internal::get_shared_vector_data(
8751  destination,
8752  this->dof_access_index ==
8754  this->active_fe_index,
8755  this->dof_info);
8756 
8757  if (this->data->data.front().fe_degree > 0 &&
8758  fast_evaluation_supported(this->data->data.front().fe_degree,
8759  this->data->data.front().n_q_points_1d) &&
8761  dim,
8762  typename VectorType::value_type,
8763  VectorizedArrayType>::
8764  supports(integration_flag,
8765  *this->data,
8766  internal::get_beginning<typename VectorType::value_type>(
8767  destination),
8768  this->dof_info->index_storage_variants[this->dof_access_index]
8769  [this->cell]))
8770  {
8771  if (fe_degree > -1)
8772  {
8774  dim,
8775  typename VectorType::value_type,
8776  VectorizedArrayType>::template run<fe_degree,
8777  n_q_points_1d>(
8778  n_components,
8779  integration_flag,
8780  internal::get_beginning<typename VectorType::value_type>(
8781  destination),
8782  shared_vector_data,
8783  *this);
8784  }
8785  else
8786  {
8788  dim,
8789  typename VectorType::value_type,
8790  VectorizedArrayType>::integrate(n_components,
8791  integration_flag,
8792  internal::get_beginning<
8793  typename VectorType::value_type>(
8794  destination),
8795  shared_vector_data,
8796  *this);
8797  }
8798  }
8799  else
8800  {
8801  integrate(integration_flag);
8802  this->distribute_local_to_global(destination);
8803  }
8804 }
8805 
8806 
8807 
8808 template <int dim,
8809  int fe_degree,
8810  int n_q_points_1d,
8811  int n_components_,
8812  typename Number,
8813  typename VectorizedArrayType>
8815 FEFaceEvaluation<dim,
8816  fe_degree,
8817  n_q_points_1d,
8818  n_components_,
8819  Number,
8820  VectorizedArrayType>::dof_indices() const
8821 {
8822  return {0U, dofs_per_cell};
8823 }
8824 
8825 
8826 
8827 template <int dim,
8828  int fe_degree,
8829  int n_q_points_1d,
8830  int n_components_,
8831  typename Number,
8832  typename VectorizedArrayType>
8833 bool
8834 FEEvaluation<dim,
8835  fe_degree,
8836  n_q_points_1d,
8837  n_components_,
8838  Number,
8839  VectorizedArrayType>::
8840  fast_evaluation_supported(const unsigned int given_degree,
8841  const unsigned int give_n_q_points_1d)
8842 {
8843  return fe_degree == -1 ?
8845  fast_evaluation_supported(given_degree, give_n_q_points_1d) :
8846  true;
8847 }
8848 
8849 
8850 
8851 template <int dim,
8852  int fe_degree,
8853  int n_q_points_1d,
8854  int n_components_,
8855  typename Number,
8856  typename VectorizedArrayType>
8857 bool
8858 FEFaceEvaluation<dim,
8859  fe_degree,
8860  n_q_points_1d,
8861  n_components_,
8862  Number,
8863  VectorizedArrayType>::
8864  fast_evaluation_supported(const unsigned int given_degree,
8865  const unsigned int give_n_q_points_1d)
8866 {
8867  return fe_degree == -1 ?
8869  fast_evaluation_supported(given_degree, give_n_q_points_1d) :
8870  true;
8871 }
8872 
8873 
8874 
8875 /*------------------------- end FEFaceEvaluation ------------------------- */
8876 
8877 
8878 #endif // ifndef DOXYGEN
8879 
8880 
8882 
8883 #endif
hessian_type get_hessian(unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
value_type get_laplacian(const unsigned int q_point) const
void submit_value(const gradient_type val_in, const unsigned int q_point)
value_type get_normal_derivative(const unsigned int q_point) const
void submit_normal_derivative(const gradient_type grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
void submit_gradient(const value_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< 1 > &mapping, const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< 1, VectorizedArrayType, is_face > *other)
value_type get_divergence(const unsigned int q_point) const
value_type get_dof_value(const unsigned int dof) const
FEEvaluationAccess(const MatrixFree< 1, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
value_type get_dof_value(const unsigned int dof) const
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
value_type get_laplacian(const unsigned int q_point) const
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
value_type get_normal_derivative(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 1, dim, Tensor< 1, dim, VectorizedArrayType >> grad_in, const unsigned int q_point)
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
Tensor< 3, dim, VectorizedArrayType > get_hessian(const unsigned int q_point) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
VectorizedArrayType get_divergence(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_value(const Tensor< 1, dim, VectorizedArrayType > val_in, const unsigned int q_point)
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
static constexpr unsigned int dimension
Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > gradient_type
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const FEEvaluationAccess &other)
static constexpr unsigned int n_components
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
value_type get_dof_value(const unsigned int dof) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
static constexpr unsigned int dimension
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type >> *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask, const bool apply_constraints=true) const
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type >> *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask) const
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
VectorizedArrayType get_divergence(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
static constexpr unsigned int n_components
gradient_type get_hessian_diagonal(const unsigned int q_point) const
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > get_hessian(const unsigned int q_point) const
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
FEEvaluationData & operator=(const FEEvaluationData &other)
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const DoFInfo * dof_info
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
void integrate(const bool integrate_values, const bool integrate_gradients)
void reinit(const std::array< unsigned int, VectorizedArrayType::size()> &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
FEEvaluation & operator=(const FEEvaluation &other)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access >> &cell)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
FEEvaluation(const FEEvaluation &other)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
static constexpr unsigned int tensor_dofs_per_cell
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array)
static constexpr unsigned int static_dofs_per_component
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int give_n_q_points_1d)
static constexpr unsigned int static_n_q_points
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate(const bool integrate_values, const bool integrate_gradients)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
Definition: mapping.h:311
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const Number * constraint_pool_begin(const unsigned int pool_index) const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
const Number * constraint_pool_end(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
bool indices_initialized() const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
unsigned int n_components() const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
Definition: tensor.h:503
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
UpdateFlags
Point< 2 > second
Definition: grid_out.cc:4604
Point< 2 > first
Definition: grid_out.cc:4603
unsigned int cell_index
Definition: grid_tools.cc:1129
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
std::string to_string(const T &t)
Definition: patterns.h:2403
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
EvaluationFlags
The EvaluationFlags enum.
static const char U
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:493
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm &mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
void run(const Iterator &begin, const typename identity< Iterator >::type &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition: work_stream.h:474
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
static const unsigned int invalid_unsigned_int
Definition: types.h:201
boost::integer_range< IncrementableType > iota_view
Definition: iota_view.h:46
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &eval, const bool sum_into_values_array=false)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition: dof_info.h:718
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
Definition: dof_info.h:554
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
Definition: dof_info.h:592
std::vector< unsigned int > component_to_base_index
Definition: dof_info.h:705
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
Definition: dof_info.h:581
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition: dof_info.h:517
std::vector< QuadratureDescriptor > descriptor
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
Definition: task_info.h:499
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)