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