Reference documentation for deal.II version 9.4.1
\(\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\}}\)
Loading...
Searching...
No Matches
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 " +
2910 std::to_string(
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 " +
2916 std::to_string(
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
5750 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5751 const VectorizedArrayType inv_det =
5752 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5753 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5754 this->jacobian[0][2][2];
5755
5756 // J * u * det(J^-1)
5757 for (unsigned int comp = 0; comp < n_components; ++comp)
5758 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5759 jac[comp][comp] * inv_det;
5760 }
5761 else
5762 {
5763 // Affine or general cell
5764 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5765 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5766 this->jacobian[q_point] :
5767 this->jacobian[0];
5769 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5770 transpose(invert(inv_t_jac)) :
5771 this->jacobian[1];
5772
5773 // Derivatives are reordered for faces. Need to take this into account
5774 const VectorizedArrayType inv_det =
5775 (is_face && dim == 2 && this->get_face_no() < 2) ?
5776 -determinant(inv_t_jac) :
5777 determinant(inv_t_jac);
5778 // J * u * det(J^-1)
5779 for (unsigned int comp = 0; comp < n_components; ++comp)
5780 {
5781 value_out[comp] =
5782 this->values_quad[q_point] * jac[comp][0] * inv_det;
5783 for (unsigned int e = 1; e < dim; ++e)
5784 value_out[comp] +=
5785 this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
5786 }
5787 }
5788 return value_out;
5789 }
5790 else
5791 {
5792 // No Piola needed
5793 return BaseClass::get_value(q_point);
5794 }
5795}
5796
5797template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5800 get_gradient(const unsigned int q_point) const
5801{
5802 if (this->data->element_type ==
5804 {
5805 // Piola transform is required
5806# ifdef DEBUG
5807 Assert(this->gradients_quad_initialized == true,
5809# endif
5810
5811 AssertIndexRange(q_point, this->n_quadrature_points);
5812 Assert(this->jacobian != nullptr,
5814 "update_gradients"));
5815 const std::size_t nqp = this->n_quadrature_points;
5817
5818 if (!is_face &&
5820 {
5821 // Cartesian cell
5822 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5823 this->jacobian[0];
5824 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5825 const VectorizedArrayType inv_det =
5826 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5827 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5828 this->jacobian[0][2][2];
5829
5830 // J * grad_quad * J^-1 * det(J^-1)
5831 for (unsigned int d = 0; d < dim; ++d)
5832 for (unsigned int comp = 0; comp < n_components; ++comp)
5833 grad_out[comp][d] =
5834 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
5835 inv_t_jac[d][d] * jac[comp][comp] * inv_det;
5836 }
5837 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5838 {
5839 // Affine cell
5840 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5841 this->jacobian[0];
5842 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5843
5844 // Derivatives are reordered for faces. Need to take this into account
5845 const VectorizedArrayType inv_det =
5846 (is_face && dim == 2 && this->get_face_no() < 2) ?
5847 -determinant(inv_t_jac) :
5848 determinant(inv_t_jac);
5849
5850 VectorizedArrayType tmp;
5851 // J * grad_quad * J^-1 * det(J^-1)
5852 for (unsigned int comp = 0; comp < n_components; ++comp)
5853 for (unsigned int d = 0; d < dim; ++d)
5854 {
5855 tmp = 0;
5856 for (unsigned int f = 0; f < dim; ++f)
5857 for (unsigned int e = 0; e < dim; ++e)
5858 tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
5859 this->gradients_quad[(f * dim + e) * nqp + q_point];
5860
5861 grad_out[comp][d] = tmp;
5862 }
5863 }
5864 else
5865 {
5866 // General cell
5867 // Here we need the jacobian gradient and not the inverse which is
5868 // stored in this->jacobian_gradients
5870 }
5871 return grad_out;
5872 }
5873 else
5874 {
5875 return BaseClass::get_gradient(q_point);
5876 }
5877}
5878
5879
5880
5881template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5882inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5884 get_divergence(const unsigned int q_point) const
5885{
5886# ifdef DEBUG
5887 Assert(this->gradients_quad_initialized == true,
5889# endif
5890 AssertIndexRange(q_point, this->n_quadrature_points);
5891 Assert(this->jacobian != nullptr,
5893 "update_gradients"));
5894
5895 VectorizedArrayType divergence;
5896 const std::size_t nqp = this->n_quadrature_points;
5897
5898 if (this->data->element_type ==
5900 {
5901 if (!is_face &&
5903 {
5904 // Cartesian cell
5905 const VectorizedArrayType inv_det =
5906 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5907 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5908 this->jacobian[0][2][2];
5909
5910 // div * det(J^-1)
5911 divergence = this->gradients_quad[q_point] * inv_det;
5912 for (unsigned int d = 1; d < dim; ++d)
5913 divergence +=
5914 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5915 }
5916 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5917 {
5918 // Affine cell
5919 // Derivatives are reordered for faces. Need to take this into account
5920 const VectorizedArrayType inv_det =
5921 (is_face && dim == 2 && this->get_face_no() < 2) ?
5922 -determinant(this->jacobian[0]) :
5923 determinant(this->jacobian[0]);
5924
5925 // div * det(J^-1)
5926 divergence = this->gradients_quad[q_point] * inv_det;
5927 for (unsigned int d = 1; d < dim; ++d)
5928 divergence +=
5929 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
5930 }
5931 else
5932 {
5933 // General cell
5934 Assert(false, ExcNotImplemented());
5935 }
5936 }
5937 else
5938 {
5939 if (!is_face &&
5941 {
5942 // Cartesian cell
5943 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
5944 for (unsigned int d = 1; d < dim; ++d)
5945 divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
5946 this->jacobian[0][d][d];
5947 }
5948 else
5949 {
5950 // cell with general/constant Jacobian
5952 this->cell_type == internal::MatrixFreeFunctions::general ?
5953 this->jacobian[q_point] :
5954 this->jacobian[0];
5955 divergence = jac[0][0] * this->gradients_quad[q_point];
5956 for (unsigned int e = 1; e < dim; ++e)
5957 divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
5958 for (unsigned int d = 1; d < dim; ++d)
5959 for (unsigned int e = 0; e < dim; ++e)
5960 divergence +=
5961 jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
5962 }
5963 }
5964 return divergence;
5965}
5966
5967
5968
5969template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5972 get_symmetric_gradient(const unsigned int q_point) const
5973{
5974 // copy from generic function into dim-specialization function
5975 const auto grad = get_gradient(q_point);
5976 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
5977 VectorizedArrayType half = Number(0.5);
5978 for (unsigned int d = 0; d < dim; ++d)
5979 symmetrized[d] = grad[d][d];
5980 switch (dim)
5981 {
5982 case 1:
5983 break;
5984 case 2:
5985 symmetrized[2] = grad[0][1] + grad[1][0];
5986 symmetrized[2] *= half;
5987 break;
5988 case 3:
5989 symmetrized[3] = grad[0][1] + grad[1][0];
5990 symmetrized[3] *= half;
5991 symmetrized[4] = grad[0][2] + grad[2][0];
5992 symmetrized[4] *= half;
5993 symmetrized[5] = grad[1][2] + grad[2][1];
5994 symmetrized[5] *= half;
5995 break;
5996 default:
5997 Assert(false, ExcNotImplemented());
5998 }
6000}
6001
6002
6003
6004template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6006 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6008 const unsigned int q_point) const
6009{
6010 // copy from generic function into dim-specialization function
6011 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6012 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6013 switch (dim)
6014 {
6015 case 1:
6016 Assert(false,
6017 ExcMessage(
6018 "Computing the curl in 1d is not a useful operation"));
6019 break;
6020 case 2:
6021 curl[0] = grad[1][0] - grad[0][1];
6022 break;
6023 case 3:
6024 curl[0] = grad[2][1] - grad[1][2];
6025 curl[1] = grad[0][2] - grad[2][0];
6026 curl[2] = grad[1][0] - grad[0][1];
6027 break;
6028 default:
6029 Assert(false, ExcNotImplemented());
6030 }
6031 return curl;
6032}
6033
6034
6035
6036template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6039 get_hessian_diagonal(const unsigned int q_point) const
6040{
6041 return BaseClass::get_hessian_diagonal(q_point);
6042}
6043
6044
6045
6046template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6049 const unsigned int q_point) const
6050{
6051# ifdef DEBUG
6052 Assert(this->hessians_quad_initialized == true,
6054# endif
6055 AssertIndexRange(q_point, this->n_quadrature_points);
6056 return BaseClass::get_hessian(q_point);
6057}
6058
6059
6060template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6061inline DEAL_II_ALWAYS_INLINE void
6064 const unsigned int q_point)
6065{
6066 if (this->data->element_type ==
6068 {
6069 // Piola transform is required
6070 AssertIndexRange(q_point, this->n_quadrature_points);
6071 Assert(this->J_value != nullptr,
6073 "update_value"));
6074# ifdef DEBUG
6075 Assert(this->is_reinitialized, ExcNotInitialized());
6076 this->values_quad_submitted = true;
6077# endif
6078
6079 const std::size_t nqp = this->n_quadrature_points;
6080 if (!is_face &&
6082 {
6083 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6084 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6085
6086 for (unsigned int comp = 0; comp < n_components; ++comp)
6087 this->values_quad[comp * nqp + q_point] =
6088 val_in[comp] * weight * jac[comp][comp];
6089 }
6090 else
6091 {
6092 // Affine or general cell
6093 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6094 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6095 this->jacobian[q_point] :
6096 this->jacobian[0];
6098 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6099 transpose(invert(inv_t_jac)) :
6100 this->jacobian[1];
6101
6102 // Derivatives are reordered for faces. Need to take this into account
6103 // and 1/inv_det != J_value for faces
6104 const VectorizedArrayType fac =
6105 (!is_face) ?
6106 this->quadrature_weights[q_point] :
6107 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6108 this->J_value[q_point] :
6109 this->J_value[0] * this->quadrature_weights[q_point]) *
6110 ((dim == 2 && this->get_face_no() < 2) ?
6111 -determinant(inv_t_jac) :
6112 determinant(inv_t_jac)));
6113
6114 // J^T * u * factor
6115 for (unsigned int comp = 0; comp < n_components; ++comp)
6116 {
6117 this->values_quad[comp * nqp + q_point] =
6118 val_in[0] * jac[0][comp] * fac;
6119 for (unsigned int e = 1; e < dim; ++e)
6120 this->values_quad[comp * nqp + q_point] +=
6121 val_in[e] * jac[e][comp] * fac;
6122 }
6123 }
6124 }
6125 else
6126 {
6127 // No Piola transform
6128 BaseClass::submit_value(val_in, q_point);
6129 }
6130}
6131
6132template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6133inline DEAL_II_ALWAYS_INLINE void
6136 const unsigned int q_point)
6137{
6138 if (this->data->element_type ==
6140 {
6141 // Piola transform is required
6142
6143# ifdef DEBUG
6144 Assert(this->is_reinitialized, ExcNotInitialized());
6145# endif
6146 AssertIndexRange(q_point, this->n_quadrature_points);
6147 Assert(this->J_value != nullptr,
6149 "update_gradients"));
6150 Assert(this->jacobian != nullptr,
6152 "update_gradients"));
6153# ifdef DEBUG
6154 this->gradients_quad_submitted = true;
6155# endif
6156
6157 const std::size_t nqp = this->n_quadrature_points;
6158 if (!is_face &&
6160 {
6161 // Cartesian cell
6162 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6163 this->jacobian[0];
6164 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6165 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6166 for (unsigned int d = 0; d < dim; ++d)
6167 for (unsigned int comp = 0; comp < n_components; ++comp)
6168 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6169 grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
6170 }
6171 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6172 {
6173 // Affine cell
6174 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6175 this->jacobian[0];
6176 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6177
6178 // Derivatives are reordered for faces. Need to take this into account
6179 // and 1/inv_det != J_value for faces
6180 const VectorizedArrayType fac =
6181 (!is_face) ? this->quadrature_weights[q_point] :
6182 this->J_value[0] * this->quadrature_weights[q_point] *
6183 ((dim == 2 && this->get_face_no() < 2) ?
6184 -determinant(inv_t_jac) :
6185 determinant(inv_t_jac));
6186
6187 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6188 for (unsigned int comp = 0; comp < n_components; ++comp)
6189 for (unsigned int d = 0; d < dim; ++d)
6190 {
6191 VectorizedArrayType tmp = 0;
6192 for (unsigned int f = 0; f < dim; ++f)
6193 for (unsigned int e = 0; e < dim; ++e)
6194 tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac;
6195
6196 this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
6197 }
6198 }
6199 else
6200 {
6201 // General cell
6203 }
6204 }
6205 else
6206 {
6207 BaseClass::submit_gradient(grad_in, q_point);
6208 }
6209}
6210
6211
6212
6213template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6214inline DEAL_II_ALWAYS_INLINE void
6217 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6218 const unsigned int q_point)
6219{
6220 if (this->data->element_type ==
6222 {
6223 // Piola transform is required
6224 const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
6226 submit_gradient(grad, q_point);
6227 }
6228 else
6229 {
6230 BaseClass::submit_gradient(grad_in, q_point);
6231 }
6232}
6233
6234
6235
6236template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6237inline DEAL_II_ALWAYS_INLINE void
6239 submit_divergence(const VectorizedArrayType div_in,
6240 const unsigned int q_point)
6241{
6242# ifdef DEBUG
6243 Assert(this->is_reinitialized, ExcNotInitialized());
6244# endif
6245 AssertIndexRange(q_point, this->n_quadrature_points);
6246 Assert(this->J_value != nullptr,
6248 "update_gradients"));
6249 Assert(this->jacobian != nullptr,
6251 "update_gradients"));
6252# ifdef DEBUG
6253 this->gradients_quad_submitted = true;
6254# endif
6255
6256 const std::size_t nqp = this->n_quadrature_points;
6257 if (this->data->element_type ==
6259 {
6260 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6261 {
6262 // Affine cell
6263
6264 // Derivatives are reordered for faces. Need to take this into account
6265 // and 1/inv_det != J_value for faces
6266 const VectorizedArrayType fac =
6267 ((!is_face) ?
6268 1 :
6269 this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ?
6270 -determinant(this->jacobian[0]) :
6271 determinant(this->jacobian[0]))) *
6272 this->quadrature_weights[q_point] * div_in;
6273
6274 for (unsigned int d = 0; d < dim; ++d)
6275 {
6276 this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
6277 for (unsigned int e = d + 1; e < dim; ++e)
6278 {
6279 this->gradients_quad[(dim * d + e) * nqp + q_point] =
6280 VectorizedArrayType();
6281 this->gradients_quad[(dim * e + d) * nqp + q_point] =
6282 VectorizedArrayType();
6283 }
6284 }
6285 }
6286 else
6287 {
6288 // General cell
6290 }
6291 }
6292 else
6293 {
6294 if (!is_face &&
6296 {
6297 const VectorizedArrayType fac =
6298 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6299 for (unsigned int d = 0; d < dim; ++d)
6300 {
6301 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6302 (fac * this->jacobian[0][d][d]);
6303 for (unsigned int e = d + 1; e < dim; ++e)
6304 {
6305 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6306 VectorizedArrayType();
6307 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6308 VectorizedArrayType();
6309 }
6310 }
6311 }
6312 else
6313 {
6315 this->cell_type == internal::MatrixFreeFunctions::general ?
6316 this->jacobian[q_point] :
6317 this->jacobian[0];
6318 const VectorizedArrayType fac =
6319 (this->cell_type == internal::MatrixFreeFunctions::general ?
6320 this->J_value[q_point] :
6321 this->J_value[0] * this->quadrature_weights[q_point]) *
6322 div_in;
6323 for (unsigned int d = 0; d < dim; ++d)
6324 {
6325 for (unsigned int e = 0; e < dim; ++e)
6326 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6327 jac[d][e] * fac;
6328 }
6329 }
6330 }
6331}
6332
6333
6334
6335template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6336inline DEAL_II_ALWAYS_INLINE void
6340 const unsigned int q_point)
6341{
6343 this->data->element_type !=
6346
6347 // could have used base class operator, but that involves some overhead
6348 // which is inefficient. it is nice to have the symmetric tensor because
6349 // that saves some operations
6350# ifdef DEBUG
6351 Assert(this->is_reinitialized, ExcNotInitialized());
6352# endif
6353 AssertIndexRange(q_point, this->n_quadrature_points);
6354 Assert(this->J_value != nullptr,
6356 "update_gradients"));
6357 Assert(this->jacobian != nullptr,
6359 "update_gradients"));
6360# ifdef DEBUG
6361 this->gradients_quad_submitted = true;
6362# endif
6363
6364 const std::size_t nqp = this->n_quadrature_points;
6365 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6366 {
6367 const VectorizedArrayType JxW =
6368 this->J_value[0] * this->quadrature_weights[q_point];
6369 for (unsigned int d = 0; d < dim; ++d)
6370 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6371 (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6372 for (unsigned int e = 0, counter = dim; e < dim; ++e)
6373 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6374 {
6375 const VectorizedArrayType value =
6376 sym_grad.access_raw_entry(counter) * JxW;
6377 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6378 value * this->jacobian[0][d][d];
6379 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6380 value * this->jacobian[0][e][e];
6381 }
6382 }
6383 // general/affine cell type
6384 else
6385 {
6386 const VectorizedArrayType JxW =
6387 this->cell_type == internal::MatrixFreeFunctions::general ?
6388 this->J_value[q_point] :
6389 this->J_value[0] * this->quadrature_weights[q_point];
6391 this->cell_type == internal::MatrixFreeFunctions::general ?
6392 this->jacobian[q_point] :
6393 this->jacobian[0];
6394 VectorizedArrayType weighted[dim][dim];
6395 for (unsigned int i = 0; i < dim; ++i)
6396 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6397 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6398 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6399 {
6400 const VectorizedArrayType value =
6401 sym_grad.access_raw_entry(counter) * JxW;
6402 weighted[i][j] = value;
6403 weighted[j][i] = value;
6404 }
6405 for (unsigned int comp = 0; comp < dim; ++comp)
6406 for (unsigned int d = 0; d < dim; ++d)
6407 {
6408 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6409 for (unsigned int e = 1; e < dim; ++e)
6410 new_val += jac[e][d] * weighted[comp][e];
6411 this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6412 }
6413 }
6414}
6415
6416
6417
6418template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6419inline DEAL_II_ALWAYS_INLINE void
6422 const unsigned int q_point)
6423{
6425 switch (dim)
6426 {
6427 case 1:
6428 Assert(false,
6429 ExcMessage(
6430 "Testing by the curl in 1d is not a useful operation"));
6431 break;
6432 case 2:
6433 grad[1][0] = curl[0];
6434 grad[0][1] = -curl[0];
6435 break;
6436 case 3:
6437 grad[2][1] = curl[0];
6438 grad[1][2] = -curl[0];
6439 grad[0][2] = curl[1];
6440 grad[2][0] = -curl[1];
6441 grad[1][0] = curl[2];
6442 grad[0][1] = -curl[2];
6443 break;
6444 default:
6445 Assert(false, ExcNotImplemented());
6446 }
6447 submit_gradient(grad, q_point);
6448}
6449
6450
6451/*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6452
6453
6454template <typename Number, bool is_face, typename VectorizedArrayType>
6458 const unsigned int dof_no,
6459 const unsigned int first_selected_component,
6460 const unsigned int quad_no,
6461 const unsigned int fe_degree,
6462 const unsigned int n_q_points,
6463 const bool is_interior_face,
6464 const unsigned int active_fe_index,
6465 const unsigned int active_quad_index,
6466 const unsigned int face_type)
6467 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6468 matrix_free,
6469 dof_no,
6470 first_selected_component,
6471 quad_no,
6472 fe_degree,
6473 n_q_points,
6474 is_interior_face,
6475 active_fe_index,
6476 active_quad_index,
6477 face_type)
6478{}
6479
6480
6481
6482template <typename Number, bool is_face, typename VectorizedArrayType>
6485 const Mapping<1> & mapping,
6486 const FiniteElement<1> &fe,
6487 const Quadrature<1> & quadrature,
6488 const UpdateFlags update_flags,
6489 const unsigned int first_selected_component,
6491 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6492 mapping,
6493 fe,
6494 quadrature,
6495 update_flags,
6496 first_selected_component,
6497 other)
6498{}
6499
6500
6501
6502template <typename Number, bool is_face, typename VectorizedArrayType>
6506 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6507{}
6508
6509
6510
6511template <typename Number, bool is_face, typename VectorizedArrayType>
6515{
6517 other);
6518 return *this;
6519}
6520
6521
6522
6523template <typename Number, bool is_face, typename VectorizedArrayType>
6524inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6526 const unsigned int dof) const
6527{
6528 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6529 return this->values_dofs[dof];
6530}
6531
6532
6533
6534template <typename Number, bool is_face, typename VectorizedArrayType>
6535inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6537 const unsigned int q_point) const
6538{
6539# ifdef DEBUG
6540 Assert(this->values_quad_initialized == true,
6542# endif
6543 AssertIndexRange(q_point, this->n_quadrature_points);
6544 return this->values_quad[q_point];
6545}
6546
6547
6548
6549template <typename Number, bool is_face, typename VectorizedArrayType>
6552 const unsigned int q_point) const
6553{
6554 // could use the base class gradient, but that involves too many inefficient
6555 // initialization operations on tensors
6556
6557# ifdef DEBUG
6558 Assert(this->gradients_quad_initialized == true,
6560# endif
6561 AssertIndexRange(q_point, this->n_quadrature_points);
6562
6564 this->cell_type == internal::MatrixFreeFunctions::general ?
6565 this->jacobian[q_point] :
6566 this->jacobian[0];
6567
6569 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6570
6571 return grad_out;
6572}
6573
6574
6575
6576template <typename Number, bool is_face, typename VectorizedArrayType>
6577inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6579 const unsigned int q_point) const
6580{
6581 return get_gradient(q_point)[0];
6582}
6583
6584
6585
6586template <typename Number, bool is_face, typename VectorizedArrayType>
6587inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6589 get_normal_derivative(const unsigned int q_point) const
6590{
6591 return BaseClass::get_normal_derivative(q_point)[0];
6592}
6593
6594
6595
6596template <typename Number, bool is_face, typename VectorizedArrayType>
6599 const unsigned int q_point) const
6600{
6601 return BaseClass::get_hessian(q_point)[0];
6602}
6603
6604
6605
6606template <typename Number, bool is_face, typename VectorizedArrayType>
6609 get_hessian_diagonal(const unsigned int q_point) const
6610{
6611 return BaseClass::get_hessian_diagonal(q_point)[0];
6612}
6613
6614
6615
6616template <typename Number, bool is_face, typename VectorizedArrayType>
6617inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6619 const unsigned int q_point) const
6620{
6621 return BaseClass::get_laplacian(q_point)[0];
6622}
6623
6624
6625
6626template <typename Number, bool is_face, typename VectorizedArrayType>
6629 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6630{
6631# ifdef DEBUG
6632 this->dof_values_initialized = true;
6633 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6634# endif
6635 this->values_dofs[dof] = val_in;
6636}
6637
6638
6639
6640template <typename Number, bool is_face, typename VectorizedArrayType>
6641inline DEAL_II_ALWAYS_INLINE void
6643 const VectorizedArrayType val_in,
6644 const unsigned int q_point)
6645{
6646# ifdef DEBUG
6647 Assert(this->is_reinitialized, ExcNotInitialized());
6648# endif
6649 AssertIndexRange(q_point, this->n_quadrature_points);
6650# ifdef DEBUG
6651 this->values_quad_submitted = true;
6652# endif
6653
6654 if (this->cell_type == internal::MatrixFreeFunctions::general)
6655 {
6656 const VectorizedArrayType JxW = this->J_value[q_point];
6657 this->values_quad[q_point] = val_in * JxW;
6658 }
6659 else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6660 {
6661 const VectorizedArrayType JxW =
6662 this->J_value[0] * this->quadrature_weights[q_point];
6663 this->values_quad[q_point] = val_in * JxW;
6664 }
6665}
6666
6667
6668
6669template <typename Number, bool is_face, typename VectorizedArrayType>
6670inline DEAL_II_ALWAYS_INLINE void
6673 const unsigned int q_point)
6674{
6675 submit_value(val_in[0], q_point);
6676}
6677
6678
6679
6680template <typename Number, bool is_face, typename VectorizedArrayType>
6681inline DEAL_II_ALWAYS_INLINE void
6684 const unsigned int q_point)
6685{
6686 submit_gradient(grad_in[0], q_point);
6687}
6688
6689
6690
6691template <typename Number, bool is_face, typename VectorizedArrayType>
6692inline DEAL_II_ALWAYS_INLINE void
6694 const VectorizedArrayType grad_in,
6695 const unsigned int q_point)
6696{
6697# ifdef DEBUG
6698 Assert(this->is_reinitialized, ExcNotInitialized());
6699# endif
6700 AssertIndexRange(q_point, this->n_quadrature_points);
6701# ifdef DEBUG
6702 this->gradients_quad_submitted = true;
6703# endif
6704
6706 this->cell_type == internal::MatrixFreeFunctions::general ?
6707 this->jacobian[q_point] :
6708 this->jacobian[0];
6709 const VectorizedArrayType JxW =
6710 this->cell_type == internal::MatrixFreeFunctions::general ?
6711 this->J_value[q_point] :
6712 this->J_value[0] * this->quadrature_weights[q_point];
6713
6714 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
6715}
6716
6717
6718
6719template <typename Number, bool is_face, typename VectorizedArrayType>
6720inline DEAL_II_ALWAYS_INLINE void
6723 const unsigned int q_point)
6724{
6725 submit_gradient(grad_in[0][0], q_point);
6726}
6727
6728
6729
6730template <typename Number, bool is_face, typename VectorizedArrayType>
6731inline DEAL_II_ALWAYS_INLINE void
6733 submit_normal_derivative(const VectorizedArrayType grad_in,
6734 const unsigned int q_point)
6735{
6737 grad[0] = grad_in;
6738 BaseClass::submit_normal_derivative(grad, q_point);
6739}
6740
6741
6742
6743template <typename Number, bool is_face, typename VectorizedArrayType>
6744inline DEAL_II_ALWAYS_INLINE void
6747 const unsigned int q_point)
6748{
6749 BaseClass::submit_normal_derivative(grad_in, q_point);
6750}
6751
6752
6753template <typename Number, bool is_face, typename VectorizedArrayType>
6754inline DEAL_II_ALWAYS_INLINE void
6756 const Tensor<2, 1, VectorizedArrayType> hessian_in,
6757 const unsigned int q_point)
6758{
6760 hessian[0] = hessian_in;
6761 BaseClass::submit_hessian(hessian, q_point);
6762}
6763
6764
6765template <typename Number, bool is_face, typename VectorizedArrayType>
6766inline VectorizedArrayType
6768 integrate_value() const
6769{
6770 return BaseClass::integrate_value()[0];
6771}
6772
6773
6774
6775/*-------------------------- FEEvaluation -----------------------------------*/
6776
6777
6778template <int dim,
6779 int fe_degree,
6780 int n_q_points_1d,
6781 int n_components_,
6782 typename Number,
6783 typename VectorizedArrayType>
6784inline FEEvaluation<dim,
6785 fe_degree,
6786 n_q_points_1d,
6787 n_components_,
6788 Number,
6789 VectorizedArrayType>::
6790 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6791 const unsigned int fe_no,
6792 const unsigned int quad_no,
6793 const unsigned int first_selected_component,
6794 const unsigned int active_fe_index,
6795 const unsigned int active_quad_index)
6796 : BaseClass(matrix_free,
6797 fe_no,
6798 first_selected_component,
6799 quad_no,
6800 fe_degree,
6801 static_n_q_points,
6802 true /*note: this is not a face*/,
6803 active_fe_index,
6804 active_quad_index)
6805 , dofs_per_component(this->data->dofs_per_component_on_cell)
6806 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6807 , n_q_points(this->data->n_q_points)
6808{
6809 check_template_arguments(fe_no, 0);
6810}
6811
6812
6813
6814template <int dim,
6815 int fe_degree,
6816 int n_q_points_1d,
6817 int n_components_,
6818 typename Number,
6819 typename VectorizedArrayType>
6820inline FEEvaluation<dim,
6821 fe_degree,
6822 n_q_points_1d,
6823 n_components_,
6824 Number,
6825 VectorizedArrayType>::
6826 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
6827 const std::pair<unsigned int, unsigned int> & range,
6828 const unsigned int dof_no,
6829 const unsigned int quad_no,
6830 const unsigned int first_selected_component)
6831 : FEEvaluation(matrix_free,
6832 dof_no,
6833 quad_no,
6834 first_selected_component,
6835 matrix_free.get_cell_active_fe_index(range))
6836{}
6837
6838
6839
6840template <int dim,
6841 int fe_degree,
6842 int n_q_points_1d,
6843 int n_components_,
6844 typename Number,
6845 typename VectorizedArrayType>
6846inline FEEvaluation<dim,
6847 fe_degree,
6848 n_q_points_1d,
6849 n_components_,
6850 Number,
6851 VectorizedArrayType>::
6852 FEEvaluation(const Mapping<dim> & mapping,
6853 const FiniteElement<dim> &fe,
6854 const Quadrature<1> & quadrature,
6855 const UpdateFlags update_flags,
6856 const unsigned int first_selected_component)
6857 : BaseClass(mapping,
6858 fe,
6859 quadrature,
6860 update_flags,
6861 first_selected_component,
6862 nullptr)
6863 , dofs_per_component(this->data->dofs_per_component_on_cell)
6864 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6865 , n_q_points(this->data->n_q_points)
6866{
6867 check_template_arguments(numbers::invalid_unsigned_int, 0);
6868}
6869
6870
6871
6872template <int dim,
6873 int fe_degree,
6874 int n_q_points_1d,
6875 int n_components_,
6876 typename Number,
6877 typename VectorizedArrayType>
6878inline FEEvaluation<dim,
6879 fe_degree,
6880 n_q_points_1d,
6881 n_components_,
6882 Number,
6883 VectorizedArrayType>::
6884 FEEvaluation(const FiniteElement<dim> &fe,
6885 const Quadrature<1> & quadrature,
6886 const UpdateFlags update_flags,
6887 const unsigned int first_selected_component)
6888 : BaseClass(StaticMappingQ1<dim>::mapping,
6889 fe,
6890 quadrature,
6891 update_flags,
6892 first_selected_component,
6893 nullptr)
6894 , dofs_per_component(this->data->dofs_per_component_on_cell)
6895 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6896 , n_q_points(this->data->n_q_points)
6897{
6898 check_template_arguments(numbers::invalid_unsigned_int, 0);
6899}
6900
6901
6902
6903template <int dim,
6904 int fe_degree,
6905 int n_q_points_1d,
6906 int n_components_,
6907 typename Number,
6908 typename VectorizedArrayType>
6909inline FEEvaluation<dim,
6910 fe_degree,
6911 n_q_points_1d,
6912 n_components_,
6913 Number,
6914 VectorizedArrayType>::
6915 FEEvaluation(const FiniteElement<dim> & fe,
6917 const unsigned int first_selected_component)
6918 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
6919 fe,
6920 other.mapped_geometry->get_quadrature(),
6921 other.mapped_geometry->get_fe_values().get_update_flags(),
6922 first_selected_component,
6923 &other)
6924 , dofs_per_component(this->data->dofs_per_component_on_cell)
6925 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6926 , n_q_points(this->data->n_q_points)
6927{
6928 check_template_arguments(numbers::invalid_unsigned_int, 0);
6929}
6930
6931
6932
6933template <int dim,
6934 int fe_degree,
6935 int n_q_points_1d,
6936 int n_components_,
6937 typename Number,
6938 typename VectorizedArrayType>
6939inline FEEvaluation<dim,
6940 fe_degree,
6941 n_q_points_1d,
6942 n_components_,
6943 Number,
6944 VectorizedArrayType>::FEEvaluation(const FEEvaluation
6945 &other)
6946 : BaseClass(other)
6947 , dofs_per_component(this->data->dofs_per_component_on_cell)
6948 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
6949 , n_q_points(this->data->n_q_points)
6950{
6951 check_template_arguments(numbers::invalid_unsigned_int, 0);
6952}
6953
6954
6955
6956template <int dim,
6957 int fe_degree,
6958 int n_q_points_1d,
6959 int n_components_,
6960 typename Number,
6961 typename VectorizedArrayType>
6962inline FEEvaluation<dim,
6963 fe_degree,
6964 n_q_points_1d,
6965 n_components_,
6966 Number,
6967 VectorizedArrayType> &
6968FEEvaluation<dim,
6969 fe_degree,
6970 n_q_points_1d,
6971 n_components_,
6972 Number,
6973 VectorizedArrayType>::operator=(const FEEvaluation &other)
6974{
6975 BaseClass::operator=(other);
6976 check_template_arguments(numbers::invalid_unsigned_int, 0);
6977 return *this;
6978}
6979
6980
6981
6982template <int dim,
6983 int fe_degree,
6984 int n_q_points_1d,
6985 int n_components_,
6986 typename Number,
6987 typename VectorizedArrayType>
6988inline void
6989FEEvaluation<dim,
6990 fe_degree,
6991 n_q_points_1d,
6992 n_components_,
6993 Number,
6994 VectorizedArrayType>::
6995 check_template_arguments(const unsigned int dof_no,
6996 const unsigned int first_selected_component)
6997{
6998 (void)dof_no;
6999 (void)first_selected_component;
7000
7001 Assert(
7002 this->data->dofs_per_component_on_cell > 0,
7003 ExcMessage(
7004 "There is nothing useful you can do with an FEEvaluation object with "
7005 "FE_Nothing, i.e., without DoFs! If you have passed to "
7006 "MatrixFree::reinit() a collection of finite elements also containing "
7007 "FE_Nothing, please check - before creating FEEvaluation - the category "
7008 "of the current range by calling either "
7009 "MatrixFree::get_cell_range_category(range) or "
7010 "MatrixFree::get_face_range_category(range). The returned category "
7011 "is the index of the active FE, which you can use to exclude "
7012 "FE_Nothing."));
7013
7014# ifdef DEBUG
7015 // print error message when the dimensions do not match. Propose a possible
7016 // fix
7017 if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7018 static_cast<unsigned int>(fe_degree) !=
7019 this->data->data.front().fe_degree) ||
7020 n_q_points != this->n_quadrature_points)
7021 {
7022 std::string message =
7023 "-------------------------------------------------------\n";
7024 message += "Illegal arguments in constructor/wrong template arguments!\n";
7025 message += " Called --> FEEvaluation<dim,";
7026 message += Utilities::int_to_string(fe_degree) + ",";
7027 message += Utilities::int_to_string(n_q_points_1d);
7028 message += "," + Utilities::int_to_string(n_components);
7029 message += ",Number>(data";
7030 if (first_selected_component != numbers::invalid_unsigned_int)
7031 {
7032 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7033 message += Utilities::int_to_string(this->quad_no) + ", ";
7034 message += Utilities::int_to_string(first_selected_component);
7035 }
7036 message += ")\n";
7037
7038 // check whether some other vector component has the correct number of
7039 // points
7040 unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7041 proposed_fe_comp = numbers::invalid_unsigned_int,
7042 proposed_quad_comp = numbers::invalid_unsigned_int;
7043 if (dof_no != numbers::invalid_unsigned_int)
7044 {
7045 if (static_cast<unsigned int>(fe_degree) ==
7046 this->data->data.front().fe_degree)
7047 {
7048 proposed_dof_comp = dof_no;
7049 proposed_fe_comp = first_selected_component;
7050 }
7051 else
7052 for (unsigned int no = 0; no < this->matrix_free->n_components();
7053 ++no)
7054 for (unsigned int nf = 0;
7055 nf < this->matrix_free->n_base_elements(no);
7056 ++nf)
7057 if (this->matrix_free
7058 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7059 .data.front()
7060 .fe_degree == static_cast<unsigned int>(fe_degree))
7061 {
7062 proposed_dof_comp = no;
7063 proposed_fe_comp = nf;
7064 break;
7065 }
7066 if (n_q_points ==
7067 this->mapping_data->descriptor[this->active_quad_index]
7068 .n_q_points)
7069 proposed_quad_comp = this->quad_no;
7070 else
7071 for (unsigned int no = 0;
7072 no < this->matrix_free->get_mapping_info().cell_data.size();
7073 ++no)
7074 if (this->matrix_free->get_mapping_info()
7075 .cell_data[no]
7076 .descriptor[this->active_quad_index]
7077 .n_q_points == n_q_points)
7078 {
7079 proposed_quad_comp = no;
7080 break;
7081 }
7082 }
7083 if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7084 proposed_quad_comp != numbers::invalid_unsigned_int)
7085 {
7086 if (proposed_dof_comp != first_selected_component)
7087 message += "Wrong vector component selection:\n";
7088 else
7089 message += "Wrong quadrature formula selection:\n";
7090 message += " Did you mean FEEvaluation<dim,";
7091 message += Utilities::int_to_string(fe_degree) + ",";
7092 message += Utilities::int_to_string(n_q_points_1d);
7093 message += "," + Utilities::int_to_string(n_components);
7094 message += ",Number>(data";
7095 if (dof_no != numbers::invalid_unsigned_int)
7096 {
7097 message +=
7098 ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7099 message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7100 message += Utilities::int_to_string(proposed_fe_comp);
7101 }
7102 message += ")?\n";
7103 std::string correct_pos;
7104 if (proposed_dof_comp != dof_no)
7105 correct_pos = " ^ ";
7106 else
7107 correct_pos = " ";
7108 if (proposed_quad_comp != this->quad_no)
7109 correct_pos += " ^ ";
7110 else
7111 correct_pos += " ";
7112 if (proposed_fe_comp != first_selected_component)
7113 correct_pos += " ^\n";
7114 else
7115 correct_pos += " \n";
7116 message += " " +
7117 correct_pos;
7118 }
7119 // ok, did not find the numbers specified by the template arguments in
7120 // the given list. Suggest correct template arguments
7121 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7122 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7123 message += "Wrong template arguments:\n";
7124 message += " Did you mean FEEvaluation<dim,";
7125 message +=
7126 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7127 message += Utilities::int_to_string(proposed_n_q_points_1d);
7128 message += "," + Utilities::int_to_string(n_components);
7129 message += ",Number>(data";
7130 if (dof_no != numbers::invalid_unsigned_int)
7131 {
7132 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7133 message += Utilities::int_to_string(this->quad_no);
7134 message += ", " + Utilities::int_to_string(first_selected_component);
7135 }
7136 message += ")?\n";
7137 std::string correct_pos;
7138 if (this->data->data.front().fe_degree !=
7139 static_cast<unsigned int>(fe_degree))
7140 correct_pos = " ^";
7141 else
7142 correct_pos = " ";
7143 if (proposed_n_q_points_1d != n_q_points_1d)
7144 correct_pos += " ^\n";
7145 else
7146 correct_pos += " \n";
7147 message += " " + correct_pos;
7148
7149 Assert(static_cast<unsigned int>(fe_degree) ==
7150 this->data->data.front().fe_degree &&
7151 n_q_points == this->n_quadrature_points,
7152 ExcMessage(message));
7153 }
7154 if (dof_no != numbers::invalid_unsigned_int)
7156 n_q_points,
7157 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7158# endif
7159}
7160
7161
7162
7163template <int dim,
7164 int fe_degree,
7165 int n_q_points_1d,
7166 int n_components_,
7167 typename Number,
7168 typename VectorizedArrayType>
7169inline void
7170FEEvaluation<dim,
7171 fe_degree,
7172 n_q_points_1d,
7173 n_components_,
7174 Number,
7175 VectorizedArrayType>::reinit(const unsigned int cell_index)
7176{
7177 Assert(this->mapped_geometry == nullptr,
7178 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7179 " Integer indexing is not possible"));
7180 if (this->mapped_geometry != nullptr)
7181 return;
7182
7183 Assert(this->dof_info != nullptr, ExcNotInitialized());
7184 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7185 this->cell = cell_index;
7186 this->cell_type =
7187 this->matrix_free->get_mapping_info().get_cell_type(cell_index);
7188
7189 const unsigned int offsets =
7190 this->mapping_data->data_index_offsets[cell_index];
7191 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7192 this->J_value = &this->mapping_data->JxW_values[offsets];
7193 this->jacobian_gradients =
7194 this->mapping_data->jacobian_gradients[0].data() + offsets;
7195
7196 unsigned int i = 0;
7197 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7198 ++i)
7199 this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7200 for (; i < VectorizedArrayType::size(); ++i)
7201 this->cell_ids[i] = numbers::invalid_unsigned_int;
7202
7203 if (this->mapping_data->quadrature_points.empty() == false)
7204 this->quadrature_points =
7205 &this->mapping_data->quadrature_points
7206 [this->mapping_data->quadrature_point_offsets[this->cell]];
7207
7208# ifdef DEBUG
7209 this->is_reinitialized = true;
7210 this->dof_values_initialized = false;
7211 this->values_quad_initialized = false;
7212 this->gradients_quad_initialized = false;
7213 this->hessians_quad_initialized = false;
7214# endif
7215}
7216
7217
7218
7219template <int dim,
7220 int fe_degree,
7221 int n_q_points_1d,
7222 int n_components_,
7223 typename Number,
7224 typename VectorizedArrayType>
7225inline void
7226FEEvaluation<dim,
7227 fe_degree,
7228 n_q_points_1d,
7229 n_components_,
7230 Number,
7231 VectorizedArrayType>::
7232 reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
7233{
7234 Assert(this->dof_info != nullptr, ExcNotInitialized());
7235 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7236
7237 this->cell = numbers::invalid_unsigned_int;
7238 this->cell_ids = cell_ids;
7239
7240 // determine type of cell batch
7242
7243 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7244 {
7245 const unsigned int cell_index = cell_ids[v];
7246
7248 continue;
7249
7250 this->cell_type =
7251 std::max(this->cell_type,
7252 this->matrix_free->get_mapping_info().get_cell_type(
7253 cell_index / VectorizedArrayType::size()));
7254 }
7255
7256 // allocate memory for internal data storage
7257 if (this->mapped_geometry == nullptr)
7258 this->mapped_geometry =
7259 std::make_shared<internal::MatrixFreeFunctions::
7260 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7261
7262 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7263
7264 auto &this_jacobian_data = mapping_storage.jacobians[0];
7265 auto &this_J_value_data = mapping_storage.JxW_values;
7266 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7267 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7268
7270 {
7271 if (this->mapping_data->jacobians[0].size() > 0)
7272 this_jacobian_data.resize_fast(2);
7273
7274 if (this->mapping_data->JxW_values.size() > 0)
7275 this_J_value_data.resize_fast(1);
7276
7277 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7278 this_jacobian_gradients_data.resize_fast(1);
7279
7280 if (this->mapping_data->quadrature_points.size() > 0)
7281 this_quadrature_points_data.resize_fast(1);
7282 }
7283 else
7284 {
7285 if (this->mapping_data->jacobians[0].size() > 0)
7286 this_jacobian_data.resize_fast(this->n_quadrature_points);
7287
7288 if (this->mapping_data->JxW_values.size() > 0)
7289 this_J_value_data.resize_fast(this->n_quadrature_points);
7290
7291 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7292 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7293
7294 if (this->mapping_data->quadrature_points.size() > 0)
7295 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7296 }
7297
7298 // set pointers to internal data storage
7299 this->jacobian = this_jacobian_data.data();
7300 this->J_value = this_J_value_data.data();
7301 this->jacobian_gradients = this_jacobian_gradients_data.data();
7302 this->quadrature_points = this_quadrature_points_data.data();
7303
7304 // fill internal data storage lane by lane
7305 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7306 {
7307 const unsigned int cell_index = cell_ids[v];
7308
7310 continue;
7311
7312 const unsigned int cell_batch_index =
7313 cell_index / VectorizedArrayType::size();
7314 const unsigned int offsets =
7315 this->mapping_data->data_index_offsets[cell_batch_index];
7316 const unsigned int lane = cell_index % VectorizedArrayType::size();
7317
7318 if (this->cell_type <=
7320 {
7321 // case that all cells are Cartesian or affine
7322 const unsigned int q = 0;
7323
7324 if (this->mapping_data->JxW_values.size() > 0)
7325 this_J_value_data[q][v] =
7326 this->mapping_data->JxW_values[offsets + q][lane];
7327
7328 if (this->mapping_data->jacobians[0].size() > 0)
7329 for (unsigned int q = 0; q < 2; ++q)
7330 for (unsigned int i = 0; i < dim; ++i)
7331 for (unsigned int j = 0; j < dim; ++j)
7332 this_jacobian_data[q][i][j][v] =
7333 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7334
7335 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7336 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7337 for (unsigned int j = 0; j < dim; ++j)
7338 this_jacobian_gradients_data[q][i][j][v] =
7339 this->mapping_data
7340 ->jacobian_gradients[0][offsets + q][i][j][lane];
7341
7342 if (this->mapping_data->quadrature_points.size() > 0)
7343 for (unsigned int i = 0; i < dim; ++i)
7344 this_quadrature_points_data[q][i][v] =
7345 this->mapping_data->quadrature_points
7346 [this->mapping_data
7347 ->quadrature_point_offsets[cell_batch_index] +
7348 q][i][lane];
7349 }
7350 else
7351 {
7352 // general case that at least one cell is not Cartesian or affine
7353 const auto cell_type =
7354 this->matrix_free->get_mapping_info().get_cell_type(
7355 cell_batch_index);
7356
7357 for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
7358 {
7359 const unsigned int q_src =
7360 (cell_type <=
7362 0 :
7363 q;
7364
7365 if (this->mapping_data->JxW_values.size() > 0)
7366 this_J_value_data[q][v] =
7367 this->mapping_data->JxW_values[offsets + q_src][lane];
7368
7369 if (this->mapping_data->jacobians[0].size() > 0)
7370 for (unsigned int i = 0; i < dim; ++i)
7371 for (unsigned int j = 0; j < dim; ++j)
7372 this_jacobian_data[q][i][j][v] =
7373 this->mapping_data
7374 ->jacobians[0][offsets + q_src][i][j][lane];
7375
7376 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7377 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7378 for (unsigned int j = 0; j < dim; ++j)
7379 this_jacobian_gradients_data[q][i][j][v] =
7380 this->mapping_data
7381 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7382
7383 if (this->mapping_data->quadrature_points.size() > 0)
7384 {
7385 if (cell_type <=
7387 {
7388 // affine case: quadrature points are not available but
7389 // have to be computed from the corner point and the
7390 // Jacobian
7392 this->mapping_data->quadrature_points
7393 [this->mapping_data
7394 ->quadrature_point_offsets[cell_batch_index] +
7395 0];
7396
7397 const