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