Reference documentation for deal.II version Git b43ef918fe 2022-01-29 10:49:23 +0100
\(\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 - 2021 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 
43 
44 #include <type_traits>
45 
46 
48 
49 
50 
88 template <int dim,
89  int n_components_,
90  typename Number,
91  bool is_face,
92  typename VectorizedArrayType>
94  : public FEEvaluationData<dim, VectorizedArrayType, is_face>
95 {
96 public:
97  using number_type = Number;
99  using gradient_type =
101  using hessian_type =
103  static constexpr unsigned int dimension = dim;
104  static constexpr unsigned int n_components = n_components_;
105 
142  template <typename VectorType>
143  void
144  read_dof_values(const VectorType & src,
145  const unsigned int first_index = 0,
146  const std::bitset<VectorizedArrayType::size()> &mask =
147  std::bitset<VectorizedArrayType::size()>().flip());
148 
177  template <typename VectorType>
178  void
179  read_dof_values_plain(const VectorType & src,
180  const unsigned int first_index = 0,
181  const std::bitset<VectorizedArrayType::size()> &mask =
182  std::bitset<VectorizedArrayType::size()>().flip());
183 
215  template <typename VectorType>
216  void
218  VectorType & dst,
219  const unsigned int first_index = 0,
220  const std::bitset<VectorizedArrayType::size()> &mask =
221  std::bitset<VectorizedArrayType::size()>().flip()) const;
222 
261  template <typename VectorType>
262  void
264  const unsigned int first_index = 0,
265  const std::bitset<VectorizedArrayType::size()> &mask =
266  std::bitset<VectorizedArrayType::size()>().flip()) const;
267 
271  template <typename VectorType>
272  void
274  VectorType & dst,
275  const unsigned int first_index = 0,
276  const std::bitset<VectorizedArrayType::size()> &mask =
277  std::bitset<VectorizedArrayType::size()>().flip()) const;
278 
280 
301  value_type
302  get_dof_value(const unsigned int dof) const;
303 
314  void
315  submit_dof_value(const value_type val_in, const unsigned int dof);
316 
329  value_type
330  get_value(const unsigned int q_point) const;
331 
344  void
345  submit_value(const value_type val_in, const unsigned int q_point);
346 
358  get_gradient(const unsigned int q_point) const;
359 
374  value_type
375  get_normal_derivative(const unsigned int q_point) const;
376 
389  void
390  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
391 
410  void
411  submit_normal_derivative(const value_type grad_in,
412  const unsigned int q_point);
413 
426  void
427  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
428 
441  get_hessian(const unsigned int q_point) const;
442 
453  get_hessian_diagonal(const unsigned int q_point) const;
454 
466  value_type
467  get_laplacian(const unsigned int q_point) const;
468 
469 #ifdef DOXYGEN
470  // doxygen does not anyhow mention functions coming from partial template
471  // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
472  // For now, hack in those functions manually only to fix documentation:
473 
481  get_divergence(const unsigned int q_point) const;
482 
492  get_symmetric_gradient(const unsigned int q_point) const;
493 
501  get_curl(const unsigned int q_point) const;
502 
518  void
520  const unsigned int q_point);
521 
538  void
541  const unsigned int q_point);
542 
555  void
557  const unsigned int q_point);
558 
559 #endif
560 
577  value_type
578  integrate_value() const;
579 
581 
590 
597  void
599  const VectorizedArrayType & value) const;
600 
605  template <typename T>
606  std::array<T, VectorizedArrayType::size()>
607  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
608  &array) const;
609 
614  template <typename T>
615  void
617  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
618  const std::array<T, VectorizedArrayType::size()> & value) const;
619 
624  get_matrix_free() const;
625 
626 protected:
637  const unsigned int dof_no,
638  const unsigned int first_selected_component,
639  const unsigned int quad_no,
640  const unsigned int fe_degree,
641  const unsigned int n_q_points,
642  const bool is_interior_face,
643  const unsigned int active_fe_index,
644  const unsigned int active_quad_index,
645  const unsigned int face_type);
646 
684  const Mapping<dim> & mapping,
685  const FiniteElement<dim> &fe,
686  const Quadrature<1> & quadrature,
687  const UpdateFlags update_flags,
688  const unsigned int first_selected_component,
690 
697  FEEvaluationBase(const FEEvaluationBase &other);
698 
706  operator=(const FEEvaluationBase &other);
707 
712 
719  template <typename VectorType, typename VectorOperation>
720  void
722  const VectorOperation & operation,
723  const std::array<VectorType *, n_components_> &vectors,
724  const std::array<
726  n_components_> & vectors_sm,
727  const std::bitset<VectorizedArrayType::size()> &mask,
728  const bool apply_constraints = true) const;
729 
737  template <typename VectorType, typename VectorOperation>
738  void
740  const VectorOperation & operation,
741  const std::array<VectorType *, n_components_> &vectors,
742  const std::array<
744  n_components_> & vectors_sm,
745  const std::bitset<VectorizedArrayType::size()> &mask) const;
746 
754  template <typename VectorType, typename VectorOperation>
755  void
757  const VectorOperation & operation,
758  const std::array<VectorType *, n_components_> &vectors) const;
759 
763  void
764  apply_hanging_node_constraints(const bool transpose) const;
765 
770 
775 
780  mutable std::vector<types::global_dof_index> local_dof_indices;
781 };
782 
783 
784 
792 template <int dim,
793  int n_components_,
794  typename Number,
795  bool is_face,
798  n_components_,
799  Number,
800  is_face,
801  VectorizedArrayType>
802 {
803  static_assert(
804  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
805  "Type of Number and of VectorizedArrayType do not match.");
806 
807 public:
808  using number_type = Number;
810  using gradient_type =
812  static constexpr unsigned int dimension = dim;
813  static constexpr unsigned int n_components = n_components_;
814  using BaseClass =
816 
817 protected:
827  const unsigned int dof_no,
828  const unsigned int first_selected_component,
829  const unsigned int quad_no,
830  const unsigned int fe_degree,
831  const unsigned int n_q_points,
832  const bool is_interior_face = true,
833  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
835  const unsigned int face_type = numbers::invalid_unsigned_int);
836 
842  const Mapping<dim> & mapping,
843  const FiniteElement<dim> &fe,
844  const Quadrature<1> & quadrature,
845  const UpdateFlags update_flags,
846  const unsigned int first_selected_component,
848 
853 
858  operator=(const FEEvaluationAccess &other);
859 };
860 
861 
862 
871 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
872 class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
873  : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
874 {
875  static_assert(
876  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
877  "Type of Number and of VectorizedArrayType do not match.");
878 
879 public:
880  using number_type = Number;
881  using value_type = VectorizedArrayType;
884  static constexpr unsigned int dimension = dim;
885  using BaseClass =
887 
891  value_type
892  get_dof_value(const unsigned int dof) const;
893 
897  void
898  submit_dof_value(const value_type val_in, const unsigned int dof);
899 
903  value_type
904  get_value(const unsigned int q_point) const;
905 
909  void
910  submit_value(const value_type val_in, const unsigned int q_point);
911 
915  void
917  const unsigned int q_point);
918 
923  get_gradient(const unsigned int q_point) const;
924 
928  value_type
929  get_normal_derivative(const unsigned int q_point) const;
930 
934  void
935  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
936 
940  void
941  submit_normal_derivative(const value_type grad_in,
942  const unsigned int q_point);
943 
948  get_hessian(unsigned int q_point) const;
949 
954  get_hessian_diagonal(const unsigned int q_point) const;
955 
959  void
960  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
961 
965  value_type
966  get_laplacian(const unsigned int q_point) const;
967 
971  value_type
972  integrate_value() const;
973 
974 protected:
984  const unsigned int dof_no,
985  const unsigned int first_selected_component,
986  const unsigned int quad_no,
987  const unsigned int fe_degree,
988  const unsigned int n_q_points,
989  const bool is_interior_face = true,
990  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
992  const unsigned int face_type = numbers::invalid_unsigned_int);
993 
999  const Mapping<dim> & mapping,
1000  const FiniteElement<dim> &fe,
1001  const Quadrature<1> & quadrature,
1002  const UpdateFlags update_flags,
1003  const unsigned int first_selected_component,
1005 
1009  FEEvaluationAccess(const FEEvaluationAccess &other);
1010 
1015  operator=(const FEEvaluationAccess &other);
1016 };
1017 
1018 
1019 
1029 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
1030 class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
1031  : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
1032 {
1033  static_assert(
1034  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1035  "Type of Number and of VectorizedArrayType do not match.");
1036 
1037 public:
1038  using number_type = Number;
1041  static constexpr unsigned int dimension = dim;
1042  static constexpr unsigned int n_components = dim;
1043  using BaseClass =
1045 
1050  get_gradient(const unsigned int q_point) const;
1051 
1056  VectorizedArrayType
1057  get_divergence(const unsigned int q_point) const;
1058 
1066  get_symmetric_gradient(const unsigned int q_point) const;
1067 
1073  get_curl(const unsigned int q_point) const;
1074 
1079  get_hessian(const unsigned int q_point) const;
1080 
1085  get_hessian_diagonal(const unsigned int q_point) const;
1086 
1090  void
1091  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1092 
1101  void
1103  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1104  const unsigned int q_point);
1105 
1114  void
1115  submit_divergence(const VectorizedArrayType div_in,
1116  const unsigned int q_point);
1117 
1126  void
1129  const unsigned int q_point);
1130 
1135  void
1137  const unsigned int q_point);
1138 
1139 protected:
1149  const unsigned int dof_no,
1150  const unsigned int first_selected_component,
1151  const unsigned int quad_no,
1152  const unsigned int dofs_per_cell,
1153  const unsigned int n_q_points,
1154  const bool is_interior_face = true,
1155  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1156  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1157  const unsigned int face_type = numbers::invalid_unsigned_int);
1158 
1164  const Mapping<dim> & mapping,
1165  const FiniteElement<dim> &fe,
1166  const Quadrature<1> & quadrature,
1167  const UpdateFlags update_flags,
1168  const unsigned int first_selected_component,
1170 
1174  FEEvaluationAccess(const FEEvaluationAccess &other);
1175 
1180  operator=(const FEEvaluationAccess &other);
1181 };
1182 
1183 
1192 template <typename Number, bool is_face, typename VectorizedArrayType>
1193 class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1194  : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1195 {
1196  static_assert(
1197  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1198  "Type of Number and of VectorizedArrayType do not match.");
1199 
1200 public:
1201  using number_type = Number;
1202  using value_type = VectorizedArrayType;
1205  static constexpr unsigned int dimension = 1;
1206  using BaseClass =
1208 
1212  value_type
1213  get_dof_value(const unsigned int dof) const;
1214 
1218  void
1219  submit_dof_value(const value_type val_in, const unsigned int dof);
1220 
1224  value_type
1225  get_value(const unsigned int q_point) const;
1226 
1230  void
1231  submit_value(const value_type val_in, const unsigned int q_point);
1232 
1236  void
1237  submit_value(const gradient_type val_in, const unsigned int q_point);
1238 
1243  get_gradient(const unsigned int q_point) const;
1244 
1248  value_type
1249  get_divergence(const unsigned int q_point) const;
1250 
1254  value_type
1255  get_normal_derivative(const unsigned int q_point) const;
1256 
1260  void
1261  submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1262 
1266  void
1267  submit_gradient(const value_type grad_in, const unsigned int q_point);
1268 
1272  void
1274  const unsigned int q_point);
1275 
1279  void
1280  submit_normal_derivative(const value_type grad_in,
1281  const unsigned int q_point);
1282 
1286  void
1288  const unsigned int q_point);
1289 
1293  hessian_type
1294  get_hessian(unsigned int q_point) const;
1295 
1300  get_hessian_diagonal(const unsigned int q_point) const;
1301 
1305  void
1306  submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1307 
1311  value_type
1312  get_laplacian(const unsigned int q_point) const;
1313 
1317  value_type
1318  integrate_value() const;
1319 
1320 protected:
1330  const unsigned int dof_no,
1331  const unsigned int first_selected_component,
1332  const unsigned int quad_no,
1333  const unsigned int fe_degree,
1334  const unsigned int n_q_points,
1335  const bool is_interior_face = true,
1336  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
1337  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
1338  const unsigned int face_type = numbers::invalid_unsigned_int);
1339 
1345  const Mapping<1> & mapping,
1346  const FiniteElement<1> &fe,
1347  const Quadrature<1> & quadrature,
1348  const UpdateFlags update_flags,
1349  const unsigned int first_selected_component,
1351 
1355  FEEvaluationAccess(const FEEvaluationAccess &other);
1356 
1361  operator=(const FEEvaluationAccess &other);
1362 };
1363 
1364 
1365 
1920 template <int dim,
1921  int fe_degree,
1922  int n_q_points_1d,
1923  int n_components_,
1924  typename Number,
1925  typename VectorizedArrayType>
1927  n_components_,
1928  Number,
1929  false,
1930  VectorizedArrayType>
1931 {
1932  static_assert(
1933  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1934  "Type of Number and of VectorizedArrayType do not match.");
1935 
1936 public:
1940  using BaseClass =
1942 
1946  using number_type = Number;
1947 
1954 
1961 
1965  static constexpr unsigned int dimension = dim;
1966 
1971  static constexpr unsigned int n_components = n_components_;
1972 
1979  static constexpr unsigned int static_n_q_points =
1980  Utilities::pow(n_q_points_1d, dim);
1981 
1989  static constexpr unsigned int static_dofs_per_component =
1990  Utilities::pow(fe_degree + 1, dim);
1991 
1999  static constexpr unsigned int tensor_dofs_per_cell =
2000  static_dofs_per_component * n_components;
2001 
2009  static constexpr unsigned int static_dofs_per_cell =
2010  static_dofs_per_component * n_components;
2011 
2046  FEEvaluation(
2048  const unsigned int dof_no = 0,
2049  const unsigned int quad_no = 0,
2050  const unsigned int first_selected_component = 0,
2051  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2052  const unsigned int active_quad_index = numbers::invalid_unsigned_int);
2053 
2062  const std::pair<unsigned int, unsigned int> & range,
2063  const unsigned int dof_no = 0,
2064  const unsigned int quad_no = 0,
2065  const unsigned int first_selected_component = 0);
2066 
2095  FEEvaluation(const Mapping<dim> & mapping,
2096  const FiniteElement<dim> &fe,
2097  const Quadrature<1> & quadrature,
2098  const UpdateFlags update_flags,
2099  const unsigned int first_selected_component = 0);
2100 
2106  FEEvaluation(const FiniteElement<dim> &fe,
2107  const Quadrature<1> & quadrature,
2108  const UpdateFlags update_flags,
2109  const unsigned int first_selected_component = 0);
2110 
2121  FEEvaluation(const FiniteElement<dim> & fe,
2123  const unsigned int first_selected_component = 0);
2124 
2131  FEEvaluation(const FEEvaluation &other);
2132 
2139  FEEvaluation &
2140  operator=(const FEEvaluation &other);
2141 
2150  void
2151  reinit(const unsigned int cell_batch_index);
2152 
2165  template <bool level_dof_access>
2166  void
2168 
2179  void
2181 
2185  static bool
2186  fast_evaluation_supported(const unsigned int given_degree,
2187  const unsigned int give_n_q_points_1d);
2188 
2198  void
2199  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2200 
2205  DEAL_II_DEPRECATED_EARLY void
2206  evaluate(const bool evaluate_values,
2207  const bool evaluate_gradients,
2208  const bool evaluate_hessians = false);
2209 
2222  void
2223  evaluate(const VectorizedArrayType * values_array,
2224  const EvaluationFlags::EvaluationFlags evaluation_flag);
2225 
2230  DEAL_II_DEPRECATED_EARLY void
2231  evaluate(const VectorizedArrayType *values_array,
2232  const bool evaluate_values,
2233  const bool evaluate_gradients,
2234  const bool evaluate_hessians = false);
2235 
2249  template <typename VectorType>
2250  void
2251  gather_evaluate(const VectorType & input_vector,
2252  const EvaluationFlags::EvaluationFlags evaluation_flag);
2253 
2257  template <typename VectorType>
2258  DEAL_II_DEPRECATED_EARLY void
2259  gather_evaluate(const VectorType &input_vector,
2260  const bool evaluate_values,
2261  const bool evaluate_gradients,
2262  const bool evaluate_hessians = false);
2263 
2274  void
2275  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2276 
2280  DEAL_II_DEPRECATED_EARLY void
2281  integrate(const bool integrate_values, const bool integrate_gradients);
2282 
2294  void
2295  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2296  VectorizedArrayType * values_array,
2297  const bool sum_into_values = false);
2298 
2302  DEAL_II_DEPRECATED_EARLY void
2303  integrate(const bool integrate_values,
2304  const bool integrate_gradients,
2305  VectorizedArrayType *values_array);
2306 
2320  template <typename VectorType>
2321  void
2322  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
2323  VectorType & output_vector);
2324 
2328  template <typename VectorType>
2329  DEAL_II_DEPRECATED_EARLY void
2330  integrate_scatter(const bool integrate_values,
2331  const bool integrate_gradients,
2332  VectorType &output_vector);
2333 
2340  const unsigned int dofs_per_component;
2341 
2348  const unsigned int dofs_per_cell;
2349 
2357  const unsigned int n_q_points;
2358 
2359 private:
2364  void
2365  check_template_arguments(const unsigned int fe_no,
2366  const unsigned int first_selected_component);
2367 };
2368 
2369 
2370 
2406 template <int dim,
2407  int fe_degree,
2408  int n_q_points_1d = fe_degree + 1,
2409  int n_components_ = 1,
2410  typename Number = double,
2411  typename VectorizedArrayType = VectorizedArray<Number>>
2413  n_components_,
2414  Number,
2415  true,
2416  VectorizedArrayType>
2417 {
2418  static_assert(
2419  std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2420  "Type of Number and of VectorizedArrayType do not match.");
2421 
2422 public:
2426  using BaseClass =
2428 
2432  using number_type = Number;
2433 
2440 
2447 
2451  static constexpr unsigned int dimension = dim;
2452 
2457  static constexpr unsigned int n_components = n_components_;
2458 
2466  static constexpr unsigned int static_n_q_points =
2467  Utilities::pow(n_q_points_1d, dim - 1);
2468 
2475  static constexpr unsigned int static_n_q_points_cell =
2476  Utilities::pow(n_q_points_1d, dim);
2477 
2484  static constexpr unsigned int static_dofs_per_component =
2485  Utilities::pow(fe_degree + 1, dim);
2486 
2493  static constexpr unsigned int tensor_dofs_per_cell =
2494  static_dofs_per_component * n_components;
2495 
2502  static constexpr unsigned int static_dofs_per_cell =
2503  static_dofs_per_component * n_components;
2504 
2548  const bool is_interior_face = true,
2549  const unsigned int dof_no = 0,
2550  const unsigned int quad_no = 0,
2551  const unsigned int first_selected_component = 0,
2552  const unsigned int active_fe_index = numbers::invalid_unsigned_int,
2553  const unsigned int active_quad_index = numbers::invalid_unsigned_int,
2554  const unsigned int face_type = numbers::invalid_unsigned_int);
2555 
2565  const std::pair<unsigned int, unsigned int> & range,
2566  const bool is_interior_face = true,
2567  const unsigned int dof_no = 0,
2568  const unsigned int quad_no = 0,
2569  const unsigned int first_selected_component = 0);
2570 
2581  void
2582  reinit(const unsigned int face_batch_number);
2583 
2591  void
2592  reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2593 
2597  static bool
2598  fast_evaluation_supported(const unsigned int given_degree,
2599  const unsigned int give_n_q_points_1d);
2600 
2611  void
2612  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag);
2613 
2617  DEAL_II_DEPRECATED_EARLY void
2618  evaluate(const bool evaluate_values, const bool evaluate_gradients);
2619 
2632  void
2633  evaluate(const VectorizedArrayType * values_array,
2634  const EvaluationFlags::EvaluationFlags evaluation_flag);
2635 
2639  DEAL_II_DEPRECATED_EARLY void
2640  evaluate(const VectorizedArrayType *values_array,
2641  const bool evaluate_values,
2642  const bool evaluate_gradients);
2643 
2655  template <typename VectorType>
2656  void
2657  gather_evaluate(const VectorType & input_vector,
2658  const EvaluationFlags::EvaluationFlags evaluation_flag);
2659 
2663  template <typename VectorType>
2664  DEAL_II_DEPRECATED_EARLY void
2665  gather_evaluate(const VectorType &input_vector,
2666  const bool evaluate_values,
2667  const bool evaluate_gradients);
2668 
2678  void
2679  integrate(const EvaluationFlags::EvaluationFlags integration_flag);
2680 
2684  DEAL_II_DEPRECATED_EARLY void
2685  integrate(const bool integrate_values, const bool integrate_gradients);
2686 
2695  void
2696  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
2697  VectorizedArrayType * values_array);
2698 
2702  DEAL_II_DEPRECATED_EARLY void
2703  integrate(const bool integrate_values,
2704  const bool integrate_gradients,
2705  VectorizedArrayType *values_array);
2706 
2718  template <typename VectorType>
2719  void
2720  integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
2721  VectorType & output_vector);
2722 
2726  template <typename VectorType>
2727  void
2728  integrate_scatter(const bool integrate_values,
2729  const bool integrate_gradients,
2730  VectorType &output_vector);
2731 
2738  const unsigned int dofs_per_component;
2739 
2746  const unsigned int dofs_per_cell;
2747 
2755  const unsigned int n_q_points;
2756 };
2757 
2758 
2759 
2760 namespace internal
2761 {
2762  namespace MatrixFreeFunctions
2763  {
2764  // a helper function to compute the number of DoFs of a DGP element at
2765  // compile time, depending on the degree
2766  template <int dim, int degree>
2768  {
2769  // this division is always without remainder
2770  static constexpr unsigned int value =
2771  (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2772  };
2773 
2774  // base specialization: 1d elements have 'degree+1' degrees of freedom
2775  template <int degree>
2776  struct DGP_dofs_per_component<1, degree>
2777  {
2778  static constexpr unsigned int value = degree + 1;
2779  };
2780  } // namespace MatrixFreeFunctions
2781 } // namespace internal
2782 
2783 
2784 /*----------------------- Inline functions ----------------------------------*/
2785 
2786 #ifndef DOXYGEN
2787 
2788 
2789 namespace internal
2790 {
2791  // Extract all internal data pointers and indices in a single function that
2792  // get passed on to the constructor of FEEvaluationData, avoiding to look
2793  // things up multiple times
2794  template <bool is_face,
2795  int dim,
2796  typename Number,
2797  typename VectorizedArrayType>
2799  InitializationData
2800  extract_initialization_data(
2802  const unsigned int dof_no,
2803  const unsigned int first_selected_component,
2804  const unsigned int quad_no,
2805  const unsigned int fe_degree,
2806  const unsigned int n_q_points,
2807  const unsigned int active_fe_index_given,
2808  const unsigned int active_quad_index_given,
2809  const unsigned int face_type)
2810  {
2812  InitializationData init_data;
2813 
2814  init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2815  init_data.mapping_data =
2816  &internal::MatrixFreeFunctions::
2817  MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2818  matrix_free.get_mapping_info(), quad_no);
2819 
2820  init_data.active_fe_index =
2821  fe_degree != numbers::invalid_unsigned_int ?
2822  init_data.dof_info->fe_index_from_degree(first_selected_component,
2823  fe_degree) :
2824  (active_fe_index_given != numbers::invalid_unsigned_int ?
2825  active_fe_index_given :
2826  0);
2827  init_data.active_quad_index =
2828  fe_degree == numbers::invalid_unsigned_int ?
2829  (active_quad_index_given != numbers::invalid_unsigned_int ?
2830  active_quad_index_given :
2831  std::min<unsigned int>(init_data.active_fe_index,
2832  init_data.mapping_data->descriptor.size() -
2833  1)) :
2834  init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2835 
2836  init_data.shape_info = &matrix_free.get_shape_info(
2837  dof_no,
2838  quad_no,
2839  init_data.dof_info->component_to_base_index[first_selected_component],
2840  init_data.active_fe_index,
2841  init_data.active_quad_index);
2842  init_data.descriptor =
2843  &init_data.mapping_data->descriptor
2844  [is_face ?
2845  (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2846  (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2847  init_data.active_quad_index];
2848 
2849  return init_data;
2850  }
2851 } // namespace internal
2852 
2853 
2854 
2855 /*----------------------- FEEvaluationBase ----------------------------------*/
2856 
2857 template <int dim,
2858  int n_components_,
2859  typename Number,
2860  bool is_face,
2861  typename VectorizedArrayType>
2862 inline FEEvaluationBase<dim,
2863  n_components_,
2864  Number,
2865  is_face,
2866  VectorizedArrayType>::
2869  const unsigned int dof_no,
2870  const unsigned int first_selected_component,
2871  const unsigned int quad_no,
2872  const unsigned int fe_degree,
2873  const unsigned int n_q_points,
2874  const bool is_interior_face,
2875  const unsigned int active_fe_index,
2876  const unsigned int active_quad_index,
2877  const unsigned int face_type)
2879  internal::extract_initialization_data<is_face>(matrix_free,
2880  dof_no,
2882  quad_no,
2883  fe_degree,
2884  n_q_points,
2887  face_type),
2888  is_interior_face,
2889  quad_no,
2890  first_selected_component)
2891  , scratch_data_array(matrix_free.acquire_scratch_data())
2892  , matrix_free(&matrix_free)
2893 {
2894  this->set_data_pointers(scratch_data_array, n_components_);
2895  Assert(
2896  this->dof_info->start_components.back() == 1 ||
2897  static_cast<int>(n_components_) <=
2898  static_cast<int>(
2899  this->dof_info->start_components
2900  [this->dof_info->component_to_base_index[first_selected_component] +
2901  1]) -
2903  ExcMessage(
2904  "You tried to construct a vector-valued evaluator with " +
2905  std::to_string(n_components) +
2906  " components. However, "
2907  "the current base element has only " +
2909  this->dof_info->start_components
2910  [this->dof_info->component_to_base_index[first_selected_component] +
2911  1] -
2912  first_selected_component) +
2913  " components left when starting from local element index " +
2915  first_selected_component -
2916  this->dof_info->start_components
2917  [this->dof_info->component_to_base_index[first_selected_component]]) +
2918  " (global index " + std::to_string(first_selected_component) + ")"));
2919 
2920  // do not check for correct dimensions of data fields here, should be done
2921  // in derived classes
2922 }
2923 
2924 
2925 
2926 template <int dim,
2927  int n_components_,
2928  typename Number,
2929  bool is_face,
2930  typename VectorizedArrayType>
2931 inline FEEvaluationBase<dim,
2932  n_components_,
2933  Number,
2934  is_face,
2935  VectorizedArrayType>::
2937  const Mapping<dim> & mapping,
2938  const FiniteElement<dim> &fe,
2939  const Quadrature<1> & quadrature,
2940  const UpdateFlags update_flags,
2941  const unsigned int first_selected_component,
2944  other != nullptr &&
2945  other->mapped_geometry->get_quadrature() == quadrature ?
2946  other->mapped_geometry :
2947  std::make_shared<internal::MatrixFreeFunctions::
2948  MappingDataOnTheFly<dim, VectorizedArrayType>>(
2949  mapping,
2950  quadrature,
2951  update_flags),
2952  n_components_,
2953  first_selected_component)
2955  , matrix_free(nullptr)
2956 {
2957  const unsigned int base_element_number =
2958  fe.component_to_base_index(first_selected_component).first;
2959  Assert(fe.element_multiplicity(base_element_number) == 1 ||
2960  fe.element_multiplicity(base_element_number) -
2961  first_selected_component >=
2962  n_components_,
2963  ExcMessage("The underlying element must at least contain as many "
2964  "components as requested by this class"));
2965  (void)base_element_number;
2966 
2967  Assert(this->data == nullptr, ExcInternalError());
2968  this->data =
2970  Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
2971  fe,
2972  fe.component_to_base_index(first_selected_component).first);
2973 
2974  this->set_data_pointers(scratch_data_array, n_components_);
2975 }
2976 
2977 
2978 
2979 template <int dim,
2980  int n_components_,
2981  typename Number,
2982  bool is_face,
2983  typename VectorizedArrayType>
2984 inline FEEvaluationBase<dim,
2985  n_components_,
2986  Number,
2987  is_face,
2988  VectorizedArrayType>::
2990  n_components_,
2991  Number,
2992  is_face,
2993  VectorizedArrayType> &other)
2995  , scratch_data_array(other.matrix_free == nullptr ?
2997  other.matrix_free->acquire_scratch_data())
2998  , matrix_free(other.matrix_free)
2999 {
3000  if (other.matrix_free == nullptr)
3001  {
3002  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3003  this->data =
3005  *other.data);
3006 
3007  // Create deep copy of mapped geometry for use in parallel
3008  this->mapped_geometry =
3009  std::make_shared<internal::MatrixFreeFunctions::
3010  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3011  other.mapped_geometry->get_fe_values().get_mapping(),
3012  other.mapped_geometry->get_quadrature(),
3013  other.mapped_geometry->get_fe_values().get_update_flags());
3014  this->mapping_data = &this->mapped_geometry->get_data_storage();
3015  this->cell = 0;
3016 
3017  this->jacobian =
3018  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3019  this->J_value =
3020  this->mapped_geometry->get_data_storage().JxW_values.begin();
3021  this->jacobian_gradients =
3022  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3023  this->quadrature_points =
3024  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3025  }
3026 
3027  this->set_data_pointers(scratch_data_array, n_components_);
3028 }
3029 
3030 
3031 
3032 template <int dim,
3033  int n_components_,
3034  typename Number,
3035  bool is_face,
3036  typename VectorizedArrayType>
3037 inline FEEvaluationBase<dim,
3038  n_components_,
3039  Number,
3040  is_face,
3041  VectorizedArrayType> &
3043 operator=(const FEEvaluationBase<dim,
3044  n_components_,
3045  Number,
3046  is_face,
3047  VectorizedArrayType> &other)
3048 {
3049  // release old memory
3050  if (matrix_free == nullptr)
3051  {
3052  delete this->data;
3053  delete scratch_data_array;
3054  }
3055  else
3056  {
3057  matrix_free->release_scratch_data(scratch_data_array);
3058  }
3059 
3061 
3062  matrix_free = other.matrix_free;
3063 
3064  if (other.matrix_free == nullptr)
3065  {
3066  Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3067  this->data =
3069  *other.data);
3071 
3072  // Create deep copy of mapped geometry for use in parallel
3073  this->mapped_geometry =
3074  std::make_shared<internal::MatrixFreeFunctions::
3075  MappingDataOnTheFly<dim, VectorizedArrayType>>(
3076  other.mapped_geometry->get_fe_values().get_mapping(),
3077  other.mapped_geometry->get_quadrature(),
3078  other.mapped_geometry->get_fe_values().get_update_flags());
3079  this->cell = 0;
3080  this->mapping_data = &this->mapped_geometry->get_data_storage();
3081  this->jacobian =
3082  this->mapped_geometry->get_data_storage().jacobians[0].begin();
3083  this->J_value =
3084  this->mapped_geometry->get_data_storage().JxW_values.begin();
3085  this->jacobian_gradients =
3086  this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3087  this->quadrature_points =
3088  this->mapped_geometry->get_data_storage().quadrature_points.begin();
3089  }
3090  else
3091  {
3092  scratch_data_array = matrix_free->acquire_scratch_data();
3093  }
3094 
3095  this->set_data_pointers(scratch_data_array, n_components_);
3096 
3097  return *this;
3098 }
3099 
3100 
3101 
3102 template <int dim,
3103  int n_components_,
3104  typename Number,
3105  bool is_face,
3106  typename VectorizedArrayType>
3107 inline FEEvaluationBase<dim,
3108  n_components_,
3109  Number,
3110  is_face,
3111  VectorizedArrayType>::~FEEvaluationBase()
3112 {
3113  if (matrix_free != nullptr)
3114  {
3115  try
3116  {
3117  matrix_free->release_scratch_data(scratch_data_array);
3118  }
3119  catch (...)
3120  {}
3121  }
3122  else
3123  {
3124  delete scratch_data_array;
3125  delete this->data;
3126  }
3127 }
3128 
3129 
3130 
3131 template <int dim,
3132  int n_components_,
3133  typename Number,
3134  bool is_face,
3135  typename VectorizedArrayType>
3138  get_matrix_free() const
3139 {
3140  Assert(matrix_free != nullptr,
3141  ExcMessage(
3142  "FEEvaluation was not initialized with a MatrixFree object!"));
3143  return *matrix_free;
3144 }
3145 
3146 
3147 namespace internal
3148 {
3149  template <int dim,
3150  int n_components,
3151  typename Number,
3152  bool is_face,
3153  typename VectorizedArrayType,
3154  typename VectorizedArrayType2,
3155  typename GlobalVectorType,
3156  typename FU>
3157  inline void
3158  process_cell_data(
3159  const FEEvaluationBase<dim,
3160  n_components,
3161  Number,
3162  is_face,
3163  VectorizedArrayType> & fe_eval,
3165  GlobalVectorType & array,
3166  VectorizedArrayType2 & out,
3167  const FU & fu)
3168  {
3169  (void)matrix_free;
3170  Assert(matrix_free != nullptr, ExcNotImplemented());
3171  AssertDimension(array.size(),
3172  matrix_free->get_task_info().cell_partition_data.back());
3173 
3174  // 1) collect ids of cell
3175  const auto cells = fe_eval.get_cell_ids();
3176 
3177  // 2) actually gather values
3178  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
3179  if (cells[i] != numbers::invalid_unsigned_int)
3180  fu(out[i],
3181  array[cells[i] / VectorizedArrayType::size()]
3182  [cells[i] % VectorizedArrayType::size()]);
3183  }
3184 } // namespace internal
3185 
3186 
3187 
3188 template <int dim,
3189  int n_components_,
3190  typename Number,
3191  bool is_face,
3192  typename VectorizedArrayType>
3193 inline VectorizedArrayType
3196 {
3197  VectorizedArrayType out = Number(1.);
3198  internal::process_cell_data(
3199  *this, this->matrix_free, array, out, [](auto &local, const auto &global) {
3200  local = global;
3201  });
3202  return out;
3203 }
3204 
3205 
3206 
3207 template <int dim,
3208  int n_components_,
3209  typename Number,
3210  bool is_face,
3211  typename VectorizedArrayType>
3212 inline void
3215  const VectorizedArrayType & in) const
3216 {
3217  internal::process_cell_data(
3218  *this, this->matrix_free, array, in, [](const auto &local, auto &global) {
3219  global = local;
3220  });
3221 }
3222 
3223 
3224 
3225 template <int dim,
3226  int n_components_,
3227  typename Number,
3228  bool is_face,
3229  typename VectorizedArrayType>
3230 template <typename T>
3231 inline std::array<T, VectorizedArrayType::size()>
3233  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
3234  &array) const
3235 {
3236  std::array<T, VectorizedArrayType::size()> out;
3237  internal::process_cell_data(
3238  *this, this->matrix_free, array, out, [](auto &local, const auto &global) {
3239  local = global;
3240  });
3241  return out;
3242 }
3243 
3244 
3245 
3246 template <int dim,
3247  int n_components_,
3248  typename Number,
3249  bool is_face,
3250  typename VectorizedArrayType>
3251 template <typename T>
3252 inline void
3255  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
3256  const std::array<T, VectorizedArrayType::size()> & in) const
3257 {
3258  internal::process_cell_data(
3259  *this, this->matrix_free, array, in, [](const auto &local, auto &global) {
3260  global = local;
3261  });
3262 }
3263 
3264 
3265 
3266 namespace internal
3267 {
3268  // given a block vector return the underlying vector type
3269  // including constness (specified by bool)
3270  template <typename VectorType, bool>
3271  struct ConstBlockVectorSelector;
3272 
3273  template <typename VectorType>
3274  struct ConstBlockVectorSelector<VectorType, true>
3275  {
3276  using BaseVectorType = const typename VectorType::BlockType;
3277  };
3278 
3279  template <typename VectorType>
3280  struct ConstBlockVectorSelector<VectorType, false>
3281  {
3282  using BaseVectorType = typename VectorType::BlockType;
3283  };
3284 
3285  // allows to select between block vectors and non-block vectors, which
3286  // allows to use a unified interface for extracting blocks on block vectors
3287  // and doing nothing on usual vectors
3288  template <typename VectorType, bool>
3289  struct BlockVectorSelector;
3290 
3291  template <typename VectorType>
3292  struct BlockVectorSelector<VectorType, true>
3293  {
3294  using BaseVectorType = typename ConstBlockVectorSelector<
3295  VectorType,
3296  std::is_const<VectorType>::value>::BaseVectorType;
3297 
3298  static BaseVectorType *
3299  get_vector_component(VectorType &vec, const unsigned int component)
3300  {
3301  AssertIndexRange(component, vec.n_blocks());
3302  return &vec.block(component);
3303  }
3304  };
3305 
3306  template <typename VectorType>
3307  struct BlockVectorSelector<VectorType, false>
3308  {
3309  using BaseVectorType = VectorType;
3310 
3311  static BaseVectorType *
3312  get_vector_component(VectorType &vec, const unsigned int component)
3313  {
3314  // FEEvaluation allows to combine several vectors from a scalar
3315  // FiniteElement into a "vector-valued" FEEvaluation object with
3316  // multiple components. These components can be extracted with the other
3317  // get_vector_component functions. If we do not get a vector of vectors
3318  // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3319  // must make sure that we do not duplicate the components in input
3320  // and/or duplicate the resulting integrals. In such a case, we should
3321  // only get the zeroth component in the vector contained set nullptr for
3322  // the others which allows us to catch unintended use in
3323  // read_write_operation.
3324  if (component == 0)
3325  return &vec;
3326  else
3327  return nullptr;
3328  }
3329  };
3330 
3331  template <typename VectorType>
3332  struct BlockVectorSelector<std::vector<VectorType>, false>
3333  {
3334  using BaseVectorType = VectorType;
3335 
3336  static BaseVectorType *
3337  get_vector_component(std::vector<VectorType> &vec,
3338  const unsigned int component)
3339  {
3340  AssertIndexRange(component, vec.size());
3341  return &vec[component];
3342  }
3343  };
3344 
3345  template <typename VectorType>
3346  struct BlockVectorSelector<const std::vector<VectorType>, false>
3347  {
3348  using BaseVectorType = const VectorType;
3349 
3350  static const BaseVectorType *
3351  get_vector_component(const std::vector<VectorType> &vec,
3352  const unsigned int component)
3353  {
3354  AssertIndexRange(component, vec.size());
3355  return &vec[component];
3356  }
3357  };
3358 
3359  template <typename VectorType>
3360  struct BlockVectorSelector<std::vector<VectorType *>, false>
3361  {
3362  using BaseVectorType = VectorType;
3363 
3364  static BaseVectorType *
3365  get_vector_component(std::vector<VectorType *> &vec,
3366  const unsigned int component)
3367  {
3368  AssertIndexRange(component, vec.size());
3369  return vec[component];
3370  }
3371  };
3372 
3373  template <typename VectorType>
3374  struct BlockVectorSelector<const std::vector<VectorType *>, false>
3375  {
3376  using BaseVectorType = const VectorType;
3377 
3378  static const BaseVectorType *
3379  get_vector_component(const std::vector<VectorType *> &vec,
3380  const unsigned int component)
3381  {
3382  AssertIndexRange(component, vec.size());
3383  return vec[component];
3384  }
3385  };
3386 } // namespace internal
3387 
3388 
3389 
3390 template <int dim,
3391  int n_components_,
3392  typename Number,
3393  bool is_face,
3394  typename VectorizedArrayType>
3395 template <typename VectorType, typename VectorOperation>
3396 inline void
3399  const VectorOperation & operation,
3400  const std::array<VectorType *, n_components_> &src,
3401  const std::array<
3403  n_components_> & src_sm,
3404  const std::bitset<VectorizedArrayType::size()> &mask,
3405  const bool apply_constraints) const
3406 {
3407  // Case 1: No MatrixFree object given, simple case because we do not need to
3408  // process constraints and need not care about vectorization -> go to
3409  // separate function
3410  if (this->matrix_free == nullptr)
3411  {
3412  read_write_operation_global(operation, src);
3413  return;
3414  }
3415 
3416  Assert(this->dof_info != nullptr, ExcNotInitialized());
3417  Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3418  if (this->n_fe_components == 1)
3419  for (unsigned int comp = 0; comp < n_components; ++comp)
3420  {
3421  Assert(src[comp] != nullptr,
3422  ExcMessage("The finite element underlying this FEEvaluation "
3423  "object is scalar, but you requested " +
3424  std::to_string(n_components) +
3425  " components via the template argument in "
3426  "FEEvaluation. In that case, you must pass an "
3427  "std::vector<VectorType> or a BlockVector to " +
3428  "read_dof_values and distribute_local_to_global."));
3429  internal::check_vector_compatibility(*src[comp], *this->dof_info);
3430  }
3431  else
3432  {
3434  }
3435 
3436  // Case 2: contiguous indices which use reduced storage of indices and can
3437  // use vectorized load/store operations -> go to separate function
3438  if (this->cell != numbers::invalid_unsigned_int)
3439  {
3441  this->cell,
3442  this->dof_info->index_storage_variants[this->dof_access_index].size());
3443  if (this->dof_info->index_storage_variants
3444  [is_face ? this->dof_access_index :
3446  [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::
3447  IndexStorageVariants::contiguous)
3448  {
3449  read_write_operation_contiguous(operation, src, src_sm, mask);
3450  return;
3451  }
3452  }
3453 
3454  // Case 3: standard operation with one index per degree of freedom -> go on
3455  // here
3456  constexpr unsigned int n_lanes = VectorizedArrayType::size();
3457  Assert(mask.count() == n_lanes,
3458  ExcNotImplemented("Masking currently not implemented for "
3459  "non-contiguous DoF storage"));
3460 
3461  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3462  this->get_cell_ids();
3463 
3464  bool has_hn_constraints = false;
3465 
3466  if (is_face == false)
3467  {
3468  for (unsigned int v = 0; v < n_lanes; ++v)
3469  if (cells[v] != numbers::invalid_unsigned_int &&
3470  this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3472  [cells[v] * this->n_fe_components +
3473  this->first_selected_component] !=
3475  has_hn_constraints = true;
3476  }
3477 
3478  std::integral_constant<bool,
3480  vector_selector;
3481 
3482  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3483  std::array<VectorizedArrayType *, n_components> values_dofs;
3484  for (unsigned int c = 0; c < n_components; ++c)
3485  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3486  c * dofs_per_component;
3487 
3488  if (this->dof_info->index_storage_variants
3489  [is_face ? this->dof_access_index :
3491  [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3492  IndexStorageVariants::interleaved &&
3493  (has_hn_constraints == false))
3494  {
3495  const unsigned int *dof_indices =
3496  this->dof_info->dof_indices_interleaved.data() +
3497  this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes]
3498  .first +
3499  this->dof_info
3501  [this->first_selected_component] *
3502  n_lanes;
3503  if (n_components == 1 || this->n_fe_components == 1)
3504  for (unsigned int i = 0; i < dofs_per_component;
3505  ++i, dof_indices += n_lanes)
3506  for (unsigned int comp = 0; comp < n_components; ++comp)
3507  operation.process_dof_gather(dof_indices,
3508  *src[comp],
3509  0,
3510  values_dofs[comp][i],
3511  vector_selector);
3512  else
3513  for (unsigned int comp = 0; comp < n_components; ++comp)
3514  for (unsigned int i = 0; i < dofs_per_component;
3515  ++i, dof_indices += n_lanes)
3516  operation.process_dof_gather(
3517  dof_indices, *src[0], 0, values_dofs[comp][i], vector_selector);
3518  return;
3519  }
3520 
3521  // Allocate pointers, then initialize all of them to nullptrs and
3522  // below overwrite the ones we actually use:
3523  std::array<const unsigned int *, n_lanes> dof_indices;
3524  dof_indices.fill(nullptr);
3525 
3526  // Assign the appropriate cell ids for face/cell case and get the pointers
3527  // to the dof indices of the cells on all lanes
3528 
3529  bool has_constraints = false;
3530  const unsigned int n_components_read =
3531  this->n_fe_components > 1 ? n_components : 1;
3532 
3533  if (is_face)
3534  {
3535  for (unsigned int v = 0; v < n_lanes; ++v)
3536  {
3537  if (cells[v] == numbers::invalid_unsigned_int)
3538  continue;
3539 
3540  Assert(cells[v] < this->dof_info->row_starts.size() - 1,
3541  ExcInternalError());
3542  const std::pair<unsigned int, unsigned int> *my_index_start =
3543  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3544  this->first_selected_component];
3545 
3546  // check whether any of the SIMD lanes has constraints, i.e., the
3547  // constraint indicator which is the second entry of row_starts
3548  // increments on this cell
3549  if (my_index_start[n_components_read].second !=
3550  my_index_start[0].second)
3551  has_constraints = true;
3552 
3553  dof_indices[v] =
3554  this->dof_info->dof_indices.data() + my_index_start[0].first;
3555  }
3556  }
3557  else
3558  {
3559  for (unsigned int v = 0; v < n_lanes; ++v)
3560  {
3561  if (cells[v] == numbers::invalid_unsigned_int)
3562  continue;
3563 
3564  const std::pair<unsigned int, unsigned int> *my_index_start =
3565  &this->dof_info->row_starts[cells[v] * this->n_fe_components +
3566  this->first_selected_component];
3567  if (my_index_start[n_components_read].second !=
3568  my_index_start[0].second)
3569  has_constraints = true;
3570 
3571  if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3573  [cells[v] * this->n_fe_components +
3574  this->first_selected_component] !=
3576  has_hn_constraints = true;
3577 
3578  Assert(my_index_start[n_components_read].first ==
3579  my_index_start[0].first ||
3580  my_index_start[0].first < this->dof_info->dof_indices.size(),
3581  ExcIndexRange(0,
3582  my_index_start[0].first,
3583  this->dof_info->dof_indices.size()));
3584  dof_indices[v] =
3585  this->dof_info->dof_indices.data() + my_index_start[0].first;
3586  }
3587  }
3588 
3589  if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3590  return i != numbers::invalid_unsigned_int;
3591  }) < n_lanes)
3592  for (unsigned int comp = 0; comp < n_components; ++comp)
3593  for (unsigned int i = 0; i < dofs_per_component; ++i)
3594  operation.process_empty(values_dofs[comp][i]);
3595 
3596  // Case where we have no constraints throughout the whole cell: Can go
3597  // through the list of DoFs directly
3598  if (!has_constraints && apply_constraints)
3599  {
3600  if (n_components == 1 || this->n_fe_components == 1)
3601  {
3602  for (unsigned int v = 0; v < n_lanes; ++v)
3603  {
3604  if (cells[v] == numbers::invalid_unsigned_int)
3605  continue;
3606 
3607  for (unsigned int i = 0; i < dofs_per_component; ++i)
3608  for (unsigned int comp = 0; comp < n_components; ++comp)
3609  operation.process_dof(dof_indices[v][i],
3610  *src[comp],
3611  values_dofs[comp][i][v]);
3612  }
3613  }
3614  else
3615  {
3616  for (unsigned int comp = 0; comp < n_components; ++comp)
3617  for (unsigned int v = 0; v < n_lanes; ++v)
3618  {
3619  if (cells[v] == numbers::invalid_unsigned_int)
3620  continue;
3621 
3622  for (unsigned int i = 0; i < dofs_per_component; ++i)
3623  operation.process_dof(
3624  dof_indices[v][comp * dofs_per_component + i],
3625  *src[0],
3626  values_dofs[comp][i][v]);
3627  }
3628  }
3629  return;
3630  }
3631 
3632  // In the case where there are some constraints to be resolved, loop over
3633  // all vector components that are filled and then over local dofs. ind_local
3634  // holds local number on cell, index iterates over the elements of
3635  // index_local_to_global and dof_indices points to the global indices stored
3636  // in index_local_to_global
3637 
3638  for (unsigned int v = 0; v < n_lanes; ++v)
3639  {
3640  if (cells[v] == numbers::invalid_unsigned_int)
3641  continue;
3642 
3643  const unsigned int cell_index = cells[v];
3644  const unsigned int cell_dof_index =
3645  cell_index * this->n_fe_components + this->first_selected_component;
3646  const unsigned int n_components_read =
3647  this->n_fe_components > 1 ? n_components : 1;
3648  unsigned int index_indicators =
3649  this->dof_info->row_starts[cell_dof_index].second;
3650  unsigned int next_index_indicators =
3651  this->dof_info->row_starts[cell_dof_index + 1].second;
3652 
3653  // For read_dof_values_plain, redirect the dof_indices field to the
3654  // unconstrained indices
3655  if (apply_constraints == false &&
3656  (this->dof_info->row_starts[cell_dof_index].second !=
3657  this->dof_info->row_starts[cell_dof_index + n_components_read]
3658  .second ||
3659  (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
3660  this->dof_info->hanging_node_constraint_masks[cell_dof_index] !=
3662  {
3663  Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
3665  ExcNotInitialized());
3666  dof_indices[v] =
3667  this->dof_info->plain_dof_indices.data() +
3668  this->dof_info
3670  [this->first_selected_component] +
3672  next_index_indicators = index_indicators;
3673  }
3674 
3675  if (n_components == 1 || this->n_fe_components == 1)
3676  {
3677  unsigned int ind_local = 0;
3678  for (; index_indicators != next_index_indicators; ++index_indicators)
3679  {
3680  const std::pair<unsigned short, unsigned short> indicator =
3681  this->dof_info->constraint_indicator[index_indicators];
3682  // run through values up to next constraint
3683  for (unsigned int j = 0; j < indicator.first; ++j)
3684  for (unsigned int comp = 0; comp < n_components; ++comp)
3685  operation.process_dof(dof_indices[v][j],
3686  *src[comp],
3687  values_dofs[comp][ind_local + j][v]);
3688 
3689  ind_local += indicator.first;
3690  dof_indices[v] += indicator.first;
3691 
3692  // constrained case: build the local value as a linear
3693  // combination of the global value according to constraints
3694  Number value[n_components];
3695  for (unsigned int comp = 0; comp < n_components; ++comp)
3696  operation.pre_constraints(values_dofs[comp][ind_local][v],
3697  value[comp]);
3698 
3699  const Number *data_val =
3700  this->matrix_free->constraint_pool_begin(indicator.second);
3701  const Number *end_pool =
3702  this->matrix_free->constraint_pool_end(indicator.second);
3703  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3704  for (unsigned int comp = 0; comp < n_components; ++comp)
3705  operation.process_constraint(*dof_indices[v],
3706  *data_val,
3707  *src[comp],
3708  value[comp]);
3709 
3710  for (unsigned int comp = 0; comp < n_components; ++comp)
3711  operation.post_constraints(value[comp],
3712  values_dofs[comp][ind_local][v]);
3713  ind_local++;
3714  }
3715 
3716  AssertIndexRange(ind_local, dofs_per_component + 1);
3717 
3718  for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3719  for (unsigned int comp = 0; comp < n_components; ++comp)
3720  operation.process_dof(*dof_indices[v],
3721  *src[comp],
3722  values_dofs[comp][ind_local][v]);
3723  }
3724  else
3725  {
3726  // case with vector-valued finite elements where all components are
3727  // included in one single vector. Assumption: first come all entries
3728  // to the first component, then all entries to the second one, and
3729  // so on. This is ensured by the way MatrixFree reads out the
3730  // indices.
3731  for (unsigned int comp = 0; comp < n_components; ++comp)
3732  {
3733  unsigned int ind_local = 0;
3734 
3735  // check whether there is any constraint on the current cell
3736  for (; index_indicators != next_index_indicators;
3737  ++index_indicators)
3738  {
3739  const std::pair<unsigned short, unsigned short> indicator =
3740  this->dof_info->constraint_indicator[index_indicators];
3741 
3742  // run through values up to next constraint
3743  for (unsigned int j = 0; j < indicator.first; ++j)
3744  operation.process_dof(dof_indices[v][j],
3745  *src[0],
3746  values_dofs[comp][ind_local + j][v]);
3747  ind_local += indicator.first;
3748  dof_indices[v] += indicator.first;
3749 
3750  // constrained case: build the local value as a linear
3751  // combination of the global value according to constraints
3752  Number value;
3753  operation.pre_constraints(values_dofs[comp][ind_local][v],
3754  value);
3755 
3756  const Number *data_val =
3757  this->matrix_free->constraint_pool_begin(indicator.second);
3758  const Number *end_pool =
3759  this->matrix_free->constraint_pool_end(indicator.second);
3760 
3761  for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3762  operation.process_constraint(*dof_indices[v],
3763  *data_val,
3764  *src[0],
3765  value);
3766 
3767  operation.post_constraints(value,
3768  values_dofs[comp][ind_local][v]);
3769  ind_local++;
3770  }
3771 
3772  AssertIndexRange(ind_local, dofs_per_component + 1);
3773 
3774  // get the dof values past the last constraint
3775  for (; ind_local < dofs_per_component;
3776  ++dof_indices[v], ++ind_local)
3777  {
3778  AssertIndexRange(*dof_indices[v], src[0]->size());
3779  operation.process_dof(*dof_indices[v],
3780  *src[0],
3781  values_dofs[comp][ind_local][v]);
3782  }
3783 
3784  if (apply_constraints == true && comp + 1 < n_components)
3785  next_index_indicators =
3786  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
3787  }
3788  }
3789  }
3790 }
3791 
3792 
3793 
3794 template <int dim,
3795  int n_components_,
3796  typename Number,
3797  bool is_face,
3798  typename VectorizedArrayType>
3799 template <typename VectorType, typename VectorOperation>
3800 inline void
3803  const VectorOperation & operation,
3804  const std::array<VectorType *, n_components_> &src) const
3805 {
3807 
3808  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3809  unsigned int index = this->first_selected_component * dofs_per_component;
3810  for (unsigned int comp = 0; comp < n_components; ++comp)
3811  {
3812  for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3813  {
3814  operation.process_empty(
3815  this->values_dofs[comp * dofs_per_component + i]);
3816  operation.process_dof_global(
3818  *src[0],
3819  this->values_dofs[comp * dofs_per_component + i][0]);
3820  }
3821  }
3822 }
3823 
3824 
3825 
3826 template <int dim,
3827  int n_components_,
3828  typename Number,
3829  bool is_face,
3830  typename VectorizedArrayType>
3831 template <typename VectorType, typename VectorOperation>
3832 inline void
3835  const VectorOperation & operation,
3836  const std::array<VectorType *, n_components_> &src,
3837  const std::array<
3839  n_components_> & vectors_sm,
3840  const std::bitset<VectorizedArrayType::size()> &mask) const
3841 {
3842  // This functions processes the functions read_dof_values,
3843  // distribute_local_to_global, and set_dof_values with the same code for
3844  // contiguous cell indices (DG case). The distinction between these three
3845  // cases is made by the input VectorOperation that either reads values from
3846  // a vector and puts the data into the local data field or write local data
3847  // into the vector. Certain operations are no-ops for the given use case.
3848 
3849  std::integral_constant<bool,
3851  vector_selector;
3853  is_face ? this->dof_access_index :
3855  const unsigned int n_lanes = mask.count();
3856 
3857  const std::vector<unsigned int> &dof_indices_cont =
3858  this->dof_info->dof_indices_contiguous[ind];
3859 
3860  const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3861  std::array<VectorizedArrayType *, n_components> values_dofs;
3862  for (unsigned int c = 0; c < n_components; ++c)
3863  values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3864  c * dofs_per_component;
3865 
3867 
3868  // Simple case: We have contiguous storage, so we can simply copy out the
3869  // data
3870  if ((this->dof_info->index_storage_variants[ind][this->cell] ==
3872  interleaved_contiguous &&
3873  n_lanes == VectorizedArrayType::size()) &&
3874  !(is_face &&
3875  this->dof_access_index ==
3877  this->is_interior_face == false) &&
3878  !(!is_face && !this->is_interior_face))
3879  {
3880  const unsigned int dof_index =
3881  dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3882  this->dof_info
3884  [this->first_selected_component] *
3885  VectorizedArrayType::size();
3886  if (n_components == 1 || this->n_fe_components == 1)
3887  for (unsigned int comp = 0; comp < n_components; ++comp)
3888  operation.process_dofs_vectorized(dofs_per_component,
3889  dof_index,
3890  *src[comp],
3891  values_dofs[comp],
3892  vector_selector);
3893  else
3894  operation.process_dofs_vectorized(dofs_per_component * n_components,
3895  dof_index,
3896  *src[0],
3897  values_dofs[0],
3898  vector_selector);
3899  return;
3900  }
3901 
3902  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3903  this->get_cell_or_face_ids();
3904 
3905  // More general case: Must go through the components one by one and apply
3906  // some transformations
3907  const unsigned int n_filled_lanes =
3908  this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
3909 
3910  const bool is_ecl =
3911  (this->dof_access_index ==
3913  this->is_interior_face == false) ||
3914  (!is_face && !this->is_interior_face);
3915 
3916  if (vectors_sm[0] != nullptr)
3917  {
3918  const auto compute_vector_ptrs = [&](const unsigned int comp) {
3919  std::array<typename VectorType::value_type *,
3920  VectorizedArrayType::size()>
3921  vector_ptrs = {};
3922 
3923  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3924  {
3925  if (mask[v] == false)
3926  {
3927  vector_ptrs[v] = nullptr;
3928  continue;
3929  }
3930 
3932  ExcNotImplemented());
3933  Assert(ind < this->dof_info->dof_indices_contiguous_sm.size(),
3934  ExcIndexRange(
3935  ind, 0, this->dof_info->dof_indices_contiguous_sm.size()));
3936  Assert(cells[v] <
3937  this->dof_info->dof_indices_contiguous_sm[ind].size(),
3938  ExcIndexRange(
3939  cells[v],
3940  0,
3941  this->dof_info->dof_indices_contiguous_sm[ind].size()));
3942 
3943  const auto &temp =
3944  this->dof_info->dof_indices_contiguous_sm[ind][cells[v]];
3945 
3946  if (temp.first != numbers::invalid_unsigned_int)
3947  vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3948  vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3950  [this->active_fe_index][this->first_selected_component]);
3951  else
3952  vector_ptrs[v] = nullptr;
3953  }
3954  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3955  ++v)
3956  vector_ptrs[v] = nullptr;
3957 
3958  return vector_ptrs;
3959  };
3960 
3961  if (n_filled_lanes == VectorizedArrayType::size() &&
3962  n_lanes == VectorizedArrayType::size() && !is_ecl)
3963  {
3964  if (n_components == 1 || this->n_fe_components == 1)
3965  {
3966  for (unsigned int comp = 0; comp < n_components; ++comp)
3967  {
3968  auto vector_ptrs = compute_vector_ptrs(comp);
3969  operation.process_dofs_vectorized_transpose(
3970  dofs_per_component,
3971  vector_ptrs,
3972  values_dofs[comp],
3973  vector_selector);
3974  }
3975  }
3976  else
3977  {
3978  auto vector_ptrs = compute_vector_ptrs(0);
3979  operation.process_dofs_vectorized_transpose(dofs_per_component *
3980  n_components,
3981  vector_ptrs,
3982  &values_dofs[0][0],
3983  vector_selector);
3984  }
3985  }
3986  else
3987  for (unsigned int comp = 0; comp < n_components; ++comp)
3988  {
3989  auto vector_ptrs = compute_vector_ptrs(
3990  (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3991 
3992  for (unsigned int i = 0; i < dofs_per_component; ++i)
3993  operation.process_empty(values_dofs[comp][i]);
3994 
3995  if (n_components == 1 || this->n_fe_components == 1)
3996  {
3997  for (unsigned int v = 0; v < n_filled_lanes; ++v)
3998  if (mask[v] == true)
3999  for (unsigned int i = 0; i < dofs_per_component; ++i)
4000  operation.process_dof(vector_ptrs[v][i],
4001  values_dofs[comp][i][v]);
4002  }
4003  else
4004  {
4005  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4006  if (mask[v] == true)
4007  for (unsigned int i = 0; i < dofs_per_component; ++i)
4008  operation.process_dof(
4009  vector_ptrs[v][i + comp * dofs_per_component],
4010  values_dofs[comp][i][v]);
4011  }
4012  }
4013  return;
4014  }
4015 
4016  unsigned int dof_indices[VectorizedArrayType::size()];
4017 
4018  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4019  {
4021  dof_indices[v] =
4022  dof_indices_cont[cells[v]] +
4023  this->dof_info
4025  [this->first_selected_component] *
4026  this->dof_info->dof_indices_interleave_strides[ind][cells[v]];
4027  }
4028 
4029  for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4030  dof_indices[v] = numbers::invalid_unsigned_int;
4031 
4032  // In the case with contiguous cell indices, we know that there are no
4033  // constraints and that the indices within each element are contiguous
4034  if (n_filled_lanes == VectorizedArrayType::size() &&
4035  n_lanes == VectorizedArrayType::size() && !is_ecl)
4036  {
4037  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4039  contiguous)
4040  {
4041  if (n_components == 1 || this->n_fe_components == 1)
4042  for (unsigned int comp = 0; comp < n_components; ++comp)
4043  operation.process_dofs_vectorized_transpose(dofs_per_component,
4044  dof_indices,
4045  *src[comp],
4046  values_dofs[comp],
4047  vector_selector);
4048  else
4049  operation.process_dofs_vectorized_transpose(dofs_per_component *
4050  n_components,
4051  dof_indices,
4052  *src[0],
4053  &values_dofs[0][0],
4054  vector_selector);
4055  }
4056  else if (this->dof_info->index_storage_variants[ind][this->cell] ==
4058  interleaved_contiguous_strided)
4059  {
4060  if (n_components == 1 || this->n_fe_components == 1)
4061  for (unsigned int i = 0; i < dofs_per_component; ++i)
4062  {
4063  for (unsigned int comp = 0; comp < n_components; ++comp)
4064  operation.process_dof_gather(dof_indices,
4065  *src[comp],
4066  i * VectorizedArrayType::size(),
4067  values_dofs[comp][i],
4068  vector_selector);
4069  }
4070  else
4071  for (unsigned int comp = 0; comp < n_components; ++comp)
4072  for (unsigned int i = 0; i < dofs_per_component; ++i)
4073  {
4074  operation.process_dof_gather(dof_indices,
4075  *src[0],
4076  (comp * dofs_per_component + i) *
4077  VectorizedArrayType::size(),
4078  values_dofs[comp][i],
4079  vector_selector);
4080  }
4081  }
4082  else
4083  {
4084  Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
4086  IndexStorageVariants::interleaved_contiguous_mixed_strides,
4087  ExcNotImplemented());
4088  const unsigned int *offsets =
4090  [ind][VectorizedArrayType::size() * this->cell];
4091  if (n_components == 1 || this->n_fe_components == 1)
4092  for (unsigned int i = 0; i < dofs_per_component; ++i)
4093  {
4094  for (unsigned int comp = 0; comp < n_components; ++comp)
4095  operation.process_dof_gather(dof_indices,
4096  *src[comp],
4097  0,
4098  values_dofs[comp][i],
4099  vector_selector);
4101  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4102  dof_indices[v] += offsets[v];
4103  }
4104  else
4105  for (unsigned int comp = 0; comp < n_components; ++comp)
4106  for (unsigned int i = 0; i < dofs_per_component; ++i)
4107  {
4108  operation.process_dof_gather(dof_indices,
4109  *src[0],
4110  0,
4111  values_dofs[comp][i],
4112  vector_selector);
4114  for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4115  dof_indices[v] += offsets[v];
4116  }
4117  }
4118  }
4119  else
4120  for (unsigned int comp = 0; comp < n_components; ++comp)
4121  {
4122  for (unsigned int i = 0; i < dofs_per_component; ++i)
4123  operation.process_empty(values_dofs[comp][i]);
4124  if (this->dof_info->index_storage_variants[ind][this->cell] ==
4126  contiguous)
4127  {
4128  if (n_components == 1 || this->n_fe_components == 1)
4129  {
4130  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4131  if (mask[v] == true)
4132  for (unsigned int i = 0; i < dofs_per_component; ++i)
4133  operation.process_dof(dof_indices[v] + i,
4134  *src[comp],
4135  values_dofs[comp][i][v]);
4136  }
4137  else
4138  {
4139  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4140  if (mask[v] == true)
4141  for (unsigned int i = 0; i < dofs_per_component; ++i)
4142  operation.process_dof(dof_indices[v] + i +
4143  comp * dofs_per_component,
4144  *src[0],
4145  values_dofs[comp][i][v]);
4146  }
4147  }
4148  else
4149  {
4150  const unsigned int *offsets =
4152  [ind][VectorizedArrayType::size() * this->cell];
4153  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4154  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4155  if (n_components == 1 || this->n_fe_components == 1)
4156  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4157  {
4158  if (mask[v] == true)
4159  for (unsigned int i = 0; i < dofs_per_component; ++i)
4160  operation.process_dof(dof_indices[v] + i * offsets[v],
4161  *src[comp],
4162  values_dofs[comp][i][v]);
4163  }
4164  else
4165  {
4166  for (unsigned int v = 0; v < n_filled_lanes; ++v)
4167  if (mask[v] == true)
4168  for (unsigned int i = 0; i < dofs_per_component; ++i)
4169  operation.process_dof(dof_indices[v] +
4170  (i + comp * dofs_per_component) *
4171  offsets[v],
4172  *src[0],
4173  values_dofs[comp][i][v]);
4174  }
4175  }
4176  }
4177 }
4178 
4179 namespace internal
4180 {
4181  template <typename Number,
4182  typename VectorType,
4183  typename std::enable_if<!IsBlockVector<VectorType>::value,
4184  VectorType>::type * = nullptr>
4185  decltype(std::declval<VectorType>().begin())
4186  get_beginning(VectorType &vec)
4187  {
4188  return vec.begin();
4189  }
4190 
4191  template <typename Number,
4192  typename VectorType,
4193  typename std::enable_if<IsBlockVector<VectorType>::value,
4194  VectorType>::type * = nullptr>
4195  typename VectorType::value_type *
4196  get_beginning(VectorType &)
4197  {
4198  return nullptr;
4199  }
4200 
4201  template <typename VectorType,
4202  typename std::enable_if<has_shared_vector_data<VectorType>::value,
4203  VectorType>::type * = nullptr>
4204  const std::vector<ArrayView<const typename VectorType::value_type>> *
4205  get_shared_vector_data(VectorType & vec,
4206  const bool is_valid_mode_for_sm,
4207  const unsigned int active_fe_index,
4209  {
4210  // note: no hp is supported
4211  if (is_valid_mode_for_sm &&
4212  dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4213  .size() > 0 &&
4214  active_fe_index == 0)
4215  return &vec.shared_vector_data();
4216  else
4217  return nullptr;
4218  }
4219 
4220  template <typename VectorType,
4221  typename std::enable_if<!has_shared_vector_data<VectorType>::value,
4222  VectorType>::type * = nullptr>
4223  const std::vector<ArrayView<const typename VectorType::value_type>> *
4224  get_shared_vector_data(VectorType &,
4225  const bool,
4226  const unsigned int,
4228  {
4229  return nullptr;
4230  }
4231 
4232  template <int n_components, typename VectorType>
4233  std::pair<
4234  std::array<typename internal::BlockVectorSelector<
4235  VectorType,
4236  IsBlockVector<VectorType>::value>::BaseVectorType *,
4237  n_components>,
4238  std::array<
4239  const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4240  VectorType,
4241  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4242  n_components>>
4243  get_vector_data(VectorType & src,
4244  const unsigned int first_index,
4245  const bool is_valid_mode_for_sm,
4246  const unsigned int active_fe_index,
4248  {
4249  // select between block vectors and non-block vectors. Note that the number
4250  // of components is checked in the internal data
4251  std::pair<
4252  std::array<typename internal::BlockVectorSelector<
4253  VectorType,
4254  IsBlockVector<VectorType>::value>::BaseVectorType *,
4255  n_components>,
4256  std::array<
4257  const std::vector<
4258  ArrayView<const typename internal::BlockVectorSelector<
4259  VectorType,
4260  IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4261  n_components>>
4262  src_data;
4263 
4264  for (unsigned int d = 0; d < n_components; ++d)
4265  src_data.first[d] = internal::BlockVectorSelector<
4266  VectorType,
4267  IsBlockVector<VectorType>::value>::get_vector_component(src,
4268  d +
4269  first_index);
4270 
4271  for (unsigned int d = 0; d < n_components; ++d)
4272  src_data.second[d] = get_shared_vector_data(
4273  const_cast<typename internal::BlockVectorSelector<
4274  typename std::remove_const<VectorType>::type,
4276  BaseVectorType &>(*src_data.first[d]),
4277  is_valid_mode_for_sm,
4278  active_fe_index,
4279  dof_info);
4280 
4281  return src_data;
4282  }
4283 } // namespace internal
4284 
4285 
4286 
4287 template <int dim,
4288  int n_components_,
4289  typename Number,
4290  bool is_face,
4291  typename VectorizedArrayType>
4292 inline void
4295 {
4296  if (this->dof_info == nullptr ||
4297  this->dof_info->hanging_node_constraint_masks.size() == 0)
4298  return; // nothing to do with faces
4299 
4300  constexpr unsigned int n_lanes = VectorizedArrayType::size();
4301  std::array<internal::MatrixFreeFunctions::ConstraintKinds, n_lanes>
4302  constraint_mask;
4303 
4304  bool hn_available = false;
4305 
4306  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
4307  this->get_cell_ids();
4308 
4309  for (unsigned int v = 0; v < n_lanes; ++v)
4310  {
4311  if (cells[v] == numbers::invalid_unsigned_int)
4312  {
4313  constraint_mask[v] =
4315  continue;
4316  }
4317 
4318  const unsigned int cell_index = cells[v];
4319  const unsigned int cell_dof_index =
4320  cell_index * this->n_fe_components + this->first_selected_component;
4321 
4322  const auto mask =
4323  this->dof_info->hanging_node_constraint_masks[cell_dof_index];
4324  constraint_mask[v] = mask;
4325 
4326  hn_available |=
4328  }
4329 
4330  if (hn_available == false)
4331  return; // no hanging node on cell batch -> nothing to do
4332 
4334  apply(n_components,
4335  this->data->data.front().fe_degree,
4336  this->get_shape_info(),
4337  transpose,
4338  constraint_mask,
4339  this->values_dofs);
4340 }
4341 
4342 
4343 
4344 template <int dim,
4345  int n_components_,
4346  typename Number,
4347  bool is_face,
4348  typename VectorizedArrayType>
4349 template <typename VectorType>
4350 inline void
4352  read_dof_values(const VectorType & src,
4353  const unsigned int first_index,
4354  const std::bitset<VectorizedArrayType::size()> &mask)
4355 {
4356  const auto src_data = internal::get_vector_data<n_components_>(
4357  src,
4358  first_index,
4359  this->dof_access_index ==
4361  this->active_fe_index,
4362  this->dof_info);
4363 
4365  read_write_operation(reader, src_data.first, src_data.second, mask, true);
4366 
4368 
4369 # ifdef DEBUG
4370  this->dof_values_initialized = true;
4371 # endif
4372 }
4373 
4374 
4375 
4376 template <int dim,
4377  int n_components_,
4378  typename Number,
4379  bool is_face,
4380  typename VectorizedArrayType>
4381 template <typename VectorType>
4382 inline void
4385  const unsigned int first_index,
4386  const std::bitset<VectorizedArrayType::size()> &mask)
4387 {
4388  const auto src_data = internal::get_vector_data<n_components_>(
4389  src,
4390  first_index,
4391  this->dof_access_index ==
4393  this->active_fe_index,
4394  this->dof_info);
4395 
4397  read_write_operation(reader, src_data.first, src_data.second, mask, false);
4398 
4399 # ifdef DEBUG
4400  this->dof_values_initialized = true;
4401 # endif
4402 }
4403 
4404 
4405 
4406 template <int dim,
4407  int n_components_,
4408  typename Number,
4409  bool is_face,
4410  typename VectorizedArrayType>
4411 template <typename VectorType>
4412 inline void
4415  VectorType & dst,
4416  const unsigned int first_index,
4417  const std::bitset<VectorizedArrayType::size()> &mask) const
4418 {
4419 # ifdef DEBUG
4420  Assert(this->dof_values_initialized == true,
4422 # endif
4423 
4425 
4426  const auto dst_data = internal::get_vector_data<n_components_>(
4427  dst,
4428  first_index,
4429  this->dof_access_index ==
4431  this->active_fe_index,
4432  this->dof_info);
4433 
4435  distributor;
4436  read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4437 }
4438 
4439 
4440 
4441 template <int dim,
4442  int n_components_,
4443  typename Number,
4444  bool is_face,
4445  typename VectorizedArrayType>
4446 template <typename VectorType>
4447 inline void
4450  const unsigned int first_index,
4451  const std::bitset<VectorizedArrayType::size()> &mask) const
4452 {
4453 # ifdef DEBUG
4454  Assert(this->dof_values_initialized == true,
4456 # endif
4457 
4458  const auto dst_data = internal::get_vector_data<n_components_>(
4459  dst,
4460  first_index,
4461  this->dof_access_index ==
4463  this->active_fe_index,
4464  this->dof_info);
4465 
4467  read_write_operation(setter, dst_data.first, dst_data.second, mask);
4468 }
4469 
4470 
4471 
4472 template <int dim,
4473  int n_components_,
4474  typename Number,
4475  bool is_face,
4476  typename VectorizedArrayType>
4477 template <typename VectorType>
4478 inline void
4481  VectorType & dst,
4482  const unsigned int first_index,
4483  const std::bitset<VectorizedArrayType::size()> &mask) const
4484 {
4485 # ifdef DEBUG
4486  Assert(this->dof_values_initialized == true,
4488 # endif
4489 
4490  const auto dst_data = internal::get_vector_data<n_components_>(
4491  dst,
4492  first_index,
4493  this->dof_access_index ==
4495  this->active_fe_index,
4496  this->dof_info);
4497 
4499  read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4500 }
4501 
4502 
4503 
4504 /*------------------------------ access to data fields ----------------------*/
4505 
4506 
4507 
4508 template <int dim,
4509  int n_components_,
4510  typename Number,
4511  bool is_face,
4512  typename VectorizedArrayType>
4515  get_dof_value(const unsigned int dof) const
4516 {
4518  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4520  for (unsigned int comp = 0; comp < n_components; ++comp)
4521  return_value[comp] = this->values_dofs[comp * dofs + dof];
4522  return return_value;
4523 }
4524 
4525 
4526 
4527 template <int dim,
4528  int n_components_,
4529  typename Number,
4530  bool is_face,
4531  typename VectorizedArrayType>
4534  get_value(const unsigned int q_point) const
4535 {
4536 # ifdef DEBUG
4537  Assert(this->values_quad_initialized == true,
4539 # endif
4540 
4541  AssertIndexRange(q_point, this->n_quadrature_points);
4542  const std::size_t nqp = this->n_quadrature_points;
4544  for (unsigned int comp = 0; comp < n_components; ++comp)
4545  return_value[comp] = this->values_quad[comp * nqp + q_point];
4546  return return_value;
4547 }
4548 
4549 
4550 
4551 template <int dim,
4552  int n_components_,
4553  typename Number,
4554  bool is_face,
4555  typename VectorizedArrayType>
4556 inline DEAL_II_ALWAYS_INLINE
4559  get_gradient(const unsigned int q_point) const
4560 {
4561 # ifdef DEBUG
4562  Assert(this->gradients_quad_initialized == true,
4564 # endif
4565 
4566  AssertIndexRange(q_point, this->n_quadrature_points);
4567  Assert(this->jacobian != nullptr,
4569  "update_gradients"));
4570  const std::size_t nqp = this->n_quadrature_points;
4572 
4573  // Cartesian cell
4574  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4575  {
4576  for (unsigned int d = 0; d < dim; ++d)
4577  for (unsigned int comp = 0; comp < n_components; ++comp)
4578  grad_out[comp][d] =
4579  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4580  this->jacobian[0][d][d];
4581  }
4582  // cell with general/affine Jacobian
4583  else
4584  {
4587  q_point :
4588  0];
4589  for (unsigned int comp = 0; comp < n_components; ++comp)
4590  for (unsigned int d = 0; d < dim; ++d)
4591  {
4592  grad_out[comp][d] =
4593  jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4594  for (unsigned int e = 1; e < dim; ++e)
4595  grad_out[comp][d] +=
4596  jac[d][e] *
4597  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4598  }
4599  }
4600  return grad_out;
4601 }
4602 
4603 
4604 
4605 template <int dim,
4606  int n_components_,
4607  typename Number,
4608  bool is_face,
4609  typename VectorizedArrayType>
4612  get_normal_derivative(const unsigned int q_point) const
4613 {
4614  AssertIndexRange(q_point, this->n_quadrature_points);
4615 # ifdef DEBUG
4616  Assert(this->gradients_quad_initialized == true,
4618 # endif
4619 
4620  Assert(this->normal_x_jacobian != nullptr,
4622  "update_gradients"));
4623 
4624  const std::size_t nqp = this->n_quadrature_points;
4626 
4628  for (unsigned int comp = 0; comp < n_components; ++comp)
4629  grad_out[comp] =
4630  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4631  (this->normal_x_jacobian[0][dim - 1]);
4632  else
4633  {
4634  const std::size_t index =
4635  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4636  for (unsigned int comp = 0; comp < n_components; ++comp)
4637  {
4638  grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4639  this->normal_x_jacobian[index][0];
4640  for (unsigned int d = 1; d < dim; ++d)
4641  grad_out[comp] +=
4642  this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4643  this->normal_x_jacobian[index][d];
4644  }
4645  }
4646  return grad_out;
4647 }
4648 
4649 
4650 
4651 namespace internal
4652 {
4653  // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4654  // store the lower diagonal because of symmetry
4655  template <typename VectorizedArrayType>
4656  inline void
4657  hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4658  const VectorizedArrayType *const hessians,
4659  const unsigned int,
4660  VectorizedArrayType (&tmp)[1][1])
4661  {
4662  tmp[0][0] = jac[0][0] * hessians[0];
4663  }
4664 
4665  template <typename VectorizedArrayType>
4666  inline void
4667  hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4668  const VectorizedArrayType *const hessians,
4669  const unsigned int nqp,
4670  VectorizedArrayType (&tmp)[2][2])
4671  {
4672  for (unsigned int d = 0; d < 2; ++d)
4673  {
4674  tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4675  tmp[1][d] =
4676  (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4677  }
4678  }
4679 
4680  template <typename VectorizedArrayType>
4681  inline void
4682  hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4683  const VectorizedArrayType *const hessians,
4684  const unsigned int nqp,
4685  VectorizedArrayType (&tmp)[3][3])
4686  {
4687  for (unsigned int d = 0; d < 3; ++d)
4688  {
4689  tmp[0][d] =
4690  (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4691  jac[d][2] * hessians[4 * nqp]);
4692  tmp[1][d] =
4693  (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4694  jac[d][2] * hessians[5 * nqp]);
4695  tmp[2][d] =
4696  (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4697  jac[d][2] * hessians[2 * nqp]);
4698  }
4699  }
4700 } // namespace internal
4701 
4702 
4703 
4704 template <int dim,
4705  int n_components_,
4706  typename Number,
4707  bool is_face,
4708  typename VectorizedArrayType>
4711  get_hessian(const unsigned int q_point) const
4712 {
4713 # ifdef DEBUG
4714  Assert(this->hessians_quad_initialized == true,
4716 # endif
4717  AssertIndexRange(q_point, this->n_quadrature_points);
4718 
4719  Assert(this->jacobian != nullptr,
4721  "update_hessian"));
4724  0 :
4725  q_point];
4726 
4728 
4729  const std::size_t nqp = this->n_quadrature_points;
4730  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4731 
4732  // Cartesian cell
4733  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4734  {
4735  for (unsigned int comp = 0; comp < n_components; ++comp)
4736  {
4737  for (unsigned int d = 0; d < dim; ++d)
4738  hessian_out[comp][d][d] =
4739  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4740  (jac[d][d] * jac[d][d]);
4741  switch (dim)
4742  {
4743  case 1:
4744  break;
4745  case 2:
4746  hessian_out[comp][0][1] =
4747  this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4748  (jac[0][0] * jac[1][1]);
4749  break;
4750  case 3:
4751  hessian_out[comp][0][1] =
4752  this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4753  (jac[0][0] * jac[1][1]);
4754  hessian_out[comp][0][2] =
4755  this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4756  (jac[0][0] * jac[2][2]);
4757  hessian_out[comp][1][2] =
4758  this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4759  (jac[1][1] * jac[2][2]);
4760  break;
4761  default:
4762  Assert(false, ExcNotImplemented());
4763  }
4764  for (unsigned int d = 0; d < dim; ++d)
4765  for (unsigned int e = d + 1; e < dim; ++e)
4766  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4767  }
4768  }
4769  // cell with general Jacobian, but constant within the cell
4770  else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4771  {
4772  for (unsigned int comp = 0; comp < n_components; ++comp)
4773  {
4774  VectorizedArrayType tmp[dim][dim];
4775  internal::hessian_unit_times_jac(
4776  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4777 
4778  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4779  for (unsigned int d = 0; d < dim; ++d)
4780  for (unsigned int e = d; e < dim; ++e)
4781  {
4782  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4783  for (unsigned int f = 1; f < dim; ++f)
4784  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4785  }
4786 
4787  // no J' * grad(u) part here because the Jacobian is constant
4788  // throughout the cell and hence, its derivative is zero
4789 
4790  // take symmetric part
4791  for (unsigned int d = 0; d < dim; ++d)
4792  for (unsigned int e = d + 1; e < dim; ++e)
4793  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4794  }
4795  }
4796  // cell with general Jacobian
4797  else
4798  {
4799  const auto &jac_grad = this->jacobian_gradients[q_point];
4800  for (unsigned int comp = 0; comp < n_components; ++comp)
4801  {
4802  VectorizedArrayType tmp[dim][dim];
4803  internal::hessian_unit_times_jac(
4804  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4805 
4806  // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4807  for (unsigned int d = 0; d < dim; ++d)
4808  for (unsigned int e = d; e < dim; ++e)
4809  {
4810  hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4811  for (unsigned int f = 1; f < dim; ++f)
4812  hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4813  }
4814 
4815  // add diagonal part of J' * grad(u)
4816  for (unsigned int d = 0; d < dim; ++d)
4817  for (unsigned int e = 0; e < dim; ++e)
4818  hessian_out[comp][d][d] +=
4819  jac_grad[d][e] *
4820  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4821 
4822  // add off-diagonal part of J' * grad(u)
4823  for (unsigned int d = 0, count = dim; d < dim; ++d)
4824  for (unsigned int e = d + 1; e < dim; ++e, ++count)
4825  for (unsigned int f = 0; f < dim; ++f)
4826  hessian_out[comp][d][e] +=
4827  jac_grad[count][f] *
4828  this->gradients_quad[(comp * dim + f) * nqp + q_point];
4829 
4830  // take symmetric part
4831  for (unsigned int d = 0; d < dim; ++d)
4832  for (unsigned int e = d + 1; e < dim; ++e)
4833  hessian_out[comp][e][d] = hessian_out[comp][d][e];
4834  }
4835  }
4836  return hessian_out;
4837 }
4838 
4839 
4840 
4841 template <int dim,
4842  int n_components_,
4843  typename Number,
4844  bool is_face,
4845  typename VectorizedArrayType>
4848  get_hessian_diagonal(const unsigned int q_point) const
4849 {
4850  Assert(!is_face, ExcNotImplemented());
4851 # ifdef DEBUG
4852  Assert(this->hessians_quad_initialized == true,
4854 # endif
4855  AssertIndexRange(q_point, this->n_quadrature_points);
4856 
4857  Assert(this->jacobian != nullptr, ExcNotImplemented());
4860  0 :
4861  q_point];
4862 
4863  const std::size_t nqp = this->n_quadrature_points;
4864  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4866 
4867  // Cartesian cell
4868  if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4869  {
4870  for (unsigned int comp = 0; comp < n_components; ++comp)
4871  for (unsigned int d = 0; d < dim; ++d)
4872  hessian_out[comp][d] =
4873  this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4874  (jac[d][d] * jac[d][d]);
4875  }
4876  // cell with general Jacobian, but constant within the cell
4877  else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4878  {
4879  for (unsigned int comp = 0; comp < n_components; ++comp)
4880  {
4881  // compute laplacian before the gradient because it needs to access
4882  // unscaled gradient data
4883  VectorizedArrayType tmp[dim][dim];
4884  internal::hessian_unit_times_jac(
4885  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4886 
4887  // compute only the trace part of hessian, J * tmp = J *
4888  // hess_unit(u) * J^T
4889  for (unsigned int d = 0; d < dim; ++d)
4890  {
4891  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4892  for (unsigned int f = 1; f < dim; ++f)
4893  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4894  }
4895  }
4896  }
4897  // cell with general Jacobian
4898  else
4899  {
4900  const auto &jac_grad = this->jacobian_gradients[q_point];
4901  for (unsigned int comp = 0; comp < n_components; ++comp)
4902  {
4903  // compute laplacian before the gradient because it needs to access
4904  // unscaled gradient data
4905  VectorizedArrayType tmp[dim][dim];
4906  internal::hessian_unit_times_jac(
4907  jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4908 
4909  // compute only the trace part of hessian, J * tmp = J *
4910  // hess_unit(u) * J^T
4911  for (unsigned int d = 0; d < dim; ++d)
4912  {
4913  hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4914  for (unsigned int f = 1; f < dim; ++f)
4915  hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4916  }
4917 
4918  for (unsigned int d = 0; d < dim; ++d)
4919  for (unsigned int e = 0; e < dim; ++e)
4920  hessian_out[comp][d] +=
4921  jac_grad[d][e] *
4922  this->gradients_quad[(comp * dim + e) * nqp + q_point];
4923  }
4924  }
4925  return hessian_out;
4926 }
4927 
4928 
4929 
4930 template <int dim,
4931  int n_components_,
4932  typename Number,
4933  bool is_face,
4934  typename VectorizedArrayType>
4937  get_laplacian(const unsigned int q_point) const
4938 {
4939  Assert(is_face == false, ExcNotImplemented());
4940 # ifdef DEBUG
4941  Assert(this->hessians_quad_initialized == true,
4943 # endif
4944  AssertIndexRange(q_point, this->n_quadrature_points);
4945 
4947  const auto hess_diag = get_hessian_diagonal(q_point);
4948  for (unsigned int comp = 0; comp < n_components; ++comp)
4949  {
4950  laplacian_out[comp] = hess_diag[comp][0];
4951  for (unsigned int d = 1; d < dim; ++d)
4952  laplacian_out[comp] += hess_diag[comp][d];
4953  }
4954  return laplacian_out;
4955 }
4956 
4957 
4958 
4959 template <int dim,
4960  int n_components_,
4961  typename Number,
4962  bool is_face,
4963  typename VectorizedArrayType>
4964 inline DEAL_II_ALWAYS_INLINE void
4967  const unsigned int dof)
4968 {
4969 # ifdef DEBUG
4970  this->dof_values_initialized = true;
4971 # endif
4972  const std::size_t dofs = this->data->dofs_per_component_on_cell;
4974  for (unsigned int comp = 0; comp < n_components; ++comp)
4975  this->values_dofs[comp * dofs + dof] = val_in[comp];
4976 }
4977 
4978 
4979 
4980 template <int dim,
4981  int n_components_,
4982  typename Number,
4983  bool is_face,
4984  typename VectorizedArrayType>
4985 inline DEAL_II_ALWAYS_INLINE void
4988  const unsigned int q_point)
4989 {
4990 # ifdef DEBUG
4992 # endif
4993  AssertIndexRange(q_point, this->n_quadrature_points);
4994  Assert(this->J_value != nullptr,
4996  "update_values"));
4997 # ifdef DEBUG
4998  this->values_quad_submitted = true;
4999 # endif
5000 
5001  const std::size_t nqp = this->n_quadrature_points;
5003  {
5004  const VectorizedArrayType JxW =
5005  this->J_value[0] * this->quadrature_weights[q_point];
5006  for (unsigned int comp = 0; comp < n_components; ++comp)
5007  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5008  }
5009  else
5010  {
5011  const VectorizedArrayType JxW = this->J_value[q_point];
5012  for (unsigned int comp = 0; comp < n_components; ++comp)
5013  this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5014  }
5015 }
5016 
5017 
5018 
5019 template <int dim,
5020  int n_components_,
5021  typename Number,
5022  bool is_face,
5023  typename VectorizedArrayType>
5024 inline DEAL_II_ALWAYS_INLINE void
5027  const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5028  const unsigned int q_point)
5029 {
5030 # ifdef DEBUG
5032 # endif
5033  AssertIndexRange(q_point, this->n_quadrature_points);
5034  Assert(this->J_value != nullptr,
5036  "update_gradients"));
5037  Assert(this->jacobian != nullptr,
5039  "update_gradients"));
5040 # ifdef DEBUG
5041  this->gradients_quad_submitted = true;
5042 # endif
5043 
5044  const std::size_t nqp = this->n_quadrature_points;
5045  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5046  {
5047  const VectorizedArrayType JxW =
5048  this->J_value[0] * this->quadrature_weights[q_point];
5049  for (unsigned int d = 0; d < dim; ++d)
5050  {
5051  const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
5052  for (unsigned int comp = 0; comp < n_components; ++comp)
5053  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5054  grad_in[comp][d] * factor;
5055  }
5056  }
5057  else
5058  {
5061  this->jacobian[q_point] :
5062  this->jacobian[0];
5063  const VectorizedArrayType JxW =
5065  this->J_value[q_point] :
5066  this->J_value[0] * this->quadrature_weights[q_point];
5067  for (unsigned int comp = 0; comp < n_components; ++comp)
5068  for (unsigned int d = 0; d < dim; ++d)
5069  {
5070  VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5071  for (unsigned int e = 1; e < dim; ++e)
5072  new_val += (jac[e][d] * grad_in[comp][e]);
5073  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5074  new_val * JxW;
5075  }
5076  }
5077 }
5078 
5079 
5080 
5081 template <int dim,
5082  int n_components_,
5083  typename Number,
5084  bool is_face,
5085  typename VectorizedArrayType>
5086 inline DEAL_II_ALWAYS_INLINE void
5090  const unsigned int q_point)
5091 {
5092  AssertIndexRange(q_point, this->n_quadrature_points);
5093  Assert(this->normal_x_jacobian != nullptr,
5095  "update_gradients"));
5096 # ifdef DEBUG
5097  this->gradients_quad_submitted = true;
5098 # endif
5099 
5100  const std::size_t nqp = this->n_quadrature_points;
5102  for (unsigned int comp = 0; comp < n_components; ++comp)
5103  {
5104  for (unsigned int d = 0; d < dim - 1; ++d)
5105  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5106  VectorizedArrayType();
5107  this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5108  grad_in[comp] *
5109  (this->normal_x_jacobian[0][dim - 1] * this->J_value[0] *
5110  this->quadrature_weights[q_point]);
5111  }
5112  else
5113  {
5114  const unsigned int index =
5115  this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5117  this->normal_x_jacobian[index];
5118  for (unsigned int comp = 0; comp < n_components; ++comp)
5119  {
5120  VectorizedArrayType factor = grad_in[comp] * this->J_value[index];
5122  factor = factor * this->quadrature_weights[q_point];
5123  for (unsigned int d = 0; d < dim; ++d)
5124  this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5125  factor * jac[d];
5126  }
5127  }
5128 }
5129 
5130 
5131 
5132 template <int dim,
5133  int n_components_,
5134  typename Number,
5135  bool is_face,
5136  typename VectorizedArrayType>
5137 inline DEAL_II_ALWAYS_INLINE void
5140  const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5141  hessian_in,
5142  const unsigned int q_point)
5143 {
5144 # ifdef DEBUG
5146 # endif
5147  AssertIndexRange(q_point, this->n_quadrature_points);
5148  Assert(this->J_value != nullptr,
5150  "update_hessians"));
5151  Assert(this->jacobian != nullptr,
5153  "update_hessians"));
5154 # ifdef DEBUG
5155  this->hessians_quad_submitted = true;
5156 # endif
5157 
5158  // compute hessian_unit = J^T * hessian_in(u) * J
5159  const std::size_t nqp = this->n_quadrature_points;
5160  constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5161  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5162  {
5163  const VectorizedArrayType JxW =
5164  this->J_value[0] * this->quadrature_weights[q_point];
5165 
5166  // diagonal part
5167  for (unsigned int d = 0; d < dim; ++d)
5168  {
5169  const auto jac_d = this->jacobian[0][d][d];
5170  const VectorizedArrayType factor = jac_d * jac_d * JxW;
5171  for (unsigned int comp = 0; comp < n_components; ++comp)
5172  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5173  hessian_in[comp][d][d] * factor;
5174  }
5175 
5176  // off diagonal part
5177  for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5178  for (unsigned int e = 0; e < d; ++e, ++off_dia)
5179  {
5180  const auto jac_d = this->jacobian[0][d][d];
5181  const auto jac_e = this->jacobian[0][e][e];
5182  const VectorizedArrayType factor = jac_d * jac_e * JxW;
5183  for (unsigned int comp = 0; comp < n_components; ++comp)
5184  this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5185  (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5186  }
5187  }
5188  // cell with general Jacobian, but constant within the cell
5190  {
5191  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5192  const VectorizedArrayType JxW =
5193  this->J_value[0] * this->quadrature_weights[q_point];
5194  for (unsigned int comp = 0; comp < n_components; ++comp)
5195  {
5196  // 1. tmp = hessian_in(u) * J
5197  VectorizedArrayType tmp[dim][dim];
5198  for (unsigned int i = 0; i < dim; ++i)
5199  for (unsigned int j = 0; j < dim; ++j)
5200  {
5201  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5202  for (unsigned int k = 1; k < dim; ++k)
5203  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5204  }
5205 
5206  // 2. hessian_unit = J^T * tmp
5207  VectorizedArrayType tmp2[dim][dim];
5208  for (unsigned int i = 0; i < dim; ++i)
5209  for (unsigned int j = 0; j < dim; ++j)
5210  {
5211  tmp2[i][j] = jac[0][i] * tmp[0][j];
5212  for (unsigned int k = 1; k < dim; ++k)
5213  tmp2[i][j] += jac[k][i] * tmp[k][j];
5214  }
5215 
5216  // diagonal part
5217  for (unsigned int d = 0; d < dim; ++d)
5218  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5219  tmp2[d][d] * JxW;
5220 
5221  // off diagonal part
5222  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5223  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5224  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5225  (tmp2[d][e] + tmp2[e][d]) * JxW;
5226  }
5227  }
5228  else
5229  {
5230  const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5231  const VectorizedArrayType JxW = this->J_value[q_point];
5232  const auto &jac_grad = this->jacobian_gradients[q_point];
5233  for (unsigned int comp = 0; comp < n_components; ++comp)
5234  {
5235  // 1. tmp = hessian_in(u) * J
5236  VectorizedArrayType tmp[dim][dim];
5237  for (unsigned int i = 0; i < dim; ++i)
5238  for (unsigned int j = 0; j < dim; ++j)
5239  {
5240  tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5241  for (unsigned int k = 1; k < dim; ++k)
5242  tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5243  }
5244 
5245  // 2. hessian_unit = J^T * tmp
5246  VectorizedArrayType tmp2[dim][dim];
5247  for (unsigned int i = 0; i < dim; ++i)
5248  for (unsigned int j = 0; j < dim; ++j)
5249  {
5250  tmp2[i][j] = jac[0][i] * tmp[0][j];
5251  for (unsigned int k = 1; k < dim; ++k)
5252  tmp2[i][j] += jac[k][i] * tmp[k][j];
5253  }
5254 
5255  // diagonal part
5256  for (unsigned int d = 0; d < dim; ++d)
5257  this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5258  tmp2[d][d] * JxW;
5259 
5260  // off diagonal part
5261  for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5262  for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5263  this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5264  (tmp2[d][e] + tmp2[e][d]) * JxW;
5265 
5266  // 3. gradient_unit = J' ** hessian_in
5267  for (unsigned int d = 0; d < dim; ++d)
5268  {
5269  VectorizedArrayType sum = 0;
5270  for (unsigned int e = 0; e < dim; ++e)
5271  sum += hessian_in[comp][e][e] * jac_grad[e][d];
5272  for (unsigned int e = 0, count = dim; e < dim; ++e)
5273  for (unsigned int f = e + 1; f < dim; ++f, ++count)
5274  sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5275  jac_grad[count][d];
5276  this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
5277  q_point] = sum * JxW;
5278  }
5279  }
5280  }
5281 }
5282 
5283 
5284 
5285 template <int dim,
5286  int n_components_,
5287  typename Number,
5288  bool is_face,
5289  typename VectorizedArrayType>
5292  integrate_value() const
5293 {
5294 # ifdef DEBUG
5296  Assert(this->values_quad_submitted == true,
5298 # endif
5299 
5301  const std::size_t nqp = this->n_quadrature_points;
5302  for (unsigned int q = 0; q < nqp; ++q)
5303  for (unsigned int comp = 0; comp < n_components; ++comp)
5304  return_value[comp] += this->values_quad[comp * nqp + q];
5305  return (return_value);
5306 }
5307 
5308 
5309 
5310 /*----------------------- FEEvaluationAccess --------------------------------*/
5311 
5312 
5313 template <int dim,
5314  int n_components_,
5315  typename Number,
5316  bool is_face,
5317  typename VectorizedArrayType>
5318 inline FEEvaluationAccess<dim,
5319  n_components_,
5320  Number,
5321  is_face,
5322  VectorizedArrayType>::
5323  FEEvaluationAccess(
5325  const unsigned int dof_no,
5326  const unsigned int first_selected_component,
5327  const unsigned int quad_no,
5328  const unsigned int fe_degree,
5329  const unsigned int n_q_points,
5330  const bool is_interior_face,
5331  const unsigned int active_fe_index,
5332  const unsigned int active_quad_index,
5333  const unsigned int face_type)
5335  matrix_free,
5336  dof_no,
5338  quad_no,
5339  fe_degree,
5340  n_q_points,
5344  face_type)
5345 {}
5346 
5347 
5348 
5349 template <int dim,
5350  int n_components_,
5351  typename Number,
5352  bool is_face,
5353  typename VectorizedArrayType>
5354 inline FEEvaluationAccess<dim,
5355  n_components_,
5356  Number,
5357  is_face,
5358  VectorizedArrayType>::
5359  FEEvaluationAccess(
5360  const Mapping<dim> & mapping,
5361  const FiniteElement<dim> &fe,
5362  const Quadrature<1> & quadrature,
5363  const UpdateFlags update_flags,
5364  const unsigned int first_selected_component,
5367  mapping,
5368  fe,
5369  quadrature,
5370  update_flags,
5372  other)
5373 {}
5374 
5375 
5376 
5377 template <int dim,
5378  int n_components_,
5379  typename Number,
5380  bool is_face,
5381  typename VectorizedArrayType>
5382 inline FEEvaluationAccess<dim,
5383  n_components_,
5384  Number,
5385  is_face,
5386  VectorizedArrayType>::
5387  FEEvaluationAccess(const FEEvaluationAccess<dim,
5388  n_components_,
5389  Number,
5390  is_face,
5391  VectorizedArrayType> &other)
5393  other)
5394 {}
5395 
5396 
5397 
5398 template <int dim,
5399  int n_components_,
5400  typename Number,
5401  bool is_face,
5402  typename VectorizedArrayType>
5403 inline FEEvaluationAccess<dim,
5404  n_components_,
5405  Number,
5406  is_face,
5407  VectorizedArrayType> &
5409 operator=(const FEEvaluationAccess<dim,
5410  n_components_,
5411  Number,
5412  is_face,
5413  VectorizedArrayType> &other)
5414 {
5415  this->FEEvaluationBase<dim,
5416  n_components_,
5417  Number,
5418  is_face,
5419  VectorizedArrayType>::operator=(other);
5420  return *this;
5421 }
5422 
5423 
5424 
5425 /*-------------------- FEEvaluationAccess scalar ----------------------------*/
5426 
5427 
5428 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5432  const unsigned int dof_no,
5433  const unsigned int first_selected_component,
5434  const unsigned int quad_no,
5435  const unsigned int fe_degree,
5436  const unsigned int n_q_points,
5437  const bool is_interior_face,
5438  const unsigned int active_fe_index,
5439  const unsigned int active_quad_index,
5440  const unsigned int face_type)
5442  matrix_free,
5443  dof_no,
5445  quad_no,
5446  fe_degree,
5447  n_q_points,
5451  face_type)
5452 {}
5453 
5454 
5455 
5456 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5459  const Mapping<dim> & mapping,
5460  const FiniteElement<dim> &fe,
5461  const Quadrature<1> & quadrature,
5462  const UpdateFlags update_flags,
5463  const unsigned int first_selected_component,
5466  mapping,
5467  fe,
5468  quadrature,
5469  update_flags,
5471  other)
5472 {}
5473 
5474 
5475 
5476 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5480  &other)
5482 {}
5483 
5484 
5485 
5486 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5490 {
5491  this
5492  ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5493  other);
5494  return *this;
5495 }
5496 
5497 
5498 
5499 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5500 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5502  const unsigned int dof) const
5503 {
5505  return this->values_dofs[dof];
5506 }
5507 
5508 
5509 
5510 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5511 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5513  const unsigned int q_point) const
5514 {
5515 # ifdef DEBUG
5516  Assert(this->values_quad_initialized == true,
5518 # endif
5519  AssertIndexRange(q_point, this->n_quadrature_points);
5520  return this->values_quad[q_point];
5521 }
5522 
5523 
5524 
5525 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5526 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5528  get_normal_derivative(const unsigned int q_point) const
5529 {
5530  return BaseClass::get_normal_derivative(q_point)[0];
5531 }
5532 
5533 
5534 
5535 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5538  const unsigned int q_point) const
5539 {
5540  // could use the base class gradient, but that involves too many expensive
5541  // initialization operations on tensors
5542 
5543 # ifdef DEBUG
5544  Assert(this->gradients_quad_initialized == true,
5546 # endif
5547  AssertIndexRange(q_point, this->n_quadrature_points);
5548 
5549  Assert(this->jacobian != nullptr,
5551  "update_gradients"));
5552 
5554 
5555  const std::size_t nqp = this->n_quadrature_points;
5556  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5557  {
5558  for (unsigned int d = 0; d < dim; ++d)
5559  grad_out[d] =
5560  this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5561  }
5562  // cell with general/affine Jacobian
5563  else
5564  {
5567  q_point :
5568  0];
5569  for (unsigned int d = 0; d < dim; ++d)
5570  {
5571  grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5572  for (unsigned int e = 1; e < dim; ++e)
5573  grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5574  }
5575  }
5576  return grad_out;
5577 }
5578 
5579 
5580 
5581 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5584  const unsigned int q_point) const
5585 {
5586  return BaseClass::get_hessian(q_point)[0];
5587 }
5588 
5589 
5590 
5591 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5594  get_hessian_diagonal(const unsigned int q_point) const
5595 {
5596  return BaseClass::get_hessian_diagonal(q_point)[0];
5597 }
5598 
5599 
5600 
5601 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5602 inline VectorizedArrayType
5604  const unsigned int q_point) const
5605 {
5606  return BaseClass::get_laplacian(q_point)[0];
5607 }
5608 
5609 
5610 
5611 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5612 inline void DEAL_II_ALWAYS_INLINE
5614  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5615 {
5616 # ifdef DEBUG
5617  this->dof_values_initialized = true;
5619 # endif
5620  this->values_dofs[dof] = val_in;
5621 }
5622 
5623 
5624 
5625 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5626 inline void DEAL_II_ALWAYS_INLINE
5628  const VectorizedArrayType val_in,
5629  const unsigned int q_point)
5630 {
5631 # ifdef DEBUG
5633 # endif
5634  AssertIndexRange(q_point, this->n_quadrature_points);
5635  Assert(this->J_value != nullptr,
5637  "update_value"));
5638 # ifdef DEBUG
5639  this->values_quad_submitted = true;
5640 # endif
5641 
5643  {
5644  const VectorizedArrayType JxW =
5645  this->J_value[0] * this->quadrature_weights[q_point];
5646  this->values_quad[q_point] = val_in * JxW;
5647  }
5648  else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5649  {
5650  this->values_quad[q_point] = val_in * this->J_value[q_point];
5651  }
5652 }
5653 
5654 
5655 
5656 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5657 inline DEAL_II_ALWAYS_INLINE void
5659  const Tensor<1, 1, VectorizedArrayType> val_in,
5660  const unsigned int q_point)
5661 {
5662  submit_value(val_in[0], q_point);
5663 }
5664 
5665 
5666 
5667 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5668 inline DEAL_II_ALWAYS_INLINE void
5670  submit_normal_derivative(const VectorizedArrayType grad_in,
5671  const unsigned int q_point)
5672 {
5674  grad[0] = grad_in;
5675  BaseClass::submit_normal_derivative(grad, q_point);
5676 }
5677 
5678 
5679 
5680 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5681 inline DEAL_II_ALWAYS_INLINE void
5684  const unsigned int q_point)
5685 {
5686 # ifdef DEBUG
5688 # endif
5689  AssertIndexRange(q_point, this->n_quadrature_points);
5690  Assert(this->J_value != nullptr,
5692  "update_gradients"));
5693  Assert(this->jacobian != nullptr,
5695  "update_gradients"));
5696 # ifdef DEBUG
5697  this->gradients_quad_submitted = true;
5698 # endif
5699 
5700  const std::size_t nqp = this->n_quadrature_points;
5701  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5702  {
5703  const VectorizedArrayType JxW =
5704  this->J_value[0] * this->quadrature_weights[q_point];
5705  for (unsigned int d = 0; d < dim; ++d)
5706  this->gradients_quad[d * nqp + q_point] =
5707  (grad_in[d] * this->jacobian[0][d][d] * JxW);
5708  }
5709  // general/affine cell type
5710  else
5711  {
5714  this->jacobian[q_point] :
5715  this->jacobian[0];
5716  const VectorizedArrayType JxW =
5718  this->J_value[q_point] :
5719  this->J_value[0] * this->quadrature_weights[q_point];
5720  for (unsigned int d = 0; d < dim; ++d)
5721  {
5722  VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5723  for (unsigned int e = 1; e < dim; ++e)
5724  new_val += jac[e][d] * grad_in[e];
5725  this->gradients_quad[d * nqp + q_point] = new_val * JxW;
5726  }
5727  }
5728 }
5729 
5730 
5731 
5732 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5733 inline DEAL_II_ALWAYS_INLINE void
5736  const unsigned int q_point)
5737 {
5739  hessian[0] = hessian_in;
5740  BaseClass::submit_hessian(hessian, q_point);
5741 }
5742 
5743 
5744 
5745 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5746 inline VectorizedArrayType
5748  integrate_value() const
5749 {
5750  return BaseClass::integrate_value()[0];
5751 }
5752 
5753 
5754 
5755 /*----------------- FEEvaluationAccess vector-valued ------------------------*/
5756 
5757 
5758 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5762  const unsigned int dof_no,
5763  const unsigned int first_selected_component,
5764  const unsigned int quad_no,
5765  const unsigned int fe_degree,
5766  const unsigned int n_q_points,
5767  const bool is_interior_face,
5768  const unsigned int active_fe_index,
5769  const unsigned int active_quad_index,
5770  const unsigned int face_type)
5772  matrix_free,
5773  dof_no,
5775  quad_no,
5776  fe_degree,
5777  n_q_points,
5781  face_type)
5782 {}
5783 
5784 
5785 
5786 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5789  const Mapping<dim> & mapping,
5790  const FiniteElement<dim> &fe,
5791  const Quadrature<1> & quadrature,
5792  const UpdateFlags update_flags,
5793  const unsigned int first_selected_component,
5796  mapping,
5797  fe,
5798  quadrature,
5799  update_flags,
5801  other)
5802 {}
5803 
5804 
5805 
5806 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5810  &other)
5812 {}
5813 
5814 
5815 
5816 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5820  &other)
5821 {
5823  operator=(other);
5824  return *this;
5825 }
5826 
5827 
5828 
5829 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5832  get_gradient(const unsigned int q_point) const
5833 {
5834  return BaseClass::get_gradient(q_point);
5835 }
5836 
5837 
5838 
5839 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5840 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5842  get_divergence(const unsigned int q_point) const
5843 {
5844 # ifdef DEBUG
5845  Assert(this->gradients_quad_initialized == true,
5847 # endif
5848  AssertIndexRange(q_point, this->n_quadrature_points);
5849  Assert(this->jacobian != nullptr,
5851  "update_gradients"));
5852 
5853  VectorizedArrayType divergence;
5854  const std::size_t nqp = this->n_quadrature_points;
5855 
5856  // Cartesian cell
5857  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5858  {
5859  divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
5860  for (unsigned int d = 1; d < dim; ++d)
5861  divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
5862  this->jacobian[0][d][d];
5863  }
5864  // cell with general/constant Jacobian
5865  else
5866  {
5869  this->jacobian[q_point] :
5870  this->jacobian[0];
5871  divergence = jac[0][0] * this->gradients_quad[q_point];
5872  for (unsigned int e = 1; e < dim; ++e)
5873  divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
5874  for (unsigned int d = 1; d < dim; ++d)
5875  for (unsigned int e = 0; e < dim; ++e)
5876  divergence +=
5877  jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
5878  }
5879  return divergence;
5880 }
5881 
5882 
5883 
5884 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5887  get_symmetric_gradient(const unsigned int q_point) const
5888 {
5889  // copy from generic function into dim-specialization function
5890  const auto grad = get_gradient(q_point);
5891  VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5892  VectorizedArrayType half = Number(0.5);
5893  for (unsigned int d = 0; d < dim; ++d)
5894  symmetrized[d] = grad[d][d];
5895  switch (dim)
5896  {
5897  case 1:
5898  break;
5899  case 2:
5900  symmetrized[2] = grad[0][1] + grad[1][0];
5901  symmetrized[2] *= half;
5902  break;
5903  case 3:
5904  symmetrized[3] = grad[0][1] + grad[1][0];
5905  symmetrized[3] *= half;
5906  symmetrized[4] = grad[0][2] + grad[2][0];
5907  symmetrized[4] *= half;
5908  symmetrized[5] = grad[1][2] + grad[2][1];
5909  symmetrized[5] *= half;
5910  break;
5911  default:
5912  Assert(false, ExcNotImplemented());
5913  }
5915 }
5916 
5917 
5918 
5919 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5920 inline DEAL_II_ALWAYS_INLINE
5923  const unsigned int q_point) const
5924 {
5925  // copy from generic function into dim-specialization function
5926  const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
5928  switch (dim)
5929  {
5930  case 1:
5931  Assert(false,
5932  ExcMessage(
5933  "Computing the curl in 1d is not a useful operation"));
5934  break;
5935  case 2:
5936  curl[0] = grad[1][0] - grad[0][1];
5937  break;
5938  case 3:
5939  curl[0] = grad[2][1] - grad[1][2];
5940  curl[1] = grad[0][2] - grad[2][0];
5941  curl[2] = grad[1][0] - grad[0][1];
5942  break;
5943  default:
5944  Assert(false, ExcNotImplemented());
5945  }
5946  return curl;
5947 }
5948 
5949 
5950 
5951 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5954  get_hessian_diagonal(const unsigned int q_point) const
5955 {
5956  return BaseClass::get_hessian_diagonal(q_point);
5957 }
5958 
5959 
5960 
5961 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5964  const unsigned int q_point) const
5965 {
5966 # ifdef DEBUG
5967  Assert(this->hessians_quad_initialized == true,
5969 # endif
5970  AssertIndexRange(q_point, this->n_quadrature_points);
5971  return BaseClass::get_hessian(q_point);
5972 }
5973 
5974 
5975 
5976 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5977 inline DEAL_II_ALWAYS_INLINE void
5980  const unsigned int q_point)
5981 {
5982  BaseClass::submit_gradient(grad_in, q_point);
5983 }
5984 
5985 
5986 
5987 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5988 inline DEAL_II_ALWAYS_INLINE void
5991  const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
5992  const unsigned int q_point)
5993 {
5994  BaseClass::submit_gradient(grad_in, q_point);
5995 }
5996 
5997 
5998 
5999 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6000 inline DEAL_II_ALWAYS_INLINE void
6002  submit_divergence(const VectorizedArrayType div_in,
6003  const unsigned int q_point)
6004 {
6005 # ifdef DEBUG
6007 # endif
6008  AssertIndexRange(q_point, this->n_quadrature_points);
6009  Assert(this->J_value != nullptr,
6011  "update_gradients"));
6012  Assert(this->jacobian != nullptr,
6014  "update_gradients"));
6015 # ifdef DEBUG
6016  this->gradients_quad_submitted = true;
6017 # endif
6018 
6019  const std::size_t nqp = this->n_quadrature_points;
6020  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6021  {
6022  const VectorizedArrayType fac =
6023  this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6024  for (unsigned int d = 0; d < dim; ++d)
6025  {
6026  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6027  (fac * this->jacobian[0][d][d]);
6028  for (unsigned int e = d + 1; e < dim; ++e)
6029  {
6030  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6031  VectorizedArrayType();
6032  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6033  VectorizedArrayType();
6034  }
6035  }
6036  }
6037  else
6038  {
6041  this->jacobian[q_point] :
6042  this->jacobian[0];
6043  const VectorizedArrayType fac =
6045  this->J_value[q_point] :
6046  this->J_value[0] * this->quadrature_weights[q_point]) *
6047  div_in;
6048  for (unsigned int d = 0; d < dim; ++d)
6049  {
6050  for (unsigned int e = 0; e < dim; ++e)
6051  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6052  jac[d][e] * fac;
6053  }
6054  }
6055 }
6056 
6057 
6058 
6059 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6060 inline DEAL_II_ALWAYS_INLINE void
6064  const unsigned int q_point)
6065 {
6066  // could have used base class operator, but that involves some overhead
6067  // which is inefficient. it is nice to have the symmetric tensor because
6068  // that saves some operations
6069 # ifdef DEBUG
6071 # endif
6072  AssertIndexRange(q_point, this->n_quadrature_points);
6073  Assert(this->J_value != nullptr,
6075  "update_gradients"));
6076  Assert(this->jacobian != nullptr,
6078  "update_gradients"));
6079 # ifdef DEBUG
6080  this->gradients_quad_submitted = true;
6081 # endif
6082 
6083  const std::size_t nqp = this->n_quadrature_points;
6084  if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6085  {
6086  const VectorizedArrayType JxW =
6087  this->J_value[0] * this->quadrature_weights[q_point];
6088  for (unsigned int d = 0; d < dim; ++d)
6089  this->gradients_quad[(d * dim + d) * nqp + q_point] =
6090  (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6091  for (unsigned int e = 0, counter = dim; e < dim; ++e)
6092  for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6093  {
6094  const VectorizedArrayType value =
6095  sym_grad.access_raw_entry(counter) * JxW;
6096  this->gradients_quad[(e * dim + d) * nqp + q_point] =
6097  value * this->jacobian[0][d][d];
6098  this->gradients_quad[(d * dim + e) * nqp + q_point] =
6099  value * this->jacobian[0][e][e];
6100  }
6101  }
6102  // general/affine cell type
6103  else
6104  {
6105  const VectorizedArrayType JxW =
6107  this->J_value[q_point] :
6108  this->J_value[0] * this->quadrature_weights[q_point];
6111  this->jacobian[q_point] :
6112  this->jacobian[0];
6113  VectorizedArrayType weighted[dim][dim];
6114  for (unsigned int i = 0; i < dim; ++i)
6115  weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6116  for (unsigned int i = 0, counter = dim; i < dim; ++i)
6117  for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6118  {
6119  const VectorizedArrayType value =
6120  sym_grad.access_raw_entry(counter) * JxW;
6121  weighted[i][j] = value;
6122  weighted[j][i] = value;
6123  }
6124  for (unsigned int comp = 0; comp < dim; ++comp)
6125  for (unsigned int d = 0; d < dim; ++d)
6126  {
6127  VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6128  for (unsigned int e = 1; e < dim; ++e)
6129  new_val += jac[e][d] * weighted[comp][e];
6130  this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6131  }
6132  }
6133 }
6134 
6135 
6136 
6137 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6138 inline DEAL_II_ALWAYS_INLINE void
6141  const unsigned int q_point)
6142 {
6144  switch (dim)
6145  {
6146  case 1:
6147  Assert(false,
6148  ExcMessage(
6149  "Testing by the curl in 1d is not a useful operation"));
6150  break;
6151  case 2:
6152  grad[1][0] = curl[0];
6153  grad[0][1] = -curl[0];
6154  break;
6155  case 3:
6156  grad[2][1] = curl[0];
6157  grad[1][2] = -curl[0];
6158  grad[0][2] = curl[1];
6159  grad[2][0] = -curl[1];
6160  grad[1][0] = curl[2];
6161  grad[0][1] = -curl[2];
6162  break;
6163  default:
6164  Assert(false, ExcNotImplemented());
6165  }
6166  submit_gradient(grad, q_point);
6167 }
6168 
6169 
6170 /*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6171 
6172 
6173 template <typename Number, bool is_face, typename VectorizedArrayType>
6176  const MatrixFree<1, Number, VectorizedArrayType> &matrix_free,
6177  const unsigned int dof_no,
6178  const unsigned int first_selected_component,
6179  const unsigned int quad_no,
6180  const unsigned int fe_degree,
6181  const unsigned int n_q_points,
6182  const bool is_interior_face,
6183  const unsigned int active_fe_index,
6184  const unsigned int active_quad_index,
6185  const unsigned int face_type)
6187  matrix_free,
6188  dof_no,
6190  quad_no,
6191  fe_degree,
6192  n_q_points,
6196  face_type)
6197 {}
6198 
6199 
6200 
6201 template <typename Number, bool is_face, typename VectorizedArrayType>
6204  const Mapping<1> & mapping,
6205  const FiniteElement<1> &fe,
6206  const Quadrature<1> & quadrature,
6207  const UpdateFlags update_flags,
6208  const unsigned int first_selected_component,
6211  mapping,
6212  fe,
6213  quadrature,
6214  update_flags,
6216  other)
6217 {}
6218 
6219 
6220 
6221 template <typename Number, bool is_face, typename VectorizedArrayType>
6226 {}
6227 
6228 
6229 
6230 template <typename Number, bool is_face, typename VectorizedArrayType>
6234 {
6236  other);
6237  return *this;
6238 }
6239 
6240 
6241 
6242 template <typename Number, bool is_face, typename VectorizedArrayType>
6243 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6245  const unsigned int dof) const
6246 {
6248  return this->values_dofs[dof];
6249 }
6250 
6251 
6252 
6253 template <typename Number, bool is_face, typename VectorizedArrayType>
6254 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6256  const unsigned int q_point) const
6257 {
6258 # ifdef DEBUG
6259  Assert(this->values_quad_initialized == true,
6261 # endif
6262  AssertIndexRange(q_point, this->n_quadrature_points);
6263  return this->values_quad[q_point];
6264 }
6265 
6266 
6267 
6268 template <typename Number, bool is_face, typename VectorizedArrayType>
6271  const unsigned int q_point) const
6272 {
6273  // could use the base class gradient, but that involves too many inefficient
6274  // initialization operations on tensors
6275 
6276 # ifdef DEBUG
6277  Assert(this->gradients_quad_initialized == true,
6279 # endif
6280  AssertIndexRange(q_point, this->n_quadrature_points);
6281 
6284  this->jacobian[q_point] :
6285  this->jacobian[0];
6286 
6288  grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6289 
6290  return grad_out;
6291 }
6292 
6293 
6294 
6295 template <typename Number, bool is_face, typename VectorizedArrayType>
6296 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6298  const unsigned int q_point) const
6299 {
6300  return get_gradient(q_point)[0];
6301 }
6302 
6303 
6304 
6305 template <typename Number, bool is_face, typename VectorizedArrayType>
6306 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6308  get_normal_derivative(const unsigned int q_point) const
6309 {
6310  return BaseClass::get_normal_derivative(q_point)[0];
6311 }
6312 
6313 
6314 
6315 template <typename Number, bool is_face, typename VectorizedArrayType>
6318  const unsigned int q_point) const
6319 {
6320  return BaseClass::get_hessian(q_point)[0];
6321 }
6322 
6323 
6324 
6325 template <typename Number, bool is_face, typename VectorizedArrayType>
6328  get_hessian_diagonal(const unsigned int q_point) const
6329 {
6330  return BaseClass::get_hessian_diagonal(q_point)[0];
6331 }
6332 
6333 
6334 
6335 template <typename Number, bool is_face, typename VectorizedArrayType>
6336 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6338  const unsigned int q_point) const
6339 {
6340  return BaseClass::get_laplacian(q_point)[0];
6341 }
6342 
6343 
6344 
6345 template <typename Number, bool is_face, typename VectorizedArrayType>
6348  submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6349 {
6350 # ifdef DEBUG
6351  this->dof_values_initialized = true;
6353 # endif
6354  this->values_dofs[dof] = val_in;
6355 }
6356 
6357 
6358 
6359 template <typename Number, bool is_face, typename VectorizedArrayType>
6360 inline DEAL_II_ALWAYS_INLINE void
6362  const VectorizedArrayType val_in,
6363  const unsigned int q_point)
6364 {
6365 # ifdef DEBUG
6367 # endif
6368  AssertIndexRange(q_point, this->n_quadrature_points);
6369 # ifdef DEBUG
6370  this->values_quad_submitted = true;
6371 # endif
6372 
6374  {
6375  const VectorizedArrayType JxW = this->J_value[q_point];
6376  this->values_quad[q_point] = val_in * JxW;
6377  }
6378  else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6379  {
6380  const VectorizedArrayType JxW =
6381  this->J_value[0] * this->quadrature_weights[q_point];
6382  this->values_quad[q_point] = val_in * JxW;
6383  }
6384 }
6385 
6386 
6387 
6388 template <typename Number, bool is_face, typename VectorizedArrayType>
6389 inline DEAL_II_ALWAYS_INLINE void
6391  const Tensor<1, 1, VectorizedArrayType> val_in,
6392  const unsigned int q_point)
6393 {
6394  submit_value(val_in[0], q_point);
6395 }
6396 
6397 
6398 
6399 template <typename Number, bool is_face, typename VectorizedArrayType>
6400 inline DEAL_II_ALWAYS_INLINE void
6402  const Tensor<1, 1, VectorizedArrayType> grad_in,
6403  const unsigned int q_point)
6404 {
6405  submit_gradient(grad_in[0], q_point);
6406 }
6407 
6408 
6409 
6410 template <typename Number, bool is_face, typename VectorizedArrayType>
6411 inline DEAL_II_ALWAYS_INLINE void
6413  const VectorizedArrayType grad_in,
6414  const unsigned int q_point)
6415 {
6416 # ifdef DEBUG
6418 # endif
6419  AssertIndexRange(q_point, this->n_quadrature_points);
6420 # ifdef DEBUG
6421  this->gradients_quad_submitted = true;
6422 # endif
6423 
6426  this->jacobian[q_point] :
6427  this->jacobian[0];
6428  const VectorizedArrayType JxW =
6430  this->J_value[q_point] :
6431  this->J_value[0] * this->quadrature_weights[q_point];
6432 
6433  this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6434 }
6435 
6436 
6437 
6438 template <typename Number, bool is_face, typename VectorizedArrayType>
6439 inline DEAL_II_ALWAYS_INLINE void
6441  const Tensor<2, 1, VectorizedArrayType> grad_in,
6442  const unsigned int q_point)
6443 {
6444  submit_gradient(grad_in[0][0], q_point);
6445 }
6446 
6447 
6448 
6449 template <typename Number, bool is_face, typename VectorizedArrayType>
6450 inline DEAL_II_ALWAYS_INLINE void
6452  submit_normal_derivative(const VectorizedArrayType grad_in,
6453  const unsigned int q_point)
6454 {
6456  grad[0] = grad_in;
6457  BaseClass::submit_normal_derivative(grad, q_point);
6458 }
6459 
6460 
6461 
6462 template <typename Number, bool is_face, typename VectorizedArrayType>
6463 inline DEAL_II_ALWAYS_INLINE void
6466  const unsigned int q_point)
6467 {
6468  BaseClass::submit_normal_derivative(grad_in, q_point);
6469 }
6470 
6471 
6472 template <typename Number, bool is_face, typename VectorizedArrayType>
6473 inline DEAL_II_ALWAYS_INLINE void
6475  const Tensor<2, 1, VectorizedArrayType> hessian_in,
6476  const unsigned int q_point)
6477 {
6479  hessian[0] = hessian_in;
6480  BaseClass::submit_hessian(hessian, q_point);
6481 }
6482 
6483 
6484 template <typename Number, bool is_face, typename VectorizedArrayType>
6485 inline VectorizedArrayType
6487  integrate_value() const
6488 {
6489  return BaseClass::integrate_value()[0];
6490 }
6491 
6492 
6493 
6494 /*-------------------------- FEEvaluation -----------------------------------*/
6495 
6496 
6497 template <int dim,
6498  int fe_degree,
6499  int n_q_points_1d,
6500  int n_components_,
6501  typename Number,
6502  typename VectorizedArrayType>
6503 inline FEEvaluation<dim,
6504  fe_degree,
6505  n_q_points_1d,
6506  n_components_,
6507  Number,
6508  VectorizedArrayType>::
6510  const unsigned int fe_no,
6511  const unsigned int quad_no,
6512  const unsigned int first_selected_component,
6513  const unsigned int active_fe_index,
6514  const unsigned int active_quad_index)
6515  : BaseClass(matrix_free,
6516  fe_no,
6517  first_selected_component,
6518  quad_no,
6519  fe_degree,
6520  static_n_q_points,
6521  true /*note: this is not a face*/,
6522  active_fe_index,
6523  active_quad_index)
6524  , dofs_per_component(this->data->dofs_per_component_on_cell)
6525  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6526  , n_q_points(this->data->n_q_points)
6527 {
6528  check_template_arguments(fe_no, 0);
6529 }
6530 
6531 
6532 
6533 template <int dim,
6534  int fe_degree,
6535  int n_q_points_1d,
6536  int n_components_,
6537  typename Number,
6538  typename VectorizedArrayType>
6539 inline FEEvaluation<dim,
6540  fe_degree,
6541  n_q_points_1d,
6542  n_components_,
6543  Number,
6544  VectorizedArrayType>::
6546  const std::pair<unsigned int, unsigned int> & range,
6547  const unsigned int dof_no,
6548  const unsigned int quad_no,
6549  const unsigned int first_selected_component)
6550  : FEEvaluation(matrix_free,
6551  dof_no,
6552  quad_no,
6553  first_selected_component,
6554  matrix_free.get_cell_active_fe_index(range))
6555 {}
6556 
6557 
6558 
6559 template <int dim,
6560  int fe_degree,
6561  int n_q_points_1d,
6562  int n_components_,
6563  typename Number,
6564  typename VectorizedArrayType>
6565 inline FEEvaluation<dim,
6566  fe_degree,
6567  n_q_points_1d,
6568  n_components_,
6569  Number,
6570  VectorizedArrayType>::
6571  FEEvaluation(const Mapping<dim> & mapping,
6572  const FiniteElement<dim> &fe,
6573  const Quadrature<1> & quadrature,
6574  const UpdateFlags update_flags,
6575  const unsigned int first_selected_component)
6576  : BaseClass(mapping,
6577  fe,
6578  quadrature,
6579  update_flags,
6580  first_selected_component,
6581  nullptr)
6582  , dofs_per_component(this->data->dofs_per_component_on_cell)
6583  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6584  , n_q_points(this->data->n_q_points)
6585 {
6586  check_template_arguments(numbers::invalid_unsigned_int, 0);
6587 }
6588 
6589 
6590 
6591 template <int dim,
6592  int fe_degree,
6593  int n_q_points_1d,
6594  int n_components_,
6595  typename Number,
6596  typename VectorizedArrayType>
6597 inline FEEvaluation<dim,
6598  fe_degree,
6599  n_q_points_1d,
6600  n_components_,
6601  Number,
6602  VectorizedArrayType>::
6604  const Quadrature<1> & quadrature,
6605  const UpdateFlags update_flags,
6606  const unsigned int first_selected_component)
6607  : BaseClass(StaticMappingQ1<dim>::mapping,
6608  fe,
6609  quadrature,
6610  update_flags,
6611  first_selected_component,
6612  nullptr)
6613  , dofs_per_component(this->data->dofs_per_component_on_cell)
6614  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6615  , n_q_points(this->data->n_q_points)
6616 {
6617  check_template_arguments(numbers::invalid_unsigned_int, 0);
6618 }
6619 
6620 
6621 
6622 template <int dim,
6623  int fe_degree,
6624  int n_q_points_1d,
6625  int n_components_,
6626  typename Number,
6627  typename VectorizedArrayType>
6628 inline FEEvaluation<dim,
6629  fe_degree,
6630  n_q_points_1d,
6631  n_components_,
6632  Number,
6633  VectorizedArrayType>::
6636  const unsigned int first_selected_component)
6637  : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6638  fe,
6639  other.mapped_geometry->get_quadrature(),
6640  other.mapped_geometry->get_fe_values().get_update_flags(),
6642  &other)
6643  , dofs_per_component(this->data->dofs_per_component_on_cell)
6644  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6645  , n_q_points(this->data->n_q_points)
6646 {
6647  check_template_arguments(numbers::invalid_unsigned_int, 0);
6648 }
6649 
6650 
6651 
6652 template <int dim,
6653  int fe_degree,
6654  int n_q_points_1d,
6655  int n_components_,
6656  typename Number,
6657  typename VectorizedArrayType>
6658 inline FEEvaluation<dim,
6659  fe_degree,
6660  n_q_points_1d,
6661  n_components_,
6662  Number,
6663  VectorizedArrayType>::FEEvaluation(const FEEvaluation
6664  &other)
6665  : BaseClass(other)
6666  , dofs_per_component(this->data->dofs_per_component_on_cell)
6667  , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6668  , n_q_points(this->data->n_q_points)
6669 {
6670  check_template_arguments(numbers::invalid_unsigned_int, 0);
6671 }
6672 
6673 
6674 
6675 template <int dim,
6676  int fe_degree,
6677  int n_q_points_1d,
6678  int n_components_,
6679  typename Number,
6680  typename VectorizedArrayType>
6681 inline FEEvaluation<dim,
6682  fe_degree,
6683  n_q_points_1d,
6684  n_components_,
6685  Number,
6686  VectorizedArrayType> &
6687 FEEvaluation<dim,
6688  fe_degree,
6689  n_q_points_1d,
6690  n_components_,
6691  Number,
6692  VectorizedArrayType>::operator=(const FEEvaluation &other)
6693 {
6694  BaseClass::operator=(other);
6695  check_template_arguments(numbers::invalid_unsigned_int, 0);
6696  return *this;
6697 }
6698 
6699 
6700 
6701 template <int dim,
6702  int fe_degree,
6703  int n_q_points_1d,
6704  int n_components_,
6705  typename Number,
6706  typename VectorizedArrayType>
6707 inline void
6708 FEEvaluation<dim,
6709  fe_degree,
6710  n_q_points_1d,
6711  n_components_,
6712  Number,
6713  VectorizedArrayType>::
6714  check_template_arguments(const unsigned int dof_no,
6715  const unsigned int first_selected_component)
6716 {
6717  (void)dof_no;
6718  (void)first_selected_component;
6719 
6720 # ifdef DEBUG
6721  // print error message when the dimensions do not match. Propose a possible
6722  // fix
6723  if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
6724  static_cast<unsigned int>(fe_degree) !=
6725  this->data->data.front().fe_degree) ||
6726  n_q_points != this->n_quadrature_points)
6727  {
6728  std::string message =
6729  "-------------------------------------------------------\n";
6730  message += "Illegal arguments in constructor/wrong template arguments!\n";
6731  message += " Called --> FEEvaluation<dim,";
6732  message += Utilities::int_to_string(fe_degree) + ",";
6733  message += Utilities::int_to_string(n_q_points_1d);
6734  message += "," + Utilities::int_to_string(n_components);
6735  message += ",Number>(data";
6736  if (first_selected_component != numbers::invalid_unsigned_int)
6737  {
6738  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6739  message += Utilities::int_to_string(this->quad_no) + ", ";
6740  message += Utilities::int_to_string(first_selected_component);
6741  }
6742  message += ")\n";
6743 
6744  // check whether some other vector component has the correct number of
6745  // points
6746  unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
6747  proposed_fe_comp = numbers::invalid_unsigned_int,
6748  proposed_quad_comp = numbers::invalid_unsigned_int;
6749  if (dof_no != numbers::invalid_unsigned_int)
6750  {
6751  if (static_cast<unsigned int>(fe_degree) ==
6752  this->data->data.front().fe_degree)
6753  {
6754  proposed_dof_comp = dof_no;
6755  proposed_fe_comp = first_selected_component;
6756  }
6757  else
6758  for (unsigned int no = 0; no < this->matrix_free->n_components();
6759  ++no)
6760  for (unsigned int nf = 0;
6761  nf < this->matrix_free->n_base_elements(no);
6762  ++nf)
6763  if (this->matrix_free
6764  ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
6765  .data.front()
6766  .fe_degree == static_cast<unsigned int>(fe_degree))
6767  {
6768  proposed_dof_comp = no;
6769  proposed_fe_comp = nf;
6770  break;
6771  }
6772  if (n_q_points ==
6773  this->mapping_data->descriptor[this->active_quad_index]
6774  .n_q_points)
6775  proposed_quad_comp = this->quad_no;
6776  else
6777  for (unsigned int no = 0;
6778  no < this->matrix_free->get_mapping_info().cell_data.size();
6779  ++no)
6780  if (this->matrix_free->get_mapping_info()
6781  .cell_data[no]
6782  .descriptor[this->active_quad_index]
6783  .n_q_points == n_q_points)
6784  {
6785  proposed_quad_comp = no;
6786  break;
6787  }
6788  }
6789  if (proposed_dof_comp != numbers::invalid_unsigned_int &&
6790  proposed_quad_comp != numbers::invalid_unsigned_int)
6791  {
6792  if (proposed_dof_comp != first_selected_component)
6793  message += "Wrong vector component selection:\n";
6794  else
6795  message += "Wrong quadrature formula selection:\n";
6796  message += " Did you mean FEEvaluation<dim,";
6797  message += Utilities::int_to_string(fe_degree) + ",";
6798  message += Utilities::int_to_string(n_q_points_1d);
6799  message += "," + Utilities::int_to_string(n_components);
6800  message += ",Number>(data";
6801  if (dof_no != numbers::invalid_unsigned_int)
6802  {
6803  message +=
6804  ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
6805  message += Utilities::int_to_string(proposed_quad_comp) + ", ";
6806  message += Utilities::int_to_string(proposed_fe_comp);
6807  }
6808  message += ")?\n";
6809  std::string correct_pos;
6810  if (proposed_dof_comp != dof_no)
6811  correct_pos = " ^ ";
6812  else
6813  correct_pos = " ";
6814  if (proposed_quad_comp != this->quad_no)
6815  correct_pos += " ^ ";
6816  else
6817  correct_pos += " ";
6818  if (proposed_fe_comp != first_selected_component)
6819  correct_pos += " ^\n";
6820  else
6821  correct_pos += " \n";
6822  message += " " +
6823  correct_pos;
6824  }
6825  // ok, did not find the numbers specified by the template arguments in
6826  // the given list. Suggest correct template arguments
6827  const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
6828  std::pow(1.001 * this->n_quadrature_points, 1. / dim));
6829  message += "Wrong template arguments:\n";
6830  message += " Did you mean FEEvaluation<dim,";
6831  message +=
6832  Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
6833  message += Utilities::int_to_string(proposed_n_q_points_1d);
6834  message += "," + Utilities::int_to_string(n_components);
6835  message += ",Number>(data";
6836  if (dof_no != numbers::invalid_unsigned_int)
6837  {
6838  message += ", " + Utilities::int_to_string(dof_no) + ", ";
6839  message += Utilities::int_to_string(this->quad_no);
6840  message += ", " + Utilities::int_to_string(first_selected_component);
6841  }
6842  message += ")?\n";
6843  std::string correct_pos;
6844  if (this->data->data.front().fe_degree !=
6845  static_cast<unsigned int>(fe_degree))
6846  correct_pos = " ^";
6847  else
6848  correct_pos = " ";
6849  if (proposed_n_q_points_1d != n_q_points_1d)
6850  correct_pos += " ^\n";
6851  else
6852  correct_pos += " \n";
6853  message += " " + correct_pos;
6854 
6855  Assert(static_cast<unsigned int>(fe_degree) ==
6856  this->data->data.front().fe_degree &&
6857  n_q_points == this->n_quadrature_points,
6858  ExcMessage(message));
6859  }
6860  if (dof_no != numbers::invalid_unsigned_int)
6862  n_q_points,
6863  this->mapping_data->descriptor[this->active_quad_index].n_q_points);
6864 # endif
6865 }
6866 
6867 
6868 
6869 template <int dim,
6870  int fe_degree,
6871  int n_q_points_1d,
6872  int n_components_,
6873  typename Number,
6874  typename VectorizedArrayType>
6875 inline void
6876 FEEvaluation<dim,
6877  fe_degree,
6878  n_q_points_1d,
6879  n_components_,
6880  Number,
6881  VectorizedArrayType>::reinit(const unsigned int cell_index)
6882 {
6883  Assert(this->mapped_geometry == nullptr,
6884  ExcMessage("FEEvaluation was initialized without a matrix-free object."
6885  " Integer indexing is not possible"));
6886  if (this->mapped_geometry != nullptr)
6887  return;
6888 
6889  Assert(this->dof_info != nullptr, ExcNotInitialized());
6890  Assert(this->mapping_data != nullptr, ExcNotInitialized());
6891  this->cell = cell_index;
6892  this->cell_type =
6893  this->matrix_free->get_mapping_info().get_cell_type(cell_index);
6894 
6895  const unsigned int offsets =
6897  this->jacobian = &this->mapping_data->jacobians[0][offsets];
6898  this->J_value = &this->mapping_data->JxW_values[offsets];
6899  this->jacobian_gradients =
6900  this->mapping_data->jacobian_gradients[0].data() + offsets;
6901 
6902  unsigned int i = 0;
6903  for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
6904  ++i)
6905  this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
6906  for (; i < VectorizedArrayType::size(); ++i)
6908 
6909  if (this->mapping_data->quadrature_points.empty() == false)
6910  this->quadrature_points =
6912  [this->mapping_data->quadrature_point_offsets[this->cell]];
6913 
6914 # ifdef DEBUG
6915  this->is_reinitialized = true;
6916  this->dof_values_initialized = false;
6917  this->values_quad_initialized = false;
6918  this->gradients_quad_initialized = false;
6919  this->hessians_quad_initialized = false;
6920 # endif
6921 }
6922 
6923 
6924 
6925 template <int dim,
6926  int fe_degree,
6927  int n_q_points_1d,
6928  int n_components_,
6929  typename Number,
6930  typename VectorizedArrayType>
6931 template <bool level_dof_access>
6932 inline void
6933 FEEvaluation<dim,
6934  fe_degree,
6935  n_q_points_1d,
6936  n_components_,
6937  Number,
6938  VectorizedArrayType>::
6940 {
6941  Assert(this->matrix_free == nullptr,
6942  ExcMessage("Cannot use initialization from cell iterator if "
6943  "initialized from MatrixFree object. Use variant for "
6944  "on the fly computation with arguments as for FEValues "
6945  "instead"));
6946  Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
6947  this->mapped_geometry->reinit(
6948  static_cast<typename Triangulation<dim>::cell_iterator>(cell));
6949  this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
6950  if (level_dof_access)
6951  cell->get_mg_dof_indices(this->local_dof_indices);
6952  else
6953  cell->get_dof_indices(this->local_dof_indices);
6954 
6955 # ifdef DEBUG
6956  this->is_reinitialized = true;
6957 # endif
6958 }
6959 
6960 
6961 
6962 template <int dim,
6963  int fe_degree,
6964  int n_q_points_1d,
6965  int n_components_,
6966  typename Number,
6967  typename VectorizedArrayType>
6968 inline void
6969 FEEvaluation<dim,
6970  fe_degree,
6971  n_q_points_1d,
6972  n_components_,
6973  Number,
6974  VectorizedArrayType>::
6976 {
6977  Assert(this->matrix_free == 0,
6978  ExcMessage("Cannot use initialization from cell iterator if "
6979  "initialized from MatrixFree object. Use variant for "
6980  "on the fly computation with arguments as for FEValues "
6981  "instead"));
6982  Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
6983  this->mapped_geometry->reinit(cell);
6984 
6985 # ifdef DEBUG
6986  this->is_reinitialized = true;
6987 # endif
6988 }
6989 
6990 
6991 
6992 template <int dim,
6993  int fe_degree,
6994  int n_q_points_1d,
6995  int n_components_,
6996  typename Number,
6997  typename VectorizedArrayType>
6998 inline void
6999 FEEvaluation<dim,
7000  fe_degree,
7001  n_q_points_1d,
7002  n_components_,
7003  Number,
7004  VectorizedArrayType>::evaluate(const bool evaluate_values,
7005  const bool evaluate_gradients,
7006  const bool evaluate_hessians)
7007 {
7008 # ifdef DEBUG
7009  Assert(this->dof_values_initialized == true,
7011 # endif
7012  evaluate(this->values_dofs,
7013  evaluate_values,
7014  evaluate_gradients,
7015  evaluate_hessians);
7016 }
7017 
7018 
7019 template <int dim,
7020  int fe_degree,
7021  int n_q_points_1d,
7022  int n_components_,
7023  typename Number,
7024  typename VectorizedArrayType>
7025 inline void
7026 FEEvaluation<dim,
7027  fe_degree,
7028  n_q_points_1d,
7029  n_components_,
7030  Number,
7031  VectorizedArrayType>::
7032  evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7033 {
7034 # ifdef DEBUG
7035  Assert(this->dof_values_initialized == true,
7037 # endif
7038  evaluate(this->values_dofs, evaluation_flags);
7039 }
7040 
7041 
7042 
7043 template <int dim,
7044  int fe_degree,
7045  int n_q_points_1d,
7046  int n_components_,
7047  typename Number,
7048  typename VectorizedArrayType>
7049 inline void
7050 FEEvaluation<dim,
7051  fe_degree,
7052  n_q_points_1d,
7053  n_components_,
7054  Number,
7055  VectorizedArrayType>::evaluate(const VectorizedArrayType
7056  * values_array,
7057  const bool evaluate_values,
7058  const bool evaluate_gradients,
7059  const bool evaluate_hessians)
7060 {
7062  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7063  ((evaluate_gradients) ? EvaluationFlags::gradients :
7064  EvaluationFlags::nothing) |
7065  ((evaluate_hessians) ? EvaluationFlags::hessians :
7066  EvaluationFlags::nothing);
7067 
7068  evaluate(values_array, flag);
7069 }
7070 
7071 
7072 
7073 template <int dim,
7074  int fe_degree,
7075  int n_q_points_1d,
7076  int n_components_,
7077  typename Number,
7078  typename VectorizedArrayType>
7079 inline void
7080 FEEvaluation<dim,
7081  fe_degree,
7082  n_q_points_1d,
7083  n_components_,
7084  Number,
7085  VectorizedArrayType>::
7086  evaluate(const VectorizedArrayType * values_array,
7087  const EvaluationFlags::EvaluationFlags evaluation_flag)
7088 {
7089  const bool hessians_on_general_cells =
7090  evaluation_flag & EvaluationFlags::hessians &&
7092  EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7093  if (hessians_on_general_cells)
7094  evaluation_flag_actual |= EvaluationFlags::gradients;
7095 
7096  if (fe_degree > -1)
7097  {
7099  evaluate(n_components, evaluation_flag_actual, values_array, *this);
7100  }
7101  else
7102  {
7104  n_components,
7105  evaluation_flag_actual,
7106  const_cast<VectorizedArrayType *>(values_array),
7107  *this);
7108  }
7109 
7110 # ifdef DEBUG
7111  if ((evaluation_flag_actual & EvaluationFlags::values) != 0u)
7112  this->values_quad_initialized = true;
7113  if ((evaluation_flag_actual & EvaluationFlags::gradients) != 0u)
7114  this->gradients_quad_initialized = true;
7115  if ((evaluation_flag_actual & EvaluationFlags::hessians) != 0u)
7116  this->hessians_quad_initialized = true;
7117 # endif
7118 }
7119 
7120 
7121 
7122 template <int dim,
7123  int fe_degree,
7124  int n_q_points_1d,
7125  int n_components_,
7126  typename Number,
7127  typename VectorizedArrayType>
7128 template <typename VectorType>
7129 inline void
7130 FEEvaluation<
7131  dim,
7132  fe_degree,
7133  n_q_points_1d,
7134  n_components_,
7135  Number,
7136  VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7137  const bool evaluate_values,
7138  const bool evaluate_gradients,
7139  const bool evaluate_hessians)
7140 {
7142  ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7143  ((evaluate_gradients) ? EvaluationFlags::gradients :
7144  EvaluationFlags::nothing) |
7145  ((evaluate_hessians) ? EvaluationFlags::hessians :
7146  EvaluationFlags::nothing);
7147 
7148  gather_evaluate(input_vector, flag);
7149 }
7150 
7151 
7152 namespace internal
7153 {
7157  template <typename Number,
7158  typename VectorizedArrayType,
7159  typename VectorType,
7160  typename EvaluatorType,
7161  typename std::enable_if<
7163  std::is_same<decltype(std::declval<VectorType>().begin()),
7164  Number *>::value,
7165  VectorType>::type * = nullptr>
7166  VectorizedArrayType *
7167  check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
7168  {
7169  const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
7170  const auto & dof_info = fe_eval.get_dof_info();
7171 
7172  // If the index storage is interleaved and contiguous and the vector
7173  // storage has the correct alignment, we can directly pass the pointer
7174  // into the vector to the evaluate() and integrate() calls, without
7175  // reading the vector entries into a separate data field. This saves some
7176  // operations.
7177  if (std::is_same<typename VectorType::value_type, Number>::value &&
7181  interleaved_contiguous &&
7182  reinterpret_cast<std::size_t>(
7183  vector.begin() +
7186  [cell * VectorizedArrayType::size()]) %
7187  sizeof(VectorizedArrayType) ==
7188  0)
7189  {
7190  return reinterpret_cast<VectorizedArrayType *>(
7191  vector.begin() +
7194  [cell * VectorizedArrayType::size()] +
7196  [fe_eval.get_active_fe_index()]
7197  [fe_eval.get_first_selected_component()] *
7198  VectorizedArrayType::size());
7199  }
7200  else
7201  return nullptr;
7202  }
7203 
7207  template <typename Number,
7208  typename VectorizedArrayType,
7209  typename VectorType,
7210  typename EvaluatorType,
7211  typename std::enable_if<
7213  !std::is_same<decltype(std::declval<VectorType>().begin()),
7214  Number *>::value,
7215  VectorType>::type * = nullptr>
7216  VectorizedArrayType *
7217  check_vector_access_inplace(const EvaluatorType &, VectorType &)
7218  {
7219  return nullptr;
7220  }
7221 } // namespace internal
7222 
7223 
7224 
7225 template <int dim,
7226  int fe_degree,
7227  int n_q_points_1d,
7228  int n_components_,
7229  typename Number,
7230  typename VectorizedArrayType>
7231 template <typename VectorType>
7232 inline void
7233 FEEvaluation<dim,
7234  fe_degree,
7235  n_q_points_1d,
7236  n_components_,
7237  Number,
7238  VectorizedArrayType>::
7239  gather_evaluate(const VectorType & input_vector,
7240  const EvaluationFlags::EvaluationFlags evaluation_flag)
7241 {
7242  const VectorizedArrayType *src_ptr =
7243  internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
7244  *this, input_vector);
7245  if (src_ptr != nullptr)
7246  evaluate(src_ptr, evaluation_flag);
7247  else
7248  {
7249  this->read_dof_values(input_vector);
7250  evaluate(this->begin_dof_values(), evaluation_flag);
7251  }
7252 }
7253 
7254 
7255 
7256 template <int dim,
7257  int fe_degree,
7258  int n_q_points_1d,
7259  int n_components_,
7260  typename Number,
7261  typename VectorizedArrayType>
7262 inline void
7263 FEEvaluation<dim,
7264  fe_degree,
7265  n_q_points_1d,
7266  n_components_,
7267  Number,
7268  VectorizedArrayType>::integrate(const bool integrate_values,
7269  const bool integrate_gradients)
7270 {
7271  integrate(integrate_values, integrate_gradients, this->values_dofs);
7272 
7273 # ifdef DEBUG
7274  this->dof_values_initialized = true;
7275 # endif
7276 }
7277 
7278 
7279 
7280 template <int dim,
7281  int fe_degree,
7282  int n_q_points_1d,
7283  int n_components_,
7284  typename Number,
7285  typename VectorizedArrayType>
7286 inline void
7287 FEEvaluation<dim,
7288  fe_degree,
7289  n_q_points_1d,
7290  n_components_,
7291  Number,
7292  VectorizedArrayType>::
7293  integrate(const EvaluationFlags::EvaluationFlags integration_flag)
7294 {
7295  integrate(integration_flag, this->values_dofs);
7296 
7297 # ifdef DEBUG
7298  this->dof_values_initialized = true;
7299 # endif
7300 }
7301 
7302 
7303 
7304 template <int dim,
7305  int fe_degree,
7306  int n_q_points_1d,
7307  int n_components_,
7308  typename Number,
7309  typename VectorizedArrayType>
7310 inline void
7311 FEEvaluation<dim,
7312  fe_degree,
7313  n_q_points_1d,
7314  n_components_,
7315  Number,
7316  VectorizedArrayType>::integrate(const bool integrate_values,
7317  const bool integrate_gradients,
7318  VectorizedArrayType *values_array)
7319 {
7321  (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
7322  (integrate_gradients ? EvaluationFlags::gradients :
7324  integrate(flag, values_array);
7325 }
7326 
7327 
7328 
7329 template <int dim,
7330  int fe_degree,
7331  int n_q_points_1d,
7332  int n_components_,
7333  typename Number,
7334  typename VectorizedArrayType>
7335 inline void
7336 FEEvaluation<dim,
7337  fe_degree,
7338  n_q_points_1d,
7339  n_components_,
7340  Number,
7341  VectorizedArrayType>::
7342  integrate(const EvaluationFlags::EvaluationFlags integration_flag,
7343  VectorizedArrayType * values_array,
7344  const bool sum_into_values_array)
7345 {
7346 # ifdef DEBUG
7347  if ((integration_flag & EvaluationFlags::values) != 0u)
7348  Assert(this->values_quad_submitted == true,
7350  if ((integration_flag & EvaluationFlags::gradients) != 0u)
7351  Assert(this->gradients_quad_submitted == true,
7353  if ((integration_flag & EvaluationFlags::hessians) != 0u)
7354  Assert(this->hessians_quad_submitted == true,
7356 # endif
7357  Assert(this->matrix_free != nullptr ||
7358  this->mapped_geometry->is_initialized(),
7359  ExcNotInitialized());
7360 
7361  Assert(
7362  (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
7363  EvaluationFlags::hessians)) == 0,
7364  ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
7365  "EvaluationFlags::hessians are supported."));
7366 
7367  EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
7368  if (integration_flag & EvaluationFlags::hessians &&
7370  {
7371  unsigned int size = n_components * dim * n_q_points;
7372  if ((integration_flag & EvaluationFlags::gradients) != 0u)
7373  {
7374  for (unsigned int i = 0; i < size; ++i)
7375  this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
7376  }
7377  else
7378  {
7379  for (unsigned int i = 0; i < size; ++i)
7380  this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
7381  integration_flag_actual |= EvaluationFlags::gradients;
7382  }
7383  }
7384 
7385  if (fe_degree > -1)
7386  {
7388  integrate(n_components,
7389  integration_flag_actual,
7390  values_array,
7391  *this,
7392  sum_into_values_array);
7393  }
7394  else
7395  {
7397  n_components,
7398  integration_flag_actual,
7399  values_array,
7400  *this,
7401  sum_into_values_array);
7402  }
7403 
7404 # ifdef DEBUG
7405  this->dof_values_initialized = true;
7406 # endif
7407 }
7408 
7409 
7410 
7411 template <int dim,
7412  int fe_degree,
7413  int n_q_points_1d,
7414  int n_components_,
7415  typename Number,
7416  typename VectorizedArrayType>
7417 template <typename VectorType>
7418 inline void
7419 FEEvaluation<
7420  dim,
7421  fe_degree,
7422  n_q_points_1d,
7423  n_components_,
7424  Number,
7425  VectorizedArrayType>::integrate_scatter(const bool integrate_values,
7426  const bool integrate_gradients,
7427  VectorType &destination)
7428 {
7430  ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7431  ((integrate_gradients) ? EvaluationFlags::gradients :
7432  EvaluationFlags::nothing);
7433 
7434  integrate_scatter(flag, destination);
7435 }
7436 
7437