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