Reference documentation for deal.II version GIT relicensing-14-gd371fb467a 2024-02-26 20:50:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2012 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_matrix_free_fe_evaluation_h
17#define dealii_matrix_free_fe_evaluation_h
18
19
20#include <deal.II/base/config.h>
21
28
30
43
44#include <type_traits>
45
46
48
49
50
88template <int dim,
89 int n_components_,
90 typename Number,
91 bool is_face,
92 typename VectorizedArrayType>
94 : public FEEvaluationData<dim, VectorizedArrayType, is_face>
95{
96public:
97 using number_type = Number;
103 static constexpr unsigned int dimension = dim;
104 static constexpr unsigned int n_components = n_components_;
105 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
106
143 template <typename VectorType>
144 void
146 const VectorType &src,
147 const unsigned int first_index = 0,
148 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
149
178 template <typename VectorType>
179 void
181 const VectorType &src,
182 const unsigned int first_index = 0,
183 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip());
184
216 template <typename VectorType>
217 void
219 VectorType &dst,
220 const unsigned int first_index = 0,
221 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
222
261 template <typename VectorType>
262 void
264 VectorType &dst,
265 const unsigned int first_index = 0,
266 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
267
271 template <typename VectorType>
272 void
274 VectorType &dst,
275 const unsigned int first_index = 0,
276 const std::bitset<n_lanes> &mask = std::bitset<n_lanes>().flip()) const;
277
301 get_dof_value(const unsigned int dof) const;
302
313 void
314 submit_dof_value(const value_type val_in, const unsigned int dof);
315
329 get_value(const unsigned int q_point) const;
330
343 void
344 submit_value(const value_type val_in, const unsigned int q_point);
345
357 get_gradient(const unsigned int q_point) const;
358
374 get_normal_derivative(const unsigned int q_point) const;
375
388 void
389 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
390
409 void
411 const unsigned int q_point);
412
425 void
426 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
427
440 get_hessian(const unsigned int q_point) const;
441
452 get_hessian_diagonal(const unsigned int q_point) const;
453
466 get_laplacian(const unsigned int q_point) const;
467
468#ifdef DOXYGEN
469 // doxygen does not anyhow mention functions coming from partial template
470 // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
471 // For now, hack in those functions manually only to fix documentation:
472
479 VectorizedArrayType
480 get_divergence(const unsigned int q_point) const;
481
491 get_symmetric_gradient(const unsigned int q_point) const;
492
499 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
500 get_curl(const unsigned int q_point) const;
501
517 void
518 submit_divergence(const VectorizedArrayType div_in,
519 const unsigned int q_point);
520
537 void
540 const unsigned int q_point);
541
554 void
556 const unsigned int q_point);
557
558#endif
559
578
586
587protected:
598 const unsigned int dof_no,
599 const unsigned int first_selected_component,
600 const unsigned int quad_no,
601 const unsigned int fe_degree,
602 const unsigned int n_q_points,
603 const bool is_interior_face,
604 const unsigned int active_fe_index,
605 const unsigned int active_quad_index,
606 const unsigned int face_type);
607
645 const Mapping<dim> &mapping,
646 const FiniteElement<dim> &fe,
647 const Quadrature<1> &quadrature,
648 const UpdateFlags update_flags,
649 const unsigned int first_selected_component,
651
659
668
673
680 template <typename VectorType, typename VectorOperation>
681 void
683 const VectorOperation &operation,
684 const std::array<VectorType *, n_components_> &vectors,
685 const std::array<
687 n_components_> &vectors_sm,
688 const std::bitset<n_lanes> &mask,
689 const bool apply_constraints = true) const;
690
698 template <typename VectorType, typename VectorOperation>
699 void
701 const VectorOperation &operation,
702 const std::array<VectorType *, n_components_> &vectors,
703 const std::array<
705 n_components_> &vectors_sm,
706 const std::bitset<n_lanes> &mask) const;
707
715 template <typename VectorType, typename VectorOperation>
716 void
718 const VectorOperation &operation,
719 const std::array<VectorType *, n_components_> &vectors) const;
720
724 void
726
731
736
741 mutable std::vector<types::global_dof_index> local_dof_indices;
742};
743
744
745
753template <int dim,
754 int n_components_,
755 typename Number,
756 bool is_face,
757 typename VectorizedArrayType = VectorizedArray<Number>>
759 n_components_,
760 Number,
761 is_face,
762 VectorizedArrayType>
763{
764 static_assert(
765 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
766 "Type of Number and of VectorizedArrayType do not match.");
767
768public:
769 using number_type = Number;
773 static constexpr unsigned int dimension = dim;
774 static constexpr unsigned int n_components = n_components_;
775 using BaseClass =
777
778protected:
788 const unsigned int dof_no,
789 const unsigned int first_selected_component,
790 const unsigned int quad_no,
791 const unsigned int fe_degree,
792 const unsigned int n_q_points,
793 const bool is_interior_face = true,
796 const unsigned int face_type = numbers::invalid_unsigned_int);
797
803 const Mapping<dim> &mapping,
804 const FiniteElement<dim> &fe,
805 const Quadrature<1> &quadrature,
806 const UpdateFlags update_flags,
807 const unsigned int first_selected_component,
809
814
820};
821
822
823
832template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
833class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
834 : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
835{
836 static_assert(
837 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
838 "Type of Number and of VectorizedArrayType do not match.");
839
840public:
841 using number_type = Number;
842 using value_type = VectorizedArrayType;
845 static constexpr unsigned int dimension = dim;
846 using BaseClass =
848
853 get_dof_value(const unsigned int dof) const;
854
858 void
859 submit_dof_value(const value_type val_in, const unsigned int dof);
860
865 get_value(const unsigned int q_point) const;
866
870 void
871 submit_value(const value_type val_in, const unsigned int q_point);
872
876 void
878 const unsigned int q_point);
879
884 get_gradient(const unsigned int q_point) const;
885
890 get_normal_derivative(const unsigned int q_point) const;
891
895 void
896 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
897
901 void
903 const unsigned int q_point);
904
909 get_hessian(unsigned int q_point) const;
910
915 get_hessian_diagonal(const unsigned int q_point) const;
916
920 void
921 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
922
927 get_laplacian(const unsigned int q_point) const;
928
934
935protected:
945 const unsigned int dof_no,
946 const unsigned int first_selected_component,
947 const unsigned int quad_no,
948 const unsigned int fe_degree,
949 const unsigned int n_q_points,
950 const bool is_interior_face = true,
953 const unsigned int face_type = numbers::invalid_unsigned_int);
954
960 const Mapping<dim> &mapping,
961 const FiniteElement<dim> &fe,
962 const Quadrature<1> &quadrature,
963 const UpdateFlags update_flags,
964 const unsigned int first_selected_component,
966
971
977};
978
979
980
990template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
991class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
992 : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
993{
994 static_assert(
995 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
996 "Type of Number and of VectorizedArrayType do not match.");
997
998public:
999 using number_type = Number;
1002 static constexpr unsigned int dimension = dim;
1003 static constexpr unsigned int n_components = dim;
1006
1011 get_value(const unsigned int q_point) const;
1012
1017 get_gradient(const unsigned int q_point) const;
1018
1023 VectorizedArrayType
1024 get_divergence(const unsigned int q_point) const;
1025
1033 get_symmetric_gradient(const unsigned int q_point) const;
1034
1039 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1040 get_curl(const unsigned int q_point) const;
1041
1046 get_hessian(const unsigned int q_point) const;
1047
1052 get_hessian_diagonal(const unsigned int q_point) const;
1053
1057 void
1059 const unsigned int q_point);
1060
1064 void
1065 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1066
1075 void
1077 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1078 const unsigned int q_point);
1079
1088 void
1089 submit_divergence(const VectorizedArrayType div_in,
1090 const unsigned int q_point);
1091
1100 void
1103 const unsigned int q_point);
1104
1109 void
1111 const unsigned int q_point);
1112
1113protected:
1123 const unsigned int dof_no,
1124 const unsigned int first_selected_component,
1125 const unsigned int quad_no,
1126 const unsigned int dofs_per_cell,
1127 const unsigned int n_q_points,
1128 const bool is_interior_face = true,
1131 const unsigned int face_type = numbers::invalid_unsigned_int);
1132
1138 const Mapping<dim> &mapping,
1139 const FiniteElement<dim> &fe,
1140 const Quadrature<1> &quadrature,
1141 const UpdateFlags update_flags,
1142 const unsigned int first_selected_component,
1144
1149
1155};
1156
1157
1166template <typename Number, bool is_face, typename VectorizedArrayType>
1167class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1168 : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1169{
1170 static_assert(
1171 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1172 "Type of Number and of VectorizedArrayType do not match.");
1173
1174public:
1175 using number_type = Number;
1176 using value_type = VectorizedArrayType;
1179 static constexpr unsigned int dimension = 1;
1182
1187 get_dof_value(const unsigned int dof) const;
1188
1192 void
1193 submit_dof_value(const value_type val_in, const unsigned int dof);
1194
1199 get_value(const unsigned int q_point) const;
1200
1204 void
1205 submit_value(const value_type val_in, const unsigned int q_point);
1206
1210 void
1211 submit_value(const gradient_type val_in, const unsigned int q_point);
1212
1217 get_gradient(const unsigned int q_point) const;
1218
1223 get_divergence(const unsigned int q_point) const;
1224
1229 get_normal_derivative(const unsigned int q_point) const;
1230
1234 void
1235 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1236
1240 void
1241 submit_gradient(const value_type grad_in, const unsigned int q_point);
1242
1246 void
1248 const unsigned int q_point);
1249
1253 void
1255 const unsigned int q_point);
1256
1260 void
1262 const unsigned int q_point);
1263
1268 get_hessian(unsigned int q_point) const;
1269
1274 get_hessian_diagonal(const unsigned int q_point) const;
1275
1279 void
1280 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1281
1286 get_laplacian(const unsigned int q_point) const;
1287
1293
1294protected:
1304 const unsigned int dof_no,
1305 const unsigned int first_selected_component,
1306 const unsigned int quad_no,
1307 const unsigned int fe_degree,
1308 const unsigned int n_q_points,
1309 const bool is_interior_face = true,
1312 const unsigned int face_type = numbers::invalid_unsigned_int);
1313
1319 const Mapping<1> &mapping,
1320 const FiniteElement<1> &fe,
1321 const Quadrature<1> &quadrature,
1322 const UpdateFlags update_flags,
1323 const unsigned int first_selected_component,
1325
1330
1336};
1337
1338
1339
1903template <int dim,
1904 int fe_degree,
1905 int n_q_points_1d,
1906 int n_components_,
1907 typename Number,
1908 typename VectorizedArrayType>
1910 n_components_,
1911 Number,
1912 false,
1913 VectorizedArrayType>
1914{
1915 static_assert(
1916 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
1917 "Type of Number and of VectorizedArrayType do not match.");
1918
1919public:
1925
1929 using number_type = Number;
1930
1937
1944
1948 static constexpr unsigned int dimension = dim;
1949
1954 static constexpr unsigned int n_components = n_components_;
1955
1959 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
1960
1969 static constexpr unsigned int static_n_q_points =
1970 Utilities::pow(n_q_points_1d, dim);
1971
1981 static constexpr unsigned int static_dofs_per_component =
1982 Utilities::pow(fe_degree + 1, dim);
1983
1993 static constexpr unsigned int tensor_dofs_per_cell =
1995
2005 static constexpr unsigned int static_dofs_per_cell =
2007
2044 const unsigned int dof_no = 0,
2045 const unsigned int quad_no = 0,
2046 const unsigned int first_selected_component = 0,
2049
2058 const std::pair<unsigned int, unsigned int> &range,
2059 const unsigned int dof_no = 0,
2060 const unsigned int quad_no = 0,
2061 const unsigned int first_selected_component = 0);
2062
2092 const FiniteElement<dim> &fe,
2093 const Quadrature<1> &quadrature,
2094 const UpdateFlags update_flags,
2095 const unsigned int first_selected_component = 0);
2096
2103 const Quadrature<1> &quadrature,
2104 const UpdateFlags update_flags,
2105 const unsigned int first_selected_component = 0);
2106
2119 const unsigned int first_selected_component = 0);
2120
2128
2135 FEEvaluation &
2137
2146 void
2147 reinit(const unsigned int cell_batch_index);
2148
2155 void
2156 reinit(const std::array<unsigned int, n_lanes> &cell_ids);
2157
2170 template <bool level_dof_access>
2171 void
2173
2184 void
2186
2190 static bool
2191 fast_evaluation_supported(const unsigned int given_degree,
2192 const unsigned int given_n_q_points_1d);
2193
2203 void
2205
2218 void
2219 evaluate(const VectorizedArrayType *values_array,
2220 const EvaluationFlags::EvaluationFlags evaluation_flag);
2221
2235 template <typename VectorType>
2236 void
2237 gather_evaluate(const VectorType &input_vector,
2238 const EvaluationFlags::EvaluationFlags evaluation_flag);
2239
2250 void
2252
2264 void
2266 VectorizedArrayType *values_array,
2267 const bool sum_into_values = false);
2268
2282 template <typename VectorType>
2283 void
2285 VectorType &output_vector);
2286
2294
2301 const unsigned int dofs_per_component;
2302
2309 const unsigned int dofs_per_cell;
2310
2318 const unsigned int n_q_points;
2319
2320private:
2325 void
2326 check_template_arguments(const unsigned int fe_no,
2327 const unsigned int first_selected_component);
2328};
2329
2330
2331
2367template <int dim,
2368 int fe_degree,
2369 int n_q_points_1d = fe_degree + 1,
2370 int n_components_ = 1,
2371 typename Number = double,
2372 typename VectorizedArrayType = VectorizedArray<Number>>
2374 n_components_,
2375 Number,
2376 true,
2377 VectorizedArrayType>
2378{
2379 static_assert(
2380 std::is_same_v<Number, typename VectorizedArrayType::value_type>,
2381 "Type of Number and of VectorizedArrayType do not match.");
2382
2383public:
2389
2393 using number_type = Number;
2394
2401
2408
2412 static constexpr unsigned int dimension = dim;
2413
2418 static constexpr unsigned int n_components = n_components_;
2419
2423 static constexpr unsigned int n_lanes = VectorizedArrayType::size();
2424
2434 static constexpr unsigned int static_n_q_points =
2435 Utilities::pow(n_q_points_1d, dim - 1);
2436
2445 static constexpr unsigned int static_n_q_points_cell =
2446 Utilities::pow(n_q_points_1d, dim);
2447
2456 static constexpr unsigned int static_dofs_per_component =
2457 Utilities::pow(fe_degree + 1, dim);
2458
2467 static constexpr unsigned int tensor_dofs_per_cell =
2469
2478 static constexpr unsigned int static_dofs_per_cell =
2480
2524 const bool is_interior_face = true,
2525 const unsigned int dof_no = 0,
2526 const unsigned int quad_no = 0,
2527 const unsigned int first_selected_component = 0,
2530 const unsigned int face_type = numbers::invalid_unsigned_int);
2531
2541 const std::pair<unsigned int, unsigned int> &range,
2542 const bool is_interior_face = true,
2543 const unsigned int dof_no = 0,
2544 const unsigned int quad_no = 0,
2545 const unsigned int first_selected_component = 0);
2546
2557 void
2558 reinit(const unsigned int face_batch_number);
2559
2567 void
2568 reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2569
2573 static bool
2574 fast_evaluation_supported(const unsigned int given_degree,
2575 const unsigned int given_n_q_points_1d);
2576
2587 void
2589
2602 void
2603 evaluate(const VectorizedArrayType *values_array,
2604 const EvaluationFlags::EvaluationFlags evaluation_flag);
2605
2610 void
2612
2617 void
2618 project_to_face(const VectorizedArrayType *values_array,
2619 const EvaluationFlags::EvaluationFlags evaluation_flag);
2620
2625 void
2627
2639 template <typename VectorType>
2640 void
2641 gather_evaluate(const VectorType &input_vector,
2642 const EvaluationFlags::EvaluationFlags evaluation_flag);
2643
2653 void
2655 const bool sum_into_values = false);
2656
2665 void
2667 VectorizedArrayType *values_array,
2668 const bool sum_into_values = false);
2669
2676 void
2678
2683 void
2685 const bool sum_into_values = false);
2686
2691 void
2693 VectorizedArrayType *values_array,
2694 const bool sum_into_values = false);
2695
2707 template <typename VectorType>
2708 void
2710 VectorType &output_vector);
2711
2715 template <typename VectorType>
2716 void
2717 integrate_scatter(const bool integrate_values,
2718 const bool integrate_gradients,
2719 VectorType &output_vector);
2720
2728
2733 bool
2735
2750
2755 unsigned int
2757
2762 unsigned int
2764
2771 const unsigned int dofs_per_component;
2772
2779 const unsigned int dofs_per_cell;
2780
2788 const unsigned int n_q_points;
2789};
2790
2791
2792
2793namespace internal
2794{
2795 namespace MatrixFreeFunctions
2796 {
2797 // a helper function to compute the number of DoFs of a DGP element at
2798 // compile time, depending on the degree
2799 template <int dim, int degree>
2801 {
2802 // this division is always without remainder
2803 static constexpr unsigned int value =
2804 (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2805 };
2806
2807 // base specialization: 1d elements have 'degree+1' degrees of freedom
2808 template <int degree>
2809 struct DGP_dofs_per_component<1, degree>
2810 {
2811 static constexpr unsigned int value = degree + 1;
2812 };
2813 } // namespace MatrixFreeFunctions
2814} // namespace internal
2815
2816
2817/*----------------------- Inline functions ----------------------------------*/
2818
2819#ifndef DOXYGEN
2820
2821
2822namespace internal
2823{
2824 // Extract all internal data pointers and indices in a single function that
2825 // get passed on to the constructor of FEEvaluationData, avoiding to look
2826 // things up multiple times
2827 template <bool is_face,
2828 int dim,
2829 typename Number,
2830 typename VectorizedArrayType>
2832 InitializationData
2833 extract_initialization_data(
2835 const unsigned int dof_no,
2836 const unsigned int first_selected_component,
2837 const unsigned int quad_no,
2838 const unsigned int fe_degree,
2839 const unsigned int n_q_points,
2840 const unsigned int active_fe_index_given,
2841 const unsigned int active_quad_index_given,
2842 const unsigned int face_type)
2843 {
2845 InitializationData init_data;
2846
2847 init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2848 init_data.mapping_data =
2849 &internal::MatrixFreeFunctions::
2850 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2851 matrix_free.get_mapping_info(), quad_no);
2852
2853 init_data.active_fe_index =
2854 fe_degree != numbers::invalid_unsigned_int ?
2855 init_data.dof_info->fe_index_from_degree(first_selected_component,
2856 fe_degree) :
2857 (active_fe_index_given != numbers::invalid_unsigned_int ?
2858 active_fe_index_given :
2859 0);
2860
2861 init_data.active_quad_index =
2862 fe_degree == numbers::invalid_unsigned_int ?
2863 (active_quad_index_given != numbers::invalid_unsigned_int ?
2864 active_quad_index_given :
2865 std::min<unsigned int>(
2866 init_data.active_fe_index,
2867 init_data.mapping_data->descriptor.size() /
2868 (is_face ? std::max<unsigned int>(1, dim - 1) : 1) -
2869 1)) :
2870 init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2871
2872 init_data.shape_info = &matrix_free.get_shape_info(
2873 dof_no,
2874 quad_no,
2875 init_data.dof_info->component_to_base_index[first_selected_component],
2876 init_data.active_fe_index,
2877 init_data.active_quad_index);
2878 init_data.descriptor =
2879 &init_data.mapping_data->descriptor
2880 [is_face ?
2881 (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2882 (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2883 init_data.active_quad_index];
2884
2885 return init_data;
2886 }
2887} // namespace internal
2888
2889
2890
2891/*----------------------- FEEvaluationBase ----------------------------------*/
2892
2893template <int dim,
2894 int n_components_,
2895 typename Number,
2896 bool is_face,
2897 typename VectorizedArrayType>
2898inline FEEvaluationBase<dim,
2899 n_components_,
2900 Number,
2901 is_face,
2902 VectorizedArrayType>::
2903 FEEvaluationBase(
2905 const unsigned int dof_no,
2906 const unsigned int first_selected_component,
2907 const unsigned int quad_no,
2908 const unsigned int fe_degree,
2909 const unsigned int n_q_points,
2910 const bool is_interior_face,
2911 const unsigned int active_fe_index,
2912 const unsigned int active_quad_index,
2913 const unsigned int face_type)
2914 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2915 internal::extract_initialization_data<is_face>(matrix_free,
2916 dof_no,
2917 first_selected_component,
2918 quad_no,
2919 fe_degree,
2920 n_q_points,
2921 active_fe_index,
2922 active_quad_index,
2923 face_type),
2924 is_interior_face,
2925 quad_no,
2926 first_selected_component)
2927 , scratch_data_array(matrix_free.acquire_scratch_data())
2928 , matrix_free(&matrix_free)
2929{
2930 this->set_data_pointers(scratch_data_array, n_components_);
2931 Assert(
2932 this->dof_info->start_components.back() == 1 ||
2933 static_cast<int>(n_components_) <=
2934 static_cast<int>(
2935 this->dof_info->start_components
2936 [this->dof_info->component_to_base_index[first_selected_component] +
2937 1]) -
2938 first_selected_component,
2939 ExcMessage(
2940 "You tried to construct a vector-valued evaluator with " +
2941 std::to_string(n_components) +
2942 " components. However, "
2943 "the current base element has only " +
2944 std::to_string(
2945 this->dof_info->start_components
2946 [this->dof_info->component_to_base_index[first_selected_component] +
2947 1] -
2948 first_selected_component) +
2949 " components left when starting from local element index " +
2950 std::to_string(
2951 first_selected_component -
2952 this->dof_info->start_components
2953 [this->dof_info->component_to_base_index[first_selected_component]]) +
2954 " (global index " + std::to_string(first_selected_component) + ")"));
2955
2956 // do not check for correct dimensions of data fields here, should be done
2957 // in derived classes
2958}
2959
2960
2961
2962template <int dim,
2963 int n_components_,
2964 typename Number,
2965 bool is_face,
2966 typename VectorizedArrayType>
2967inline FEEvaluationBase<dim,
2968 n_components_,
2969 Number,
2970 is_face,
2971 VectorizedArrayType>::
2972 FEEvaluationBase(
2973 const Mapping<dim> &mapping,
2974 const FiniteElement<dim> &fe,
2975 const Quadrature<1> &quadrature,
2976 const UpdateFlags update_flags,
2977 const unsigned int first_selected_component,
2979 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2980 other != nullptr &&
2981 other->mapped_geometry->get_quadrature() == quadrature ?
2982 other->mapped_geometry :
2983 std::make_shared<internal::MatrixFreeFunctions::
2984 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2985 mapping,
2986 quadrature,
2987 update_flags),
2988 n_components_,
2989 first_selected_component)
2990 , scratch_data_array(new AlignedVector<VectorizedArrayType>())
2991 , matrix_free(nullptr)
2992{
2993 const unsigned int base_element_number =
2994 fe.component_to_base_index(first_selected_component).first;
2995 Assert(fe.element_multiplicity(base_element_number) == 1 ||
2996 fe.element_multiplicity(base_element_number) -
2997 first_selected_component >=
2998 n_components_,
2999 ExcMessage("The underlying element must at least contain as many "
3000 "components as requested by this class"));
3001 (void)base_element_number;
3002
3003 Assert(this->data == nullptr, ExcInternalError());
3004 this->data =
3006 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
3007 fe,
3008 fe.component_to_base_index(first_selected_component).first);
3009
3010 this->set_data_pointers(scratch_data_array, n_components_);
3011}
3012
3013
3014
3015template <int dim,
3016 int n_components_,
3017 typename Number,
3018 bool is_face,
3019 typename VectorizedArrayType>
3020inline FEEvaluationBase<dim,
3021 n_components_,
3022 Number,
3023 is_face,
3024 VectorizedArrayType>::
3025 FEEvaluationBase(const FEEvaluationBase<dim,
3026 n_components_,
3027 Number,
3028 is_face,
3029 VectorizedArrayType> &other)
3030 : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
3031 , scratch_data_array(other.matrix_free == nullptr ?
3032 new AlignedVector<VectorizedArrayType>() :
3033 other.matrix_free->acquire_scratch_data())
3034 , matrix_free(other.matrix_free)
3035{
3036 if (other.matrix_free == nullptr)
3037 {
3038 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3039 this->data =
3041 *other.data);
3042
3043 // Create deep copy of mapped geometry for use in parallel
3044 this->mapped_geometry =
3045 std::make_shared<internal::MatrixFreeFunctions::
3046 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3047 other.mapped_geometry->get_fe_values().get_mapping(),
3048 other.mapped_geometry->get_quadrature(),
3049 other.mapped_geometry->get_fe_values().get_update_flags());
3050 this->mapping_data = &this->mapped_geometry->get_data_storage();
3051 this->cell = 0;
3052
3053 this->jacobian =
3054 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3055 this->J_value =
3056 this->mapped_geometry->get_data_storage().JxW_values.begin();
3057 this->jacobian_gradients =
3058 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3059 this->jacobian_gradients_non_inverse =
3060 this->mapped_geometry->get_data_storage()
3061 .jacobian_gradients_non_inverse[0]
3062 .begin();
3063 this->quadrature_points =
3064 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3065 }
3066
3067 this->set_data_pointers(scratch_data_array, n_components_);
3068}
3069
3070
3071
3072template <int dim,
3073 int n_components_,
3074 typename Number,
3075 bool is_face,
3076 typename VectorizedArrayType>
3077inline FEEvaluationBase<dim,
3078 n_components_,
3079 Number,
3080 is_face,
3081 VectorizedArrayType> &
3083operator=(const FEEvaluationBase<dim,
3084 n_components_,
3085 Number,
3086 is_face,
3087 VectorizedArrayType> &other)
3088{
3089 // release old memory
3090 if (matrix_free == nullptr)
3091 {
3092 delete this->data;
3093 delete scratch_data_array;
3094 }
3095 else
3096 {
3097 matrix_free->release_scratch_data(scratch_data_array);
3098 }
3099
3101
3102 matrix_free = other.matrix_free;
3103
3104 if (other.matrix_free == nullptr)
3105 {
3106 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3107 this->data =
3109 *other.data);
3110 scratch_data_array = new AlignedVector<VectorizedArrayType>();
3111
3112 // Create deep copy of mapped geometry for use in parallel
3113 this->mapped_geometry =
3114 std::make_shared<internal::MatrixFreeFunctions::
3115 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3116 other.mapped_geometry->get_fe_values().get_mapping(),
3117 other.mapped_geometry->get_quadrature(),
3118 other.mapped_geometry->get_fe_values().get_update_flags());
3119 this->cell = 0;
3120 this->mapping_data = &this->mapped_geometry->get_data_storage();
3121 this->jacobian =
3122 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3123 this->J_value =
3124 this->mapped_geometry->get_data_storage().JxW_values.begin();
3125 this->jacobian_gradients =
3126 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3127 this->jacobian_gradients_non_inverse =
3128 this->mapped_geometry->get_data_storage()
3129 .jacobian_gradients_non_inverse[0]
3130 .begin();
3131 this->quadrature_points =
3132 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3133 }
3134 else
3135 {
3136 scratch_data_array = matrix_free->acquire_scratch_data();
3137 }
3138
3139 this->set_data_pointers(scratch_data_array, n_components_);
3140
3141 return *this;
3142}
3143
3144
3145
3146template <int dim,
3147 int n_components_,
3148 typename Number,
3149 bool is_face,
3150 typename VectorizedArrayType>
3151inline FEEvaluationBase<dim,
3152 n_components_,
3153 Number,
3154 is_face,
3155 VectorizedArrayType>::~FEEvaluationBase()
3156{
3157 if (matrix_free != nullptr)
3158 {
3159 try
3160 {
3161 matrix_free->release_scratch_data(scratch_data_array);
3162 }
3163 catch (...)
3164 {}
3165 }
3166 else
3167 {
3168 delete scratch_data_array;
3169 delete this->data;
3170 }
3171}
3172
3173
3174
3175template <int dim,
3176 int n_components_,
3177 typename Number,
3178 bool is_face,
3179 typename VectorizedArrayType>
3182 get_matrix_free() const
3183{
3184 Assert(matrix_free != nullptr,
3185 ExcMessage(
3186 "FEEvaluation was not initialized with a MatrixFree object!"));
3187 return *matrix_free;
3188}
3189
3190
3191
3192namespace internal
3193{
3194 // given a block vector return the underlying vector type
3195 // including constness (specified by bool)
3196 template <typename VectorType, bool>
3197 struct ConstBlockVectorSelector;
3198
3199 template <typename VectorType>
3200 struct ConstBlockVectorSelector<VectorType, true>
3201 {
3202 using BaseVectorType = const typename VectorType::BlockType;
3203 };
3204
3205 template <typename VectorType>
3206 struct ConstBlockVectorSelector<VectorType, false>
3207 {
3208 using BaseVectorType = typename VectorType::BlockType;
3209 };
3210
3211 // allows to select between block vectors and non-block vectors, which
3212 // allows to use a unified interface for extracting blocks on block vectors
3213 // and doing nothing on usual vectors
3214 template <typename VectorType, bool>
3215 struct BlockVectorSelector;
3216
3217 template <typename VectorType>
3218 struct BlockVectorSelector<VectorType, true>
3219 {
3220 using BaseVectorType = typename ConstBlockVectorSelector<
3221 VectorType,
3222 std::is_const_v<VectorType>>::BaseVectorType;
3223
3224 static BaseVectorType *
3225 get_vector_component(VectorType &vec, const unsigned int component)
3226 {
3227 AssertIndexRange(component, vec.n_blocks());
3228 return &vec.block(component);
3229 }
3230 };
3231
3232 template <typename VectorType>
3233 struct BlockVectorSelector<VectorType, false>
3234 {
3235 using BaseVectorType = VectorType;
3236
3237 static BaseVectorType *
3238 get_vector_component(VectorType &vec, const unsigned int component)
3239 {
3240 // FEEvaluation allows to combine several vectors from a scalar
3241 // FiniteElement into a "vector-valued" FEEvaluation object with
3242 // multiple components. These components can be extracted with the other
3243 // get_vector_component functions. If we do not get a vector of vectors
3244 // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3245 // must make sure that we do not duplicate the components in input
3246 // and/or duplicate the resulting integrals. In such a case, we should
3247 // only get the zeroth component in the vector contained set nullptr for
3248 // the others which allows us to catch unintended use in
3249 // read_write_operation.
3250 if (component == 0)
3251 return &vec;
3252 else
3253 return nullptr;
3254 }
3255 };
3256
3257 template <typename VectorType>
3258 struct BlockVectorSelector<std::vector<VectorType>, false>
3259 {
3260 using BaseVectorType = VectorType;
3261
3262 static BaseVectorType *
3263 get_vector_component(std::vector<VectorType> &vec,
3264 const unsigned int component)
3265 {
3266 AssertIndexRange(component, vec.size());
3267 return &vec[component];
3268 }
3269 };
3270
3271 template <typename VectorType>
3272 struct BlockVectorSelector<const std::vector<VectorType>, false>
3273 {
3274 using BaseVectorType = const VectorType;
3275
3276 static const BaseVectorType *
3277 get_vector_component(const std::vector<VectorType> &vec,
3278 const unsigned int component)
3279 {
3280 AssertIndexRange(component, vec.size());
3281 return &vec[component];
3282 }
3283 };
3284
3285 template <typename VectorType>
3286 struct BlockVectorSelector<std::vector<VectorType *>, false>
3287 {
3288 using BaseVectorType = VectorType;
3289
3290 static BaseVectorType *
3291 get_vector_component(std::vector<VectorType *> &vec,
3292 const unsigned int component)
3293 {
3294 AssertIndexRange(component, vec.size());
3295 return vec[component];
3296 }
3297 };
3298
3299 template <typename VectorType>
3300 struct BlockVectorSelector<const std::vector<VectorType *>, false>
3301 {
3302 using BaseVectorType = const VectorType;
3303
3304 static const BaseVectorType *
3305 get_vector_component(const std::vector<VectorType *> &vec,
3306 const unsigned int component)
3307 {
3308 AssertIndexRange(component, vec.size());
3309 return vec[component];
3310 }
3311 };
3312} // namespace internal
3313
3314
3315
3316template <int dim,
3317 int n_components_,
3318 typename Number,
3319 bool is_face,
3320 typename VectorizedArrayType>
3321template <typename VectorType, typename VectorOperation>
3322inline void
3325 const VectorOperation &operation,
3326 const std::array<VectorType *, n_components_> &src,
3327 const std::array<
3329 n_components_> &src_sm,
3330 const std::bitset<n_lanes> &mask,
3331 const bool apply_constraints) const
3332{
3333 // Case 1: No MatrixFree object given, simple case because we do not need to
3334 // process constraints and need not care about vectorization -> go to
3335 // separate function
3336 if (this->matrix_free == nullptr)
3337 {
3338 read_write_operation_global(operation, src);
3339 return;
3340 }
3341
3342 Assert(this->dof_info != nullptr, ExcNotInitialized());
3343 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3344 Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3345 if (this->n_fe_components == 1)
3346 for (unsigned int comp = 0; comp < n_components; ++comp)
3347 {
3348 Assert(src[comp] != nullptr,
3349 ExcMessage("The finite element underlying this FEEvaluation "
3350 "object is scalar, but you requested " +
3351 std::to_string(n_components) +
3352 " components via the template argument in "
3353 "FEEvaluation. In that case, you must pass an "
3354 "std::vector<VectorType> or a BlockVector to " +
3355 "read_dof_values and distribute_local_to_global."));
3357 *this->matrix_free,
3358 *this->dof_info);
3359 }
3360 else
3361 {
3363 *this->matrix_free,
3364 *this->dof_info);
3365 }
3366
3367 const bool accesses_exterior_dofs =
3368 this->dof_access_index ==
3370 this->is_interior_face() == false;
3371
3372 // Case 2: contiguous indices which use reduced storage of indices and can
3373 // use vectorized load/store operations -> go to separate function
3374 if (this->cell != numbers::invalid_unsigned_int)
3375 {
3377 this->cell,
3378 dof_info.index_storage_variants[this->dof_access_index].size());
3379
3380 bool is_contiguous = true;
3381 // check if exterior cells are not contiguous (ECL case)
3382 if (accesses_exterior_dofs)
3383 {
3384 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
3385 const unsigned int n_filled_lanes =
3388 [this->cell];
3389 // we have to check all filled lanes which are active in the mask
3390 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3391 if (mask[v] == true &&
3392 dof_info.index_storage_variants
3394 [cells[v] / n_lanes] <
3396 contiguous)
3397 is_contiguous = false;
3398 } // or if cell/face batch is not contiguous
3399 else if (dof_info.index_storage_variants
3400 [is_face ?
3401 this->dof_access_index :
3402 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
3403 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
3404 IndexStorageVariants::contiguous)
3405 {
3406 is_contiguous = false;
3407 }
3408
3409 if (is_contiguous)
3410 {
3411 read_write_operation_contiguous(operation, src, src_sm, mask);
3412 return;
3413 }
3414 }
3415
3416 // Case 3: standard operation with one index per degree of freedom -> go on
3417 // here
3418 std::array<unsigned int, n_lanes> cells = this->get_cell_ids();
3419
3420 const bool masking_is_active = mask.count() < n_lanes;
3421 if (masking_is_active)
3422 for (unsigned int v = 0; v < n_lanes; ++v)
3423 if (mask[v] == false)
3425
3426 bool has_hn_constraints = false;
3427
3428 if (is_face == false)
3429 {
3430 if (!dof_info.hanging_node_constraint_masks.empty() &&
3431 !dof_info.hanging_node_constraint_masks_comp.empty() &&
3432 dof_info
3433 .hanging_node_constraint_masks_comp[this->active_fe_index]
3434 [this->first_selected_component])
3435 for (unsigned int v = 0; v < n_lanes; ++v)
3436 if (cells[v] != numbers::invalid_unsigned_int &&
3437 dof_info.hanging_node_constraint_masks[cells[v]] !=
3440 has_hn_constraints = true;
3441 }
3442
3443 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3444 vector_selector;
3445
3446 const bool use_vectorized_path =
3447 !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
3448
3449 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3450 std::array<VectorizedArrayType *, n_components> values_dofs;
3451 for (unsigned int c = 0; c < n_components; ++c)
3452 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3453 c * dofs_per_component;
3454
3455 if (this->cell != numbers::invalid_unsigned_int &&
3456 dof_info.index_storage_variants
3457 [is_face ? this->dof_access_index :
3458 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
3459 [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3460 IndexStorageVariants::interleaved &&
3461 use_vectorized_path)
3462 {
3463 const unsigned int *dof_indices =
3464 dof_info.dof_indices_interleaved.data() +
3465 dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
3466 .first +
3467 this->dof_info
3468 ->component_dof_indices_offset[this->active_fe_index]
3469 [this->first_selected_component] *
3470 n_lanes;
3471
3472 std::array<typename VectorType::value_type *, n_components> src_ptrs;
3473 if (n_components == 1 || this->n_fe_components == 1)
3474 for (unsigned int comp = 0; comp < n_components; ++comp)
3475 src_ptrs[comp] =
3476 const_cast<typename VectorType::value_type *>(src[comp]->begin());
3477 else
3478 src_ptrs[0] =
3479 const_cast<typename VectorType::value_type *>(src[0]->begin());
3480
3481 if (n_components == 1 || this->n_fe_components == 1)
3482 for (unsigned int i = 0; i < dofs_per_component;
3483 ++i, dof_indices += n_lanes)
3484 for (unsigned int comp = 0; comp < n_components; ++comp)
3485 operation.process_dof_gather(dof_indices,
3486 *src[comp],
3487 0,
3488 src_ptrs[comp],
3489 values_dofs[comp][i],
3490 vector_selector);
3491 else
3492 for (unsigned int comp = 0; comp < n_components; ++comp)
3493 for (unsigned int i = 0; i < dofs_per_component;
3494 ++i, dof_indices += n_lanes)
3495 operation.process_dof_gather(dof_indices,
3496 *src[0],
3497 0,
3498 src_ptrs[0],
3499 values_dofs[comp][i],
3500 vector_selector);
3501 return;
3502 }
3503
3504 // Allocate pointers, then initialize all of them to nullptrs and
3505 // below overwrite the ones we actually use:
3506 std::array<const unsigned int *, n_lanes> dof_indices;
3507 dof_indices.fill(nullptr);
3508
3509 // Assign the appropriate cell ids for face/cell case and get the pointers
3510 // to the dof indices of the cells on all lanes
3511
3512 bool has_constraints = false;
3513 const unsigned int n_components_read =
3514 this->n_fe_components > 1 ? n_components : 1;
3515
3516 if (is_face)
3517 {
3518 for (unsigned int v = 0; v < n_lanes; ++v)
3519 {
3520 if (cells[v] == numbers::invalid_unsigned_int)
3521 continue;
3522
3523 Assert(cells[v] < dof_info.row_starts.size() - 1, ExcInternalError());
3524 const std::pair<unsigned int, unsigned int> *my_index_start =
3525 &dof_info.row_starts[cells[v] * this->n_fe_components +
3526 this->first_selected_component];
3527
3528 // check whether any of the SIMD lanes has constraints, i.e., the
3529 // constraint indicator which is the second entry of row_starts
3530 // increments on this cell
3531 if (my_index_start[n_components_read].second !=
3532 my_index_start[0].second)
3533 has_constraints = true;
3534
3535 dof_indices[v] =
3536 dof_info.dof_indices.data() + my_index_start[0].first;
3537 }
3538 }
3539 else
3540 {
3541 for (unsigned int v = 0; v < n_lanes; ++v)
3542 {
3543 if (cells[v] == numbers::invalid_unsigned_int)
3544 continue;
3545
3546 const std::pair<unsigned int, unsigned int> *my_index_start =
3547 &dof_info.row_starts[cells[v] * this->n_fe_components +
3548 this->first_selected_component];
3549 if (my_index_start[n_components_read].second !=
3550 my_index_start[0].second)
3551 has_constraints = true;
3552
3553 if (dof_info.hanging_node_constraint_masks.size() > 0 &&
3554 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3555 dof_info.hanging_node_constraint_masks[cells[v]] !=
3558 dof_info.hanging_node_constraint_masks_comp
3559 [this->active_fe_index][this->first_selected_component])
3560 has_hn_constraints = true;
3561
3562 Assert(my_index_start[n_components_read].first ==
3563 my_index_start[0].first ||
3564 my_index_start[0].first < dof_info.dof_indices.size(),
3565 ExcIndexRange(0,
3566 my_index_start[0].first,
3567 dof_info.dof_indices.size()));
3568 dof_indices[v] =
3569 dof_info.dof_indices.data() + my_index_start[0].first;
3570 }
3571 }
3572
3573 if (std::count_if(cells.begin(), cells.end(), [](const auto i) {
3574 return i != numbers::invalid_unsigned_int;
3575 }) < n_lanes)
3576 for (unsigned int comp = 0; comp < n_components; ++comp)
3577 for (unsigned int i = 0; i < dofs_per_component; ++i)
3578 operation.process_empty(values_dofs[comp][i]);
3579
3580 // Case where we have no constraints throughout the whole cell: Can go
3581 // through the list of DoFs directly
3582 if (!has_constraints && apply_constraints)
3583 {
3584 if (n_components == 1 || this->n_fe_components == 1)
3585 {
3586 for (unsigned int v = 0; v < n_lanes; ++v)
3587 {
3588 if (cells[v] == numbers::invalid_unsigned_int)
3589 continue;
3590
3591 for (unsigned int i = 0; i < dofs_per_component; ++i)
3592 for (unsigned int comp = 0; comp < n_components; ++comp)
3593 operation.process_dof(dof_indices[v][i],
3594 *src[comp],
3595 values_dofs[comp][i][v]);
3596 }
3597 }
3598 else
3599 {
3600 for (unsigned int comp = 0; comp < n_components; ++comp)
3601 for (unsigned int v = 0; v < n_lanes; ++v)
3602 {
3603 if (cells[v] == numbers::invalid_unsigned_int)
3604 continue;
3605
3606 for (unsigned int i = 0; i < dofs_per_component; ++i)
3607 operation.process_dof(
3608 dof_indices[v][comp * dofs_per_component + i],
3609 *src[0],
3610 values_dofs[comp][i][v]);
3611 }
3612 }
3613 return;
3614 }
3615
3616 // In the case where there are some constraints to be resolved, loop over
3617 // all vector components that are filled and then over local dofs. ind_local
3618 // holds local number on cell, index iterates over the elements of
3619 // index_local_to_global and dof_indices points to the global indices stored
3620 // in index_local_to_global
3621
3622 for (unsigned int v = 0; v < n_lanes; ++v)
3623 {
3624 if (cells[v] == numbers::invalid_unsigned_int)
3625 continue;
3626
3627 const unsigned int cell_index = cells[v];
3628 const unsigned int cell_dof_index =
3629 cell_index * this->n_fe_components + this->first_selected_component;
3630 const unsigned int n_components_read =
3631 this->n_fe_components > 1 ? n_components : 1;
3632 unsigned int index_indicators =
3633 dof_info.row_starts[cell_dof_index].second;
3634 unsigned int next_index_indicators =
3635 dof_info.row_starts[cell_dof_index + 1].second;
3636
3637 // For read_dof_values_plain, redirect the dof_indices field to the
3638 // unconstrained indices
3639 if (apply_constraints == false &&
3640 (dof_info.row_starts[cell_dof_index].second !=
3641 dof_info.row_starts[cell_dof_index + n_components_read].second ||
3642 ((dof_info.hanging_node_constraint_masks.size() > 0 &&
3643 dof_info.hanging_node_constraint_masks_comp.size() > 0 &&
3647 dof_info.hanging_node_constraint_masks_comp
3648 [this->active_fe_index][this->first_selected_component])))
3649 {
3653 dof_indices[v] =
3654 dof_info.plain_dof_indices.data() +
3655 this->dof_info
3656 ->component_dof_indices_offset[this->active_fe_index]
3657 [this->first_selected_component] +
3659 next_index_indicators = index_indicators;
3660 }
3661
3662 if (n_components == 1 || this->n_fe_components == 1)
3663 {
3664 unsigned int ind_local = 0;
3665 for (; index_indicators != next_index_indicators; ++index_indicators)
3666 {
3667 const std::pair<unsigned short, unsigned short> indicator =
3668 dof_info.constraint_indicator[index_indicators];
3669 // run through values up to next constraint
3670 for (unsigned int j = 0; j < indicator.first; ++j)
3671 for (unsigned int comp = 0; comp < n_components; ++comp)
3672 operation.process_dof(dof_indices[v][j],
3673 *src[comp],
3674 values_dofs[comp][ind_local + j][v]);
3675
3676 ind_local += indicator.first;
3677 dof_indices[v] += indicator.first;
3678
3679 // constrained case: build the local value as a linear
3680 // combination of the global value according to constraints
3681 Number value[n_components];
3682 for (unsigned int comp = 0; comp < n_components; ++comp)
3683 operation.pre_constraints(values_dofs[comp][ind_local][v],
3684 value[comp]);
3685
3686 const Number *data_val =
3687 this->matrix_free->constraint_pool_begin(indicator.second);
3688 const Number *end_pool =
3689 this->matrix_free->constraint_pool_end(indicator.second);
3690 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3691 for (unsigned int comp = 0; comp < n_components; ++comp)
3692 operation.process_constraint(*dof_indices[v],
3693 *data_val,
3694 *src[comp],
3695 value[comp]);
3696
3697 for (unsigned int comp = 0; comp < n_components; ++comp)
3698 operation.post_constraints(value[comp],
3699 values_dofs[comp][ind_local][v]);
3700 ++ind_local;
3701 }
3702
3703 AssertIndexRange(ind_local, dofs_per_component + 1);
3704
3705 for (; ind_local < dofs_per_component; ++dof_indices[v], ++ind_local)
3706 for (unsigned int comp = 0; comp < n_components; ++comp)
3707 operation.process_dof(*dof_indices[v],
3708 *src[comp],
3709 values_dofs[comp][ind_local][v]);
3710 }
3711 else
3712 {
3713 // case with vector-valued finite elements where all components are
3714 // included in one single vector. Assumption: first come all entries
3715 // to the first component, then all entries to the second one, and
3716 // so on. This is ensured by the way MatrixFree reads out the
3717 // indices.
3718 for (unsigned int comp = 0; comp < n_components; ++comp)
3719 {
3720 unsigned int ind_local = 0;
3721
3722 // check whether there is any constraint on the current cell
3723 for (; index_indicators != next_index_indicators;
3724 ++index_indicators)
3725 {
3726 const std::pair<unsigned short, unsigned short> indicator =
3727 dof_info.constraint_indicator[index_indicators];
3728
3729 // run through values up to next constraint
3730 for (unsigned int j = 0; j < indicator.first; ++j)
3731 operation.process_dof(dof_indices[v][j],
3732 *src[0],
3733 values_dofs[comp][ind_local + j][v]);
3734 ind_local += indicator.first;
3735 dof_indices[v] += indicator.first;
3736
3737 // constrained case: build the local value as a linear
3738 // combination of the global value according to constraints
3739 Number value;
3740 operation.pre_constraints(values_dofs[comp][ind_local][v],
3741 value);
3742
3743 const Number *data_val =
3744 this->matrix_free->constraint_pool_begin(indicator.second);
3745 const Number *end_pool =
3746 this->matrix_free->constraint_pool_end(indicator.second);
3747
3748 for (; data_val != end_pool; ++data_val, ++dof_indices[v])
3749 operation.process_constraint(*dof_indices[v],
3750 *data_val,
3751 *src[0],
3752 value);
3753
3754 operation.post_constraints(value,
3755 values_dofs[comp][ind_local][v]);
3756 ++ind_local;
3757 }
3758
3759 AssertIndexRange(ind_local, dofs_per_component + 1);
3760
3761 // get the dof values past the last constraint
3762 for (; ind_local < dofs_per_component;
3763 ++dof_indices[v], ++ind_local)
3764 {
3765 AssertIndexRange(*dof_indices[v], src[0]->size());
3766 operation.process_dof(*dof_indices[v],
3767 *src[0],
3768 values_dofs[comp][ind_local][v]);
3769 }
3770
3771 if (apply_constraints == true && comp + 1 < n_components)
3772 next_index_indicators =
3773 dof_info.row_starts[cell_dof_index + comp + 2].second;
3774 }
3775 }
3776 }
3777}
3778
3779
3780
3781template <int dim,
3782 int n_components_,
3783 typename Number,
3784 bool is_face,
3785 typename VectorizedArrayType>
3786template <typename VectorType, typename VectorOperation>
3787inline void
3790 const VectorOperation &operation,
3791 const std::array<VectorType *, n_components_> &src) const
3792{
3793 Assert(!local_dof_indices.empty(), ExcNotInitialized());
3794
3795 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3796 unsigned int index = this->first_selected_component * dofs_per_component;
3797 for (unsigned int comp = 0; comp < n_components; ++comp)
3798 {
3799 for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
3800 {
3801 operation.process_empty(
3802 this->values_dofs[comp * dofs_per_component + i]);
3803 operation.process_dof_global(
3804 local_dof_indices[this->data->lexicographic_numbering[index]],
3805 *src[0],
3806 this->values_dofs[comp * dofs_per_component + i][0]);
3807 }
3808 }
3809}
3810
3811
3812
3813template <int dim,
3814 int n_components_,
3815 typename Number,
3816 bool is_face,
3817 typename VectorizedArrayType>
3818template <typename VectorType, typename VectorOperation>
3819inline void
3822 const VectorOperation &operation,
3823 const std::array<VectorType *, n_components_> &src,
3824 const std::array<
3826 n_components_> &vectors_sm,
3827 const std::bitset<n_lanes> &mask) const
3828{
3829 // This functions processes the functions read_dof_values,
3830 // distribute_local_to_global, and set_dof_values with the same code for
3831 // contiguous cell indices (DG case). The distinction between these three
3832 // cases is made by the input VectorOperation that either reads values from
3833 // a vector and puts the data into the local data field or write local data
3834 // into the vector. Certain operations are no-ops for the given use case.
3835
3836 std::bool_constant<internal::is_vectorizable<VectorType, Number>::value>
3837 vector_selector;
3839 is_face ? this->dof_access_index :
3841 const unsigned int n_active_lanes = mask.count();
3842
3843 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3844 const std::vector<unsigned int> &dof_indices_cont =
3845 dof_info.dof_indices_contiguous[ind];
3846
3847 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3848 std::array<VectorizedArrayType *, n_components> values_dofs;
3849 for (unsigned int c = 0; c < n_components; ++c)
3850 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3851 c * dofs_per_component;
3852
3854
3855 const bool accesses_exterior_dofs =
3856 this->dof_access_index ==
3858 this->is_interior_face() == false;
3859
3860 // Simple case: We have contiguous storage, so we can simply copy out the
3861 // data
3862 if (dof_info.index_storage_variants[ind][this->cell] ==
3864 interleaved_contiguous &&
3865 n_active_lanes == n_lanes && !accesses_exterior_dofs)
3866 {
3867 const unsigned int dof_index =
3868 dof_indices_cont[this->cell * n_lanes] +
3869 this->dof_info
3870 ->component_dof_indices_offset[this->active_fe_index]
3871 [this->first_selected_component] *
3872 n_lanes;
3873 if (n_components == 1 || this->n_fe_components == 1)
3874 for (unsigned int comp = 0; comp < n_components; ++comp)
3875 operation.process_dofs_vectorized(dofs_per_component,
3876 dof_index,
3877 *src[comp],
3878 values_dofs[comp],
3879 vector_selector);
3880 else
3881 operation.process_dofs_vectorized(dofs_per_component * n_components,
3882 dof_index,
3883 *src[0],
3884 values_dofs[0],
3885 vector_selector);
3886 return;
3887 }
3888
3889 const std::array<unsigned int, n_lanes> &cells = this->get_cell_or_face_ids();
3890
3891 // More general case: Must go through the components one by one and apply
3892 // some transformations
3893 const unsigned int n_filled_lanes =
3894 dof_info.n_vectorization_lanes_filled[ind][this->cell];
3895
3896 const bool use_vectorized_path = n_filled_lanes == n_lanes &&
3897 n_active_lanes == n_lanes &&
3898 !accesses_exterior_dofs;
3899
3900 if (vectors_sm[0] != nullptr)
3901 {
3902 const auto compute_vector_ptrs = [&](const unsigned int comp) {
3903 std::array<typename VectorType::value_type *, n_lanes> vector_ptrs = {};
3904
3905 const auto upper_bound =
3906 std::min<unsigned int>(n_filled_lanes, n_lanes);
3907 for (unsigned int v = 0; v < upper_bound; ++v)
3908 {
3909 if (mask[v] == false)
3910 {
3911 vector_ptrs[v] = nullptr;
3912 continue;
3913 }
3914
3917 Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3918 ExcIndexRange(ind,
3919 0,
3920 dof_info.dof_indices_contiguous_sm.size()));
3921 Assert(
3922 cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3923 ExcIndexRange(cells[v],
3924 0,
3925 dof_info.dof_indices_contiguous_sm[ind].size()));
3926
3927 const auto &temp =
3928 dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3929
3930 if (temp.first != numbers::invalid_unsigned_int)
3931 vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3932 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3934 [this->active_fe_index][this->first_selected_component]);
3935 else
3936 vector_ptrs[v] = nullptr;
3937 }
3938 for (unsigned int v = n_filled_lanes; v < n_lanes; ++v)
3939 vector_ptrs[v] = nullptr;
3940
3941 return vector_ptrs;
3942 };
3943
3944 if (use_vectorized_path)
3945 {
3946 if (n_components == 1 || this->n_fe_components == 1)
3947 {
3948 for (unsigned int comp = 0; comp < n_components; ++comp)
3949 {
3950 auto vector_ptrs = compute_vector_ptrs(comp);
3951 operation.process_dofs_vectorized_transpose(
3952 dofs_per_component,
3953 vector_ptrs,
3954 values_dofs[comp],
3955 vector_selector);
3956 }
3957 }
3958 else
3959 {
3960 auto vector_ptrs = compute_vector_ptrs(0);
3961 operation.process_dofs_vectorized_transpose(dofs_per_component *
3962 n_components,
3963 vector_ptrs,
3964 &values_dofs[0][0],
3965 vector_selector);
3966 }
3967 }
3968 else
3969 for (unsigned int comp = 0; comp < n_components; ++comp)
3970 {
3971 auto vector_ptrs = compute_vector_ptrs(
3972 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3973
3974 for (unsigned int i = 0; i < dofs_per_component; ++i)
3975 operation.process_empty(values_dofs[comp][i]);
3976
3977 if (n_components == 1 || this->n_fe_components == 1)
3978 {
3979 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3980 if (mask[v] == true)
3981 for (unsigned int i = 0; i < dofs_per_component; ++i)
3982 operation.process_dof(vector_ptrs[v][i],
3983 values_dofs[comp][i][v]);
3984 }
3985 else
3986 {
3987 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3988 if (mask[v] == true)
3989 for (unsigned int i = 0; i < dofs_per_component; ++i)
3990 operation.process_dof(
3991 vector_ptrs[v][i + comp * dofs_per_component],
3992 values_dofs[comp][i][v]);
3993 }
3994 }
3995 return;
3996 }
3997
3998 std::array<unsigned int, n_lanes> dof_indices;
3999 std::fill(dof_indices.begin(),
4000 dof_indices.end(),
4002
4003 Assert(n_filled_lanes <= n_lanes, ExcInternalError());
4004 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4005 {
4006 Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
4008 if (mask[v] == true)
4009 dof_indices[v] =
4010 dof_indices_cont[cells[v]] +
4011 this->dof_info
4012 ->component_dof_indices_offset[this->active_fe_index]
4013 [this->first_selected_component] *
4014 dof_info.dof_indices_interleave_strides[ind][cells[v]];
4015 }
4016
4017 // In the case with contiguous cell indices, we know that there are no
4018 // constraints and that the indices within each element are contiguous
4019 if (use_vectorized_path)
4020 {
4021 if (dof_info.index_storage_variants[ind][this->cell] ==
4023 contiguous)
4024 {
4025 if (n_components == 1 || this->n_fe_components == 1)
4026 for (unsigned int comp = 0; comp < n_components; ++comp)
4027 operation.process_dofs_vectorized_transpose(dofs_per_component,
4028 dof_indices.data(),
4029 *src[comp],
4030 values_dofs[comp],
4031 vector_selector);
4032 else
4033 operation.process_dofs_vectorized_transpose(dofs_per_component *
4034 n_components,
4035 dof_indices.data(),
4036 *src[0],
4037 &values_dofs[0][0],
4038 vector_selector);
4039 }
4040 else if (dof_info.index_storage_variants[ind][this->cell] ==
4042 interleaved_contiguous_strided)
4043 {
4044 std::array<typename VectorType::value_type *, n_components> src_ptrs;
4045 if (n_components == 1 || this->n_fe_components == 1)
4046 for (unsigned int comp = 0; comp < n_components; ++comp)
4047 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
4048 src[comp]->begin());
4049 else
4050 src_ptrs[0] =
4051 const_cast<typename VectorType::value_type *>(src[0]->begin());
4052
4053 if (n_components == 1 || this->n_fe_components == 1)
4054 for (unsigned int i = 0; i < dofs_per_component; ++i)
4055 {
4056 for (unsigned int comp = 0; comp < n_components; ++comp)
4057 operation.process_dof_gather(dof_indices.data(),
4058 *src[comp],
4059 i * n_lanes,
4060 src_ptrs[comp] + i * n_lanes,
4061 values_dofs[comp][i],
4062 vector_selector);
4063 }
4064 else
4065 for (unsigned int comp = 0; comp < n_components; ++comp)
4066 for (unsigned int i = 0; i < dofs_per_component; ++i)
4067 {
4068 operation.process_dof_gather(
4069 dof_indices.data(),
4070 *src[0],
4071 (comp * dofs_per_component + i) * n_lanes,
4072 src_ptrs[0] + (comp * dofs_per_component + i) * n_lanes,
4073 values_dofs[comp][i],
4074 vector_selector);
4075 }
4076 }
4077 else
4078 {
4079 Assert(dof_info.index_storage_variants[ind][this->cell] ==
4081 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4083 std::array<typename VectorType::value_type *, n_components> src_ptrs;
4084 if (n_components == 1 || this->n_fe_components == 1)
4085 for (unsigned int comp = 0; comp < n_components; ++comp)
4086 src_ptrs[comp] = const_cast<typename VectorType::value_type *>(
4087 src[comp]->begin());
4088 else
4089 src_ptrs[0] =
4090 const_cast<typename VectorType::value_type *>(src[0]->begin());
4091
4092 const unsigned int *offsets =
4093 &dof_info.dof_indices_interleave_strides[ind][n_lanes * this->cell];
4094 if (n_components == 1 || this->n_fe_components == 1)
4095 for (unsigned int i = 0; i < dofs_per_component; ++i)
4096 {
4097 for (unsigned int comp = 0; comp < n_components; ++comp)
4098 operation.process_dof_gather(dof_indices.data(),
4099 *src[comp],
4100 0,
4101 src_ptrs[comp],
4102 values_dofs[comp][i],
4103 vector_selector);
4105 for (unsigned int v = 0; v < n_lanes; ++v)
4106 dof_indices[v] += offsets[v];
4107 }
4108 else
4109 for (unsigned int comp = 0; comp < n_components; ++comp)
4110 for (unsigned int i = 0; i < dofs_per_component; ++i)
4111 {
4112 operation.process_dof_gather(dof_indices.data(),
4113 *src[0],
4114 0,
4115 src_ptrs[0],
4116 values_dofs[comp][i],
4117 vector_selector);
4119 for (unsigned int v = 0; v < n_lanes; ++v)
4120 dof_indices[v] += offsets[v];
4121 }
4122 }
4123 }
4124 else
4125 for (unsigned int comp = 0; comp < n_components; ++comp)
4126 {
4127 for (unsigned int i = 0; i < dofs_per_component; ++i)
4128 operation.process_empty(values_dofs[comp][i]);
4129 if (accesses_exterior_dofs)
4130 {
4131 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4132 if (mask[v] == true)
4133 {
4134 if (dof_info.index_storage_variants
4135 [ind][cells[v] / VectorizedArrayType::size()] ==
4138 {
4139 if (n_components == 1 || this->n_fe_components == 1)
4140 {
4141 for (unsigned int i = 0; i < dofs_per_component; ++i)
4142 operation.process_dof(dof_indices[v] + i,
4143 *src[comp],
4144 values_dofs[comp][i][v]);
4145 }
4146 else
4147 {
4148 for (unsigned int i = 0; i < dofs_per_component; ++i)
4149 operation.process_dof(dof_indices[v] + i +
4150 comp * dofs_per_component,
4151 *src[0],
4152 values_dofs[comp][i][v]);
4153 }
4154 }
4155 else
4156 {
4157 const unsigned int offset =
4158 dof_info.dof_indices_interleave_strides[ind][cells[v]];
4159 AssertIndexRange(offset, VectorizedArrayType::size() + 1);
4160 if (n_components == 1 || this->n_fe_components == 1)
4161 {
4162 for (unsigned int i = 0; i < dofs_per_component; ++i)
4163 operation.process_dof(dof_indices[v] + i * offset,
4164 *src[comp],
4165 values_dofs[comp][i][v]);
4166 }
4167 else
4168 {
4169 for (unsigned int i = 0; i < dofs_per_component; ++i)
4170 operation.process_dof(
4171 dof_indices[v] +
4172 (i + comp * dofs_per_component) * offset,
4173 *src[0],
4174 values_dofs[comp][i][v]);
4175 }
4176 }
4177 }
4178 }
4179 else
4180 {
4181 if (dof_info.index_storage_variants[ind][this->cell] ==
4183 contiguous)
4184 {
4185 if (n_components == 1 || this->n_fe_components == 1)
4186 {
4187 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4188 if (mask[v] == true)
4189 for (unsigned int i = 0; i < dofs_per_component; ++i)
4190 operation.process_dof(dof_indices[v] + i,
4191 *src[comp],
4192 values_dofs[comp][i][v]);
4193 }
4194 else
4195 {
4196 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4197 if (mask[v] == true)
4198 for (unsigned int i = 0; i < dofs_per_component; ++i)
4199 operation.process_dof(dof_indices[v] + i +
4200 comp * dofs_per_component,
4201 *src[0],
4202 values_dofs[comp][i][v]);
4203 }
4204 }
4205 else
4206 {
4207 const unsigned int *offsets =
4209 [ind][VectorizedArrayType::size() * this->cell];
4210 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4211 AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4212 if (n_components == 1 || this->n_fe_components == 1)
4213 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4214 {
4215 if (mask[v] == true)
4216 for (unsigned int i = 0; i < dofs_per_component; ++i)
4217 operation.process_dof(dof_indices[v] + i * offsets[v],
4218 *src[comp],
4219 values_dofs[comp][i][v]);
4220 }
4221 else
4222 {
4223 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4224 if (mask[v] == true)
4225 for (unsigned int i = 0; i < dofs_per_component; ++i)
4226 operation.process_dof(
4227 dof_indices[v] +
4228 (i + comp * dofs_per_component) * offsets[v],
4229 *src[0],
4230 values_dofs[comp][i][v]);
4231 }
4232 }
4233 }
4234 }
4235}
4236
4237namespace internal
4238{
4239 template <
4240 typename Number,
4241 typename VectorType,
4242 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4243 decltype(std::declval<VectorType>().begin())
4244 get_beginning(VectorType &vec)
4245 {
4246 return vec.begin();
4247 }
4248
4249 template <
4250 typename Number,
4251 typename VectorType,
4252 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4253 typename VectorType::value_type *
4254 get_beginning(VectorType &)
4255 {
4256 return nullptr;
4257 }
4258
4259 template <typename VectorType,
4260 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4261 nullptr>
4262 const std::vector<ArrayView<const typename VectorType::value_type>> *
4263 get_shared_vector_data(VectorType *vec,
4264 const bool is_valid_mode_for_sm,
4265 const unsigned int active_fe_index,
4267 {
4268 // note: no hp is supported
4269 if (is_valid_mode_for_sm &&
4270 dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4271 .size() > 0 &&
4272 active_fe_index == 0)
4273 return &vec->shared_vector_data();
4274 else
4275 return nullptr;
4276 }
4277
4278 template <typename VectorType,
4279 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4280 * = nullptr>
4281 const std::vector<ArrayView<const typename VectorType::value_type>> *
4282 get_shared_vector_data(VectorType *,
4283 const bool,
4284 const unsigned int,
4286 {
4287 return nullptr;
4288 }
4289
4290 template <int n_components, typename VectorType>
4291 std::pair<
4292 std::array<typename internal::BlockVectorSelector<
4293 VectorType,
4294 IsBlockVector<VectorType>::value>::BaseVectorType *,
4295 n_components>,
4296 std::array<
4297 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4298 VectorType,
4299 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4300 n_components>>
4301 get_vector_data(VectorType &src,
4302 const unsigned int first_index,
4303 const bool is_valid_mode_for_sm,
4304 const unsigned int active_fe_index,
4306 {
4307 // select between block vectors and non-block vectors. Note that the number
4308 // of components is checked in the internal data
4309 std::pair<
4310 std::array<typename internal::BlockVectorSelector<
4311 VectorType,
4312 IsBlockVector<VectorType>::value>::BaseVectorType *,
4313 n_components>,
4314 std::array<
4315 const std::vector<
4316 ArrayView<const typename internal::BlockVectorSelector<
4317 VectorType,
4318 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4319 n_components>>
4320 src_data;
4321
4322 for (unsigned int d = 0; d < n_components; ++d)
4323 src_data.first[d] = internal::BlockVectorSelector<
4324 VectorType,
4325 IsBlockVector<VectorType>::value>::get_vector_component(src,
4326 d +
4327 first_index);
4328
4329 for (unsigned int d = 0; d < n_components; ++d)
4330 src_data.second[d] = get_shared_vector_data(
4331 const_cast<typename internal::BlockVectorSelector<
4332 std::remove_const_t<VectorType>,
4334 *>(src_data.first[d]),
4335 is_valid_mode_for_sm,
4336 active_fe_index,
4337 dof_info);
4338
4339 return src_data;
4340 }
4341} // namespace internal
4342
4343
4344
4345template <int dim,
4346 int n_components_,
4347 typename Number,
4348 bool is_face,
4349 typename VectorizedArrayType>
4350inline void
4353{
4354 if (this->dof_info == nullptr ||
4355 this->dof_info->hanging_node_constraint_masks.empty() ||
4356 this->dof_info->hanging_node_constraint_masks_comp.empty() ||
4357 this->dof_info->hanging_node_constraint_masks_comp
4358 [this->active_fe_index][this->first_selected_component] == false)
4359 return; // nothing to do with faces
4360
4361 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4362 constraint_mask;
4363
4364 bool hn_available = false;
4365
4366 const std::array<unsigned int, n_lanes> &cells = this->get_cell_ids();
4367
4368 for (unsigned int v = 0; v < n_lanes; ++v)
4369 {
4370 if (cells[v] == numbers::invalid_unsigned_int)
4371 {
4372 constraint_mask[v] = internal::MatrixFreeFunctions::
4374 continue;
4375 }
4376
4377 const unsigned int cell_index = cells[v];
4378 const auto mask =
4380 constraint_mask[v] = mask;
4381
4382 hn_available |= (mask != internal::MatrixFreeFunctions::
4384 }
4385
4386 if (hn_available == false)
4387 return; // no hanging node on cell batch -> nothing to do
4388
4390 apply(n_components,
4391 this->data->data.front().fe_degree,
4392 this->get_shape_info(),
4393 transpose,
4394 constraint_mask,
4395 this->values_dofs);
4396}
4397
4398
4399
4400template <int dim,
4401 int n_components_,
4402 typename Number,
4403 bool is_face,
4404 typename VectorizedArrayType>
4405template <typename VectorType>
4406inline void
4408 read_dof_values(const VectorType &src,
4409 const unsigned int first_index,
4410 const std::bitset<n_lanes> &mask)
4411{
4412 const auto src_data = internal::get_vector_data<n_components_>(
4413 src,
4414 first_index,
4415 this->dof_access_index ==
4417 this->active_fe_index,
4418 this->dof_info);
4419
4421 read_write_operation(reader, src_data.first, src_data.second, mask, true);
4422
4423 apply_hanging_node_constraints(false);
4424
4425# ifdef DEBUG
4426 this->dof_values_initialized = true;
4427# endif
4428}
4429
4430
4431
4432template <int dim,
4433 int n_components_,
4434 typename Number,
4435 bool is_face,
4436 typename VectorizedArrayType>
4437template <typename VectorType>
4438inline void
4440 read_dof_values_plain(const VectorType &src,
4441 const unsigned int first_index,
4442 const std::bitset<n_lanes> &mask)
4443{
4444 const auto src_data = internal::get_vector_data<n_components_>(
4445 src,
4446 first_index,
4447 this->dof_access_index ==
4449 this->active_fe_index,
4450 this->dof_info);
4451
4453 read_write_operation(reader, src_data.first, src_data.second, mask, false);
4454
4455# ifdef DEBUG
4456 this->dof_values_initialized = true;
4457# endif
4458}
4459
4460
4461
4462template <int dim,
4463 int n_components_,
4464 typename Number,
4465 bool is_face,
4466 typename VectorizedArrayType>
4467template <typename VectorType>
4468inline void
4470 distribute_local_to_global(VectorType &dst,
4471 const unsigned int first_index,
4472 const std::bitset<n_lanes> &mask) const
4473{
4474# ifdef DEBUG
4475 Assert(this->dof_values_initialized == true,
4477# endif
4478
4479 apply_hanging_node_constraints(true);
4480
4481 const auto dst_data = internal::get_vector_data<n_components_>(
4482 dst,
4483 first_index,
4484 this->dof_access_index ==
4486 this->active_fe_index,
4487 this->dof_info);
4488
4490 distributor;
4491 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4492}
4493
4494
4495
4496template <int dim,
4497 int n_components_,
4498 typename Number,
4499 bool is_face,
4500 typename VectorizedArrayType>
4501template <typename VectorType>
4502inline void
4504 set_dof_values(VectorType &dst,
4505 const unsigned int first_index,
4506 const std::bitset<n_lanes> &mask) const
4507{
4508# ifdef DEBUG
4509 Assert(this->dof_values_initialized == true,
4511# endif
4512
4513 const auto dst_data = internal::get_vector_data<n_components_>(
4514 dst,
4515 first_index,
4516 this->dof_access_index ==
4518 this->active_fe_index,
4519 this->dof_info);
4520
4522 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4523}
4524
4525
4526
4527template <int dim,
4528 int n_components_,
4529 typename Number,
4530 bool is_face,
4531 typename VectorizedArrayType>
4532template <typename VectorType>
4533inline void
4535 set_dof_values_plain(VectorType &dst,
4536 const unsigned int first_index,
4537 const std::bitset<n_lanes> &mask) const
4538{
4539# ifdef DEBUG
4540 Assert(this->dof_values_initialized == true,
4542# endif
4543
4544 const auto dst_data = internal::get_vector_data<n_components_>(
4545 dst,
4546 first_index,
4547 this->dof_access_index ==
4549 this->active_fe_index,
4550 this->dof_info);
4551
4553 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4554}
4555
4556
4557
4558/*------------------------------ access to data fields ----------------------*/
4559
4560
4561
4562template <int dim,
4563 int n_components_,
4564 typename Number,
4565 bool is_face,
4566 typename VectorizedArrayType>
4569 get_dof_value(const unsigned int dof) const
4570{
4571 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4572 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4574 for (unsigned int comp = 0; comp < n_components; ++comp)
4575 return_value[comp] = this->values_dofs[comp * dofs + dof];
4576 return return_value;
4577}
4578
4579
4580
4581template <int dim,
4582 int n_components_,
4583 typename Number,
4584 bool is_face,
4585 typename VectorizedArrayType>
4588 get_value(const unsigned int q_point) const
4589{
4590# ifdef DEBUG
4591 Assert(this->values_quad_initialized == true,
4593# endif
4594
4595 AssertIndexRange(q_point, this->n_quadrature_points);
4596 const std::size_t nqp = this->n_quadrature_points;
4598 for (unsigned int comp = 0; comp < n_components; ++comp)
4599 return_value[comp] = this->values_quad[comp * nqp + q_point];
4600 return return_value;
4601}
4602
4603
4604
4605template <int dim,
4606 int n_components_,
4607 typename Number,
4608 bool is_face,
4609 typename VectorizedArrayType>
4613 get_gradient(const unsigned int q_point) const
4614{
4615# ifdef DEBUG
4616 Assert(this->gradients_quad_initialized == true,
4618# endif
4619
4620 AssertIndexRange(q_point, this->n_quadrature_points);
4621 Assert(this->jacobian != nullptr,
4623 "update_gradients"));
4624 const std::size_t nqp = this->n_quadrature_points;
4626
4627 // Cartesian cell
4628 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4629 {
4630 for (unsigned int comp = 0; comp < n_components; ++comp)
4631 for (unsigned int d = 0; d < dim; ++d)
4632 grad_out[comp][d] =
4633 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4634 this->jacobian[0][d][d];
4635 }
4636 // cell with general/affine Jacobian
4637 else
4638 {
4640 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4641 q_point :
4642 0];
4643 for (unsigned int comp = 0; comp < n_components; ++comp)
4644 for (unsigned int d = 0; d < dim; ++d)
4645 {
4646 grad_out[comp][d] =
4647 jac[d][0] * this->gradients_quad[(comp * nqp + q_point) * dim];
4648 for (unsigned int e = 1; e < dim; ++e)
4649 grad_out[comp][d] +=
4650 jac[d][e] *
4651 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4652 }
4653 }
4654 return grad_out;
4655}
4656
4657
4658
4659template <int dim,
4660 int n_components_,
4661 typename Number,
4662 bool is_face,
4663 typename VectorizedArrayType>
4666 get_normal_derivative(const unsigned int q_point) const
4667{
4668 AssertIndexRange(q_point, this->n_quadrature_points);
4669# ifdef DEBUG
4670 Assert(this->gradients_quad_initialized == true,
4672# endif
4673
4674 Assert(this->normal_x_jacobian != nullptr,
4676 "update_gradients"));
4677
4678 const std::size_t nqp = this->n_quadrature_points;
4680
4681 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4682 for (unsigned int comp = 0; comp < n_components; ++comp)
4683 grad_out[comp] =
4684 this->gradients_quad[(comp * nqp + q_point) * dim + dim - 1] *
4685 (this->normal_x_jacobian[0][dim - 1]);
4686 else
4687 {
4688 const std::size_t index =
4689 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4690 for (unsigned int comp = 0; comp < n_components; ++comp)
4691 {
4692 grad_out[comp] = this->gradients_quad[(comp * nqp + q_point) * dim] *
4693 this->normal_x_jacobian[index][0];
4694 for (unsigned int d = 1; d < dim; ++d)
4695 grad_out[comp] +=
4696 this->gradients_quad[(comp * nqp + q_point) * dim + d] *
4697 this->normal_x_jacobian[index][d];
4698 }
4699 }
4700 return grad_out;
4701}
4702
4703
4704
4705namespace internal
4706{
4707 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4708 // store the lower diagonal because of symmetry
4709 template <typename VectorizedArrayType>
4710 inline void
4711 hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4712 const VectorizedArrayType *const hessians,
4713 const unsigned int,
4714 VectorizedArrayType (&tmp)[1][1])
4715 {
4716 tmp[0][0] = jac[0][0] * hessians[0];
4717 }
4718
4719 template <typename VectorizedArrayType>
4720 inline void
4721 hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4722 const VectorizedArrayType *const hessians,
4723 const unsigned int nqp,
4724 VectorizedArrayType (&tmp)[2][2])
4725 {
4726 for (unsigned int d = 0; d < 2; ++d)
4727 {
4728 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4729 tmp[1][d] =
4730 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4731 }
4732 }
4733
4734 template <typename VectorizedArrayType>
4735 inline void
4736 hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4737 const VectorizedArrayType *const hessians,
4738 const unsigned int nqp,
4739 VectorizedArrayType (&tmp)[3][3])
4740 {
4741 for (unsigned int d = 0; d < 3; ++d)
4742 {
4743 tmp[0][d] =
4744 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4745 jac[d][2] * hessians[4 * nqp]);
4746 tmp[1][d] =
4747 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4748 jac[d][2] * hessians[5 * nqp]);
4749 tmp[2][d] =
4750 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4751 jac[d][2] * hessians[2 * nqp]);
4752 }
4753 }
4754} // namespace internal
4755
4756
4757
4758template <int dim,
4759 int n_components_,
4760 typename Number,
4761 bool is_face,
4762 typename VectorizedArrayType>
4765 get_hessian(const unsigned int q_point) const
4766{
4767# ifdef DEBUG
4768 Assert(this->hessians_quad_initialized == true,
4770# endif
4771 AssertIndexRange(q_point, this->n_quadrature_points);
4772
4773 Assert(this->jacobian != nullptr,
4775 "update_hessian"));
4777 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4778 0 :
4779 q_point];
4780
4782
4783 const std::size_t nqp = this->n_quadrature_points;
4784 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4785
4786 // Cartesian cell
4787 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4788 {
4789 for (unsigned int comp = 0; comp < n_components; ++comp)
4790 {
4791 for (unsigned int d = 0; d < dim; ++d)
4792 hessian_out[comp][d][d] =
4793 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4794 (jac[d][d] * jac[d][d]);
4795 switch (dim)
4796 {
4797 case 1:
4798 break;
4799 case 2:
4800 hessian_out[comp][0][1] =
4801 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4802 (jac[0][0] * jac[1][1]);
4803 break;
4804 case 3:
4805 hessian_out[comp][0][1] =
4806 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4807 (jac[0][0] * jac[1][1]);
4808 hessian_out[comp][0][2] =
4809 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4810 (jac[0][0] * jac[2][2]);
4811 hessian_out[comp][1][2] =
4812 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4813 (jac[1][1] * jac[2][2]);
4814 break;
4815 default:
4817 }
4818 for (unsigned int d = 0; d < dim; ++d)
4819 for (unsigned int e = d + 1; e < dim; ++e)
4820 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4821 }
4822 }
4823 // cell with general Jacobian, but constant within the cell
4824 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4825 {
4826 for (unsigned int comp = 0; comp < n_components; ++comp)
4827 {
4828 VectorizedArrayType tmp[dim][dim];
4829 internal::hessian_unit_times_jac(
4830 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4831
4832 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4833 for (unsigned int d = 0; d < dim; ++d)
4834 for (unsigned int e = d; e < dim; ++e)
4835 {
4836 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4837 for (unsigned int f = 1; f < dim; ++f)
4838 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4839 }
4840
4841 // no J' * grad(u) part here because the Jacobian is constant
4842 // throughout the cell and hence, its derivative is zero
4843
4844 // take symmetric part
4845 for (unsigned int d = 0; d < dim; ++d)
4846 for (unsigned int e = d + 1; e < dim; ++e)
4847 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4848 }
4849 }
4850 // cell with general Jacobian
4851 else
4852 {
4853 const auto &jac_grad = this->jacobian_gradients[q_point];
4854 for (unsigned int comp = 0; comp < n_components; ++comp)
4855 {
4856 VectorizedArrayType tmp[dim][dim];
4857 internal::hessian_unit_times_jac(
4858 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4859
4860 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4861 for (unsigned int d = 0; d < dim; ++d)
4862 for (unsigned int e = d; e < dim; ++e)
4863 {
4864 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4865 for (unsigned int f = 1; f < dim; ++f)
4866 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4867 }
4868
4869 // add diagonal part of J' * grad(u)
4870 for (unsigned int d = 0; d < dim; ++d)
4871 for (unsigned int e = 0; e < dim; ++e)
4872 hessian_out[comp][d][d] +=
4873 jac_grad[d][e] *
4874 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4875
4876 // add off-diagonal part of J' * grad(u)
4877 for (unsigned int d = 0, count = dim; d < dim; ++d)
4878 for (unsigned int e = d + 1; e < dim; ++e, ++count)
4879 for (unsigned int f = 0; f < dim; ++f)
4880 hessian_out[comp][d][e] +=
4881 jac_grad[count][f] *
4882 this->gradients_quad[(comp * nqp + q_point) * dim + f];
4883
4884 // take symmetric part
4885 for (unsigned int d = 0; d < dim; ++d)
4886 for (unsigned int e = d + 1; e < dim; ++e)
4887 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4888 }
4889 }
4890 return hessian_out;
4891}
4892
4893
4894
4895template <int dim,
4896 int n_components_,
4897 typename Number,
4898 bool is_face,
4899 typename VectorizedArrayType>
4902 get_hessian_diagonal(const unsigned int q_point) const
4903{
4904 Assert(!is_face, ExcNotImplemented());
4905# ifdef DEBUG
4906 Assert(this->hessians_quad_initialized == true,
4908# endif
4909 AssertIndexRange(q_point, this->n_quadrature_points);
4910
4911 Assert(this->jacobian != nullptr, ExcNotImplemented());
4913 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4914 0 :
4915 q_point];
4916
4917 const std::size_t nqp = this->n_quadrature_points;
4918 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4920
4921 // Cartesian cell
4922 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4923 {
4924 for (unsigned int comp = 0; comp < n_components; ++comp)
4925 for (unsigned int d = 0; d < dim; ++d)
4926 hessian_out[comp][d] =
4927 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4928 (jac[d][d] * jac[d][d]);
4929 }
4930 // cell with general Jacobian, but constant within the cell
4931 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4932 {
4933 for (unsigned int comp = 0; comp < n_components; ++comp)
4934 {
4935 // compute laplacian before the gradient because it needs to access
4936 // unscaled gradient data
4937 VectorizedArrayType tmp[dim][dim];
4938 internal::hessian_unit_times_jac(
4939 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4940
4941 // compute only the trace part of hessian, J * tmp = J *
4942 // hess_unit(u) * J^T
4943 for (unsigned int d = 0; d < dim; ++d)
4944 {
4945 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4946 for (unsigned int f = 1; f < dim; ++f)
4947 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4948 }
4949 }
4950 }
4951 // cell with general Jacobian
4952 else
4953 {
4954 const auto &jac_grad = this->jacobian_gradients[q_point];
4955 for (unsigned int comp = 0; comp < n_components; ++comp)
4956 {
4957 // compute laplacian before the gradient because it needs to access
4958 // unscaled gradient data
4959 VectorizedArrayType tmp[dim][dim];
4960 internal::hessian_unit_times_jac(
4961 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4962
4963 // compute only the trace part of hessian, J * tmp = J *
4964 // hess_unit(u) * J^T
4965 for (unsigned int d = 0; d < dim; ++d)
4966 {
4967 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4968 for (unsigned int f = 1; f < dim; ++f)
4969 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4970 }
4971
4972 for (unsigned int d = 0; d < dim; ++d)
4973 for (unsigned int e = 0; e < dim; ++e)
4974 hessian_out[comp][d] +=
4975 jac_grad[d][e] *
4976 this->gradients_quad[(comp * nqp + q_point) * dim + e];
4977 }
4978 }
4979 return hessian_out;
4980}
4981
4982
4983
4984template <int dim,
4985 int n_components_,
4986 typename Number,
4987 bool is_face,
4988 typename VectorizedArrayType>
4991 get_laplacian(const unsigned int q_point) const
4992{
4993 Assert(is_face == false, ExcNotImplemented());
4994# ifdef DEBUG
4995 Assert(this->hessians_quad_initialized == true,
4997# endif
4998 AssertIndexRange(q_point, this->n_quadrature_points);
4999
5001 const auto hess_diag = get_hessian_diagonal(q_point);
5002 for (unsigned int comp = 0; comp < n_components; ++comp)
5003 {
5004 laplacian_out[comp] = hess_diag[comp][0];
5005 for (unsigned int d = 1; d < dim; ++d)
5006 laplacian_out[comp] += hess_diag[comp][d];
5007 }
5008 return laplacian_out;
5009}
5010
5011
5012
5013template <int dim,
5014 int n_components_,
5015 typename Number,
5016 bool is_face,
5017 typename VectorizedArrayType>
5018inline DEAL_II_ALWAYS_INLINE void
5021 const unsigned int dof)
5022{
5023# ifdef DEBUG
5024 this->dof_values_initialized = true;
5025# endif
5026 const std::size_t dofs = this->data->dofs_per_component_on_cell;
5027 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5028 for (unsigned int comp = 0; comp < n_components; ++comp)
5029 this->values_dofs[comp * dofs + dof] = val_in[comp];
5030}
5031
5032
5033
5034template <int dim,
5035 int n_components_,
5036 typename Number,
5037 bool is_face,
5038 typename VectorizedArrayType>
5039inline DEAL_II_ALWAYS_INLINE void
5042 const unsigned int q_point)
5043{
5044# ifdef DEBUG
5045 Assert(this->is_reinitialized, ExcNotInitialized());
5046# endif
5047 AssertIndexRange(q_point, this->n_quadrature_points);
5048 Assert(this->J_value != nullptr,
5050 "update_values"));
5051# ifdef DEBUG
5052 this->values_quad_submitted = true;
5053# endif
5054
5055 const std::size_t nqp = this->n_quadrature_points;
5056 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5057 {
5058 const VectorizedArrayType JxW =
5059 this->J_value[0] * this->quadrature_weights[q_point];
5060 for (unsigned int comp = 0; comp < n_components; ++comp)
5061 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5062 }
5063 else
5064 {
5065 const VectorizedArrayType JxW = this->J_value[q_point];
5066 for (unsigned int comp = 0; comp < n_components; ++comp)
5067 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5068 }
5069}
5070
5071
5072
5073template <int dim,
5074 int n_components_,
5075 typename Number,
5076 bool is_face,
5077 typename VectorizedArrayType>
5078inline DEAL_II_ALWAYS_INLINE void
5081 const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5082 const unsigned int q_point)
5083{
5084# ifdef DEBUG
5085 Assert(this->is_reinitialized, ExcNotInitialized());
5086# endif
5087 AssertIndexRange(q_point, this->n_quadrature_points);
5088 Assert(this->J_value != nullptr,
5090 "update_gradients"));
5091 Assert(this->jacobian != nullptr,
5093 "update_gradients"));
5094# ifdef DEBUG
5095 this->gradients_quad_submitted = true;
5096# endif
5097
5098 const std::size_t nqp_d = this->n_quadrature_points * dim;
5099 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5100
5101 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5102 {
5103 const VectorizedArrayType JxW =
5104 this->J_value[0] * this->quadrature_weights[q_point];
5105 std::array<VectorizedArrayType, dim> jac;
5106 for (unsigned int d = 0; d < dim; ++d)
5107 jac[d] = this->jacobian[0][d][d];
5108 for (unsigned int d = 0; d < dim; ++d)
5109 {
5110 const VectorizedArrayType factor = jac[d] * JxW;
5111 for (unsigned int comp = 0; comp < n_components; ++comp)
5112 gradients[comp * nqp_d + d] = grad_in[comp][d] * factor;
5113 }
5114 }
5115 else
5116 {
5118 this->cell_type > internal::MatrixFreeFunctions::affine ?
5119 this->jacobian[q_point] :
5120 this->jacobian[0];
5121 const VectorizedArrayType JxW =
5122 this->cell_type > internal::MatrixFreeFunctions::affine ?
5123 this->J_value[q_point] :
5124 this->J_value[0] * this->quadrature_weights[q_point];
5125 for (unsigned int comp = 0; comp < n_components; ++comp)
5126 for (unsigned int d = 0; d < dim; ++d)
5127 {
5128 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5129 for (unsigned int e = 1; e < dim; ++e)
5130 new_val += (jac[e][d] * grad_in[comp][e]);
5131 gradients[comp * nqp_d + d] = new_val * JxW;
5132 }
5133 }
5134}
5135
5136
5137
5138template <int dim,
5139 int n_components_,
5140 typename Number,
5141 bool is_face,
5142 typename VectorizedArrayType>
5143inline DEAL_II_ALWAYS_INLINE void
5147 const unsigned int q_point)
5148{
5149 AssertIndexRange(q_point, this->n_quadrature_points);
5150 Assert(this->normal_x_jacobian != nullptr,
5152 "update_gradients"));
5153# ifdef DEBUG
5154 this->gradients_quad_submitted = true;
5155# endif
5156
5157 const std::size_t nqp_d = this->n_quadrature_points * dim;
5158 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
5159
5160 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5161 {
5162 const VectorizedArrayType JxW_jac = this->J_value[0] *
5163 this->quadrature_weights[q_point] *
5164 this->normal_x_jacobian[0][dim - 1];
5165 for (unsigned int comp = 0; comp < n_components; ++comp)
5166 {
5167 for (unsigned int d = 0; d < dim - 1; ++d)
5168 gradients[comp * nqp_d + d] = VectorizedArrayType();
5169 gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac;
5170 }
5171 }
5172 else
5173 {
5174 const unsigned int index =
5175 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5177 this->normal_x_jacobian[index];
5178 const VectorizedArrayType JxW =
5179 (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5180 this->J_value[index] * this->quadrature_weights[q_point] :
5181 this->J_value[index];
5182 for (unsigned int comp = 0; comp < n_components; ++comp)
5183 {
5184 for (unsigned int d = 0; d < dim; ++d)
5185 gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d];
5186 }
5187 }
5188}
5189
5190
5191
5192template <int dim,
5193 int n_components_,
5194 typename Number,
5195 bool is_face,
5196 typename VectorizedArrayType>
5197inline DEAL_II_ALWAYS_INLINE void
5200 const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5201 hessian_in,
5202 const unsigned int q_point)
5203{
5204# ifdef DEBUG
5205 Assert(this->is_reinitialized, ExcNotInitialized());
5206# endif
5207 AssertIndexRange(q_point, this->n_quadrature_points);
5208 Assert(this->J_value != nullptr,
5210 "update_hessians"));
5211 Assert(this->jacobian != nullptr,
5213 "update_hessians"));
5214# ifdef DEBUG
5215 this->hessians_quad_submitted = true;
5216# endif
5217
5218 // compute hessian_unit = J^T * hessian_in(u) * J
5219 const std::size_t nqp = this->n_quadrature_points;
5220 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5221 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5222 {
5223 const VectorizedArrayType JxW =
5224 this->J_value[0] * this->quadrature_weights[q_point];
5225
5226 // diagonal part
5227 for (unsigned int d = 0; d < dim; ++d)
5228 {
5229 const auto jac_d = this->jacobian[0][d][d];
5230 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5231 for (unsigned int comp = 0; comp < n_components; ++comp)
5232 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5233 hessian_in[comp][d][d] * factor;
5234 }
5235
5236 // off diagonal part
5237 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5238 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5239 {
5240 const auto jac_d = this->jacobian[0][d][d];
5241 const auto jac_e = this->jacobian[0][e][e];
5242 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5243 for (unsigned int comp = 0; comp < n_components; ++comp)
5244 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5245 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5246 }
5247 }
5248 // cell with general Jacobian, but constant within the cell
5249 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5250 {
5251 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5252 const VectorizedArrayType JxW =
5253 this->J_value[0] * this->quadrature_weights[q_point];
5254 for (unsigned int comp = 0; comp < n_components; ++comp)
5255 {
5256 // 1. tmp = hessian_in(u) * J
5257 VectorizedArrayType tmp[dim][dim];
5258 for (unsigned int i = 0; i < dim; ++i)
5259 for (unsigned int j = 0; j < dim; ++j)
5260 {
5261 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5262 for (unsigned int k = 1; k < dim; ++k)
5263 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5264 }
5265
5266 // 2. hessian_unit = J^T * tmp
5267 VectorizedArrayType tmp2[dim][dim];
5268 for (unsigned int i = 0; i < dim; ++i)
5269 for (unsigned int j = 0; j < dim; ++j)
5270 {
5271 tmp2[i][j] = jac[0][i] * tmp[0][j];
5272 for (unsigned int k = 1; k < dim; ++k)
5273 tmp2[i][j] += jac[k][i] * tmp[k][j];
5274 }
5275
5276 // diagonal part
5277 for (unsigned int d = 0; d < dim; ++d)
5278 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5279 tmp2[d][d] * JxW;
5280
5281 // off diagonal part
5282 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5283 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5284 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5285 (tmp2[d][e] + tmp2[e][d]) * JxW;
5286 }
5287 }
5288 else
5289 {
5290 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5291 const VectorizedArrayType JxW = this->J_value[q_point];
5292 const auto &jac_grad = this->jacobian_gradients[q_point];
5293 for (unsigned int comp = 0; comp < n_components; ++comp)
5294 {
5295 // 1. tmp = hessian_in(u) * J
5296 VectorizedArrayType tmp[dim][dim];
5297 for (unsigned int i = 0; i < dim; ++i)
5298 for (unsigned int j = 0; j < dim; ++j)
5299 {
5300 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5301 for (unsigned int k = 1; k < dim; ++k)
5302 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5303 }
5304
5305 // 2. hessian_unit = J^T * tmp
5306 VectorizedArrayType tmp2[dim][dim];
5307 for (unsigned int i = 0; i < dim; ++i)
5308 for (unsigned int j = 0; j < dim; ++j)
5309 {
5310 tmp2[i][j] = jac[0][i] * tmp[0][j];
5311 for (unsigned int k = 1; k < dim; ++k)
5312 tmp2[i][j] += jac[k][i] * tmp[k][j];
5313 }
5314
5315 // diagonal part
5316 for (unsigned int d = 0; d < dim; ++d)
5317 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5318 tmp2[d][d] * JxW;
5319
5320 // off diagonal part
5321 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5322 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5323 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5324 (tmp2[d][e] + tmp2[e][d]) * JxW;
5325
5326 // 3. gradient_unit = J' ** hessian_in
5327 for (unsigned int d = 0; d < dim; ++d)
5328 {
5329 VectorizedArrayType sum = 0;
5330 for (unsigned int e = 0; e < dim; ++e)
5331 sum += hessian_in[comp][e][e] * jac_grad[e][d];
5332 for (unsigned int e = 0, count = dim; e < dim; ++e)
5333 for (unsigned int f = e + 1; f < dim; ++f, ++count)
5334 sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5335 jac_grad[count][d];
5336 this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim +
5337 d] = sum * JxW;
5338 }
5339 }
5340 }
5341}
5342
5343
5344
5345template <int dim,
5346 int n_components_,
5347 typename Number,
5348 bool is_face,
5349 typename VectorizedArrayType>
5352 integrate_value() const
5353{
5354# ifdef DEBUG
5355 Assert(this->is_reinitialized, ExcNotInitialized());
5356 Assert(this->values_quad_submitted == true,
5358# endif
5359
5361 const std::size_t nqp = this->n_quadrature_points;
5362 for (unsigned int q = 0; q < nqp; ++q)
5363 for (unsigned int comp = 0; comp < n_components; ++comp)
5364 return_value[comp] += this->values_quad[comp * nqp + q];
5365 return (return_value);
5366}
5367
5368
5369
5370/*----------------------- FEEvaluationAccess --------------------------------*/
5371
5372
5373template <int dim,
5374 int n_components_,
5375 typename Number,
5376 bool is_face,
5377 typename VectorizedArrayType>
5378inline FEEvaluationAccess<dim,
5379 n_components_,
5380 Number,
5381 is_face,
5382 VectorizedArrayType>::
5383 FEEvaluationAccess(
5385 const unsigned int dof_no,
5386 const unsigned int first_selected_component,
5387 const unsigned int quad_no,
5388 const unsigned int fe_degree,
5389 const unsigned int n_q_points,
5390 const bool is_interior_face,
5391 const unsigned int active_fe_index,
5392 const unsigned int active_quad_index,
5393 const unsigned int face_type)
5394 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5395 matrix_free,
5396 dof_no,
5397 first_selected_component,
5398 quad_no,
5399 fe_degree,
5400 n_q_points,
5401 is_interior_face,
5402 active_fe_index,
5403 active_quad_index,
5404 face_type)
5405{}
5406
5407
5408
5409template <int dim,
5410 int n_components_,
5411 typename Number,
5412 bool is_face,
5413 typename VectorizedArrayType>
5414inline FEEvaluationAccess<dim,
5415 n_components_,
5416 Number,
5417 is_face,
5418 VectorizedArrayType>::
5419 FEEvaluationAccess(
5420 const Mapping<dim> &mapping,
5421 const FiniteElement<dim> &fe,
5422 const Quadrature<1> &quadrature,
5423 const UpdateFlags update_flags,
5424 const unsigned int first_selected_component,
5426 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5427 mapping,
5428 fe,
5429 quadrature,
5430 update_flags,
5431 first_selected_component,
5432 other)
5433{}
5434
5435
5436
5437template <int dim,
5438 int n_components_,
5439 typename Number,
5440 bool is_face,
5441 typename VectorizedArrayType>
5442inline FEEvaluationAccess<dim,
5443 n_components_,
5444 Number,
5445 is_face,
5446 VectorizedArrayType>::
5447 FEEvaluationAccess(const FEEvaluationAccess<dim,
5448 n_components_,
5449 Number,
5450 is_face,
5451 VectorizedArrayType> &other)
5452 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5453 other)
5454{}
5455
5456
5457
5458template <int dim,
5459 int n_components_,
5460 typename Number,
5461 bool is_face,
5462 typename VectorizedArrayType>
5463inline FEEvaluationAccess<dim,
5464 n_components_,
5465 Number,
5466 is_face,
5467 VectorizedArrayType> &
5470 n_components_,
5471 Number,
5472 is_face,
5473 VectorizedArrayType> &other)
5474{
5475 this->FEEvaluationBase<dim,
5476 n_components_,
5477 Number,
5478 is_face,
5479 VectorizedArrayType>::operator=(other);
5480 return *this;
5481}
5482
5483
5484
5485/*-------------------- FEEvaluationAccess scalar ----------------------------*/
5486
5487
5488template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5492 const unsigned int dof_no,
5493 const unsigned int first_selected_component,
5494 const unsigned int quad_no,
5495 const unsigned int fe_degree,
5496 const unsigned int n_q_points,
5497 const bool is_interior_face,
5498 const unsigned int active_fe_index,
5499 const unsigned int active_quad_index,
5500 const unsigned int face_type)
5501 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5502 matrix_free,
5503 dof_no,
5504 first_selected_component,
5505 quad_no,
5506 fe_degree,
5507 n_q_points,
5508 is_interior_face,
5509 active_fe_index,
5510 active_quad_index,
5511 face_type)
5512{}
5513
5514
5515
5516template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5519 const Mapping<dim> &mapping,
5520 const FiniteElement<dim> &fe,
5521 const Quadrature<1> &quadrature,
5522 const UpdateFlags update_flags,
5523 const unsigned int first_selected_component,
5525 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5526 mapping,
5527 fe,
5528 quadrature,
5529 update_flags,
5530 first_selected_component,
5531 other)
5532{}
5533
5534
5535
5536template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5540 &other)
5541 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5542{}
5543
5544
5545
5546template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5550{
5551 this
5552 ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5553 other);
5554 return *this;
5555}
5556
5557
5558
5559template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5560inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5562 const unsigned int dof) const
5563{
5564 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5565 return this->values_dofs[dof];
5566}
5567
5568
5569
5570template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5571inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5573 const unsigned int q_point) const
5574{
5575# ifdef DEBUG
5576 Assert(this->values_quad_initialized == true,
5578# endif
5579 AssertIndexRange(q_point, this->n_quadrature_points);
5580 return this->values_quad[q_point];
5581}
5582
5583
5584
5585template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5586inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5588 get_normal_derivative(const unsigned int q_point) const
5589{
5590 return BaseClass::get_normal_derivative(q_point)[0];
5591}
5592
5593
5594
5595template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5598 const unsigned int q_point_in) const
5599{
5600 // could use the base class gradient, but that involves too many expensive
5601 // initialization operations on tensors
5602
5603# ifdef DEBUG
5604 Assert(this->gradients_quad_initialized == true,
5606# endif
5607 AssertIndexRange(q_point_in, this->n_quadrature_points);
5608
5609 Assert(this->jacobian != nullptr,
5611 "update_gradients"));
5612
5614
5615 const std::size_t q_point = q_point_in;
5616 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5617 {
5618 for (unsigned int d = 0; d < dim; ++d)
5619 grad_out[d] =
5620 this->gradients_quad[dim * q_point + d] * this->jacobian[0][d][d];
5621 }
5622 // cell with general/affine Jacobian
5623 else
5624 {
5626 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5627 q_point :
5628 0];
5629 for (unsigned int d = 0; d < dim; ++d)
5630 {
5631 grad_out[d] = jac[d][0] * this->gradients_quad[dim * q_point];
5632 for (unsigned int e = 1; e < dim; ++e)
5633 grad_out[d] += jac[d][e] * this->gradients_quad[dim * q_point + e];
5634 }
5635 }
5636 return grad_out;
5637}
5638
5639
5640
5641template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5644 const unsigned int q_point) const
5645{
5646 return BaseClass::get_hessian(q_point)[0];
5647}
5648
5649
5650
5651template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5654 get_hessian_diagonal(const unsigned int q_point) const
5655{
5656 return BaseClass::get_hessian_diagonal(q_point)[0];
5657}
5658
5659
5660
5661template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5662inline VectorizedArrayType
5664 const unsigned int q_point) const
5665{
5666 return BaseClass::get_laplacian(q_point)[0];
5667}
5668
5669
5670
5671template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5672inline void DEAL_II_ALWAYS_INLINE
5674 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5675{
5676# ifdef DEBUG
5677 this->dof_values_initialized = true;
5678 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5679# endif
5680 this->values_dofs[dof] = val_in;
5681}
5682
5683
5684
5685template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5686inline void DEAL_II_ALWAYS_INLINE
5688 const VectorizedArrayType val_in,
5689 const unsigned int q_point)
5690{
5691# ifdef DEBUG
5692 Assert(this->is_reinitialized, ExcNotInitialized());
5693# endif
5694 AssertIndexRange(q_point, this->n_quadrature_points);
5695 Assert(this->J_value != nullptr,
5697 "update_value"));
5698# ifdef DEBUG
5699 this->values_quad_submitted = true;
5700# endif
5701
5702 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5703 {
5704 const VectorizedArrayType JxW =
5705 this->J_value[0] * this->quadrature_weights[q_point];
5706 this->values_quad[q_point] = val_in * JxW;
5707 }
5708 else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5709 {
5710 this->values_quad[q_point] = val_in * this->J_value[q_point];
5711 }
5712}
5713
5714
5715
5716template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5717inline DEAL_II_ALWAYS_INLINE void
5720 const unsigned int q_point)
5721{
5722 submit_value(val_in[0], q_point);
5723}
5724
5725
5726
5727template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5728inline DEAL_II_ALWAYS_INLINE void
5730 submit_normal_derivative(const VectorizedArrayType grad_in,
5731 const unsigned int q_point)
5732{
5734 grad[0] = grad_in;
5735 BaseClass::submit_normal_derivative(grad, q_point);
5736}
5737
5738
5739
5740template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5741inline DEAL_II_ALWAYS_INLINE void
5744 const unsigned int q_point_in)
5745{
5746# ifdef DEBUG
5747 Assert(this->is_reinitialized, ExcNotInitialized());
5748# endif
5749 AssertIndexRange(q_point_in, this->n_quadrature_points);
5750 Assert(this->J_value != nullptr,
5752 "update_gradients"));
5753 Assert(this->jacobian != nullptr,
5755 "update_gradients"));
5756# ifdef DEBUG
5757 this->gradients_quad_submitted = true;
5758# endif
5759
5760 const std::size_t q_point = q_point_in;
5761 VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point;
5762 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5763 {
5764 const VectorizedArrayType JxW =
5765 this->J_value[0] * this->quadrature_weights[q_point];
5766
5767 // Make sure the compiler does not think 'jacobian' is aliased with
5768 // 'gradients_quad'
5769 std::array<VectorizedArrayType, dim> jac;
5770 for (unsigned int d = 0; d < dim; ++d)
5771 jac[d] = this->jacobian[0][d][d];
5772
5773 for (unsigned int d = 0; d < dim; ++d)
5774 grad_ptr[d] = grad_in[d] * jac[d] * JxW;
5775 }
5776 // general/affine cell type
5777 else
5778 {
5780 this->cell_type > internal::MatrixFreeFunctions::affine ?
5781 this->jacobian[q_point] :
5782 this->jacobian[0];
5783 const VectorizedArrayType JxW =
5784 this->cell_type > internal::MatrixFreeFunctions::affine ?
5785 this->J_value[q_point] :
5786 this->J_value[0] * this->quadrature_weights[q_point];
5787 for (unsigned int d = 0; d < dim; ++d)
5788 {
5789 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5790 for (unsigned int e = 1; e < dim; ++e)
5791 new_val += jac[e][d] * grad_in[e];
5792 grad_ptr[d] = new_val * JxW;
5793 }
5794 }
5795}
5796
5797
5798
5799template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5800inline DEAL_II_ALWAYS_INLINE void
5803 const unsigned int q_point)
5804{
5806 hessian[0] = hessian_in;
5807 BaseClass::submit_hessian(hessian, q_point);
5808}
5809
5810
5811
5812template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5813inline VectorizedArrayType
5815 integrate_value() const
5816{
5817 return BaseClass::integrate_value()[0];
5818}
5819
5820
5821
5822/*----------------- FEEvaluationAccess vector-valued ------------------------*/
5823
5824
5825template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5829 const unsigned int dof_no,
5830 const unsigned int first_selected_component,
5831 const unsigned int quad_no,
5832 const unsigned int fe_degree,
5833 const unsigned int n_q_points,
5834 const bool is_interior_face,
5835 const unsigned int active_fe_index,
5836 const unsigned int active_quad_index,
5837 const unsigned int face_type)
5838 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5839 matrix_free,
5840 dof_no,
5841 first_selected_component,
5842 quad_no,
5843 fe_degree,
5844 n_q_points,
5845 is_interior_face,
5846 active_fe_index,
5847 active_quad_index,
5848 face_type)
5849{}
5850
5851
5852
5853template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5856 const Mapping<dim> &mapping,
5857 const FiniteElement<dim> &fe,
5858 const Quadrature<1> &quadrature,
5859 const UpdateFlags update_flags,
5860 const unsigned int first_selected_component,
5862 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5863 mapping,
5864 fe,
5865 quadrature,
5866 update_flags,
5867 first_selected_component,
5868 other)
5869{}
5870
5871
5872
5873template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5877 &other)
5878 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5879{}
5880
5881
5882
5883template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5887 &other)
5888{
5890 operator=(other);
5891 return *this;
5892}
5893
5894
5895template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5898 const unsigned int q_point) const
5899{
5900 if (this->data->element_type ==
5902 {
5903 // Piola transform is required
5904# ifdef DEBUG
5905 Assert(this->values_quad_initialized == true,
5907# endif
5908
5909 AssertIndexRange(q_point, this->n_quadrature_points);
5910 Assert(this->J_value != nullptr,
5912 "update_values"));
5913 const std::size_t nqp = this->n_quadrature_points;
5915
5916 if (!is_face &&
5918 {
5919 // Cartesian cell
5920 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5921 const VectorizedArrayType inv_det =
5922 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5923 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5924 this->jacobian[0][2][2];
5925
5926 // J * u * det(J^-1)
5927 for (unsigned int comp = 0; comp < n_components; ++comp)
5928 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5929 jac[comp][comp] * inv_det;
5930 }
5931 else
5932 {
5933 // Affine or general cell
5934 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
5935 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5936 this->jacobian[q_point] :
5937 this->jacobian[0];
5939 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5940 transpose(invert(inv_t_jac)) :
5941 this->jacobian[1];
5942
5943 // Derivatives are reordered for faces. Need to take this into account
5944 const VectorizedArrayType inv_det =
5945 (is_face && dim == 2 && this->get_face_no() < 2) ?
5946 -determinant(inv_t_jac) :
5947 determinant(inv_t_jac);
5948 // J * u * det(J^-1)
5949 for (unsigned int comp = 0; comp < n_components; ++comp)
5950 {
5951 value_out[comp] = this->values_quad[q_point] * jac[comp][0];
5952 for (unsigned int e = 1; e < dim; ++e)
5953 value_out[comp] +=
5954 this->values_quad[e * nqp + q_point] * jac[comp][e];
5955 value_out[comp] *= inv_det;
5956 }
5957 }
5958 return value_out;
5959 }
5960 else
5961 {
5962 // No Piola needed
5963 return BaseClass::get_value(q_point);
5964 }
5965}
5966
5967
5968
5969template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5972 get_gradient(const unsigned int q_point) const
5973{
5974 if (this->data->element_type ==
5976 {
5977 // Piola transform is required
5978# ifdef DEBUG
5979 Assert(this->gradients_quad_initialized == true,
5981# endif
5982
5983 AssertIndexRange(q_point, this->n_quadrature_points);
5984 Assert(this->jacobian != nullptr,
5986 "update_gradients"));
5987 const std::size_t nqp = this->n_quadrature_points;
5988 const std::size_t nqp_d = nqp * dim;
5990 const VectorizedArrayType *gradients =
5991 this->gradients_quad + q_point * dim;
5992
5993
5994 if (!is_face &&
5996 {
5997 // Cartesian cell
5998 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5999 this->jacobian[0];
6000 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6001 const VectorizedArrayType inv_det =
6002 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
6003 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
6004 this->jacobian[0][2][2];
6005
6006 // J * grad_quad * J^-1 * det(J^-1)
6007 for (unsigned int d = 0; d < dim; ++d)
6008 for (unsigned int comp = 0; comp < n_components; ++comp)
6009 grad_out[comp][d] = gradients[comp * nqp_d + d] *
6010 inv_t_jac[d][d] * (jac[comp][comp] * inv_det);
6011 }
6012 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6013 {
6014 // Affine cell
6015 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6016 this->jacobian[0];
6017 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6018
6019 // Derivatives are reordered for faces. Need to take this into account
6020 const VectorizedArrayType inv_det =
6021 (is_face && dim == 2 && this->get_face_no() < 2) ?
6022 -determinant(inv_t_jac) :
6023 determinant(inv_t_jac);
6024
6025 VectorizedArrayType tmp[dim][dim];
6026 // J * grad_quad * J^-1 * det(J^-1)
6027 for (unsigned int d = 0; d < dim; ++d)
6028 for (unsigned int e = 0; e < dim; ++e)
6029 {
6030 tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0];
6031 for (unsigned int f = 1; f < dim; ++f)
6032 tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f];
6033 }
6034 for (unsigned int comp = 0; comp < n_components; ++comp)
6035 for (unsigned int d = 0; d < dim; ++d)
6036 {
6037 VectorizedArrayType res = jac[comp][0] * tmp[d][0];
6038 for (unsigned int f = 1; f < dim; ++f)
6039 res += jac[comp][f] * tmp[d][f];
6040
6041 grad_out[comp][d] = res * inv_det;
6042 }
6043 }
6044 else
6045 {
6046 // General cell
6047
6048 // This assert could be removed if we make sure that this is updated
6049 // even though update_hessians or update_jacobian_grads is not passed,
6050 // i.e make the necessary changes in
6051 // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
6052 Assert(this->jacobian_gradients_non_inverse != nullptr,
6054 "update_hessians"));
6055
6056 const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
6057 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
6058 this->jacobian[q_point];
6059
6060 // Derivatives are reordered for faces. Need to take this into account
6061 const VectorizedArrayType inv_det =
6062 (is_face && dim == 2 && this->get_face_no() < 2) ?
6063 -determinant(inv_t_jac) :
6064 determinant(inv_t_jac);
6065 const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac);
6066
6067 // (J * grad_quad) * J^-1 * det(J^-1), part in braces
6068 VectorizedArrayType tmp[dim][dim];
6069 for (unsigned int d = 0; d < dim; ++d)
6070 for (unsigned int e = 0; e < dim; ++e)
6071 {
6072 tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e];
6073 for (unsigned int f = 1; f < dim; ++f)
6074 tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e];
6075 }
6076
6077 // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms
6078 // outside braces with gradient part from above
6079 for (unsigned int d = 0; d < dim; ++d)
6080 {
6081 for (unsigned int e = 0; e < dim; ++e)
6082 tmp[e][d] +=
6083 jac_grad[e][d] * this->values_quad[e * nqp + q_point];
6084 for (unsigned int f = 0, r = dim; f < dim; ++f)
6085 for (unsigned int k = f + 1; k < dim; ++k, ++r)
6086 {
6087 tmp[k][d] +=
6088 jac_grad[r][d] * this->values_quad[f * nqp + q_point];
6089 tmp[f][d] +=
6090 jac_grad[r][d] * this->values_quad[k * nqp + q_point];
6091 }
6092 }
6093
6094 // Apply J^{-1} appearing in both terms outside braces above
6095 for (unsigned int d = 0; d < dim; ++d)
6096 for (unsigned int e = 0; e < dim; ++e)
6097 {
6098 VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0];
6099 for (unsigned int f = 1; f < dim; ++f)
6100 res += tmp[f][d] * inv_t_jac[e][f];
6101 grad_out[d][e] = res;
6102 }
6103
6104 // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})),
6105 // which can be expressed as a rank-1 update tmp[d] * tmp4[e], where
6106 // tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1})
6107 VectorizedArrayType tmp3[dim], tmp4[dim];
6108 for (unsigned int d = 0; d < dim; ++d)
6109 {
6110 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
6111 for (unsigned int e = 1; e < dim; ++e)
6112 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
6113 }
6114 for (unsigned int e = 0, k = dim; e < dim; ++e)
6115 for (unsigned int f = e + 1; f < dim; ++k, ++f)
6116 for (unsigned int d = 0; d < dim; ++d)
6117 {
6118 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
6119 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
6120 }
6121 for (unsigned int d = 0; d < dim; ++d)
6122 {
6123 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
6124 for (unsigned int e = 1; e < dim; ++e)
6125 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
6126 }
6127
6128 VectorizedArrayType tmp2[dim];
6129 for (unsigned int d = 0; d < dim; ++d)
6130 {
6131 tmp2[d] = t_jac[0][d] * this->values_quad[q_point];
6132 for (unsigned e = 1; e < dim; ++e)
6133 tmp2[d] += t_jac[e][d] * this->values_quad[e * nqp + q_point];
6134 }
6135
6136 for (unsigned int d = 0; d < dim; ++d)
6137 for (unsigned int e = 0; e < dim; ++e)
6138 {
6139 grad_out[d][e] -= tmp4[e] * tmp2[d];
6140
6141 // finally multiply by det(J^{-1}) necessary in all
6142 // contributions above
6143 grad_out[d][e] *= inv_det;
6144 }
6145 }
6146 return grad_out;
6147 }
6148 else
6149 {
6150 return BaseClass::get_gradient(q_point);
6151 }
6152}
6153
6154
6155
6156template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6157inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6159 get_divergence(const unsigned int q_point) const
6160{
6161# ifdef DEBUG
6162 Assert(this->gradients_quad_initialized == true,
6164# endif
6165 AssertIndexRange(q_point, this->n_quadrature_points);
6166 Assert(this->jacobian != nullptr,
6168 "update_gradients"));
6169
6170 VectorizedArrayType divergence;
6171 const std::size_t nqp = this->n_quadrature_points;
6172
6173 if (this->data->element_type ==
6175 {
6176 VectorizedArrayType inv_det =
6177 (!is_face &&
6178 this->cell_type == internal::MatrixFreeFunctions::cartesian) ?
6179 this->jacobian[0][0][0] *
6180 ((dim == 2) ? this->jacobian[0][1][1] :
6181 this->jacobian[0][1][1] * this->jacobian[0][2][2]) :
6182 determinant(this->jacobian[this->cell_type >
6183 internal::MatrixFreeFunctions::affine ?
6184 q_point :
6185 0]);
6186
6187 // on faces in 2d, the determinant has the wrong sign due to ordering of
6188 // derivatives
6189 if (is_face && dim == 2 && this->get_face_no() < 2)
6190 inv_det = -inv_det;
6191
6192 // div * det(J^-1)
6193 divergence = this->gradients_quad[q_point * dim];
6194 for (unsigned int d = 1; d < dim; ++d)
6195 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d];
6196 divergence *= inv_det;
6197 }
6198 else
6199 {
6200 if (!is_face &&
6202 {
6203 // Cartesian cell
6204 divergence =
6205 this->gradients_quad[q_point * dim] * this->jacobian[0][0][0];
6206 for (unsigned int d = 1; d < dim; ++d)
6207 divergence += this->gradients_quad[(d * nqp + q_point) * dim + d] *
6208 this->jacobian[0][d][d];
6209 }
6210 else
6211 {
6212 // cell with general/constant Jacobian
6214 this->cell_type == internal::MatrixFreeFunctions::general ?
6215 this->jacobian[q_point] :
6216 this->jacobian[0];
6217 divergence = jac[0][0] * this->gradients_quad[q_point * dim];
6218 for (unsigned int e = 1; e < dim; ++e)
6219 divergence += jac[0][e] * this->gradients_quad[q_point * dim + e];
6220 for (unsigned int d = 1; d < dim; ++d)
6221 for (unsigned int e = 0; e < dim; ++e)
6222 divergence +=
6223 jac[d][e] * this->gradients_quad[(d * nqp + q_point) * dim + e];
6224 }
6225 }
6226 return divergence;
6227}
6228
6229
6230
6231template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6234 get_symmetric_gradient(const unsigned int q_point) const
6235{
6236 // copy from generic function into dim-specialization function
6237 const auto grad = get_gradient(q_point);
6238 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6239 VectorizedArrayType half = Number(0.5);
6240 for (unsigned int d = 0; d < dim; ++d)
6241 symmetrized[d] = grad[d][d];
6242 switch (dim)
6243 {
6244 case 1:
6245 break;
6246 case 2:
6247 symmetrized[2] = grad[0][1] + grad[1][0];
6248 symmetrized[2] *= half;
6249 break;
6250 case 3:
6251 symmetrized[3] = grad[0][1] + grad[1][0];
6252 symmetrized[3] *= half;
6253 symmetrized[4] = grad[0][2] + grad[2][0];
6254 symmetrized[4] *= half;
6255 symmetrized[5] = grad[1][2] + grad[2][1];
6256 symmetrized[5] *= half;
6257 break;
6258 default:
6260 }
6262}
6263
6264
6265
6266template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6268 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6270 const unsigned int q_point) const
6271{
6272 // copy from generic function into dim-specialization function
6273 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6274 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6275 switch (dim)
6276 {
6277 case 1:
6278 Assert(false,
6279 ExcMessage(
6280 "Computing the curl in 1d is not a useful operation"));
6281 break;
6282 case 2:
6283 curl[0] = grad[1][0] - grad[0][1];
6284 break;
6285 case 3:
6286 curl[0] = grad[2][1] - grad[1][2];
6287 curl[1] = grad[0][2] - grad[2][0];
6288 curl[2] = grad[1][0] - grad[0][1];
6289 break;
6290 default:
6292 }
6293 return curl;
6294}
6295
6296
6297
6298template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6301 get_hessian_diagonal(const unsigned int q_point) const
6302{
6303 return BaseClass::get_hessian_diagonal(q_point);
6304}
6305
6306
6307
6308template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6311 const unsigned int q_point) const
6312{
6313# ifdef DEBUG
6314 Assert(this->hessians_quad_initialized == true,
6316# endif
6317 AssertIndexRange(q_point, this->n_quadrature_points);
6318 return BaseClass::get_hessian(q_point);
6319}
6320
6321
6322template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6323inline DEAL_II_ALWAYS_INLINE void
6326 const unsigned int q_point)
6327{
6328 if (this->data->element_type ==
6330 {
6331 // Piola transform is required
6332 AssertIndexRange(q_point, this->n_quadrature_points);
6333 Assert(this->J_value != nullptr,
6335 "update_value"));
6336# ifdef DEBUG
6337 Assert(this->is_reinitialized, ExcNotInitialized());
6338 this->values_quad_submitted = true;
6339# endif
6340
6341 VectorizedArrayType *values = this->values_quad + q_point;
6342 const std::size_t nqp = this->n_quadrature_points;
6343
6344 if (!is_face &&
6346 {
6347 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6348 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6349
6350 for (unsigned int comp = 0; comp < n_components; ++comp)
6351 values[comp * nqp] = val_in[comp] * weight * jac[comp][comp];
6352 }
6353 else
6354 {
6355 // Affine or general cell
6356 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
6357 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6358 this->jacobian[q_point] :
6359 this->jacobian[0];
6360
6361 // Derivatives are reordered for faces. Need to take this into account
6362 // and 1/inv_det != J_value for faces
6363 const VectorizedArrayType fac =
6364 (!is_face) ?
6365 this->quadrature_weights[q_point] :
6366 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6367 this->J_value[q_point] :
6368 this->J_value[0] * this->quadrature_weights[q_point]) *
6369 ((dim == 2 && this->get_face_no() < 2) ?
6370 -determinant(inv_t_jac) :
6371 determinant(inv_t_jac)));
6373 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6374 transpose(invert(inv_t_jac)) :
6375 this->jacobian[1];
6376
6377 // J^T * u * factor
6378 for (unsigned int comp = 0; comp < n_components; ++comp)
6379 {
6380 values[comp * nqp] = val_in[0] * jac[0][comp];
6381 for (unsigned int e = 1; e < dim; ++e)
6382 values[comp * nqp] += val_in[e] * jac[e][comp];
6383 values[comp * nqp] *= fac;
6384 }
6385 }
6386 }
6387 else
6388 {
6389 // No Piola transform
6390 BaseClass::submit_value(val_in, q_point);
6391 }
6392}
6393
6394
6395
6396template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6397inline DEAL_II_ALWAYS_INLINE void
6400 const unsigned int q_point)
6401{
6402 if (this->data->element_type ==
6404 {
6405 // Piola transform is required
6406
6407# ifdef DEBUG
6408 Assert(this->is_reinitialized, ExcNotInitialized());
6409# endif
6410 AssertIndexRange(q_point, this->n_quadrature_points);
6411 Assert(this->J_value != nullptr,
6413 "update_gradients"));
6414 Assert(this->jacobian != nullptr,
6416 "update_gradients"));
6417# ifdef DEBUG
6418 this->gradients_quad_submitted = true;
6419# endif
6420
6421 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
6422 VectorizedArrayType *values = this->values_from_gradients_quad + q_point;
6423 const std::size_t nqp = this->n_quadrature_points;
6424 const std::size_t nqp_d = nqp * dim;
6425
6426 if (!is_face &&
6428 {
6429 // Cartesian cell
6430 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6431 this->jacobian[0];
6432 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6433 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6434 for (unsigned int d = 0; d < dim; ++d)
6435 for (unsigned int comp = 0; comp < n_components; ++comp)
6436 gradients[comp * nqp_d + d] =
6437 grad_in[comp][d] * inv_t_jac[d][d] * (jac[comp][comp] * weight);
6438 }
6439 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6440 {
6441 // Affine cell
6442 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6443 this->jacobian[0];
6444 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6445
6446 // Derivatives are reordered for faces. Need to take this into account
6447 // and 1/inv_det != J_value for faces
6448 const VectorizedArrayType fac =
6449 (!is_face) ? this->quadrature_weights[q_point] :
6450 this->J_value[0] * this->quadrature_weights[q_point] *
6451 ((dim == 2 && this->get_face_no() < 2) ?
6452 -determinant(inv_t_jac) :
6453 determinant(inv_t_jac));
6454
6455 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6456 VectorizedArrayType tmp[dim][dim];
6457 for (unsigned int d = 0; d < dim; ++d)
6458 for (unsigned int e = 0; e < dim; ++e)
6459 {
6460 tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0];
6461 for (unsigned int f = 1; f < dim; ++f)
6462 tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f];
6463 }
6464 for (unsigned int comp = 0; comp < n_components; ++comp)
6465 for (unsigned int d = 0; d < dim; ++d)
6466 {
6467 VectorizedArrayType res = jac[0][comp] * tmp[d][0];
6468 for (unsigned int f = 1; f < dim; ++f)
6469 res += jac[f][comp] * tmp[d][f];
6470
6471 gradients[comp * nqp_d + d] = res * fac;
6472 }
6473 }
6474 else
6475 {
6476 // General cell
6477
6478 const auto jac_grad = this->jacobian_gradients_non_inverse[q_point];
6479 const Tensor<2, dim, VectorizedArrayType> inv_t_jac =
6480 this->jacobian[q_point];
6481
6482 // Derivatives are reordered for faces. Need to take this into account
6483 // and 1/inv_det != J_value for faces
6484 const VectorizedArrayType fac =
6485 (!is_face) ?
6486 this->quadrature_weights[q_point] :
6487 this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6488 -determinant(inv_t_jac) :
6489 determinant(inv_t_jac));
6490 const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac);
6491
6492 // Start evaluation for values part below to enable the compiler to
6493 // possibly re-use the same computation in get_gradient() without
6494 // interfering with stores to 'gradients'
6495 VectorizedArrayType tmp3[dim], tmp4[dim];
6496 for (unsigned int d = 0; d < dim; ++d)
6497 {
6498 tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0];
6499 for (unsigned int e = 1; e < dim; ++e)
6500 tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e];
6501 }
6502 for (unsigned int e = 0, k = dim; e < dim; ++e)
6503 for (unsigned int f = e + 1; f < dim; ++k, ++f)
6504 for (unsigned int d = 0; d < dim; ++d)
6505 {
6506 tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d];
6507 tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d];
6508 }
6509 for (unsigned int d = 0; d < dim; ++d)
6510 {
6511 tmp4[d] = tmp3[0] * inv_t_jac[d][0];
6512 for (unsigned int e = 1; e < dim; ++e)
6513 tmp4[d] += tmp3[e] * inv_t_jac[d][e];
6514 }
6515
6516 const Tensor<2, dim, VectorizedArrayType> grad_in_scaled =
6517 fac * grad_in;
6518
6519 VectorizedArrayType tmp[dim][dim];
6520
6521 // J * (J^{-1} * (grad_in * factor))
6522 for (unsigned int d = 0; d < dim; ++d)
6523 for (unsigned int e = 0; e < dim; ++e)
6524 {
6525 tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0];
6526 for (unsigned int f = 1; f < dim; ++f)
6527 tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f];
6528 }
6529
6530 for (unsigned int d = 0; d < dim; ++d)
6531 for (unsigned int e = 0; e < dim; ++e)
6532 {
6533 VectorizedArrayType res = t_jac[d][0] * tmp[e][0];
6534 for (unsigned int f = 1; f < dim; ++f)
6535 res += t_jac[d][f] * tmp[e][f];
6536
6537 gradients[d * nqp_d + e] = res;
6538 }
6539
6540 // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces
6541 // as 'tmp' from above
6542 VectorizedArrayType value[dim];
6543 for (unsigned int d = 0; d < dim; ++d)
6544 {
6545 value[d] = tmp[d][0] * jac_grad[d][0];
6546 for (unsigned int e = 1; e < dim; ++e)
6547 value[d] += tmp[d][e] * jac_grad[d][e];
6548 }
6549 for (unsigned int e = 0, k = dim; e < dim; ++e)
6550 for (unsigned int f = e + 1; f < dim; ++k, ++f)
6551 for (unsigned int d = 0; d < dim; ++d)
6552 {
6553 value[e] += tmp[f][d] * jac_grad[k][d];
6554 value[f] += tmp[e][d] * jac_grad[k][d];
6555 }
6556
6557 // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1})
6558 // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ )
6559 for (unsigned int d = 0; d < dim; ++d)
6560 {
6561 VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0];
6562 for (unsigned int e = 1; e < dim; ++e)
6563 tmp2 += grad_in_scaled[d][e] * tmp4[e];
6564 for (unsigned int e = 0; e < dim; ++e)
6565 value[e] -= t_jac[e][d] * tmp2;
6566 }
6567
6568 for (unsigned int d = 0; d < dim; ++d)
6569 values[d * nqp] = value[d];
6570 }
6571 }
6572 else
6573 {
6574 BaseClass::submit_gradient(grad_in, q_point);
6575 }
6576}
6577
6578
6579
6580template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6581inline DEAL_II_ALWAYS_INLINE void
6584 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6585 const unsigned int q_point)
6586{
6587 if (this->data->element_type ==
6589 {
6590 // Piola transform is required
6591 const Tensor<2, dim, VectorizedArrayType> grad = grad_in;
6593 submit_gradient(grad, q_point);
6594 }
6595 else
6596 {
6597 BaseClass::submit_gradient(grad_in, q_point);
6598 }
6599}
6600
6601
6602
6603template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6604inline DEAL_II_ALWAYS_INLINE void
6606 submit_divergence(const VectorizedArrayType div_in,
6607 const unsigned int q_point)
6608{
6609# ifdef DEBUG
6610 Assert(this->is_reinitialized, ExcNotInitialized());
6611# endif
6612 AssertIndexRange(q_point, this->n_quadrature_points);
6613 Assert(this->J_value != nullptr,
6615 "update_gradients"));
6616 Assert(this->jacobian != nullptr,
6618 "update_gradients"));
6619# ifdef DEBUG
6620 this->gradients_quad_submitted = true;
6621# endif
6622
6623 const std::size_t nqp_d = this->n_quadrature_points * dim;
6624 VectorizedArrayType *gradients = this->gradients_quad + q_point * dim;
6625
6626 if (this->data->element_type ==
6628 {
6629 // General cell
6630
6631 // Derivatives are reordered for faces. Need to take this into account
6632 // and 1/inv_det != J_value for faces
6633 const VectorizedArrayType fac =
6634 (!is_face) ?
6635 this->quadrature_weights[q_point] * div_in :
6636 (this->cell_type > internal::MatrixFreeFunctions::affine ?
6637 this->J_value[q_point] :
6638 this->J_value[0] * this->quadrature_weights[q_point]) *
6639 div_in *
6641 this->jacobian[this->cell_type >
6642 internal::MatrixFreeFunctions::affine ?
6643 q_point :
6644 0]) *
6645 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6646
6647 for (unsigned int d = 0; d < dim; ++d)
6648 {
6649 for (unsigned int e = 0; e < dim; ++e)
6650 gradients[d * nqp_d + e] = (d == e) ? fac : 0.;
6651 }
6652 this->divergence_is_requested = true;
6653 }
6654 else
6655 {
6656 if (!is_face &&
6658 {
6659 const VectorizedArrayType fac =
6660 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6661 for (unsigned int d = 0; d < dim; ++d)
6662 {
6663 const VectorizedArrayType jac_dd = this->jacobian[0][d][d];
6664 for (unsigned int e = 0; e < dim; ++e)
6665 gradients[d * nqp_d + e] = (d == e) ? fac * jac_dd : 0.;
6666 }
6667 }
6668 else
6669 {
6671 this->cell_type == internal::MatrixFreeFunctions::general ?
6672 this->jacobian[q_point] :
6673 this->jacobian[0];
6674 const VectorizedArrayType fac =
6675 (this->cell_type == internal::MatrixFreeFunctions::general ?
6676 this->J_value[q_point] :
6677 this->J_value[0] * this->quadrature_weights[q_point]) *
6678 div_in;
6679 for (unsigned int d = 0; d < dim; ++d)
6680 {
6681 for (unsigned int e = 0; e < dim; ++e)
6682 gradients[d * nqp_d + e] = jac[d][e] * fac;
6683 }
6684 }
6685 }
6686}
6687
6688
6689
6690template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6691inline DEAL_II_ALWAYS_INLINE void
6695 const unsigned int q_point)
6696{
6698 this->data->element_type !=
6701
6702 // could have used base class operator, but that involves some overhead
6703 // which is inefficient. it is nice to have the symmetric tensor because
6704 // that saves some operations
6705# ifdef DEBUG
6706 Assert(this->is_reinitialized, ExcNotInitialized());
6707# endif
6708 AssertIndexRange(q_point, this->n_quadrature_points);
6709 Assert(this->J_value != nullptr,
6711 "update_gradients"));
6712 Assert(this->jacobian != nullptr,
6714 "update_gradients"));
6715# ifdef DEBUG
6716 this->gradients_quad_submitted = true;
6717# endif
6718
6719 const std::size_t nqp_d = this->n_quadrature_points * dim;
6720 VectorizedArrayType *gradients = this->gradients_quad + dim * q_point;
6721 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6722 {
6723 const VectorizedArrayType JxW =
6724 this->J_value[0] * this->quadrature_weights[q_point];
6725 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
6726 for (unsigned int d = 0; d < dim; ++d)
6727 gradients[d * nqp_d + d] =
6728 (sym_grad.access_raw_entry(d) * JxW * jac[d][d]);
6729 for (unsigned int e = 0, counter = dim; e < dim; ++e)
6730 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6731 {
6732 const VectorizedArrayType value =
6733 sym_grad.access_raw_entry(counter) * JxW;
6734 gradients[e * nqp_d + d] = value * jac[d][d];
6735 gradients[d * nqp_d + e] = value * jac[e][e];
6736 }
6737 }
6738 // general/affine cell type
6739 else
6740 {
6741 const VectorizedArrayType JxW =
6742 this->cell_type == internal::MatrixFreeFunctions::general ?
6743 this->J_value[q_point] :
6744 this->J_value[0] * this->quadrature_weights[q_point];
6746 this->cell_type == internal::MatrixFreeFunctions::general ?
6747 this->jacobian[q_point] :
6748 this->jacobian[0];
6749 VectorizedArrayType weighted[dim][dim];
6750 for (unsigned int i = 0; i < dim; ++i)
6751 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6752 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6753 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6754 {
6755 const VectorizedArrayType value =
6756 sym_grad.access_raw_entry(counter) * JxW;
6757 weighted[i][j] = value;
6758 weighted[j][i] = value;
6759 }
6760 for (unsigned int comp = 0; comp < dim; ++comp)
6761 for (unsigned int d = 0; d < dim; ++d)
6762 {
6763 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6764 for (unsigned int e = 1; e < dim; ++e)
6765 new_val += jac[e][d] * weighted[comp][e];
6766 gradients[comp * nqp_d + d] = new_val;
6767 }
6768 }
6769}
6770
6771
6772
6773template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6774inline DEAL_II_ALWAYS_INLINE void
6777 const unsigned int q_point)
6778{
6780 switch (dim)
6781 {
6782 case 1:
6783 Assert(false,
6784 ExcMessage(
6785 "Testing by the curl in 1d is not a useful operation"));
6786 break;
6787 case 2:
6788 grad[1][0] = curl[0];
6789 grad[0][1] = -curl[0];
6790 break;
6791 case 3:
6792 grad[2][1] = curl[0];
6793 grad[1][2] = -curl[0];
6794 grad[0][2] = curl[1];
6795 grad[2][0] = -curl[1];
6796 grad[1][0] = curl[2];
6797 grad[0][1] = -curl[2];
6798 break;
6799 default:
6801 }
6802 submit_gradient(grad, q_point);
6803}
6804
6805
6806/*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6807
6808
6809template <typename Number, bool is_face, typename VectorizedArrayType>
6813 const unsigned int dof_no,
6814 const unsigned int first_selected_component,
6815 const unsigned int quad_no,
6816 const unsigned int fe_degree,
6817 const unsigned int n_q_points,
6818 const bool is_interior_face,
6819 const unsigned int active_fe_index,
6820 const unsigned int active_quad_index,
6821 const unsigned int face_type)
6822 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6823 matrix_free,
6824 dof_no,
6825 first_selected_component,
6826 quad_no,
6827 fe_degree,
6828 n_q_points,
6829 is_interior_face,
6830 active_fe_index,
6831 active_quad_index,
6832 face_type)
6833{}
6834
6835
6836
6837template <typename Number, bool is_face, typename VectorizedArrayType>
6840 const Mapping<1> &mapping,
6841 const FiniteElement<1> &fe,
6842 const Quadrature<1> &quadrature,
6843 const UpdateFlags update_flags,
6844 const unsigned int first_selected_component,
6846 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6847 mapping,
6848 fe,
6849 quadrature,
6850 update_flags,
6851 first_selected_component,
6852 other)
6853{}
6854
6855
6856
6857template <typename Number, bool is_face, typename VectorizedArrayType>
6861 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6862{}
6863
6864
6865
6866template <typename Number, bool is_face, typename VectorizedArrayType>
6870{
6872 other);
6873 return *this;
6874}
6875
6876
6877
6878template <typename Number, bool is_face, typename VectorizedArrayType>
6879inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6881 const unsigned int dof) const
6882{
6883 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6884 return this->values_dofs[dof];
6885}
6886
6887
6888
6889template <typename Number, bool is_face, typename VectorizedArrayType>
6890inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6892 const unsigned int q_point) const
6893{
6894# ifdef DEBUG
6895 Assert(this->values_quad_initialized == true,
6897# endif
6898 AssertIndexRange(q_point, this->n_quadrature_points);
6899 return this->values_quad[q_point];
6900}
6901
6902
6903
6904template <typename Number, bool is_face, typename VectorizedArrayType>
6907 const unsigned int q_point) const
6908{
6909 // could use the base class gradient, but that involves too many inefficient
6910 // initialization operations on tensors
6911
6912# ifdef DEBUG
6913 Assert(this->gradients_quad_initialized == true,
6915# endif
6916 AssertIndexRange(q_point, this->n_quadrature_points);
6917
6919 this->cell_type == internal::MatrixFreeFunctions::general ?
6920 this->jacobian[q_point] :
6921 this->jacobian[0];
6922
6924 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6925
6926 return grad_out;
6927}
6928
6929
6930
6931template <typename Number, bool is_face, typename VectorizedArrayType>
6932inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6934 const unsigned int q_point) const
6935{
6936 return get_gradient(q_point)[0];
6937}
6938
6939
6940
6941template <typename Number, bool is_face, typename VectorizedArrayType>
6942inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6944 get_normal_derivative(const unsigned int q_point) const
6945{
6946 return BaseClass::get_normal_derivative(q_point)[0];
6947}
6948
6949
6950
6951template <typename Number, bool is_face, typename VectorizedArrayType>
6954 const unsigned int q_point) const
6955{
6956 return BaseClass::get_hessian(q_point)[0];
6957}
6958
6959
6960
6961template <typename Number, bool is_face, typename VectorizedArrayType>
6964 get_hessian_diagonal(const unsigned int q_point) const
6965{
6966 return BaseClass::get_hessian_diagonal(q_point)[0];
6967}
6968
6969
6970
6971template <typename Number, bool is_face, typename VectorizedArrayType>
6972inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6974 const unsigned int q_point) const
6975{
6976 return BaseClass::get_laplacian(q_point)[0];
6977}
6978
6979
6980
6981template <typename Number, bool is_face, typename VectorizedArrayType>
6984 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6985{
6986# ifdef DEBUG
6987 this->dof_values_initialized = true;
6988 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6989# endif
6990 this->values_dofs[dof] = val_in;
6991}
6992
6993
6994
6995template <typename Number, bool is_face, typename VectorizedArrayType>
6996inline DEAL_II_ALWAYS_INLINE void
6998 const VectorizedArrayType val_in,
6999 const unsigned int q_point)
7000{
7001# ifdef DEBUG
7002 Assert(this->is_reinitialized, ExcNotInitialized());
7003# endif
7004 AssertIndexRange(q_point, this->n_quadrature_points);
7005# ifdef DEBUG
7006 this->values_quad_submitted = true;
7007# endif
7008
7009 if (this->cell_type == internal::MatrixFreeFunctions::general)
7010 {
7011 const VectorizedArrayType JxW = this->J_value[q_point];
7012 this->values_quad[q_point] = val_in * JxW;
7013 }
7014 else // if (this->cell_type == internal::MatrixFreeFunctions::general)
7015 {
7016 const VectorizedArrayType JxW =
7017 this->J_value[0] * this->quadrature_weights[q_point];
7018 this->values_quad[q_point] = val_in * JxW;
7019 }
7020}
7021
7022
7023
7024template <typename Number, bool is_face, typename VectorizedArrayType>
7025inline DEAL_II_ALWAYS_INLINE void
7028 const unsigned int q_point)
7029{
7030 submit_value(val_in[0], q_point);
7031}
7032
7033
7034
7035template <typename Number, bool is_face, typename VectorizedArrayType>
7036inline DEAL_II_ALWAYS_INLINE void
7039 const unsigned int q_point)
7040{
7041 submit_gradient(grad_in[0], q_point);
7042}
7043
7044
7045
7046template <typename Number, bool is_face, typename VectorizedArrayType>
7047inline DEAL_II_ALWAYS_INLINE void
7049 const VectorizedArrayType grad_in,
7050 const unsigned int q_point)
7051{
7052# ifdef DEBUG
7053 Assert(this->is_reinitialized, ExcNotInitialized());
7054# endif
7055 AssertIndexRange(q_point, this->n_quadrature_points);
7056# ifdef DEBUG
7057 this->gradients_quad_submitted = true;
7058# endif
7059
7061 this->cell_type == internal::MatrixFreeFunctions::general ?
7062 this->jacobian[q_point] :
7063 this->jacobian[0];
7064 const VectorizedArrayType JxW =
7065 this->cell_type == internal::MatrixFreeFunctions::general ?
7066 this->J_value[q_point] :
7067 this->J_value[0] * this->quadrature_weights[q_point];
7068
7069 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7070}
7071
7072
7073
7074template <typename Number, bool is_face, typename VectorizedArrayType>
7075inline DEAL_II_ALWAYS_INLINE void
7078 const unsigned int q_point)
7079{
7080 submit_gradient(grad_in[0][0], q_point);
7081}
7082
7083
7084
7085template <typename Number, bool is_face, typename VectorizedArrayType>
7086inline DEAL_II_ALWAYS_INLINE void
7088 submit_normal_derivative(const VectorizedArrayType grad_in,
7089 const unsigned int q_point)
7090{
7092 grad[0] = grad_in;
7093 BaseClass::submit_normal_derivative(grad, q_point);
7094}
7095
7096
7097
7098template <typename Number, bool is_face, typename VectorizedArrayType>
7099inline DEAL_II_ALWAYS_INLINE void
7102 const unsigned int q_point)
7103{
7104 BaseClass::submit_normal_derivative(grad_in, q_point);
7105}
7106
7107
7108template <typename Number, bool is_face, typename VectorizedArrayType>
7109inline DEAL_II_ALWAYS_INLINE void
7111 const Tensor<2, 1, VectorizedArrayType> hessian_in,
7112 const unsigned int q_point)
7113{
7115 hessian[0] = hessian_in;
7116 BaseClass::submit_hessian(hessian, q_point);
7117}
7118
7119
7120template <typename Number, bool is_face, typename VectorizedArrayType>
7121inline VectorizedArrayType
7123 integrate_value() const
7124{
7125 return BaseClass::integrate_value()[0];
7126}
7127
7128
7129
7130/*-------------------------- FEEvaluation -----------------------------------*/
7131
7132
7133template <int dim,
7134 int fe_degree,
7135 int n_q_points_1d,
7136 int n_components_,
7137 typename Number,
7138 typename VectorizedArrayType>
7139inline FEEvaluation<dim,
7140 fe_degree,
7141 n_q_points_1d,
7142 n_components_,
7143 Number,
7144 VectorizedArrayType>::
7145 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7146 const unsigned int fe_no,
7147 const unsigned int quad_no,
7148 const unsigned int first_selected_component,
7149 const unsigned int active_fe_index,
7150 const unsigned int active_quad_index)
7151 : BaseClass(matrix_free,
7152 fe_no,
7153 first_selected_component,
7154 quad_no,
7155 fe_degree,
7156 static_n_q_points,
7157 true /*note: this is not a face*/,
7158 active_fe_index,
7159 active_quad_index)
7160 , dofs_per_component(this->data->dofs_per_component_on_cell)
7161 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7162 , n_q_points(this->data->n_q_points)
7163{
7164 check_template_arguments(fe_no, 0);
7165}
7166
7167
7168
7169template <int dim,
7170 int fe_degree,
7171 int n_q_points_1d,
7172 int n_components_,
7173 typename Number,
7174 typename VectorizedArrayType>
7175inline FEEvaluation<dim,
7176 fe_degree,
7177 n_q_points_1d,
7178 n_components_,
7179 Number,
7180 VectorizedArrayType>::
7181 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7182 const std::pair<unsigned int, unsigned int> &range,
7183 const unsigned int dof_no,
7184 const unsigned int quad_no,
7185 const unsigned int first_selected_component)
7186 : FEEvaluation(matrix_free,
7187 dof_no,
7188 quad_no,
7189 first_selected_component,
7190 matrix_free.get_cell_active_fe_index(range))
7191{}
7192
7193
7194
7195template <int dim,
7196 int fe_degree,
7197 int n_q_points_1d,
7198 int n_components_,
7199 typename Number,
7200 typename VectorizedArrayType>
7201inline FEEvaluation<dim,
7202 fe_degree,
7203 n_q_points_1d,
7204 n_components_,
7205 Number,
7206 VectorizedArrayType>::
7207 FEEvaluation(const Mapping<dim> &mapping,
7208 const FiniteElement<dim> &fe,
7209 const Quadrature<1> &quadrature,
7210 const UpdateFlags update_flags,
7211 const unsigned int first_selected_component)
7212 : BaseClass(mapping,
7213 fe,
7214 quadrature,
7215 update_flags,
7216 first_selected_component,
7217 nullptr)
7218 , dofs_per_component(this->data->dofs_per_component_on_cell)
7219 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7220 , n_q_points(this->data->n_q_points)
7221{
7222 check_template_arguments(numbers::invalid_unsigned_int, 0);
7223}
7224
7225
7226
7227template <int dim,
7228 int fe_degree,
7229 int n_q_points_1d,
7230 int n_components_,
7231 typename Number,
7232 typename VectorizedArrayType>
7233inline FEEvaluation<dim,
7234 fe_degree,
7235 n_q_points_1d,
7236 n_components_,
7237 Number,
7238 VectorizedArrayType>::
7239 FEEvaluation(const FiniteElement<dim> &fe,
7240 const Quadrature<1> &quadrature,
7241 const UpdateFlags update_flags,
7242 const unsigned int first_selected_component)
7243 : BaseClass(StaticMappingQ1<dim>::mapping,
7244 fe,
7245 quadrature,
7246 update_flags,
7247 first_selected_component,
7248 nullptr)
7249 , dofs_per_component(this->data->dofs_per_component_on_cell)
7250 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7251 , n_q_points(this->data->n_q_points)
7252{
7253 check_template_arguments(numbers::invalid_unsigned_int, 0);
7254}
7255
7256
7257
7258template <int dim,
7259 int fe_degree,
7260 int n_q_points_1d,
7261 int n_components_,
7262 typename Number,
7263 typename VectorizedArrayType>
7264inline FEEvaluation<dim,
7265 fe_degree,
7266 n_q_points_1d,
7267 n_components_,
7268 Number,
7269 VectorizedArrayType>::
7270 FEEvaluation(const FiniteElement<dim> &fe,
7272 const unsigned int first_selected_component)
7273 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7274 fe,
7275 other.mapped_geometry->get_quadrature(),
7276 other.mapped_geometry->get_fe_values().get_update_flags(),
7277 first_selected_component,
7278 &other)
7279 , dofs_per_component(this->data->dofs_per_component_on_cell)
7280 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7281 , n_q_points(this->data->n_q_points)
7282{
7283 check_template_arguments(numbers::invalid_unsigned_int, 0);
7284}
7285
7286
7287
7288template <int dim,
7289 int fe_degree,
7290 int n_q_points_1d,
7291 int n_components_,
7292 typename Number,
7293 typename VectorizedArrayType>
7294inline FEEvaluation<dim,
7295 fe_degree,
7296 n