Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_evaluation.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2011 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
17#ifndef dealii_matrix_free_fe_evaluation_h
18#define dealii_matrix_free_fe_evaluation_h
19
20
21#include <deal.II/base/config.h>
22
29
31
44
45#include <type_traits>
46
47
49
50
51
89template <int dim,
90 int n_components_,
91 typename Number,
92 bool is_face,
93 typename VectorizedArrayType>
95 : public FEEvaluationData<dim, VectorizedArrayType, is_face>
96{
97public:
98 using number_type = Number;
104 static constexpr unsigned int dimension = dim;
105 static constexpr unsigned int n_components = n_components_;
106
143 template <typename VectorType>
144 void
145 read_dof_values(const VectorType & src,
146 const unsigned int first_index = 0,
147 const std::bitset<VectorizedArrayType::size()> &mask =
148 std::bitset<VectorizedArrayType::size()>().flip());
149
178 template <typename VectorType>
179 void
180 read_dof_values_plain(const VectorType & src,
181 const unsigned int first_index = 0,
182 const std::bitset<VectorizedArrayType::size()> &mask =
183 std::bitset<VectorizedArrayType::size()>().flip());
184
216 template <typename VectorType>
217 void
219 VectorType & dst,
220 const unsigned int first_index = 0,
221 const std::bitset<VectorizedArrayType::size()> &mask =
222 std::bitset<VectorizedArrayType::size()>().flip()) const;
223
262 template <typename VectorType>
263 void
264 set_dof_values(VectorType & dst,
265 const unsigned int first_index = 0,
266 const std::bitset<VectorizedArrayType::size()> &mask =
267 std::bitset<VectorizedArrayType::size()>().flip()) const;
268
272 template <typename VectorType>
273 void
275 VectorType & dst,
276 const unsigned int first_index = 0,
277 const std::bitset<VectorizedArrayType::size()> &mask =
278 std::bitset<VectorizedArrayType::size()>().flip()) const;
279
303 get_dof_value(const unsigned int dof) const;
304
315 void
316 submit_dof_value(const value_type val_in, const unsigned int dof);
317
331 get_value(const unsigned int q_point) const;
332
345 void
346 submit_value(const value_type val_in, const unsigned int q_point);
347
359 get_gradient(const unsigned int q_point) const;
360
376 get_normal_derivative(const unsigned int q_point) const;
377
390 void
391 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
392
411 void
413 const unsigned int q_point);
414
427 void
428 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
429
442 get_hessian(const unsigned int q_point) const;
443
454 get_hessian_diagonal(const unsigned int q_point) const;
455
468 get_laplacian(const unsigned int q_point) const;
469
470#ifdef DOXYGEN
471 // doxygen does not anyhow mention functions coming from partial template
472 // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
473 // For now, hack in those functions manually only to fix documentation:
474
481 VectorizedArrayType
482 get_divergence(const unsigned int q_point) const;
483
493 get_symmetric_gradient(const unsigned int q_point) const;
494
501 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
502 get_curl(const unsigned int q_point) const;
503
519 void
520 submit_divergence(const VectorizedArrayType div_in,
521 const unsigned int q_point);
522
539 void
542 const unsigned int q_point);
543
556 void
558 const unsigned int q_point);
559
560#endif
561
580
588
589protected:
600 const unsigned int dof_no,
601 const unsigned int first_selected_component,
602 const unsigned int quad_no,
603 const unsigned int fe_degree,
604 const unsigned int n_q_points,
605 const bool is_interior_face,
606 const unsigned int active_fe_index,
607 const unsigned int active_quad_index,
608 const unsigned int face_type);
609
647 const Mapping<dim> & mapping,
648 const FiniteElement<dim> &fe,
649 const Quadrature<1> & quadrature,
650 const UpdateFlags update_flags,
651 const unsigned int first_selected_component,
653
661
670
675
682 template <typename VectorType, typename VectorOperation>
683 void
685 const VectorOperation & operation,
686 const std::array<VectorType *, n_components_> &vectors,
687 const std::array<
689 n_components_> & vectors_sm,
690 const std::bitset<VectorizedArrayType::size()> &mask,
691 const bool apply_constraints = true) const;
692
700 template <typename VectorType, typename VectorOperation>
701 void
703 const VectorOperation & operation,
704 const std::array<VectorType *, n_components_> &vectors,
705 const std::array<
707 n_components_> & vectors_sm,
708 const std::bitset<VectorizedArrayType::size()> &mask) const;
709
717 template <typename VectorType, typename VectorOperation>
718 void
720 const VectorOperation & operation,
721 const std::array<VectorType *, n_components_> &vectors) const;
722
726 void
728
733
738
743 mutable std::vector<types::global_dof_index> local_dof_indices;
744};
745
746
747
755template <int dim,
756 int n_components_,
757 typename Number,
758 bool is_face,
759 typename VectorizedArrayType = VectorizedArray<Number>>
761 n_components_,
762 Number,
763 is_face,
764 VectorizedArrayType>
765{
766 static_assert(
767 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
768 "Type of Number and of VectorizedArrayType do not match.");
769
770public:
771 using number_type = Number;
775 static constexpr unsigned int dimension = dim;
776 static constexpr unsigned int n_components = n_components_;
777 using BaseClass =
779
780protected:
790 const unsigned int dof_no,
791 const unsigned int first_selected_component,
792 const unsigned int quad_no,
793 const unsigned int fe_degree,
794 const unsigned int n_q_points,
795 const bool is_interior_face = true,
798 const unsigned int face_type = numbers::invalid_unsigned_int);
799
805 const Mapping<dim> & mapping,
806 const FiniteElement<dim> &fe,
807 const Quadrature<1> & quadrature,
808 const UpdateFlags update_flags,
809 const unsigned int first_selected_component,
811
816
822};
823
824
825
834template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
835class FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>
836 : public FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>
837{
838 static_assert(
839 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
840 "Type of Number and of VectorizedArrayType do not match.");
841
842public:
843 using number_type = Number;
844 using value_type = VectorizedArrayType;
847 static constexpr unsigned int dimension = dim;
848 using BaseClass =
850
855 get_dof_value(const unsigned int dof) const;
856
860 void
861 submit_dof_value(const value_type val_in, const unsigned int dof);
862
867 get_value(const unsigned int q_point) const;
868
872 void
873 submit_value(const value_type val_in, const unsigned int q_point);
874
878 void
880 const unsigned int q_point);
881
886 get_gradient(const unsigned int q_point) const;
887
892 get_normal_derivative(const unsigned int q_point) const;
893
897 void
898 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
899
903 void
905 const unsigned int q_point);
906
911 get_hessian(unsigned int q_point) const;
912
917 get_hessian_diagonal(const unsigned int q_point) const;
918
922 void
923 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
924
929 get_laplacian(const unsigned int q_point) const;
930
936
937protected:
947 const unsigned int dof_no,
948 const unsigned int first_selected_component,
949 const unsigned int quad_no,
950 const unsigned int fe_degree,
951 const unsigned int n_q_points,
952 const bool is_interior_face = true,
955 const unsigned int face_type = numbers::invalid_unsigned_int);
956
962 const Mapping<dim> & mapping,
963 const FiniteElement<dim> &fe,
964 const Quadrature<1> & quadrature,
965 const UpdateFlags update_flags,
966 const unsigned int first_selected_component,
968
973
979};
980
981
982
992template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
993class FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>
994 : public FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>
995{
996 static_assert(
997 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
998 "Type of Number and of VectorizedArrayType do not match.");
999
1000public:
1001 using number_type = Number;
1004 static constexpr unsigned int dimension = dim;
1005 static constexpr unsigned int n_components = dim;
1008
1013 get_value(const unsigned int q_point) const;
1014
1019 get_gradient(const unsigned int q_point) const;
1020
1025 VectorizedArrayType
1026 get_divergence(const unsigned int q_point) const;
1027
1035 get_symmetric_gradient(const unsigned int q_point) const;
1036
1041 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
1042 get_curl(const unsigned int q_point) const;
1043
1048 get_hessian(const unsigned int q_point) const;
1049
1054 get_hessian_diagonal(const unsigned int q_point) const;
1055
1059 void
1061 const unsigned int q_point);
1062
1066 void
1067 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1068
1077 void
1079 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
1080 const unsigned int q_point);
1081
1090 void
1091 submit_divergence(const VectorizedArrayType div_in,
1092 const unsigned int q_point);
1093
1102 void
1105 const unsigned int q_point);
1106
1111 void
1113 const unsigned int q_point);
1114
1115protected:
1125 const unsigned int dof_no,
1126 const unsigned int first_selected_component,
1127 const unsigned int quad_no,
1128 const unsigned int dofs_per_cell,
1129 const unsigned int n_q_points,
1130 const bool is_interior_face = true,
1133 const unsigned int face_type = numbers::invalid_unsigned_int);
1134
1140 const Mapping<dim> & mapping,
1141 const FiniteElement<dim> &fe,
1142 const Quadrature<1> & quadrature,
1143 const UpdateFlags update_flags,
1144 const unsigned int first_selected_component,
1146
1151
1157};
1158
1159
1168template <typename Number, bool is_face, typename VectorizedArrayType>
1169class FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>
1170 : public FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>
1171{
1172 static_assert(
1173 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1174 "Type of Number and of VectorizedArrayType do not match.");
1175
1176public:
1177 using number_type = Number;
1178 using value_type = VectorizedArrayType;
1181 static constexpr unsigned int dimension = 1;
1184
1189 get_dof_value(const unsigned int dof) const;
1190
1194 void
1195 submit_dof_value(const value_type val_in, const unsigned int dof);
1196
1201 get_value(const unsigned int q_point) const;
1202
1206 void
1207 submit_value(const value_type val_in, const unsigned int q_point);
1208
1212 void
1213 submit_value(const gradient_type val_in, const unsigned int q_point);
1214
1219 get_gradient(const unsigned int q_point) const;
1220
1225 get_divergence(const unsigned int q_point) const;
1226
1231 get_normal_derivative(const unsigned int q_point) const;
1232
1236 void
1237 submit_gradient(const gradient_type grad_in, const unsigned int q_point);
1238
1242 void
1243 submit_gradient(const value_type grad_in, const unsigned int q_point);
1244
1248 void
1250 const unsigned int q_point);
1251
1255 void
1257 const unsigned int q_point);
1258
1262 void
1264 const unsigned int q_point);
1265
1270 get_hessian(unsigned int q_point) const;
1271
1276 get_hessian_diagonal(const unsigned int q_point) const;
1277
1281 void
1282 submit_hessian(const hessian_type hessian_in, const unsigned int q_point);
1283
1288 get_laplacian(const unsigned int q_point) const;
1289
1295
1296protected:
1306 const unsigned int dof_no,
1307 const unsigned int first_selected_component,
1308 const unsigned int quad_no,
1309 const unsigned int fe_degree,
1310 const unsigned int n_q_points,
1311 const bool is_interior_face = true,
1314 const unsigned int face_type = numbers::invalid_unsigned_int);
1315
1321 const Mapping<1> & mapping,
1322 const FiniteElement<1> &fe,
1323 const Quadrature<1> & quadrature,
1324 const UpdateFlags update_flags,
1325 const unsigned int first_selected_component,
1327
1332
1338};
1339
1340
1341
1905template <int dim,
1906 int fe_degree,
1907 int n_q_points_1d,
1908 int n_components_,
1909 typename Number,
1910 typename VectorizedArrayType>
1912 n_components_,
1913 Number,
1914 false,
1915 VectorizedArrayType>
1916{
1917 static_assert(
1918 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
1919 "Type of Number and of VectorizedArrayType do not match.");
1920
1921public:
1927
1931 using number_type = Number;
1932
1939
1946
1950 static constexpr unsigned int dimension = dim;
1951
1956 static constexpr unsigned int n_components = n_components_;
1957
1966 static constexpr unsigned int static_n_q_points =
1967 Utilities::pow(n_q_points_1d, dim);
1968
1978 static constexpr unsigned int static_dofs_per_component =
1979 Utilities::pow(fe_degree + 1, dim);
1980
1990 static constexpr unsigned int tensor_dofs_per_cell =
1992
2002 static constexpr unsigned int static_dofs_per_cell =
2004
2041 const unsigned int dof_no = 0,
2042 const unsigned int quad_no = 0,
2043 const unsigned int first_selected_component = 0,
2046
2055 const std::pair<unsigned int, unsigned int> & range,
2056 const unsigned int dof_no = 0,
2057 const unsigned int quad_no = 0,
2058 const unsigned int first_selected_component = 0);
2059
2089 const FiniteElement<dim> &fe,
2090 const Quadrature<1> & quadrature,
2091 const UpdateFlags update_flags,
2092 const unsigned int first_selected_component = 0);
2093
2100 const Quadrature<1> & quadrature,
2101 const UpdateFlags update_flags,
2102 const unsigned int first_selected_component = 0);
2103
2116 const unsigned int first_selected_component = 0);
2117
2125
2132 FEEvaluation &
2134
2143 void
2144 reinit(const unsigned int cell_batch_index);
2145
2152 void
2153 reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids);
2154
2167 template <bool level_dof_access>
2168 void
2170
2181 void
2183
2187 static bool
2188 fast_evaluation_supported(const unsigned int given_degree,
2189 const unsigned int given_n_q_points_1d);
2190
2200 void
2202
2208 evaluate(const bool evaluate_values,
2209 const bool evaluate_gradients,
2210 const bool evaluate_hessians = false);
2211
2224 void
2225 evaluate(const VectorizedArrayType * values_array,
2226 const EvaluationFlags::EvaluationFlags evaluation_flag);
2227
2233 evaluate(const VectorizedArrayType *values_array,
2234 const bool evaluate_values,
2235 const bool evaluate_gradients,
2236 const bool evaluate_hessians = false);
2237
2251 template <typename VectorType>
2252 void
2253 gather_evaluate(const VectorType & input_vector,
2254 const EvaluationFlags::EvaluationFlags evaluation_flag);
2255
2259 template <typename VectorType>
2261 gather_evaluate(const VectorType &input_vector,
2262 const bool evaluate_values,
2263 const bool evaluate_gradients,
2264 const bool evaluate_hessians = false);
2265
2276 void
2278
2283 integrate(const bool integrate_values, const bool integrate_gradients);
2284
2296 void
2298 VectorizedArrayType * values_array,
2299 const bool sum_into_values = false);
2300
2305 integrate(const bool integrate_values,
2306 const bool integrate_gradients,
2307 VectorizedArrayType *values_array);
2308
2322 template <typename VectorType>
2323 void
2325 VectorType & output_vector);
2326
2330 template <typename VectorType>
2332 integrate_scatter(const bool integrate_values,
2333 const bool integrate_gradients,
2334 VectorType &output_vector);
2335
2343
2350 const unsigned int dofs_per_component;
2351
2358 const unsigned int dofs_per_cell;
2359
2367 const unsigned int n_q_points;
2368
2369private:
2374 void
2375 check_template_arguments(const unsigned int fe_no,
2376 const unsigned int first_selected_component);
2377};
2378
2379
2380
2416template <int dim,
2417 int fe_degree,
2418 int n_q_points_1d = fe_degree + 1,
2419 int n_components_ = 1,
2420 typename Number = double,
2421 typename VectorizedArrayType = VectorizedArray<Number>>
2423 n_components_,
2424 Number,
2425 true,
2426 VectorizedArrayType>
2427{
2428 static_assert(
2429 std::is_same<Number, typename VectorizedArrayType::value_type>::value,
2430 "Type of Number and of VectorizedArrayType do not match.");
2431
2432public:
2438
2442 using number_type = Number;
2443
2450
2457
2461 static constexpr unsigned int dimension = dim;
2462
2467 static constexpr unsigned int n_components = n_components_;
2468
2478 static constexpr unsigned int static_n_q_points =
2479 Utilities::pow(n_q_points_1d, dim - 1);
2480
2489 static constexpr unsigned int static_n_q_points_cell =
2490 Utilities::pow(n_q_points_1d, dim);
2491
2500 static constexpr unsigned int static_dofs_per_component =
2501 Utilities::pow(fe_degree + 1, dim);
2502
2511 static constexpr unsigned int tensor_dofs_per_cell =
2513
2522 static constexpr unsigned int static_dofs_per_cell =
2524
2568 const bool is_interior_face = true,
2569 const unsigned int dof_no = 0,
2570 const unsigned int quad_no = 0,
2571 const unsigned int first_selected_component = 0,
2574 const unsigned int face_type = numbers::invalid_unsigned_int);
2575
2585 const std::pair<unsigned int, unsigned int> & range,
2586 const bool is_interior_face = true,
2587 const unsigned int dof_no = 0,
2588 const unsigned int quad_no = 0,
2589 const unsigned int first_selected_component = 0);
2590
2601 void
2602 reinit(const unsigned int face_batch_number);
2603
2611 void
2612 reinit(const unsigned int cell_batch_number, const unsigned int face_number);
2613
2617 static bool
2618 fast_evaluation_supported(const unsigned int given_degree,
2619 const unsigned int given_n_q_points_1d);
2620
2631 void
2633
2638 evaluate(const bool evaluate_values, const bool evaluate_gradients);
2639
2652 void
2653 evaluate(const VectorizedArrayType * values_array,
2654 const EvaluationFlags::EvaluationFlags evaluation_flag);
2655
2660 evaluate(const VectorizedArrayType *values_array,
2661 const bool evaluate_values,
2662 const bool evaluate_gradients);
2663
2675 template <typename VectorType>
2676 void
2677 gather_evaluate(const VectorType & input_vector,
2678 const EvaluationFlags::EvaluationFlags evaluation_flag);
2679
2683 template <typename VectorType>
2685 gather_evaluate(const VectorType &input_vector,
2686 const bool evaluate_values,
2687 const bool evaluate_gradients);
2688
2698 void
2700
2705 integrate(const bool integrate_values, const bool integrate_gradients);
2706
2715 void
2717 VectorizedArrayType * values_array);
2718
2723 integrate(const bool integrate_values,
2724 const bool integrate_gradients,
2725 VectorizedArrayType *values_array);
2726
2738 template <typename VectorType>
2739 void
2741 VectorType & output_vector);
2742
2746 template <typename VectorType>
2747 void
2748 integrate_scatter(const bool integrate_values,
2749 const bool integrate_gradients,
2750 VectorType &output_vector);
2751
2759
2764 bool
2766
2781
2788 const unsigned int dofs_per_component;
2789
2796 const unsigned int dofs_per_cell;
2797
2805 const unsigned int n_q_points;
2806};
2807
2808
2809
2810namespace internal
2811{
2812 namespace MatrixFreeFunctions
2813 {
2814 // a helper function to compute the number of DoFs of a DGP element at
2815 // compile time, depending on the degree
2816 template <int dim, int degree>
2818 {
2819 // this division is always without remainder
2820 static constexpr unsigned int value =
2821 (DGP_dofs_per_component<dim - 1, degree>::value * (degree + dim)) / dim;
2822 };
2823
2824 // base specialization: 1d elements have 'degree+1' degrees of freedom
2825 template <int degree>
2826 struct DGP_dofs_per_component<1, degree>
2827 {
2828 static constexpr unsigned int value = degree + 1;
2829 };
2830 } // namespace MatrixFreeFunctions
2831} // namespace internal
2832
2833
2834/*----------------------- Inline functions ----------------------------------*/
2835
2836#ifndef DOXYGEN
2837
2838
2839namespace internal
2840{
2841 // Extract all internal data pointers and indices in a single function that
2842 // get passed on to the constructor of FEEvaluationData, avoiding to look
2843 // things up multiple times
2844 template <bool is_face,
2845 int dim,
2846 typename Number,
2847 typename VectorizedArrayType>
2849 InitializationData
2850 extract_initialization_data(
2852 const unsigned int dof_no,
2853 const unsigned int first_selected_component,
2854 const unsigned int quad_no,
2855 const unsigned int fe_degree,
2856 const unsigned int n_q_points,
2857 const unsigned int active_fe_index_given,
2858 const unsigned int active_quad_index_given,
2859 const unsigned int face_type)
2860 {
2862 InitializationData init_data;
2863
2864 init_data.dof_info = &matrix_free.get_dof_info(dof_no);
2865 init_data.mapping_data =
2866 &internal::MatrixFreeFunctions::
2867 MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
2868 matrix_free.get_mapping_info(), quad_no);
2869
2870 init_data.active_fe_index =
2871 fe_degree != numbers::invalid_unsigned_int ?
2872 init_data.dof_info->fe_index_from_degree(first_selected_component,
2873 fe_degree) :
2874 (active_fe_index_given != numbers::invalid_unsigned_int ?
2875 active_fe_index_given :
2876 0);
2877 init_data.active_quad_index =
2878 fe_degree == numbers::invalid_unsigned_int ?
2879 (active_quad_index_given != numbers::invalid_unsigned_int ?
2880 active_quad_index_given :
2881 std::min<unsigned int>(init_data.active_fe_index,
2882 init_data.mapping_data->descriptor.size() -
2883 1)) :
2884 init_data.mapping_data->quad_index_from_n_q_points(n_q_points);
2885
2886 init_data.shape_info = &matrix_free.get_shape_info(
2887 dof_no,
2888 quad_no,
2889 init_data.dof_info->component_to_base_index[first_selected_component],
2890 init_data.active_fe_index,
2891 init_data.active_quad_index);
2892 init_data.descriptor =
2893 &init_data.mapping_data->descriptor
2894 [is_face ?
2895 (init_data.active_quad_index * std::max<unsigned int>(1, dim - 1) +
2896 (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
2897 init_data.active_quad_index];
2898
2899 return init_data;
2900 }
2901} // namespace internal
2902
2903
2904
2905/*----------------------- FEEvaluationBase ----------------------------------*/
2906
2907template <int dim,
2908 int n_components_,
2909 typename Number,
2910 bool is_face,
2911 typename VectorizedArrayType>
2912inline FEEvaluationBase<dim,
2913 n_components_,
2914 Number,
2915 is_face,
2916 VectorizedArrayType>::
2917 FEEvaluationBase(
2919 const unsigned int dof_no,
2920 const unsigned int first_selected_component,
2921 const unsigned int quad_no,
2922 const unsigned int fe_degree,
2923 const unsigned int n_q_points,
2924 const bool is_interior_face,
2925 const unsigned int active_fe_index,
2926 const unsigned int active_quad_index,
2927 const unsigned int face_type)
2928 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2929 internal::extract_initialization_data<is_face>(matrix_free,
2930 dof_no,
2931 first_selected_component,
2932 quad_no,
2933 fe_degree,
2934 n_q_points,
2935 active_fe_index,
2936 active_quad_index,
2937 face_type),
2938 is_interior_face,
2939 quad_no,
2940 first_selected_component)
2941 , scratch_data_array(matrix_free.acquire_scratch_data())
2942 , matrix_free(&matrix_free)
2943{
2944 this->set_data_pointers(scratch_data_array, n_components_);
2945 Assert(
2946 this->dof_info->start_components.back() == 1 ||
2947 static_cast<int>(n_components_) <=
2948 static_cast<int>(
2949 this->dof_info->start_components
2950 [this->dof_info->component_to_base_index[first_selected_component] +
2951 1]) -
2952 first_selected_component,
2953 ExcMessage(
2954 "You tried to construct a vector-valued evaluator with " +
2955 std::to_string(n_components) +
2956 " components. However, "
2957 "the current base element has only " +
2958 std::to_string(
2959 this->dof_info->start_components
2960 [this->dof_info->component_to_base_index[first_selected_component] +
2961 1] -
2962 first_selected_component) +
2963 " components left when starting from local element index " +
2964 std::to_string(
2965 first_selected_component -
2966 this->dof_info->start_components
2967 [this->dof_info->component_to_base_index[first_selected_component]]) +
2968 " (global index " + std::to_string(first_selected_component) + ")"));
2969
2970 // do not check for correct dimensions of data fields here, should be done
2971 // in derived classes
2972}
2973
2974
2975
2976template <int dim,
2977 int n_components_,
2978 typename Number,
2979 bool is_face,
2980 typename VectorizedArrayType>
2981inline FEEvaluationBase<dim,
2982 n_components_,
2983 Number,
2984 is_face,
2985 VectorizedArrayType>::
2986 FEEvaluationBase(
2987 const Mapping<dim> & mapping,
2988 const FiniteElement<dim> &fe,
2989 const Quadrature<1> & quadrature,
2990 const UpdateFlags update_flags,
2991 const unsigned int first_selected_component,
2993 : FEEvaluationData<dim, VectorizedArrayType, is_face>(
2994 other != nullptr &&
2995 other->mapped_geometry->get_quadrature() == quadrature ?
2996 other->mapped_geometry :
2997 std::make_shared<internal::MatrixFreeFunctions::
2998 MappingDataOnTheFly<dim, VectorizedArrayType>>(
2999 mapping,
3000 quadrature,
3001 update_flags),
3002 n_components_,
3003 first_selected_component)
3004 , scratch_data_array(new AlignedVector<VectorizedArrayType>())
3005 , matrix_free(nullptr)
3006{
3007 const unsigned int base_element_number =
3008 fe.component_to_base_index(first_selected_component).first;
3009 Assert(fe.element_multiplicity(base_element_number) == 1 ||
3010 fe.element_multiplicity(base_element_number) -
3011 first_selected_component >=
3012 n_components_,
3013 ExcMessage("The underlying element must at least contain as many "
3014 "components as requested by this class"));
3015 (void)base_element_number;
3016
3017 Assert(this->data == nullptr, ExcInternalError());
3018 this->data =
3020 Quadrature<(is_face ? dim - 1 : dim)>(quadrature),
3021 fe,
3022 fe.component_to_base_index(first_selected_component).first);
3023
3024 this->set_data_pointers(scratch_data_array, n_components_);
3025}
3026
3027
3028
3029template <int dim,
3030 int n_components_,
3031 typename Number,
3032 bool is_face,
3033 typename VectorizedArrayType>
3034inline FEEvaluationBase<dim,
3035 n_components_,
3036 Number,
3037 is_face,
3038 VectorizedArrayType>::
3039 FEEvaluationBase(const FEEvaluationBase<dim,
3040 n_components_,
3041 Number,
3042 is_face,
3043 VectorizedArrayType> &other)
3044 : FEEvaluationData<dim, VectorizedArrayType, is_face>(other)
3045 , scratch_data_array(other.matrix_free == nullptr ?
3046 new AlignedVector<VectorizedArrayType>() :
3047 other.matrix_free->acquire_scratch_data())
3048 , matrix_free(other.matrix_free)
3049{
3050 if (other.matrix_free == nullptr)
3051 {
3052 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3053 this->data =
3055 *other.data);
3056
3057 // Create deep copy of mapped geometry for use in parallel
3058 this->mapped_geometry =
3059 std::make_shared<internal::MatrixFreeFunctions::
3060 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3061 other.mapped_geometry->get_fe_values().get_mapping(),
3062 other.mapped_geometry->get_quadrature(),
3063 other.mapped_geometry->get_fe_values().get_update_flags());
3064 this->mapping_data = &this->mapped_geometry->get_data_storage();
3065 this->cell = 0;
3066
3067 this->jacobian =
3068 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3069 this->J_value =
3070 this->mapped_geometry->get_data_storage().JxW_values.begin();
3071 this->jacobian_gradients =
3072 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3073 this->jacobian_gradients_non_inverse =
3074 this->mapped_geometry->get_data_storage()
3075 .jacobian_gradients_non_inverse[0]
3076 .begin();
3077 this->quadrature_points =
3078 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3079 }
3080
3081 this->set_data_pointers(scratch_data_array, n_components_);
3082}
3083
3084
3085
3086template <int dim,
3087 int n_components_,
3088 typename Number,
3089 bool is_face,
3090 typename VectorizedArrayType>
3091inline FEEvaluationBase<dim,
3092 n_components_,
3093 Number,
3094 is_face,
3095 VectorizedArrayType> &
3097operator=(const FEEvaluationBase<dim,
3098 n_components_,
3099 Number,
3100 is_face,
3101 VectorizedArrayType> &other)
3102{
3103 // release old memory
3104 if (matrix_free == nullptr)
3105 {
3106 delete this->data;
3107 delete scratch_data_array;
3108 }
3109 else
3110 {
3111 matrix_free->release_scratch_data(scratch_data_array);
3112 }
3113
3115
3116 matrix_free = other.matrix_free;
3117
3118 if (other.matrix_free == nullptr)
3119 {
3120 Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
3121 this->data =
3123 *other.data);
3124 scratch_data_array = new AlignedVector<VectorizedArrayType>();
3125
3126 // Create deep copy of mapped geometry for use in parallel
3127 this->mapped_geometry =
3128 std::make_shared<internal::MatrixFreeFunctions::
3129 MappingDataOnTheFly<dim, VectorizedArrayType>>(
3130 other.mapped_geometry->get_fe_values().get_mapping(),
3131 other.mapped_geometry->get_quadrature(),
3132 other.mapped_geometry->get_fe_values().get_update_flags());
3133 this->cell = 0;
3134 this->mapping_data = &this->mapped_geometry->get_data_storage();
3135 this->jacobian =
3136 this->mapped_geometry->get_data_storage().jacobians[0].begin();
3137 this->J_value =
3138 this->mapped_geometry->get_data_storage().JxW_values.begin();
3139 this->jacobian_gradients =
3140 this->mapped_geometry->get_data_storage().jacobian_gradients[0].begin();
3141 this->jacobian_gradients_non_inverse =
3142 this->mapped_geometry->get_data_storage()
3143 .jacobian_gradients_non_inverse[0]
3144 .begin();
3145 this->quadrature_points =
3146 this->mapped_geometry->get_data_storage().quadrature_points.begin();
3147 }
3148 else
3149 {
3150 scratch_data_array = matrix_free->acquire_scratch_data();
3151 }
3152
3153 this->set_data_pointers(scratch_data_array, n_components_);
3154
3155 return *this;
3156}
3157
3158
3159
3160template <int dim,
3161 int n_components_,
3162 typename Number,
3163 bool is_face,
3164 typename VectorizedArrayType>
3165inline FEEvaluationBase<dim,
3166 n_components_,
3167 Number,
3168 is_face,
3169 VectorizedArrayType>::~FEEvaluationBase()
3170{
3171 if (matrix_free != nullptr)
3172 {
3173 try
3174 {
3175 matrix_free->release_scratch_data(scratch_data_array);
3176 }
3177 catch (...)
3178 {}
3179 }
3180 else
3181 {
3182 delete scratch_data_array;
3183 delete this->data;
3184 }
3185}
3186
3187
3188
3189template <int dim,
3190 int n_components_,
3191 typename Number,
3192 bool is_face,
3193 typename VectorizedArrayType>
3196 get_matrix_free() const
3197{
3198 Assert(matrix_free != nullptr,
3199 ExcMessage(
3200 "FEEvaluation was not initialized with a MatrixFree object!"));
3201 return *matrix_free;
3202}
3203
3204
3205
3206namespace internal
3207{
3208 // given a block vector return the underlying vector type
3209 // including constness (specified by bool)
3210 template <typename VectorType, bool>
3211 struct ConstBlockVectorSelector;
3212
3213 template <typename VectorType>
3214 struct ConstBlockVectorSelector<VectorType, true>
3215 {
3216 using BaseVectorType = const typename VectorType::BlockType;
3217 };
3218
3219 template <typename VectorType>
3220 struct ConstBlockVectorSelector<VectorType, false>
3221 {
3222 using BaseVectorType = typename VectorType::BlockType;
3223 };
3224
3225 // allows to select between block vectors and non-block vectors, which
3226 // allows to use a unified interface for extracting blocks on block vectors
3227 // and doing nothing on usual vectors
3228 template <typename VectorType, bool>
3229 struct BlockVectorSelector;
3230
3231 template <typename VectorType>
3232 struct BlockVectorSelector<VectorType, true>
3233 {
3234 using BaseVectorType = typename ConstBlockVectorSelector<
3235 VectorType,
3236 std::is_const<VectorType>::value>::BaseVectorType;
3237
3238 static BaseVectorType *
3239 get_vector_component(VectorType &vec, const unsigned int component)
3240 {
3241 AssertIndexRange(component, vec.n_blocks());
3242 return &vec.block(component);
3243 }
3244 };
3245
3246 template <typename VectorType>
3247 struct BlockVectorSelector<VectorType, false>
3248 {
3249 using BaseVectorType = VectorType;
3250
3251 static BaseVectorType *
3252 get_vector_component(VectorType &vec, const unsigned int component)
3253 {
3254 // FEEvaluation allows to combine several vectors from a scalar
3255 // FiniteElement into a "vector-valued" FEEvaluation object with
3256 // multiple components. These components can be extracted with the other
3257 // get_vector_component functions. If we do not get a vector of vectors
3258 // (std::vector<VectorType>, std::vector<VectorType*>, BlockVector), we
3259 // must make sure that we do not duplicate the components in input
3260 // and/or duplicate the resulting integrals. In such a case, we should
3261 // only get the zeroth component in the vector contained set nullptr for
3262 // the others which allows us to catch unintended use in
3263 // read_write_operation.
3264 if (component == 0)
3265 return &vec;
3266 else
3267 return nullptr;
3268 }
3269 };
3270
3271 template <typename VectorType>
3272 struct BlockVectorSelector<std::vector<VectorType>, false>
3273 {
3274 using BaseVectorType = VectorType;
3275
3276 static BaseVectorType *
3277 get_vector_component(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<const std::vector<VectorType>, false>
3287 {
3288 using BaseVectorType = const VectorType;
3289
3290 static const BaseVectorType *
3291 get_vector_component(const 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<std::vector<VectorType *>, false>
3301 {
3302 using BaseVectorType = VectorType;
3303
3304 static BaseVectorType *
3305 get_vector_component(std::vector<VectorType *> &vec,
3306 const unsigned int component)
3307 {
3308 AssertIndexRange(component, vec.size());
3309 return vec[component];
3310 }
3311 };
3312
3313 template <typename VectorType>
3314 struct BlockVectorSelector<const std::vector<VectorType *>, false>
3315 {
3316 using BaseVectorType = const VectorType;
3317
3318 static const BaseVectorType *
3319 get_vector_component(const std::vector<VectorType *> &vec,
3320 const unsigned int component)
3321 {
3322 AssertIndexRange(component, vec.size());
3323 return vec[component];
3324 }
3325 };
3326} // namespace internal
3327
3328
3329
3330template <int dim,
3331 int n_components_,
3332 typename Number,
3333 bool is_face,
3334 typename VectorizedArrayType>
3335template <typename VectorType, typename VectorOperation>
3336inline void
3339 const VectorOperation & operation,
3340 const std::array<VectorType *, n_components_> &src,
3341 const std::array<
3343 n_components_> & src_sm,
3344 const std::bitset<VectorizedArrayType::size()> &mask,
3345 const bool apply_constraints) const
3346{
3347 // Case 1: No MatrixFree object given, simple case because we do not need to
3348 // process constraints and need not care about vectorization -> go to
3349 // separate function
3350 if (this->matrix_free == nullptr)
3351 {
3352 read_write_operation_global(operation, src);
3353 return;
3354 }
3355
3356 Assert(this->dof_info != nullptr, ExcNotInitialized());
3357 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3358 Assert(this->matrix_free->indices_initialized() == true, ExcNotInitialized());
3359 if (this->n_fe_components == 1)
3360 for (unsigned int comp = 0; comp < n_components; ++comp)
3361 {
3362 Assert(src[comp] != nullptr,
3363 ExcMessage("The finite element underlying this FEEvaluation "
3364 "object is scalar, but you requested " +
3365 std::to_string(n_components) +
3366 " components via the template argument in "
3367 "FEEvaluation. In that case, you must pass an "
3368 "std::vector<VectorType> or a BlockVector to " +
3369 "read_dof_values and distribute_local_to_global."));
3371 *this->matrix_free,
3372 *this->dof_info);
3373 }
3374 else
3375 {
3377 *this->matrix_free,
3378 *this->dof_info);
3379 }
3380
3381 // Case 2: contiguous indices which use reduced storage of indices and can
3382 // use vectorized load/store operations -> go to separate function
3383 if (this->cell != numbers::invalid_unsigned_int)
3384 {
3386 this->cell,
3387 dof_info.index_storage_variants[this->dof_access_index].size());
3388
3389 bool is_contiguous = true;
3390 // check if exterior cells are not contiguous (ECL case)
3391 if (is_face && !this->interior_face &&
3392 (this->dof_access_index ==
3394 {
3395 const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3396 this->get_cell_ids();
3397 const unsigned int n_filled_lanes =
3400 [this->cell];
3401 // we have to check all filled lanes which are active in the mask
3402 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3403 if (mask[v] == true &&
3404 dof_info.index_storage_variants
3406 [cells[v] / VectorizedArrayType::size()] <
3408 contiguous)
3409 is_contiguous = false;
3410 } // or if cell/face batch is not contiguous
3411 else if (dof_info.index_storage_variants
3412 [is_face ?
3413 this->dof_access_index :
3414 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
3415 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
3416 IndexStorageVariants::contiguous)
3417 {
3418 is_contiguous = false;
3419 }
3420
3421 if (is_contiguous)
3422 {
3423 read_write_operation_contiguous(operation, src, src_sm, mask);
3424 return;
3425 }
3426 }
3427
3428 // Case 3: standard operation with one index per degree of freedom -> go on
3429 // here
3430 constexpr unsigned int n_lanes = VectorizedArrayType::size();
3431
3432 std::array<unsigned int, VectorizedArrayType::size()> cells =
3433 this->get_cell_ids();
3434
3435 const bool masking_is_active = mask.count() < n_lanes;
3436 if (masking_is_active)
3437 for (unsigned int v = 0; v < n_lanes; ++v)
3438 if (mask[v] == false)
3440
3441 bool has_hn_constraints = false;
3442
3443 if (is_face == false)
3444 {
3445 if (!dof_info.hanging_node_constraint_masks.empty() &&
3446 !dof_info.hanging_node_constraint_masks_comp.empty() &&
3447 dof_info
3448 .hanging_node_constraint_masks_comp[this->active_fe_index]
3449 [this->first_selected_component])
3450 for (unsigned int v = 0; v < n_lanes; ++v)
3451 if (cells[v] != numbers::invalid_unsigned_int &&
3452 dof_info.hanging_node_constraint_masks[cells[v]] !=
3455 has_hn_constraints = true;
3456 }
3457
3458 std::integral_constant<bool,
3459 internal::is_vectorizable<VectorType, Number>::value>
3460 vector_selector;
3461
3462 const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
3463
3464 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3465 std::array<VectorizedArrayType *, n_components> values_dofs;
3466 for (unsigned int c = 0; c < n_components; ++c)
3467 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3468 c * dofs_per_component;
3469
3470 if (this->cell != numbers::invalid_unsigned_int &&
3471 dof_info.index_storage_variants
3472 [is_face ? this->dof_access_index :
3473 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
3474 [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
3475 IndexStorageVariants::interleaved &&
3476 use_vectorized_path)
3477 {
3478 const unsigned int *dof_indices =
3479 dof_info.dof_indices_interleaved.data() +
3480 dof_info.row_starts[this->cell * this->n_fe_components * n_lanes]
3481 .first +
3482 this->dof_info
3483 ->component_dof_indices_offset[this->active_fe_index]
3484 [this->first_selected_component] *
3485 n_lanes;
3486 if (n_components == 1 || this->n_fe_components == 1)
3487 for (unsigned int i = 0; i < dofs_per_component;
3488 ++i, dof_indices += n_lanes)
3489 for (unsigned int comp = 0; comp < n_components; ++comp)
3490 operation.process_dof_gather(dof_indices,
3491 *src[comp],
3492 0,
3493 values_dofs[comp][i],
3494 vector_selector);
3495 else
3496 for (unsigned int comp = 0; comp < n_components; ++comp)
3497 for (unsigned int i = 0; i < dofs_per_component;
3498 ++i, dof_indices += n_lanes)
3499 operation.process_dof_gather(
3500 dof_indices, *src[0], 0, values_dofs[comp][i], 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<VectorizedArrayType::size()> &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::integral_constant<bool,
3837 internal::is_vectorizable<VectorType, Number>::value>
3838 vector_selector;
3840 is_face ? this->dof_access_index :
3842 const unsigned int n_lanes = mask.count();
3843
3844 const internal::MatrixFreeFunctions::DoFInfo &dof_info = *this->dof_info;
3845 const std::vector<unsigned int> & dof_indices_cont =
3846 dof_info.dof_indices_contiguous[ind];
3847
3848 const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
3849 std::array<VectorizedArrayType *, n_components> values_dofs;
3850 for (unsigned int c = 0; c < n_components; ++c)
3851 values_dofs[c] = const_cast<VectorizedArrayType *>(this->values_dofs) +
3852 c * dofs_per_component;
3853
3855
3856 // Simple case: We have contiguous storage, so we can simply copy out the
3857 // data
3858 if ((dof_info.index_storage_variants[ind][this->cell] ==
3860 interleaved_contiguous &&
3861 n_lanes == VectorizedArrayType::size()) &&
3862 !(is_face &&
3863 this->dof_access_index ==
3864 internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
3865 this->is_interior_face() == false) &&
3866 !(!is_face && !this->is_interior_face()))
3867 {
3868 const unsigned int dof_index =
3869 dof_indices_cont[this->cell * VectorizedArrayType::size()] +
3870 this->dof_info
3871 ->component_dof_indices_offset[this->active_fe_index]
3872 [this->first_selected_component] *
3873 VectorizedArrayType::size();
3874 if (n_components == 1 || this->n_fe_components == 1)
3875 for (unsigned int comp = 0; comp < n_components; ++comp)
3876 operation.process_dofs_vectorized(dofs_per_component,
3877 dof_index,
3878 *src[comp],
3879 values_dofs[comp],
3880 vector_selector);
3881 else
3882 operation.process_dofs_vectorized(dofs_per_component * n_components,
3883 dof_index,
3884 *src[0],
3885 values_dofs[0],
3886 vector_selector);
3887 return;
3888 }
3889
3890 const std::array<unsigned int, VectorizedArrayType::size()> &cells =
3891 this->get_cell_or_face_ids();
3892
3893 // More general case: Must go through the components one by one and apply
3894 // some transformations
3895 const unsigned int n_filled_lanes =
3896 dof_info.n_vectorization_lanes_filled[ind][this->cell];
3897
3898 const bool is_ecl =
3899 (this->dof_access_index ==
3901 this->is_interior_face() == false) ||
3902 (!is_face && !this->is_interior_face());
3903
3904 if (vectors_sm[0] != nullptr)
3905 {
3906 const auto compute_vector_ptrs = [&](const unsigned int comp) {
3907 std::array<typename VectorType::value_type *,
3908 VectorizedArrayType::size()>
3909 vector_ptrs = {};
3910
3911 const auto upper_bound =
3912 std::min<unsigned int>(n_filled_lanes, VectorizedArrayType::size());
3913 for (unsigned int v = 0; v < upper_bound; ++v)
3914 {
3915 if (mask[v] == false)
3916 {
3917 vector_ptrs[v] = nullptr;
3918 continue;
3919 }
3920
3923 Assert(ind < dof_info.dof_indices_contiguous_sm.size(),
3924 ExcIndexRange(ind,
3925 0,
3926 dof_info.dof_indices_contiguous_sm.size()));
3927 Assert(
3928 cells[v] < dof_info.dof_indices_contiguous_sm[ind].size(),
3929 ExcIndexRange(cells[v],
3930 0,
3931 dof_info.dof_indices_contiguous_sm[ind].size()));
3932
3933 const auto &temp =
3934 dof_info.dof_indices_contiguous_sm[ind][cells[v]];
3935
3936 if (temp.first != numbers::invalid_unsigned_int)
3937 vector_ptrs[v] = const_cast<typename VectorType::value_type *>(
3938 vectors_sm[comp]->operator[](temp.first).data() + temp.second +
3940 [this->active_fe_index][this->first_selected_component]);
3941 else
3942 vector_ptrs[v] = nullptr;
3943 }
3944 for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size();
3945 ++v)
3946 vector_ptrs[v] = nullptr;
3947
3948 return vector_ptrs;
3949 };
3950
3951 if (n_filled_lanes == VectorizedArrayType::size() &&
3952 n_lanes == VectorizedArrayType::size() && !is_ecl)
3953 {
3954 if (n_components == 1 || this->n_fe_components == 1)
3955 {
3956 for (unsigned int comp = 0; comp < n_components; ++comp)
3957 {
3958 auto vector_ptrs = compute_vector_ptrs(comp);
3959 operation.process_dofs_vectorized_transpose(
3960 dofs_per_component,
3961 vector_ptrs,
3962 values_dofs[comp],
3963 vector_selector);
3964 }
3965 }
3966 else
3967 {
3968 auto vector_ptrs = compute_vector_ptrs(0);
3969 operation.process_dofs_vectorized_transpose(dofs_per_component *
3970 n_components,
3971 vector_ptrs,
3972 &values_dofs[0][0],
3973 vector_selector);
3974 }
3975 }
3976 else
3977 for (unsigned int comp = 0; comp < n_components; ++comp)
3978 {
3979 auto vector_ptrs = compute_vector_ptrs(
3980 (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
3981
3982 for (unsigned int i = 0; i < dofs_per_component; ++i)
3983 operation.process_empty(values_dofs[comp][i]);
3984
3985 if (n_components == 1 || this->n_fe_components == 1)
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(vector_ptrs[v][i],
3991 values_dofs[comp][i][v]);
3992 }
3993 else
3994 {
3995 for (unsigned int v = 0; v < n_filled_lanes; ++v)
3996 if (mask[v] == true)
3997 for (unsigned int i = 0; i < dofs_per_component; ++i)
3998 operation.process_dof(
3999 vector_ptrs[v][i + comp * dofs_per_component],
4000 values_dofs[comp][i][v]);
4001 }
4002 }
4003 return;
4004 }
4005
4006 unsigned int dof_indices[VectorizedArrayType::size()];
4007
4008 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4009 {
4010 Assert(mask[v] == false || cells[v] != numbers::invalid_unsigned_int,
4012 if (mask[v] == true)
4013 dof_indices[v] =
4014 dof_indices_cont[cells[v]] +
4015 this->dof_info
4016 ->component_dof_indices_offset[this->active_fe_index]
4017 [this->first_selected_component] *
4018 dof_info.dof_indices_interleave_strides[ind][cells[v]];
4019 else
4020 dof_indices[v] = numbers::invalid_unsigned_int;
4021 }
4022
4023 for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
4024 dof_indices[v] = numbers::invalid_unsigned_int;
4025
4026 // In the case with contiguous cell indices, we know that there are no
4027 // constraints and that the indices within each element are contiguous
4028 if (n_filled_lanes == VectorizedArrayType::size() &&
4029 n_lanes == VectorizedArrayType::size() && !is_ecl)
4030 {
4031 if (dof_info.index_storage_variants[ind][this->cell] ==
4033 contiguous)
4034 {
4035 if (n_components == 1 || this->n_fe_components == 1)
4036 for (unsigned int comp = 0; comp < n_components; ++comp)
4037 operation.process_dofs_vectorized_transpose(dofs_per_component,
4038 dof_indices,
4039 *src[comp],
4040 values_dofs[comp],
4041 vector_selector);
4042 else
4043 operation.process_dofs_vectorized_transpose(dofs_per_component *
4044 n_components,
4045 dof_indices,
4046 *src[0],
4047 &values_dofs[0][0],
4048 vector_selector);
4049 }
4050 else if (dof_info.index_storage_variants[ind][this->cell] ==
4052 interleaved_contiguous_strided)
4053 {
4054 if (n_components == 1 || this->n_fe_components == 1)
4055 for (unsigned int i = 0; i < dofs_per_component; ++i)
4056 {
4057 for (unsigned int comp = 0; comp < n_components; ++comp)
4058 operation.process_dof_gather(dof_indices,
4059 *src[comp],
4060 i * VectorizedArrayType::size(),
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(dof_indices,
4069 *src[0],
4070 (comp * dofs_per_component + i) *
4071 VectorizedArrayType::size(),
4072 values_dofs[comp][i],
4073 vector_selector);
4074 }
4075 }
4076 else
4077 {
4078 Assert(dof_info.index_storage_variants[ind][this->cell] ==
4080 IndexStorageVariants::interleaved_contiguous_mixed_strides,
4082 const unsigned int *offsets =
4084 [ind][VectorizedArrayType::size() * this->cell];
4085 if (n_components == 1 || this->n_fe_components == 1)
4086 for (unsigned int i = 0; i < dofs_per_component; ++i)
4087 {
4088 for (unsigned int comp = 0; comp < n_components; ++comp)
4089 operation.process_dof_gather(dof_indices,
4090 *src[comp],
4091 0,
4092 values_dofs[comp][i],
4093 vector_selector);
4095 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4096 dof_indices[v] += offsets[v];
4097 }
4098 else
4099 for (unsigned int comp = 0; comp < n_components; ++comp)
4100 for (unsigned int i = 0; i < dofs_per_component; ++i)
4101 {
4102 operation.process_dof_gather(dof_indices,
4103 *src[0],
4104 0,
4105 values_dofs[comp][i],
4106 vector_selector);
4108 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
4109 dof_indices[v] += offsets[v];
4110 }
4111 }
4112 }
4113 else
4114 for (unsigned int comp = 0; comp < n_components; ++comp)
4115 {
4116 for (unsigned int i = 0; i < dofs_per_component; ++i)
4117 operation.process_empty(values_dofs[comp][i]);
4118 if (dof_info.index_storage_variants[ind][this->cell] ==
4120 contiguous)
4121 {
4122 if (n_components == 1 || this->n_fe_components == 1)
4123 {
4124 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4125 if (mask[v] == true)
4126 for (unsigned int i = 0; i < dofs_per_component; ++i)
4127 operation.process_dof(dof_indices[v] + i,
4128 *src[comp],
4129 values_dofs[comp][i][v]);
4130 }
4131 else
4132 {
4133 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4134 if (mask[v] == true)
4135 for (unsigned int i = 0; i < dofs_per_component; ++i)
4136 operation.process_dof(dof_indices[v] + i +
4137 comp * dofs_per_component,
4138 *src[0],
4139 values_dofs[comp][i][v]);
4140 }
4141 }
4142 else
4143 {
4144 const unsigned int *offsets =
4146 [ind][VectorizedArrayType::size() * this->cell];
4147 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4148 AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
4149 if (n_components == 1 || this->n_fe_components == 1)
4150 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4151 {
4152 if (mask[v] == true)
4153 for (unsigned int i = 0; i < dofs_per_component; ++i)
4154 operation.process_dof(dof_indices[v] + i * offsets[v],
4155 *src[comp],
4156 values_dofs[comp][i][v]);
4157 }
4158 else
4159 {
4160 for (unsigned int v = 0; v < n_filled_lanes; ++v)
4161 if (mask[v] == true)
4162 for (unsigned int i = 0; i < dofs_per_component; ++i)
4163 operation.process_dof(dof_indices[v] +
4164 (i + comp * dofs_per_component) *
4165 offsets[v],
4166 *src[0],
4167 values_dofs[comp][i][v]);
4168 }
4169 }
4170 }
4171}
4172
4173namespace internal
4174{
4175 template <
4176 typename Number,
4177 typename VectorType,
4178 std::enable_if_t<!IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4179 decltype(std::declval<VectorType>().begin())
4180 get_beginning(VectorType &vec)
4181 {
4182 return vec.begin();
4183 }
4184
4185 template <
4186 typename Number,
4187 typename VectorType,
4188 std::enable_if_t<IsBlockVector<VectorType>::value, VectorType> * = nullptr>
4189 typename VectorType::value_type *
4190 get_beginning(VectorType &)
4191 {
4192 return nullptr;
4193 }
4194
4195 template <typename VectorType,
4196 std::enable_if_t<has_shared_vector_data<VectorType>, VectorType> * =
4197 nullptr>
4198 const std::vector<ArrayView<const typename VectorType::value_type>> *
4199 get_shared_vector_data(VectorType * vec,
4200 const bool is_valid_mode_for_sm,
4201 const unsigned int active_fe_index,
4203 {
4204 // note: no hp is supported
4205 if (is_valid_mode_for_sm &&
4206 dof_info->dof_indices_contiguous_sm[0 /*any index (<3) should work*/]
4207 .size() > 0 &&
4208 active_fe_index == 0)
4209 return &vec->shared_vector_data();
4210 else
4211 return nullptr;
4212 }
4213
4214 template <typename VectorType,
4215 std::enable_if_t<!has_shared_vector_data<VectorType>, VectorType>
4216 * = nullptr>
4217 const std::vector<ArrayView<const typename VectorType::value_type>> *
4218 get_shared_vector_data(VectorType *,
4219 const bool,
4220 const unsigned int,
4222 {
4223 return nullptr;
4224 }
4225
4226 template <int n_components, typename VectorType>
4227 std::pair<
4228 std::array<typename internal::BlockVectorSelector<
4229 VectorType,
4230 IsBlockVector<VectorType>::value>::BaseVectorType *,
4231 n_components>,
4232 std::array<
4233 const std::vector<ArrayView<const typename internal::BlockVectorSelector<
4234 VectorType,
4235 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4236 n_components>>
4237 get_vector_data(VectorType & src,
4238 const unsigned int first_index,
4239 const bool is_valid_mode_for_sm,
4240 const unsigned int active_fe_index,
4242 {
4243 // select between block vectors and non-block vectors. Note that the number
4244 // of components is checked in the internal data
4245 std::pair<
4246 std::array<typename internal::BlockVectorSelector<
4247 VectorType,
4248 IsBlockVector<VectorType>::value>::BaseVectorType *,
4249 n_components>,
4250 std::array<
4251 const std::vector<
4252 ArrayView<const typename internal::BlockVectorSelector<
4253 VectorType,
4254 IsBlockVector<VectorType>::value>::BaseVectorType::value_type>> *,
4255 n_components>>
4256 src_data;
4257
4258 for (unsigned int d = 0; d < n_components; ++d)
4259 src_data.first[d] = internal::BlockVectorSelector<
4260 VectorType,
4261 IsBlockVector<VectorType>::value>::get_vector_component(src,
4262 d +
4263 first_index);
4264
4265 for (unsigned int d = 0; d < n_components; ++d)
4266 src_data.second[d] = get_shared_vector_data(
4267 const_cast<typename internal::BlockVectorSelector<
4268 typename std::remove_const<VectorType>::type,
4270 BaseVectorType *>(src_data.first[d]),
4271 is_valid_mode_for_sm,
4272 active_fe_index,
4273 dof_info);
4274
4275 return src_data;
4276 }
4277} // namespace internal
4278
4279
4280
4281template <int dim,
4282 int n_components_,
4283 typename Number,
4284 bool is_face,
4285 typename VectorizedArrayType>
4286inline void
4289{
4290 if (this->dof_info == nullptr ||
4291 this->dof_info->hanging_node_constraint_masks.size() == 0 ||
4292 this->dof_info->hanging_node_constraint_masks_comp.size() == 0 ||
4293 this->dof_info->hanging_node_constraint_masks_comp
4294 [this->active_fe_index][this->first_selected_component] == false)
4295 return; // nothing to do with faces
4296
4297 constexpr unsigned int n_lanes = VectorizedArrayType::size();
4298 std::array<internal::MatrixFreeFunctions::compressed_constraint_kind, n_lanes>
4299 constraint_mask;
4300
4301 bool hn_available = false;
4302
4303 const std::array<unsigned int, VectorizedArrayType::size()> &cells =
4304 this->get_cell_ids();
4305
4306 for (unsigned int v = 0; v < n_lanes; ++v)
4307 {
4308 if (cells[v] == numbers::invalid_unsigned_int)
4309 {
4310 constraint_mask[v] = internal::MatrixFreeFunctions::
4312 continue;
4313 }
4314
4315 const unsigned int cell_index = cells[v];
4316 const auto mask =
4318 constraint_mask[v] = mask;
4319
4320 hn_available |= (mask != internal::MatrixFreeFunctions::
4322 }
4323
4324 if (hn_available == false)
4325 return; // no hanging node on cell batch -> nothing to do
4326
4328 apply(n_components,
4329 this->data->data.front().fe_degree,
4330 this->get_shape_info(),
4331 transpose,
4332 constraint_mask,
4333 this->values_dofs);
4334}
4335
4336
4337
4338template <int dim,
4339 int n_components_,
4340 typename Number,
4341 bool is_face,
4342 typename VectorizedArrayType>
4343template <typename VectorType>
4344inline void
4346 read_dof_values(const VectorType & src,
4347 const unsigned int first_index,
4348 const std::bitset<VectorizedArrayType::size()> &mask)
4349{
4350 const auto src_data = internal::get_vector_data<n_components_>(
4351 src,
4352 first_index,
4353 this->dof_access_index ==
4355 this->active_fe_index,
4356 this->dof_info);
4357
4359 read_write_operation(reader, src_data.first, src_data.second, mask, true);
4360
4361 apply_hanging_node_constraints(false);
4362
4363# ifdef DEBUG
4364 this->dof_values_initialized = true;
4365# endif
4366}
4367
4368
4369
4370template <int dim,
4371 int n_components_,
4372 typename Number,
4373 bool is_face,
4374 typename VectorizedArrayType>
4375template <typename VectorType>
4376inline void
4378 read_dof_values_plain(const VectorType & src,
4379 const unsigned int first_index,
4380 const std::bitset<VectorizedArrayType::size()> &mask)
4381{
4382 const auto src_data = internal::get_vector_data<n_components_>(
4383 src,
4384 first_index,
4385 this->dof_access_index ==
4387 this->active_fe_index,
4388 this->dof_info);
4389
4391 read_write_operation(reader, src_data.first, src_data.second, mask, false);
4392
4393# ifdef DEBUG
4394 this->dof_values_initialized = true;
4395# endif
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
4409 VectorType & dst,
4410 const unsigned int first_index,
4411 const std::bitset<VectorizedArrayType::size()> &mask) const
4412{
4413# ifdef DEBUG
4414 Assert(this->dof_values_initialized == true,
4416# endif
4417
4418 apply_hanging_node_constraints(true);
4419
4420 const auto dst_data = internal::get_vector_data<n_components_>(
4421 dst,
4422 first_index,
4423 this->dof_access_index ==
4425 this->active_fe_index,
4426 this->dof_info);
4427
4429 distributor;
4430 read_write_operation(distributor, dst_data.first, dst_data.second, mask);
4431}
4432
4433
4434
4435template <int dim,
4436 int n_components_,
4437 typename Number,
4438 bool is_face,
4439 typename VectorizedArrayType>
4440template <typename VectorType>
4441inline void
4443 set_dof_values(VectorType & dst,
4444 const unsigned int first_index,
4445 const std::bitset<VectorizedArrayType::size()> &mask) const
4446{
4447# ifdef DEBUG
4448 Assert(this->dof_values_initialized == true,
4450# endif
4451
4452 const auto dst_data = internal::get_vector_data<n_components_>(
4453 dst,
4454 first_index,
4455 this->dof_access_index ==
4457 this->active_fe_index,
4458 this->dof_info);
4459
4461 read_write_operation(setter, dst_data.first, dst_data.second, mask);
4462}
4463
4464
4465
4466template <int dim,
4467 int n_components_,
4468 typename Number,
4469 bool is_face,
4470 typename VectorizedArrayType>
4471template <typename VectorType>
4472inline void
4475 VectorType & dst,
4476 const unsigned int first_index,
4477 const std::bitset<VectorizedArrayType::size()> &mask) const
4478{
4479# ifdef DEBUG
4480 Assert(this->dof_values_initialized == true,
4482# endif
4483
4484 const auto dst_data = internal::get_vector_data<n_components_>(
4485 dst,
4486 first_index,
4487 this->dof_access_index ==
4489 this->active_fe_index,
4490 this->dof_info);
4491
4493 read_write_operation(setter, dst_data.first, dst_data.second, mask, false);
4494}
4495
4496
4497
4498/*------------------------------ access to data fields ----------------------*/
4499
4500
4501
4502template <int dim,
4503 int n_components_,
4504 typename Number,
4505 bool is_face,
4506 typename VectorizedArrayType>
4509 get_dof_value(const unsigned int dof) const
4510{
4511 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4512 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4514 for (unsigned int comp = 0; comp < n_components; ++comp)
4515 return_value[comp] = this->values_dofs[comp * dofs + dof];
4516 return return_value;
4517}
4518
4519
4520
4521template <int dim,
4522 int n_components_,
4523 typename Number,
4524 bool is_face,
4525 typename VectorizedArrayType>
4528 get_value(const unsigned int q_point) const
4529{
4530# ifdef DEBUG
4531 Assert(this->values_quad_initialized == true,
4533# endif
4534
4535 AssertIndexRange(q_point, this->n_quadrature_points);
4536 const std::size_t nqp = this->n_quadrature_points;
4538 for (unsigned int comp = 0; comp < n_components; ++comp)
4539 return_value[comp] = this->values_quad[comp * nqp + q_point];
4540 return return_value;
4541}
4542
4543
4544
4545template <int dim,
4546 int n_components_,
4547 typename Number,
4548 bool is_face,
4549 typename VectorizedArrayType>
4553 get_gradient(const unsigned int q_point) const
4554{
4555# ifdef DEBUG
4556 Assert(this->gradients_quad_initialized == true,
4558# endif
4559
4560 AssertIndexRange(q_point, this->n_quadrature_points);
4561 Assert(this->jacobian != nullptr,
4563 "update_gradients"));
4564 const std::size_t nqp = this->n_quadrature_points;
4566
4567 // Cartesian cell
4568 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4569 {
4570 for (unsigned int d = 0; d < dim; ++d)
4571 for (unsigned int comp = 0; comp < n_components; ++comp)
4572 grad_out[comp][d] =
4573 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4574 this->jacobian[0][d][d];
4575 }
4576 // cell with general/affine Jacobian
4577 else
4578 {
4580 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
4581 q_point :
4582 0];
4583 for (unsigned int comp = 0; comp < n_components; ++comp)
4584 for (unsigned int d = 0; d < dim; ++d)
4585 {
4586 grad_out[comp][d] =
4587 jac[d][0] * this->gradients_quad[(comp * dim) * nqp + q_point];
4588 for (unsigned int e = 1; e < dim; ++e)
4589 grad_out[comp][d] +=
4590 jac[d][e] *
4591 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4592 }
4593 }
4594 return grad_out;
4595}
4596
4597
4598
4599template <int dim,
4600 int n_components_,
4601 typename Number,
4602 bool is_face,
4603 typename VectorizedArrayType>
4606 get_normal_derivative(const unsigned int q_point) const
4607{
4608 AssertIndexRange(q_point, this->n_quadrature_points);
4609# ifdef DEBUG
4610 Assert(this->gradients_quad_initialized == true,
4612# endif
4613
4614 Assert(this->normal_x_jacobian != nullptr,
4616 "update_gradients"));
4617
4618 const std::size_t nqp = this->n_quadrature_points;
4620
4621 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4622 for (unsigned int comp = 0; comp < n_components; ++comp)
4623 grad_out[comp] =
4624 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] *
4625 (this->normal_x_jacobian[0][dim - 1]);
4626 else
4627 {
4628 const std::size_t index =
4629 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
4630 for (unsigned int comp = 0; comp < n_components; ++comp)
4631 {
4632 grad_out[comp] = this->gradients_quad[comp * dim * nqp + q_point] *
4633 this->normal_x_jacobian[index][0];
4634 for (unsigned int d = 1; d < dim; ++d)
4635 grad_out[comp] +=
4636 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
4637 this->normal_x_jacobian[index][d];
4638 }
4639 }
4640 return grad_out;
4641}
4642
4643
4644
4645namespace internal
4646{
4647 // compute tmp = hess_unit(u) * J^T. do this manually because we do not
4648 // store the lower diagonal because of symmetry
4649 template <typename VectorizedArrayType>
4650 inline void
4651 hessian_unit_times_jac(const Tensor<2, 1, VectorizedArrayType> &jac,
4652 const VectorizedArrayType *const hessians,
4653 const unsigned int,
4654 VectorizedArrayType (&tmp)[1][1])
4655 {
4656 tmp[0][0] = jac[0][0] * hessians[0];
4657 }
4658
4659 template <typename VectorizedArrayType>
4660 inline void
4661 hessian_unit_times_jac(const Tensor<2, 2, VectorizedArrayType> &jac,
4662 const VectorizedArrayType *const hessians,
4663 const unsigned int nqp,
4664 VectorizedArrayType (&tmp)[2][2])
4665 {
4666 for (unsigned int d = 0; d < 2; ++d)
4667 {
4668 tmp[0][d] = (jac[d][0] * hessians[0] + jac[d][1] * hessians[2 * nqp]);
4669 tmp[1][d] =
4670 (jac[d][0] * hessians[2 * nqp] + jac[d][1] * hessians[1 * nqp]);
4671 }
4672 }
4673
4674 template <typename VectorizedArrayType>
4675 inline void
4676 hessian_unit_times_jac(const Tensor<2, 3, VectorizedArrayType> &jac,
4677 const VectorizedArrayType *const hessians,
4678 const unsigned int nqp,
4679 VectorizedArrayType (&tmp)[3][3])
4680 {
4681 for (unsigned int d = 0; d < 3; ++d)
4682 {
4683 tmp[0][d] =
4684 (jac[d][0] * hessians[0 * nqp] + jac[d][1] * hessians[3 * nqp] +
4685 jac[d][2] * hessians[4 * nqp]);
4686 tmp[1][d] =
4687 (jac[d][0] * hessians[3 * nqp] + jac[d][1] * hessians[1 * nqp] +
4688 jac[d][2] * hessians[5 * nqp]);
4689 tmp[2][d] =
4690 (jac[d][0] * hessians[4 * nqp] + jac[d][1] * hessians[5 * nqp] +
4691 jac[d][2] * hessians[2 * nqp]);
4692 }
4693 }
4694} // namespace internal
4695
4696
4697
4698template <int dim,
4699 int n_components_,
4700 typename Number,
4701 bool is_face,
4702 typename VectorizedArrayType>
4705 get_hessian(const unsigned int q_point) const
4706{
4707# ifdef DEBUG
4708 Assert(this->hessians_quad_initialized == true,
4710# endif
4711 AssertIndexRange(q_point, this->n_quadrature_points);
4712
4713 Assert(this->jacobian != nullptr,
4715 "update_hessian"));
4717 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4718 0 :
4719 q_point];
4720
4722
4723 const std::size_t nqp = this->n_quadrature_points;
4724 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4725
4726 // Cartesian cell
4727 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
4728 {
4729 for (unsigned int comp = 0; comp < n_components; ++comp)
4730 {
4731 for (unsigned int d = 0; d < dim; ++d)
4732 hessian_out[comp][d][d] =
4733 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4734 (jac[d][d] * jac[d][d]);
4735 switch (dim)
4736 {
4737 case 1:
4738 break;
4739 case 2:
4740 hessian_out[comp][0][1] =
4741 this->hessians_quad[(comp * hdim + 2) * nqp + q_point] *
4742 (jac[0][0] * jac[1][1]);
4743 break;
4744 case 3:
4745 hessian_out[comp][0][1] =
4746 this->hessians_quad[(comp * hdim + 3) * nqp + q_point] *
4747 (jac[0][0] * jac[1][1]);
4748 hessian_out[comp][0][2] =
4749 this->hessians_quad[(comp * hdim + 4) * nqp + q_point] *
4750 (jac[0][0] * jac[2][2]);
4751 hessian_out[comp][1][2] =
4752 this->hessians_quad[(comp * hdim + 5) * nqp + q_point] *
4753 (jac[1][1] * jac[2][2]);
4754 break;
4755 default:
4756 Assert(false, ExcNotImplemented());
4757 }
4758 for (unsigned int d = 0; d < dim; ++d)
4759 for (unsigned int e = d + 1; e < dim; ++e)
4760 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4761 }
4762 }
4763 // cell with general Jacobian, but constant within the cell
4764 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4765 {
4766 for (unsigned int comp = 0; comp < n_components; ++comp)
4767 {
4768 VectorizedArrayType tmp[dim][dim];
4769 internal::hessian_unit_times_jac(
4770 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4771
4772 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4773 for (unsigned int d = 0; d < dim; ++d)
4774 for (unsigned int e = d; e < dim; ++e)
4775 {
4776 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4777 for (unsigned int f = 1; f < dim; ++f)
4778 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4779 }
4780
4781 // no J' * grad(u) part here because the Jacobian is constant
4782 // throughout the cell and hence, its derivative is zero
4783
4784 // take symmetric part
4785 for (unsigned int d = 0; d < dim; ++d)
4786 for (unsigned int e = d + 1; e < dim; ++e)
4787 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4788 }
4789 }
4790 // cell with general Jacobian
4791 else
4792 {
4793 const auto &jac_grad = this->jacobian_gradients[q_point];
4794 for (unsigned int comp = 0; comp < n_components; ++comp)
4795 {
4796 VectorizedArrayType tmp[dim][dim];
4797 internal::hessian_unit_times_jac(
4798 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4799
4800 // compute first part of hessian, J * tmp = J * hess_unit(u) * J^T
4801 for (unsigned int d = 0; d < dim; ++d)
4802 for (unsigned int e = d; e < dim; ++e)
4803 {
4804 hessian_out[comp][d][e] = jac[d][0] * tmp[0][e];
4805 for (unsigned int f = 1; f < dim; ++f)
4806 hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
4807 }
4808
4809 // add diagonal part of J' * grad(u)
4810 for (unsigned int d = 0; d < dim; ++d)
4811 for (unsigned int e = 0; e < dim; ++e)
4812 hessian_out[comp][d][d] +=
4813 jac_grad[d][e] *
4814 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4815
4816 // add off-diagonal part of J' * grad(u)
4817 for (unsigned int d = 0, count = dim; d < dim; ++d)
4818 for (unsigned int e = d + 1; e < dim; ++e, ++count)
4819 for (unsigned int f = 0; f < dim; ++f)
4820 hessian_out[comp][d][e] +=
4821 jac_grad[count][f] *
4822 this->gradients_quad[(comp * dim + f) * nqp + q_point];
4823
4824 // take symmetric part
4825 for (unsigned int d = 0; d < dim; ++d)
4826 for (unsigned int e = d + 1; e < dim; ++e)
4827 hessian_out[comp][e][d] = hessian_out[comp][d][e];
4828 }
4829 }
4830 return hessian_out;
4831}
4832
4833
4834
4835template <int dim,
4836 int n_components_,
4837 typename Number,
4838 bool is_face,
4839 typename VectorizedArrayType>
4842 get_hessian_diagonal(const unsigned int q_point) const
4843{
4844 Assert(!is_face, ExcNotImplemented());
4845# ifdef DEBUG
4846 Assert(this->hessians_quad_initialized == true,
4848# endif
4849 AssertIndexRange(q_point, this->n_quadrature_points);
4850
4851 Assert(this->jacobian != nullptr, ExcNotImplemented());
4853 this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
4854 0 :
4855 q_point];
4856
4857 const std::size_t nqp = this->n_quadrature_points;
4858 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
4860
4861 // Cartesian cell
4862 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
4863 {
4864 for (unsigned int comp = 0; comp < n_components; ++comp)
4865 for (unsigned int d = 0; d < dim; ++d)
4866 hessian_out[comp][d] =
4867 this->hessians_quad[(comp * hdim + d) * nqp + q_point] *
4868 (jac[d][d] * jac[d][d]);
4869 }
4870 // cell with general Jacobian, but constant within the cell
4871 else if (this->cell_type == internal::MatrixFreeFunctions::affine)
4872 {
4873 for (unsigned int comp = 0; comp < n_components; ++comp)
4874 {
4875 // compute laplacian before the gradient because it needs to access
4876 // unscaled gradient data
4877 VectorizedArrayType tmp[dim][dim];
4878 internal::hessian_unit_times_jac(
4879 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4880
4881 // compute only the trace part of hessian, J * tmp = J *
4882 // hess_unit(u) * J^T
4883 for (unsigned int d = 0; d < dim; ++d)
4884 {
4885 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4886 for (unsigned int f = 1; f < dim; ++f)
4887 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4888 }
4889 }
4890 }
4891 // cell with general Jacobian
4892 else
4893 {
4894 const auto &jac_grad = this->jacobian_gradients[q_point];
4895 for (unsigned int comp = 0; comp < n_components; ++comp)
4896 {
4897 // compute laplacian before the gradient because it needs to access
4898 // unscaled gradient data
4899 VectorizedArrayType tmp[dim][dim];
4900 internal::hessian_unit_times_jac(
4901 jac, this->hessians_quad + comp * hdim * nqp + q_point, nqp, tmp);
4902
4903 // compute only the trace part of hessian, J * tmp = J *
4904 // hess_unit(u) * J^T
4905 for (unsigned int d = 0; d < dim; ++d)
4906 {
4907 hessian_out[comp][d] = jac[d][0] * tmp[0][d];
4908 for (unsigned int f = 1; f < dim; ++f)
4909 hessian_out[comp][d] += jac[d][f] * tmp[f][d];
4910 }
4911
4912 for (unsigned int d = 0; d < dim; ++d)
4913 for (unsigned int e = 0; e < dim; ++e)
4914 hessian_out[comp][d] +=
4915 jac_grad[d][e] *
4916 this->gradients_quad[(comp * dim + e) * nqp + q_point];
4917 }
4918 }
4919 return hessian_out;
4920}
4921
4922
4923
4924template <int dim,
4925 int n_components_,
4926 typename Number,
4927 bool is_face,
4928 typename VectorizedArrayType>
4931 get_laplacian(const unsigned int q_point) const
4932{
4933 Assert(is_face == false, ExcNotImplemented());
4934# ifdef DEBUG
4935 Assert(this->hessians_quad_initialized == true,
4937# endif
4938 AssertIndexRange(q_point, this->n_quadrature_points);
4939
4941 const auto hess_diag = get_hessian_diagonal(q_point);
4942 for (unsigned int comp = 0; comp < n_components; ++comp)
4943 {
4944 laplacian_out[comp] = hess_diag[comp][0];
4945 for (unsigned int d = 1; d < dim; ++d)
4946 laplacian_out[comp] += hess_diag[comp][d];
4947 }
4948 return laplacian_out;
4949}
4950
4951
4952
4953template <int dim,
4954 int n_components_,
4955 typename Number,
4956 bool is_face,
4957 typename VectorizedArrayType>
4958inline DEAL_II_ALWAYS_INLINE void
4961 const unsigned int dof)
4962{
4963# ifdef DEBUG
4964 this->dof_values_initialized = true;
4965# endif
4966 const std::size_t dofs = this->data->dofs_per_component_on_cell;
4967 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
4968 for (unsigned int comp = 0; comp < n_components; ++comp)
4969 this->values_dofs[comp * dofs + dof] = val_in[comp];
4970}
4971
4972
4973
4974template <int dim,
4975 int n_components_,
4976 typename Number,
4977 bool is_face,
4978 typename VectorizedArrayType>
4979inline DEAL_II_ALWAYS_INLINE void
4982 const unsigned int q_point)
4983{
4984# ifdef DEBUG
4985 Assert(this->is_reinitialized, ExcNotInitialized());
4986# endif
4987 AssertIndexRange(q_point, this->n_quadrature_points);
4988 Assert(this->J_value != nullptr,
4990 "update_values"));
4991# ifdef DEBUG
4992 this->values_quad_submitted = true;
4993# endif
4994
4995 const std::size_t nqp = this->n_quadrature_points;
4996 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
4997 {
4998 const VectorizedArrayType JxW =
4999 this->J_value[0] * this->quadrature_weights[q_point];
5000 for (unsigned int comp = 0; comp < n_components; ++comp)
5001 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5002 }
5003 else
5004 {
5005 const VectorizedArrayType JxW = this->J_value[q_point];
5006 for (unsigned int comp = 0; comp < n_components; ++comp)
5007 this->values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
5008 }
5009}
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 Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in,
5022 const unsigned int q_point)
5023{
5024# ifdef DEBUG
5025 Assert(this->is_reinitialized, ExcNotInitialized());
5026# endif
5027 AssertIndexRange(q_point, this->n_quadrature_points);
5028 Assert(this->J_value != nullptr,
5030 "update_gradients"));
5031 Assert(this->jacobian != nullptr,
5033 "update_gradients"));
5034# ifdef DEBUG
5035 this->gradients_quad_submitted = true;
5036# endif
5037
5038 const std::size_t nqp = this->n_quadrature_points;
5039 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5040 {
5041 const VectorizedArrayType JxW =
5042 this->J_value[0] * this->quadrature_weights[q_point];
5043 std::array<VectorizedArrayType, dim> jac;
5044 for (unsigned int d = 0; d < dim; ++d)
5045 jac[d] = this->jacobian[0][d][d];
5046 for (unsigned int d = 0; d < dim; ++d)
5047 {
5048 const VectorizedArrayType factor = jac[d] * JxW;
5049 for (unsigned int comp = 0; comp < n_components; ++comp)
5050 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5051 grad_in[comp][d] * factor;
5052 }
5053 }
5054 else
5055 {
5057 this->cell_type > internal::MatrixFreeFunctions::affine ?
5058 this->jacobian[q_point] :
5059 this->jacobian[0];
5060 const VectorizedArrayType JxW =
5061 this->cell_type > internal::MatrixFreeFunctions::affine ?
5062 this->J_value[q_point] :
5063 this->J_value[0] * this->quadrature_weights[q_point];
5064 for (unsigned int comp = 0; comp < n_components; ++comp)
5065 for (unsigned int d = 0; d < dim; ++d)
5066 {
5067 VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0];
5068 for (unsigned int e = 1; e < dim; ++e)
5069 new_val += (jac[e][d] * grad_in[comp][e]);
5070 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5071 new_val * JxW;
5072 }
5073 }
5074}
5075
5076
5077
5078template <int dim,
5079 int n_components_,
5080 typename Number,
5081 bool is_face,
5082 typename VectorizedArrayType>
5083inline DEAL_II_ALWAYS_INLINE void
5087 const unsigned int q_point)
5088{
5089 AssertIndexRange(q_point, this->n_quadrature_points);
5090 Assert(this->normal_x_jacobian != nullptr,
5092 "update_gradients"));
5093# ifdef DEBUG
5094 this->gradients_quad_submitted = true;
5095# endif
5096
5097 const std::size_t nqp = this->n_quadrature_points;
5098 if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
5099 {
5100 const VectorizedArrayType JxW_jac = this->J_value[0] *
5101 this->quadrature_weights[q_point] *
5102 this->normal_x_jacobian[0][dim - 1];
5103 for (unsigned int comp = 0; comp < n_components; ++comp)
5104 {
5105 for (unsigned int d = 0; d < dim - 1; ++d)
5106 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5107 VectorizedArrayType();
5108 this->gradients_quad[(comp * dim + dim - 1) * nqp + q_point] =
5109 grad_in[comp] * JxW_jac;
5110 }
5111 }
5112 else
5113 {
5114 const unsigned int index =
5115 this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point;
5117 this->normal_x_jacobian[index];
5118 const VectorizedArrayType JxW =
5119 (this->cell_type <= internal::MatrixFreeFunctions::affine) ?
5120 this->J_value[index] * this->quadrature_weights[q_point] :
5121 this->J_value[index];
5122 for (unsigned int comp = 0; comp < n_components; ++comp)
5123 {
5124 for (unsigned int d = 0; d < dim; ++d)
5125 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
5126 (grad_in[comp] * JxW) * jac[d];
5127 }
5128 }
5129}
5130
5131
5132
5133template <int dim,
5134 int n_components_,
5135 typename Number,
5136 bool is_face,
5137 typename VectorizedArrayType>
5138inline DEAL_II_ALWAYS_INLINE void
5141 const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>
5142 hessian_in,
5143 const unsigned int q_point)
5144{
5145# ifdef DEBUG
5146 Assert(this->is_reinitialized, ExcNotInitialized());
5147# endif
5148 AssertIndexRange(q_point, this->n_quadrature_points);
5149 Assert(this->J_value != nullptr,
5151 "update_hessians"));
5152 Assert(this->jacobian != nullptr,
5154 "update_hessians"));
5155# ifdef DEBUG
5156 this->hessians_quad_submitted = true;
5157# endif
5158
5159 // compute hessian_unit = J^T * hessian_in(u) * J
5160 const std::size_t nqp = this->n_quadrature_points;
5161 constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
5162 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5163 {
5164 const VectorizedArrayType JxW =
5165 this->J_value[0] * this->quadrature_weights[q_point];
5166
5167 // diagonal part
5168 for (unsigned int d = 0; d < dim; ++d)
5169 {
5170 const auto jac_d = this->jacobian[0][d][d];
5171 const VectorizedArrayType factor = jac_d * jac_d * JxW;
5172 for (unsigned int comp = 0; comp < n_components; ++comp)
5173 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5174 hessian_in[comp][d][d] * factor;
5175 }
5176
5177 // off diagonal part
5178 for (unsigned int d = 1, off_dia = dim; d < dim; ++d)
5179 for (unsigned int e = 0; e < d; ++e, ++off_dia)
5180 {
5181 const auto jac_d = this->jacobian[0][d][d];
5182 const auto jac_e = this->jacobian[0][e][e];
5183 const VectorizedArrayType factor = jac_d * jac_e * JxW;
5184 for (unsigned int comp = 0; comp < n_components; ++comp)
5185 this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] =
5186 (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor;
5187 }
5188 }
5189 // cell with general Jacobian, but constant within the cell
5190 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5191 {
5192 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[0];
5193 const VectorizedArrayType JxW =
5194 this->J_value[0] * this->quadrature_weights[q_point];
5195 for (unsigned int comp = 0; comp < n_components; ++comp)
5196 {
5197 // 1. tmp = hessian_in(u) * J
5198 VectorizedArrayType tmp[dim][dim];
5199 for (unsigned int i = 0; i < dim; ++i)
5200 for (unsigned int j = 0; j < dim; ++j)
5201 {
5202 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5203 for (unsigned int k = 1; k < dim; ++k)
5204 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5205 }
5206
5207 // 2. hessian_unit = J^T * tmp
5208 VectorizedArrayType tmp2[dim][dim];
5209 for (unsigned int i = 0; i < dim; ++i)
5210 for (unsigned int j = 0; j < dim; ++j)
5211 {
5212 tmp2[i][j] = jac[0][i] * tmp[0][j];
5213 for (unsigned int k = 1; k < dim; ++k)
5214 tmp2[i][j] += jac[k][i] * tmp[k][j];
5215 }
5216
5217 // diagonal part
5218 for (unsigned int d = 0; d < dim; ++d)
5219 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5220 tmp2[d][d] * JxW;
5221
5222 // off diagonal part
5223 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5224 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5225 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5226 (tmp2[d][e] + tmp2[e][d]) * JxW;
5227 }
5228 }
5229 else
5230 {
5231 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[q_point];
5232 const VectorizedArrayType JxW = this->J_value[q_point];
5233 const auto &jac_grad = this->jacobian_gradients[q_point];
5234 for (unsigned int comp = 0; comp < n_components; ++comp)
5235 {
5236 // 1. tmp = hessian_in(u) * J
5237 VectorizedArrayType tmp[dim][dim];
5238 for (unsigned int i = 0; i < dim; ++i)
5239 for (unsigned int j = 0; j < dim; ++j)
5240 {
5241 tmp[i][j] = hessian_in[comp][i][0] * jac[0][j];
5242 for (unsigned int k = 1; k < dim; ++k)
5243 tmp[i][j] += hessian_in[comp][i][k] * jac[k][j];
5244 }
5245
5246 // 2. hessian_unit = J^T * tmp
5247 VectorizedArrayType tmp2[dim][dim];
5248 for (unsigned int i = 0; i < dim; ++i)
5249 for (unsigned int j = 0; j < dim; ++j)
5250 {
5251 tmp2[i][j] = jac[0][i] * tmp[0][j];
5252 for (unsigned int k = 1; k < dim; ++k)
5253 tmp2[i][j] += jac[k][i] * tmp[k][j];
5254 }
5255
5256 // diagonal part
5257 for (unsigned int d = 0; d < dim; ++d)
5258 this->hessians_quad[(comp * hdim + d) * nqp + q_point] =
5259 tmp2[d][d] * JxW;
5260
5261 // off diagonal part
5262 for (unsigned int d = 0, off_diag = dim; d < dim; ++d)
5263 for (unsigned int e = d + 1; e < dim; ++e, ++off_diag)
5264 this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] =
5265 (tmp2[d][e] + tmp2[e][d]) * JxW;
5266
5267 // 3. gradient_unit = J' ** hessian_in
5268 for (unsigned int d = 0; d < dim; ++d)
5269 {
5270 VectorizedArrayType sum = 0;
5271 for (unsigned int e = 0; e < dim; ++e)
5272 sum += hessian_in[comp][e][e] * jac_grad[e][d];
5273 for (unsigned int e = 0, count = dim; e < dim; ++e)
5274 for (unsigned int f = e + 1; f < dim; ++f, ++count)
5275 sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) *
5276 jac_grad[count][d];
5277 this->gradients_from_hessians_quad[(comp * dim + d) * nqp +
5278 q_point] = sum * JxW;
5279 }
5280 }
5281 }
5282}
5283
5284
5285
5286template <int dim,
5287 int n_components_,
5288 typename Number,
5289 bool is_face,
5290 typename VectorizedArrayType>
5293 integrate_value() const
5294{
5295# ifdef DEBUG
5296 Assert(this->is_reinitialized, ExcNotInitialized());
5297 Assert(this->values_quad_submitted == true,
5299# endif
5300
5302 const std::size_t nqp = this->n_quadrature_points;
5303 for (unsigned int q = 0; q < nqp; ++q)
5304 for (unsigned int comp = 0; comp < n_components; ++comp)
5305 return_value[comp] += this->values_quad[comp * nqp + q];
5306 return (return_value);
5307}
5308
5309
5310
5311/*----------------------- FEEvaluationAccess --------------------------------*/
5312
5313
5314template <int dim,
5315 int n_components_,
5316 typename Number,
5317 bool is_face,
5318 typename VectorizedArrayType>
5319inline FEEvaluationAccess<dim,
5320 n_components_,
5321 Number,
5322 is_face,
5323 VectorizedArrayType>::
5324 FEEvaluationAccess(
5326 const unsigned int dof_no,
5327 const unsigned int first_selected_component,
5328 const unsigned int quad_no,
5329 const unsigned int fe_degree,
5330 const unsigned int n_q_points,
5331 const bool is_interior_face,
5332 const unsigned int active_fe_index,
5333 const unsigned int active_quad_index,
5334 const unsigned int face_type)
5335 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5336 matrix_free,
5337 dof_no,
5338 first_selected_component,
5339 quad_no,
5340 fe_degree,
5341 n_q_points,
5342 is_interior_face,
5343 active_fe_index,
5344 active_quad_index,
5345 face_type)
5346{}
5347
5348
5349
5350template <int dim,
5351 int n_components_,
5352 typename Number,
5353 bool is_face,
5354 typename VectorizedArrayType>
5355inline FEEvaluationAccess<dim,
5356 n_components_,
5357 Number,
5358 is_face,
5359 VectorizedArrayType>::
5360 FEEvaluationAccess(
5361 const Mapping<dim> & mapping,
5362 const FiniteElement<dim> &fe,
5363 const Quadrature<1> & quadrature,
5364 const UpdateFlags update_flags,
5365 const unsigned int first_selected_component,
5367 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5368 mapping,
5369 fe,
5370 quadrature,
5371 update_flags,
5372 first_selected_component,
5373 other)
5374{}
5375
5376
5377
5378template <int dim,
5379 int n_components_,
5380 typename Number,
5381 bool is_face,
5382 typename VectorizedArrayType>
5383inline FEEvaluationAccess<dim,
5384 n_components_,
5385 Number,
5386 is_face,
5387 VectorizedArrayType>::
5388 FEEvaluationAccess(const FEEvaluationAccess<dim,
5389 n_components_,
5390 Number,
5391 is_face,
5392 VectorizedArrayType> &other)
5393 : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
5394 other)
5395{}
5396
5397
5398
5399template <int dim,
5400 int n_components_,
5401 typename Number,
5402 bool is_face,
5403 typename VectorizedArrayType>
5404inline FEEvaluationAccess<dim,
5405 n_components_,
5406 Number,
5407 is_face,
5408 VectorizedArrayType> &
5411 n_components_,
5412 Number,
5413 is_face,
5414 VectorizedArrayType> &other)
5415{
5416 this->FEEvaluationBase<dim,
5417 n_components_,
5418 Number,
5419 is_face,
5420 VectorizedArrayType>::operator=(other);
5421 return *this;
5422}
5423
5424
5425
5426/*-------------------- FEEvaluationAccess scalar ----------------------------*/
5427
5428
5429template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5433 const unsigned int dof_no,
5434 const unsigned int first_selected_component,
5435 const unsigned int quad_no,
5436 const unsigned int fe_degree,
5437 const unsigned int n_q_points,
5438 const bool is_interior_face,
5439 const unsigned int active_fe_index,
5440 const unsigned int active_quad_index,
5441 const unsigned int face_type)
5442 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5443 matrix_free,
5444 dof_no,
5445 first_selected_component,
5446 quad_no,
5447 fe_degree,
5448 n_q_points,
5449 is_interior_face,
5450 active_fe_index,
5451 active_quad_index,
5452 face_type)
5453{}
5454
5455
5456
5457template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5460 const Mapping<dim> & mapping,
5461 const FiniteElement<dim> &fe,
5462 const Quadrature<1> & quadrature,
5463 const UpdateFlags update_flags,
5464 const unsigned int first_selected_component,
5466 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
5467 mapping,
5468 fe,
5469 quadrature,
5470 update_flags,
5471 first_selected_component,
5472 other)
5473{}
5474
5475
5476
5477template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5481 &other)
5482 : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(other)
5483{}
5484
5485
5486
5487template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5491{
5492 this
5493 ->FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>::operator=(
5494 other);
5495 return *this;
5496}
5497
5498
5499
5500template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5501inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5503 const unsigned int dof) const
5504{
5505 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5506 return this->values_dofs[dof];
5507}
5508
5509
5510
5511template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5512inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5514 const unsigned int q_point) const
5515{
5516# ifdef DEBUG
5517 Assert(this->values_quad_initialized == true,
5519# endif
5520 AssertIndexRange(q_point, this->n_quadrature_points);
5521 return this->values_quad[q_point];
5522}
5523
5524
5525
5526template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5527inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
5529 get_normal_derivative(const unsigned int q_point) const
5530{
5531 return BaseClass::get_normal_derivative(q_point)[0];
5532}
5533
5534
5535
5536template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5539 const unsigned int q_point) const
5540{
5541 // could use the base class gradient, but that involves too many expensive
5542 // initialization operations on tensors
5543
5544# ifdef DEBUG
5545 Assert(this->gradients_quad_initialized == true,
5547# endif
5548 AssertIndexRange(q_point, this->n_quadrature_points);
5549
5550 Assert(this->jacobian != nullptr,
5552 "update_gradients"));
5553
5555
5556 const std::size_t nqp = this->n_quadrature_points;
5557 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5558 {
5559 for (unsigned int d = 0; d < dim; ++d)
5560 grad_out[d] =
5561 this->gradients_quad[d * nqp + q_point] * this->jacobian[0][d][d];
5562 }
5563 // cell with general/affine Jacobian
5564 else
5565 {
5567 this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
5568 q_point :
5569 0];
5570 for (unsigned int d = 0; d < dim; ++d)
5571 {
5572 grad_out[d] = jac[d][0] * this->gradients_quad[q_point];
5573 for (unsigned int e = 1; e < dim; ++e)
5574 grad_out[d] += jac[d][e] * this->gradients_quad[e * nqp + q_point];
5575 }
5576 }
5577 return grad_out;
5578}
5579
5580
5581
5582template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5585 const unsigned int q_point) const
5586{
5587 return BaseClass::get_hessian(q_point)[0];
5588}
5589
5590
5591
5592template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5595 get_hessian_diagonal(const unsigned int q_point) const
5596{
5597 return BaseClass::get_hessian_diagonal(q_point)[0];
5598}
5599
5600
5601
5602template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5603inline VectorizedArrayType
5605 const unsigned int q_point) const
5606{
5607 return BaseClass::get_laplacian(q_point)[0];
5608}
5609
5610
5611
5612template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5613inline void DEAL_II_ALWAYS_INLINE
5615 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
5616{
5617# ifdef DEBUG
5618 this->dof_values_initialized = true;
5619 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
5620# endif
5621 this->values_dofs[dof] = val_in;
5622}
5623
5624
5625
5626template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5627inline void DEAL_II_ALWAYS_INLINE
5629 const VectorizedArrayType val_in,
5630 const unsigned int q_point)
5631{
5632# ifdef DEBUG
5633 Assert(this->is_reinitialized, ExcNotInitialized());
5634# endif
5635 AssertIndexRange(q_point, this->n_quadrature_points);
5636 Assert(this->J_value != nullptr,
5638 "update_value"));
5639# ifdef DEBUG
5640 this->values_quad_submitted = true;
5641# endif
5642
5643 if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5644 {
5645 const VectorizedArrayType JxW =
5646 this->J_value[0] * this->quadrature_weights[q_point];
5647 this->values_quad[q_point] = val_in * JxW;
5648 }
5649 else // if (this->cell_type < internal::MatrixFreeFunctions::general)
5650 {
5651 this->values_quad[q_point] = val_in * this->J_value[q_point];
5652 }
5653}
5654
5655
5656
5657template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5658inline DEAL_II_ALWAYS_INLINE void
5661 const unsigned int q_point)
5662{
5663 submit_value(val_in[0], q_point);
5664}
5665
5666
5667
5668template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5669inline DEAL_II_ALWAYS_INLINE void
5671 submit_normal_derivative(const VectorizedArrayType grad_in,
5672 const unsigned int q_point)
5673{
5675 grad[0] = grad_in;
5676 BaseClass::submit_normal_derivative(grad, q_point);
5677}
5678
5679
5680
5681template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5682inline DEAL_II_ALWAYS_INLINE void
5685 const unsigned int q_point)
5686{
5687# ifdef DEBUG
5688 Assert(this->is_reinitialized, ExcNotInitialized());
5689# endif
5690 AssertIndexRange(q_point, this->n_quadrature_points);
5691 Assert(this->J_value != nullptr,
5693 "update_gradients"));
5694 Assert(this->jacobian != nullptr,
5696 "update_gradients"));
5697# ifdef DEBUG
5698 this->gradients_quad_submitted = true;
5699# endif
5700
5701 const std::size_t nqp = this->n_quadrature_points;
5702 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
5703 {
5704 const VectorizedArrayType JxW =
5705 this->J_value[0] * this->quadrature_weights[q_point];
5706
5707 // Make sure the compiler does not think 'jacobian' is aliased with
5708 // 'gradients_quad'
5709 std::array<VectorizedArrayType, dim> jac;
5710 for (unsigned int d = 0; d < dim; ++d)
5711 jac[d] = this->jacobian[0][d][d];
5712
5713 for (unsigned int d = 0; d < dim; ++d)
5714 this->gradients_quad[d * nqp + q_point] = grad_in[d] * jac[d] * JxW;
5715 }
5716 // general/affine cell type
5717 else
5718 {
5720 this->cell_type > internal::MatrixFreeFunctions::affine ?
5721 this->jacobian[q_point] :
5722 this->jacobian[0];
5723 const VectorizedArrayType JxW =
5724 this->cell_type > internal::MatrixFreeFunctions::affine ?
5725 this->J_value[q_point] :
5726 this->J_value[0] * this->quadrature_weights[q_point];
5727 for (unsigned int d = 0; d < dim; ++d)
5728 {
5729 VectorizedArrayType new_val = jac[0][d] * grad_in[0];
5730 for (unsigned int e = 1; e < dim; ++e)
5731 new_val += jac[e][d] * grad_in[e];
5732 this->gradients_quad[d * nqp + q_point] = new_val * JxW;
5733 }
5734 }
5735}
5736
5737
5738
5739template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5740inline DEAL_II_ALWAYS_INLINE void
5743 const unsigned int q_point)
5744{
5746 hessian[0] = hessian_in;
5747 BaseClass::submit_hessian(hessian, q_point);
5748}
5749
5750
5751
5752template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5753inline VectorizedArrayType
5755 integrate_value() const
5756{
5757 return BaseClass::integrate_value()[0];
5758}
5759
5760
5761
5762/*----------------- FEEvaluationAccess vector-valued ------------------------*/
5763
5764
5765template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5769 const unsigned int dof_no,
5770 const unsigned int first_selected_component,
5771 const unsigned int quad_no,
5772 const unsigned int fe_degree,
5773 const unsigned int n_q_points,
5774 const bool is_interior_face,
5775 const unsigned int active_fe_index,
5776 const unsigned int active_quad_index,
5777 const unsigned int face_type)
5778 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5779 matrix_free,
5780 dof_no,
5781 first_selected_component,
5782 quad_no,
5783 fe_degree,
5784 n_q_points,
5785 is_interior_face,
5786 active_fe_index,
5787 active_quad_index,
5788 face_type)
5789{}
5790
5791
5792
5793template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5796 const Mapping<dim> & mapping,
5797 const FiniteElement<dim> &fe,
5798 const Quadrature<1> & quadrature,
5799 const UpdateFlags update_flags,
5800 const unsigned int first_selected_component,
5802 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
5803 mapping,
5804 fe,
5805 quadrature,
5806 update_flags,
5807 first_selected_component,
5808 other)
5809{}
5810
5811
5812
5813template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5817 &other)
5818 : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(other)
5819{}
5820
5821
5822
5823template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5827 &other)
5828{
5830 operator=(other);
5831 return *this;
5832}
5833
5834
5835template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5838 const unsigned int q_point) const
5839{
5840 if (this->data->element_type ==
5842 {
5843 // Piola transform is required
5844# ifdef DEBUG
5845 Assert(this->values_quad_initialized == true,
5847# endif
5848
5849 AssertIndexRange(q_point, this->n_quadrature_points);
5850 Assert(this->J_value != nullptr,
5852 "update_values"));
5853 const std::size_t nqp = this->n_quadrature_points;
5855
5856 if (!is_face &&
5858 {
5859 // Cartesian cell
5860 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
5861 const VectorizedArrayType inv_det =
5862 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5863 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5864 this->jacobian[0][2][2];
5865
5866 // J * u * det(J^-1)
5867 for (unsigned int comp = 0; comp < n_components; ++comp)
5868 value_out[comp] = this->values_quad[comp * nqp + q_point] *
5869 jac[comp][comp] * inv_det;
5870 }
5871 else
5872 {
5873 // Affine or general cell
5874 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5875 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5876 this->jacobian[q_point] :
5877 this->jacobian[0];
5879 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
5880 transpose(invert(inv_t_jac)) :
5881 this->jacobian[1];
5882
5883 // Derivatives are reordered for faces. Need to take this into account
5884 const VectorizedArrayType inv_det =
5885 (is_face && dim == 2 && this->get_face_no() < 2) ?
5886 -determinant(inv_t_jac) :
5887 determinant(inv_t_jac);
5888 // J * u * det(J^-1)
5889 for (unsigned int comp = 0; comp < n_components; ++comp)
5890 {
5891 value_out[comp] =
5892 this->values_quad[q_point] * jac[comp][0] * inv_det;
5893 for (unsigned int e = 1; e < dim; ++e)
5894 value_out[comp] +=
5895 this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det;
5896 }
5897 }
5898 return value_out;
5899 }
5900 else
5901 {
5902 // No Piola needed
5903 return BaseClass::get_value(q_point);
5904 }
5905}
5906
5907template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
5910 get_gradient(const unsigned int q_point) const
5911{
5912 if (this->data->element_type ==
5914 {
5915 // Piola transform is required
5916# ifdef DEBUG
5917 Assert(this->gradients_quad_initialized == true,
5919# endif
5920
5921 AssertIndexRange(q_point, this->n_quadrature_points);
5922 Assert(this->jacobian != nullptr,
5924 "update_gradients"));
5925 const std::size_t nqp = this->n_quadrature_points;
5927
5928 if (!is_face &&
5930 {
5931 // Cartesian cell
5932 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5933 this->jacobian[0];
5934 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5935 const VectorizedArrayType inv_det =
5936 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
5937 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
5938 this->jacobian[0][2][2];
5939
5940 // J * grad_quad * J^-1 * det(J^-1)
5941 for (unsigned int d = 0; d < dim; ++d)
5942 for (unsigned int comp = 0; comp < n_components; ++comp)
5943 grad_out[comp][d] =
5944 this->gradients_quad[(comp * dim + d) * nqp + q_point] *
5945 inv_t_jac[d][d] * jac[comp][comp] * inv_det;
5946 }
5947 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
5948 {
5949 // Affine cell
5950 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5951 this->jacobian[0];
5952 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
5953
5954 // Derivatives are reordered for faces. Need to take this into account
5955 const VectorizedArrayType inv_det =
5956 (is_face && dim == 2 && this->get_face_no() < 2) ?
5957 -determinant(inv_t_jac) :
5958 determinant(inv_t_jac);
5959
5960 VectorizedArrayType tmp;
5961 // J * grad_quad * J^-1 * det(J^-1)
5962 for (unsigned int comp = 0; comp < n_components; ++comp)
5963 for (unsigned int d = 0; d < dim; ++d)
5964 {
5965 tmp = 0;
5966 for (unsigned int f = 0; f < dim; ++f)
5967 for (unsigned int e = 0; e < dim; ++e)
5968 tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det *
5969 this->gradients_quad[(f * dim + e) * nqp + q_point];
5970
5971 grad_out[comp][d] = tmp;
5972 }
5973 }
5974 else
5975 {
5976 // General cell
5977
5978 // This assert could be removed if we make sure that this is updated
5979 // even though update_hessians or update_jacobian_grads is not passed,
5980 // i.e make the necessary changes in
5981 // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags
5982 Assert(this->jacobian_gradients_non_inverse != nullptr,
5984 "update_hessians"));
5985
5986 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
5987 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
5988 this->jacobian[q_point];
5989 const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
5990
5991 // Derivatives are reordered for faces. Need to take this into account
5992 const VectorizedArrayType inv_det =
5993 (is_face && dim == 2 && this->get_face_no() < 2) ?
5994 -determinant(inv_t_jac) :
5995 determinant(inv_t_jac);
5996
5997 VectorizedArrayType tmp;
5998 // J * grad_quad * J^-1 * det(J^-1)
5999 for (unsigned int comp = 0; comp < n_components; ++comp)
6000 for (unsigned int d = 0; d < dim; ++d)
6001 {
6002 tmp = 0;
6003 for (unsigned int f = 0; f < dim; ++f)
6004 for (unsigned int e = 0; e < dim; ++e)
6005 tmp += t_jac[f][comp] * inv_t_jac[d][e] *
6006 this->gradients_quad[(f * dim + e) * nqp + q_point];
6007
6008 grad_out[comp][d] = tmp * inv_det;
6009 }
6010
6011 // Contribution from values
6012 {
6013 // Diagonal part of jac_grad
6014
6015 // Add jac_grad * J^{-1} * values * det(J^{-1})
6016 // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
6017 for (unsigned int i = 0; i < dim; ++i)
6018 for (unsigned int j = 0; j < dim; ++j)
6019 {
6020 tmp = jac_grad[0][i] * inv_t_jac[j][0] *
6021 this->values_quad[q_point];
6022 for (unsigned int f = 1; f < dim; ++f)
6023 tmp += jac_grad[f][i] * inv_t_jac[j][f] *
6024 this->values_quad[f * nqp + q_point];
6025
6026 grad_out[i][j] += tmp * inv_det;
6027 }
6028
6029 for (unsigned int i = 0; i < dim; ++i)
6030 for (unsigned int j = 0; j < dim; ++j)
6031 {
6032 tmp = 0;
6033 for (unsigned int f = 0; f < dim; ++f)
6034 for (unsigned int n = 0; n < dim; ++n)
6035 for (unsigned int m = 0; m < dim; ++m)
6036 tmp += inv_t_jac[m][f] * jac_grad[f][m] *
6037 inv_t_jac[j][f] * t_jac[n][i] *
6038 this->values_quad[n * nqp + q_point];
6039 grad_out[i][j] -= tmp * inv_det;
6040 }
6041 }
6042
6043 {
6044 // Off-diagonal part of jac_grad
6045
6046 // Add jac_grad * J^{-1} * values * det(J^{-1})
6047 // -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1}))
6048 for (unsigned int i = 0; i < dim; ++i)
6049 for (unsigned int j = 0; j < dim; ++j)
6050 {
6051 tmp = 0;
6052 for (unsigned int r = 0, f = dim; r < dim; ++r)
6053 for (unsigned int k = r + 1; k < dim; ++k, ++f)
6054 {
6055 tmp += jac_grad[f][i] *
6056 (inv_t_jac[j][k] *
6057 this->values_quad[r * nqp + q_point] +
6058 inv_t_jac[j][r] *
6059 this->values_quad[k * nqp + q_point]);
6060 for (unsigned int n = 0; n < dim; ++n)
6061 for (unsigned int m = 0; m < dim; ++m)
6062 tmp -= jac_grad[f][m] * t_jac[n][i] *
6063 this->values_quad[n * nqp + q_point] *
6064 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6065 inv_t_jac[m][r] * inv_t_jac[j][k]);
6066 }
6067 grad_out[i][j] += tmp * inv_det;
6068 }
6069 }
6070 }
6071 return grad_out;
6072 }
6073 else
6074 {
6075 return BaseClass::get_gradient(q_point);
6076 }
6077}
6078
6079
6080
6081template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6082inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6084 get_divergence(const unsigned int q_point) const
6085{
6086# ifdef DEBUG
6087 Assert(this->gradients_quad_initialized == true,
6089# endif
6090 AssertIndexRange(q_point, this->n_quadrature_points);
6091 Assert(this->jacobian != nullptr,
6093 "update_gradients"));
6094
6095 VectorizedArrayType divergence;
6096 const std::size_t nqp = this->n_quadrature_points;
6097
6098 if (this->data->element_type ==
6100 {
6101 if (!is_face &&
6103 {
6104 // Cartesian cell
6105 const VectorizedArrayType inv_det =
6106 (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] :
6107 this->jacobian[0][0][0] * this->jacobian[0][1][1] *
6108 this->jacobian[0][2][2];
6109
6110 // div * det(J^-1)
6111 divergence = this->gradients_quad[q_point] * inv_det;
6112 for (unsigned int d = 1; d < dim; ++d)
6113 divergence +=
6114 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
6115 }
6116 else
6117 {
6118 // General cell
6119 // Derivatives are reordered for faces. Need to take this into account
6120 const VectorizedArrayType inv_det =
6122 this->jacobian[this->cell_type >
6124 q_point :
6125 0]) *
6126 Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6127
6128 // div * det(J^-1)
6129 divergence = this->gradients_quad[q_point] * inv_det;
6130 for (unsigned int d = 1; d < dim; ++d)
6131 divergence +=
6132 this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det;
6133 }
6134 }
6135 else
6136 {
6137 if (!is_face &&
6139 {
6140 // Cartesian cell
6141 divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0];
6142 for (unsigned int d = 1; d < dim; ++d)
6143 divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] *
6144 this->jacobian[0][d][d];
6145 }
6146 else
6147 {
6148 // cell with general/constant Jacobian
6150 this->cell_type == internal::MatrixFreeFunctions::general ?
6151 this->jacobian[q_point] :
6152 this->jacobian[0];
6153 divergence = jac[0][0] * this->gradients_quad[q_point];
6154 for (unsigned int e = 1; e < dim; ++e)
6155 divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point];
6156 for (unsigned int d = 1; d < dim; ++d)
6157 for (unsigned int e = 0; e < dim; ++e)
6158 divergence +=
6159 jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point];
6160 }
6161 }
6162 return divergence;
6163}
6164
6165
6166
6167template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6170 get_symmetric_gradient(const unsigned int q_point) const
6171{
6172 // copy from generic function into dim-specialization function
6173 const auto grad = get_gradient(q_point);
6174 VectorizedArrayType symmetrized[(dim * dim + dim) / 2];
6175 VectorizedArrayType half = Number(0.5);
6176 for (unsigned int d = 0; d < dim; ++d)
6177 symmetrized[d] = grad[d][d];
6178 switch (dim)
6179 {
6180 case 1:
6181 break;
6182 case 2:
6183 symmetrized[2] = grad[0][1] + grad[1][0];
6184 symmetrized[2] *= half;
6185 break;
6186 case 3:
6187 symmetrized[3] = grad[0][1] + grad[1][0];
6188 symmetrized[3] *= half;
6189 symmetrized[4] = grad[0][2] + grad[2][0];
6190 symmetrized[4] *= half;
6191 symmetrized[5] = grad[1][2] + grad[2][1];
6192 symmetrized[5] *= half;
6193 break;
6194 default:
6195 Assert(false, ExcNotImplemented());
6196 }
6198}
6199
6200
6201
6202template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6204 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType>
6206 const unsigned int q_point) const
6207{
6208 // copy from generic function into dim-specialization function
6209 const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point);
6210 Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl;
6211 switch (dim)
6212 {
6213 case 1:
6214 Assert(false,
6215 ExcMessage(
6216 "Computing the curl in 1d is not a useful operation"));
6217 break;
6218 case 2:
6219 curl[0] = grad[1][0] - grad[0][1];
6220 break;
6221 case 3:
6222 curl[0] = grad[2][1] - grad[1][2];
6223 curl[1] = grad[0][2] - grad[2][0];
6224 curl[2] = grad[1][0] - grad[0][1];
6225 break;
6226 default:
6227 Assert(false, ExcNotImplemented());
6228 }
6229 return curl;
6230}
6231
6232
6233
6234template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6237 get_hessian_diagonal(const unsigned int q_point) const
6238{
6239 return BaseClass::get_hessian_diagonal(q_point);
6240}
6241
6242
6243
6244template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6247 const unsigned int q_point) const
6248{
6249# ifdef DEBUG
6250 Assert(this->hessians_quad_initialized == true,
6252# endif
6253 AssertIndexRange(q_point, this->n_quadrature_points);
6254 return BaseClass::get_hessian(q_point);
6255}
6256
6257
6258template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6259inline DEAL_II_ALWAYS_INLINE void
6262 const unsigned int q_point)
6263{
6264 if (this->data->element_type ==
6266 {
6267 // Piola transform is required
6268 AssertIndexRange(q_point, this->n_quadrature_points);
6269 Assert(this->J_value != nullptr,
6271 "update_value"));
6272# ifdef DEBUG
6273 Assert(this->is_reinitialized, ExcNotInitialized());
6274 this->values_quad_submitted = true;
6275# endif
6276
6277 const std::size_t nqp = this->n_quadrature_points;
6278 if (!is_face &&
6280 {
6281 const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1];
6282 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6283
6284 for (unsigned int comp = 0; comp < n_components; ++comp)
6285 this->values_quad[comp * nqp + q_point] =
6286 val_in[comp] * weight * jac[comp][comp];
6287 }
6288 else
6289 {
6290 // Affine or general cell
6291 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6292 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6293 this->jacobian[q_point] :
6294 this->jacobian[0];
6296 (this->cell_type > internal::MatrixFreeFunctions::affine) ?
6297 transpose(invert(inv_t_jac)) :
6298 this->jacobian[1];
6299
6300 // Derivatives are reordered for faces. Need to take this into account
6301 // and 1/inv_det != J_value for faces
6302 const VectorizedArrayType fac =
6303 (!is_face) ?
6304 this->quadrature_weights[q_point] :
6305 (((this->cell_type > internal::MatrixFreeFunctions::affine) ?
6306 this->J_value[q_point] :
6307 this->J_value[0] * this->quadrature_weights[q_point]) *
6308 ((dim == 2 && this->get_face_no() < 2) ?
6309 -determinant(inv_t_jac) :
6310 determinant(inv_t_jac)));
6311
6312 // J^T * u * factor
6313 for (unsigned int comp = 0; comp < n_components; ++comp)
6314 {
6315 this->values_quad[comp * nqp + q_point] =
6316 val_in[0] * jac[0][comp] * fac;
6317 for (unsigned int e = 1; e < dim; ++e)
6318 this->values_quad[comp * nqp + q_point] +=
6319 val_in[e] * jac[e][comp] * fac;
6320 }
6321 }
6322 }
6323 else
6324 {
6325 // No Piola transform
6326 BaseClass::submit_value(val_in, q_point);
6327 }
6328}
6329
6330template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6331inline DEAL_II_ALWAYS_INLINE void
6334 const unsigned int q_point)
6335{
6336 if (this->data->element_type ==
6338 {
6339 // Piola transform is required
6340
6341# ifdef DEBUG
6342 Assert(this->is_reinitialized, ExcNotInitialized());
6343# endif
6344 AssertIndexRange(q_point, this->n_quadrature_points);
6345 Assert(this->J_value != nullptr,
6347 "update_gradients"));
6348 Assert(this->jacobian != nullptr,
6350 "update_gradients"));
6351# ifdef DEBUG
6352 this->gradients_quad_submitted = true;
6353# endif
6354
6355 const std::size_t nqp = this->n_quadrature_points;
6356 if (!is_face &&
6358 {
6359 // Cartesian cell
6360 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6361 this->jacobian[0];
6362 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6363 const VectorizedArrayType weight = this->quadrature_weights[q_point];
6364 for (unsigned int d = 0; d < dim; ++d)
6365 for (unsigned int comp = 0; comp < n_components; ++comp)
6366 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6367 grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight;
6368 }
6369 else if (this->cell_type <= internal::MatrixFreeFunctions::affine)
6370 {
6371 // Affine cell
6372 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6373 this->jacobian[0];
6374 const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
6375
6376 // Derivatives are reordered for faces. Need to take this into account
6377 // and 1/inv_det != J_value for faces
6378 const VectorizedArrayType fac =
6379 (!is_face) ? this->quadrature_weights[q_point] :
6380 this->J_value[0] * this->quadrature_weights[q_point] *
6381 ((dim == 2 && this->get_face_no() < 2) ?
6382 -determinant(inv_t_jac) :
6383 determinant(inv_t_jac));
6384
6385 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6386 for (unsigned int comp = 0; comp < n_components; ++comp)
6387 for (unsigned int d = 0; d < dim; ++d)
6388 {
6389 VectorizedArrayType tmp = 0;
6390 for (unsigned int f = 0; f < dim; ++f)
6391 for (unsigned int e = 0; e < dim; ++e)
6392 tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e];
6393
6394 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6395 tmp * fac;
6396 }
6397 }
6398 else
6399 {
6400 // General cell
6401
6402 const auto &jac_grad = this->jacobian_gradients_non_inverse[q_point];
6403 const Tensor<2, dim, VectorizedArrayType> &inv_t_jac =
6404 this->jacobian[q_point];
6405 const Tensor<2, dim, VectorizedArrayType> &t_jac = invert(inv_t_jac);
6406
6407 // Derivatives are reordered for faces. Need to take this into account
6408 // and 1/inv_det != J_value for faces
6409 const VectorizedArrayType fac =
6410 (!is_face) ?
6411 this->quadrature_weights[q_point] :
6412 this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ?
6413 -determinant(inv_t_jac) :
6414 determinant(inv_t_jac));
6415
6416 VectorizedArrayType tmp;
6417 // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor
6418 for (unsigned int comp = 0; comp < n_components; ++comp)
6419 for (unsigned int d = 0; d < dim; ++d)
6420 {
6421 tmp = 0;
6422 for (unsigned int f = 0; f < dim; ++f)
6423 for (unsigned int e = 0; e < dim; ++e)
6424 tmp += t_jac[comp][f] * inv_t_jac[e][d] * grad_in[f][e];
6425
6426 this->gradients_quad[(comp * dim + d) * nqp + q_point] =
6427 tmp * fac;
6428 }
6429
6430 // Contribution from values
6431 {
6432 // Diagonal part of jac_grad
6433
6434 // Add jac_grad * J^{-1} * values * factor
6435 // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
6436 for (unsigned int f = 0; f < dim; ++f)
6437 {
6438 tmp = 0;
6439 for (unsigned int i = 0; i < dim; ++i)
6440 for (unsigned int j = 0; j < dim; ++j)
6441 {
6442 tmp += inv_t_jac[j][f] * jac_grad[f][i] * grad_in[i][j];
6443 for (unsigned int m = 0; m < dim; ++m)
6444 for (unsigned int k = 0; k < dim; ++k)
6445 tmp -= inv_t_jac[m][k] * jac_grad[k][m] *
6446 inv_t_jac[j][k] * t_jac[f][i] * grad_in[i][j];
6447 }
6448 this->values_from_gradients_quad[f * nqp + q_point] = tmp * fac;
6449 }
6450 }
6451
6452 {
6453 // Off-diagonal part of jac_grad
6454
6455 // Add jac_grad * J^{-1} * values * factor
6456 for (unsigned int r = 0, f = dim; r < dim; ++r)
6457 for (unsigned int k = r + 1; k < dim; ++k, ++f)
6458 {
6459 tmp = jac_grad[f][0] * inv_t_jac[0][k] * grad_in[0][0];
6460 for (unsigned int j = 1; j < dim; ++j)
6461 tmp += jac_grad[f][0] * inv_t_jac[j][k] * grad_in[0][j];
6462 for (unsigned int i = 1; i < dim; ++i)
6463 for (unsigned int j = 0; j < dim; ++j)
6464 tmp += jac_grad[f][i] * inv_t_jac[j][k] * grad_in[i][j];
6465 this->values_from_gradients_quad[r * nqp + q_point] +=
6466 tmp * fac;
6467
6468 tmp = jac_grad[f][0] * inv_t_jac[0][r] * grad_in[0][0];
6469 for (unsigned int j = 1; j < dim; ++j)
6470 tmp += jac_grad[f][0] * inv_t_jac[j][r] * grad_in[0][j];
6471 for (unsigned int i = 1; i < dim; ++i)
6472 for (unsigned int j = 0; j < dim; ++j)
6473 tmp += jac_grad[f][i] * inv_t_jac[j][r] * grad_in[i][j];
6474 this->values_from_gradients_quad[k * nqp + q_point] +=
6475 tmp * fac;
6476 }
6477
6478 // -(J^{-T} * jac_grad * J^{-1} * J * values * factor)
6479 for (unsigned int n = 0; n < dim; ++n)
6480 {
6481 tmp = 0;
6482 for (unsigned int r = 0, f = dim; r < dim; ++r)
6483 for (unsigned int k = r + 1; k < dim; ++k, ++f)
6484 for (unsigned int i = 0; i < dim; ++i)
6485 for (unsigned int j = 0; j < dim; ++j)
6486 for (unsigned int m = 0; m < dim; ++m)
6487 tmp += jac_grad[f][m] * t_jac[n][i] * grad_in[i][j] *
6488 (inv_t_jac[m][k] * inv_t_jac[j][r] +
6489 inv_t_jac[m][r] * inv_t_jac[j][k]);
6490
6491 this->values_from_gradients_quad[n * nqp + q_point] -=
6492 tmp * fac;
6493 }
6494 }
6495 }
6496 }
6497 else
6498 {
6499 BaseClass::submit_gradient(grad_in, q_point);
6500 }
6501}
6502
6503
6504
6505template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6506inline DEAL_II_ALWAYS_INLINE void
6509 const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
6510 const unsigned int q_point)
6511{
6512 if (this->data->element_type ==
6514 {
6515 // Piola transform is required
6516 const Tensor<2, dim, VectorizedArrayType> &grad = grad_in;
6518 submit_gradient(grad, q_point);
6519 }
6520 else
6521 {
6522 BaseClass::submit_gradient(grad_in, q_point);
6523 }
6524}
6525
6526
6527
6528template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6529inline DEAL_II_ALWAYS_INLINE void
6531 submit_divergence(const VectorizedArrayType div_in,
6532 const unsigned int q_point)
6533{
6534# ifdef DEBUG
6535 Assert(this->is_reinitialized, ExcNotInitialized());
6536# endif
6537 AssertIndexRange(q_point, this->n_quadrature_points);
6538 Assert(this->J_value != nullptr,
6540 "update_gradients"));
6541 Assert(this->jacobian != nullptr,
6543 "update_gradients"));
6544# ifdef DEBUG
6545 this->gradients_quad_submitted = true;
6546# endif
6547
6548 const std::size_t nqp = this->n_quadrature_points;
6549 if (this->data->element_type ==
6551 {
6552 // General cell
6553
6554 // Derivatives are reordered for faces. Need to take this into account
6555 // and 1/inv_det != J_value for faces
6556 const VectorizedArrayType fac =
6557 (!is_face) ?
6558 this->quadrature_weights[q_point] * div_in :
6559 (this->cell_type > internal::MatrixFreeFunctions::affine ?
6560 this->J_value[q_point] :
6561 this->J_value[0] * this->quadrature_weights[q_point]) *
6562 div_in *
6564 this->jacobian[this->cell_type >
6565 internal::MatrixFreeFunctions::affine ?
6566 q_point :
6567 0]) *
6568 Number((dim == 2 && this->get_face_no() < 2) ? -1 : 1);
6569
6570 for (unsigned int d = 0; d < dim; ++d)
6571 {
6572 this->gradients_quad[(dim * d + d) * nqp + q_point] = fac;
6573 for (unsigned int e = d + 1; e < dim; ++e)
6574 {
6575 this->gradients_quad[(dim * d + e) * nqp + q_point] =
6576 VectorizedArrayType();
6577 this->gradients_quad[(dim * e + d) * nqp + q_point] =
6578 VectorizedArrayType();
6579 }
6580 }
6581 this->divergence_is_requested = true;
6582 }
6583 else
6584 {
6585 if (!is_face &&
6587 {
6588 const VectorizedArrayType fac =
6589 this->J_value[0] * this->quadrature_weights[q_point] * div_in;
6590 for (unsigned int d = 0; d < dim; ++d)
6591 {
6592 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6593 (fac * this->jacobian[0][d][d]);
6594 for (unsigned int e = d + 1; e < dim; ++e)
6595 {
6596 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6597 VectorizedArrayType();
6598 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6599 VectorizedArrayType();
6600 }
6601 }
6602 }
6603 else
6604 {
6606 this->cell_type == internal::MatrixFreeFunctions::general ?
6607 this->jacobian[q_point] :
6608 this->jacobian[0];
6609 const VectorizedArrayType fac =
6610 (this->cell_type == internal::MatrixFreeFunctions::general ?
6611 this->J_value[q_point] :
6612 this->J_value[0] * this->quadrature_weights[q_point]) *
6613 div_in;
6614 for (unsigned int d = 0; d < dim; ++d)
6615 {
6616 for (unsigned int e = 0; e < dim; ++e)
6617 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6618 jac[d][e] * fac;
6619 }
6620 }
6621 }
6622}
6623
6624
6625
6626template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6627inline DEAL_II_ALWAYS_INLINE void
6631 const unsigned int q_point)
6632{
6634 this->data->element_type !=
6637
6638 // could have used base class operator, but that involves some overhead
6639 // which is inefficient. it is nice to have the symmetric tensor because
6640 // that saves some operations
6641# ifdef DEBUG
6642 Assert(this->is_reinitialized, ExcNotInitialized());
6643# endif
6644 AssertIndexRange(q_point, this->n_quadrature_points);
6645 Assert(this->J_value != nullptr,
6647 "update_gradients"));
6648 Assert(this->jacobian != nullptr,
6650 "update_gradients"));
6651# ifdef DEBUG
6652 this->gradients_quad_submitted = true;
6653# endif
6654
6655 const std::size_t nqp = this->n_quadrature_points;
6656 if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
6657 {
6658 const VectorizedArrayType JxW =
6659 this->J_value[0] * this->quadrature_weights[q_point];
6660 for (unsigned int d = 0; d < dim; ++d)
6661 this->gradients_quad[(d * dim + d) * nqp + q_point] =
6662 (sym_grad.access_raw_entry(d) * JxW * this->jacobian[0][d][d]);
6663 for (unsigned int e = 0, counter = dim; e < dim; ++e)
6664 for (unsigned int d = e + 1; d < dim; ++d, ++counter)
6665 {
6666 const VectorizedArrayType value =
6667 sym_grad.access_raw_entry(counter) * JxW;
6668 this->gradients_quad[(e * dim + d) * nqp + q_point] =
6669 value * this->jacobian[0][d][d];
6670 this->gradients_quad[(d * dim + e) * nqp + q_point] =
6671 value * this->jacobian[0][e][e];
6672 }
6673 }
6674 // general/affine cell type
6675 else
6676 {
6677 const VectorizedArrayType JxW =
6678 this->cell_type == internal::MatrixFreeFunctions::general ?
6679 this->J_value[q_point] :
6680 this->J_value[0] * this->quadrature_weights[q_point];
6682 this->cell_type == internal::MatrixFreeFunctions::general ?
6683 this->jacobian[q_point] :
6684 this->jacobian[0];
6685 VectorizedArrayType weighted[dim][dim];
6686 for (unsigned int i = 0; i < dim; ++i)
6687 weighted[i][i] = sym_grad.access_raw_entry(i) * JxW;
6688 for (unsigned int i = 0, counter = dim; i < dim; ++i)
6689 for (unsigned int j = i + 1; j < dim; ++j, ++counter)
6690 {
6691 const VectorizedArrayType value =
6692 sym_grad.access_raw_entry(counter) * JxW;
6693 weighted[i][j] = value;
6694 weighted[j][i] = value;
6695 }
6696 for (unsigned int comp = 0; comp < dim; ++comp)
6697 for (unsigned int d = 0; d < dim; ++d)
6698 {
6699 VectorizedArrayType new_val = jac[0][d] * weighted[comp][0];
6700 for (unsigned int e = 1; e < dim; ++e)
6701 new_val += jac[e][d] * weighted[comp][e];
6702 this->gradients_quad[(comp * dim + d) * nqp + q_point] = new_val;
6703 }
6704 }
6705}
6706
6707
6708
6709template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
6710inline DEAL_II_ALWAYS_INLINE void
6713 const unsigned int q_point)
6714{
6716 switch (dim)
6717 {
6718 case 1:
6719 Assert(false,
6720 ExcMessage(
6721 "Testing by the curl in 1d is not a useful operation"));
6722 break;
6723 case 2:
6724 grad[1][0] = curl[0];
6725 grad[0][1] = -curl[0];
6726 break;
6727 case 3:
6728 grad[2][1] = curl[0];
6729 grad[1][2] = -curl[0];
6730 grad[0][2] = curl[1];
6731 grad[2][0] = -curl[1];
6732 grad[1][0] = curl[2];
6733 grad[0][1] = -curl[2];
6734 break;
6735 default:
6736 Assert(false, ExcNotImplemented());
6737 }
6738 submit_gradient(grad, q_point);
6739}
6740
6741
6742/*-------------------- FEEvaluationAccess scalar for 1d ---------------------*/
6743
6744
6745template <typename Number, bool is_face, typename VectorizedArrayType>
6749 const unsigned int dof_no,
6750 const unsigned int first_selected_component,
6751 const unsigned int quad_no,
6752 const unsigned int fe_degree,
6753 const unsigned int n_q_points,
6754 const bool is_interior_face,
6755 const unsigned int active_fe_index,
6756 const unsigned int active_quad_index,
6757 const unsigned int face_type)
6758 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6759 matrix_free,
6760 dof_no,
6761 first_selected_component,
6762 quad_no,
6763 fe_degree,
6764 n_q_points,
6765 is_interior_face,
6766 active_fe_index,
6767 active_quad_index,
6768 face_type)
6769{}
6770
6771
6772
6773template <typename Number, bool is_face, typename VectorizedArrayType>
6776 const Mapping<1> & mapping,
6777 const FiniteElement<1> &fe,
6778 const Quadrature<1> & quadrature,
6779 const UpdateFlags update_flags,
6780 const unsigned int first_selected_component,
6782 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
6783 mapping,
6784 fe,
6785 quadrature,
6786 update_flags,
6787 first_selected_component,
6788 other)
6789{}
6790
6791
6792
6793template <typename Number, bool is_face, typename VectorizedArrayType>
6797 : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(other)
6798{}
6799
6800
6801
6802template <typename Number, bool is_face, typename VectorizedArrayType>
6806{
6808 other);
6809 return *this;
6810}
6811
6812
6813
6814template <typename Number, bool is_face, typename VectorizedArrayType>
6815inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6817 const unsigned int dof) const
6818{
6819 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6820 return this->values_dofs[dof];
6821}
6822
6823
6824
6825template <typename Number, bool is_face, typename VectorizedArrayType>
6826inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6828 const unsigned int q_point) const
6829{
6830# ifdef DEBUG
6831 Assert(this->values_quad_initialized == true,
6833# endif
6834 AssertIndexRange(q_point, this->n_quadrature_points);
6835 return this->values_quad[q_point];
6836}
6837
6838
6839
6840template <typename Number, bool is_face, typename VectorizedArrayType>
6843 const unsigned int q_point) const
6844{
6845 // could use the base class gradient, but that involves too many inefficient
6846 // initialization operations on tensors
6847
6848# ifdef DEBUG
6849 Assert(this->gradients_quad_initialized == true,
6851# endif
6852 AssertIndexRange(q_point, this->n_quadrature_points);
6853
6855 this->cell_type == internal::MatrixFreeFunctions::general ?
6856 this->jacobian[q_point] :
6857 this->jacobian[0];
6858
6860 grad_out[0] = jac[0][0] * this->gradients_quad[q_point];
6861
6862 return grad_out;
6863}
6864
6865
6866
6867template <typename Number, bool is_face, typename VectorizedArrayType>
6868inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6870 const unsigned int q_point) const
6871{
6872 return get_gradient(q_point)[0];
6873}
6874
6875
6876
6877template <typename Number, bool is_face, typename VectorizedArrayType>
6878inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6880 get_normal_derivative(const unsigned int q_point) const
6881{
6882 return BaseClass::get_normal_derivative(q_point)[0];
6883}
6884
6885
6886
6887template <typename Number, bool is_face, typename VectorizedArrayType>
6890 const unsigned int q_point) const
6891{
6892 return BaseClass::get_hessian(q_point)[0];
6893}
6894
6895
6896
6897template <typename Number, bool is_face, typename VectorizedArrayType>
6900 get_hessian_diagonal(const unsigned int q_point) const
6901{
6902 return BaseClass::get_hessian_diagonal(q_point)[0];
6903}
6904
6905
6906
6907template <typename Number, bool is_face, typename VectorizedArrayType>
6908inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
6910 const unsigned int q_point) const
6911{
6912 return BaseClass::get_laplacian(q_point)[0];
6913}
6914
6915
6916
6917template <typename Number, bool is_face, typename VectorizedArrayType>
6920 submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof)
6921{
6922# ifdef DEBUG
6923 this->dof_values_initialized = true;
6924 AssertIndexRange(dof, this->data->dofs_per_component_on_cell);
6925# endif
6926 this->values_dofs[dof] = val_in;
6927}
6928
6929
6930
6931template <typename Number, bool is_face, typename VectorizedArrayType>
6932inline DEAL_II_ALWAYS_INLINE void
6934 const VectorizedArrayType val_in,
6935 const unsigned int q_point)
6936{
6937# ifdef DEBUG
6938 Assert(this->is_reinitialized, ExcNotInitialized());
6939# endif
6940 AssertIndexRange(q_point, this->n_quadrature_points);
6941# ifdef DEBUG
6942 this->values_quad_submitted = true;
6943# endif
6944
6945 if (this->cell_type == internal::MatrixFreeFunctions::general)
6946 {
6947 const VectorizedArrayType JxW = this->J_value[q_point];
6948 this->values_quad[q_point] = val_in * JxW;
6949 }
6950 else // if (this->cell_type == internal::MatrixFreeFunctions::general)
6951 {
6952 const VectorizedArrayType JxW =
6953 this->J_value[0] * this->quadrature_weights[q_point];
6954 this->values_quad[q_point] = val_in * JxW;
6955 }
6956}
6957
6958
6959
6960template <typename Number, bool is_face, typename VectorizedArrayType>
6961inline DEAL_II_ALWAYS_INLINE void
6964 const unsigned int q_point)
6965{
6966 submit_value(val_in[0], q_point);
6967}
6968
6969
6970
6971template <typename Number, bool is_face, typename VectorizedArrayType>
6972inline DEAL_II_ALWAYS_INLINE void
6975 const unsigned int q_point)
6976{
6977 submit_gradient(grad_in[0], q_point);
6978}
6979
6980
6981
6982template <typename Number, bool is_face, typename VectorizedArrayType>
6983inline DEAL_II_ALWAYS_INLINE void
6985 const VectorizedArrayType grad_in,
6986 const unsigned int q_point)
6987{
6988# ifdef DEBUG
6989 Assert(this->is_reinitialized, ExcNotInitialized());
6990# endif
6991 AssertIndexRange(q_point, this->n_quadrature_points);
6992# ifdef DEBUG
6993 this->gradients_quad_submitted = true;
6994# endif
6995
6997 this->cell_type == internal::MatrixFreeFunctions::general ?
6998 this->jacobian[q_point] :
6999 this->jacobian[0];
7000 const VectorizedArrayType JxW =
7001 this->cell_type == internal::MatrixFreeFunctions::general ?
7002 this->J_value[q_point] :
7003 this->J_value[0] * this->quadrature_weights[q_point];
7004
7005 this->gradients_quad[q_point] = jac[0][0] * grad_in * JxW;
7006}
7007
7008
7009
7010template <typename Number, bool is_face, typename VectorizedArrayType>
7011inline DEAL_II_ALWAYS_INLINE void
7014 const unsigned int q_point)
7015{
7016 submit_gradient(grad_in[0][0], q_point);
7017}
7018
7019
7020
7021template <typename Number, bool is_face, typename VectorizedArrayType>
7022inline DEAL_II_ALWAYS_INLINE void
7024 submit_normal_derivative(const VectorizedArrayType grad_in,
7025 const unsigned int q_point)
7026{
7028 grad[0] = grad_in;
7029 BaseClass::submit_normal_derivative(grad, q_point);
7030}
7031
7032
7033
7034template <typename Number, bool is_face, typename VectorizedArrayType>
7035inline DEAL_II_ALWAYS_INLINE void
7038 const unsigned int q_point)
7039{
7040 BaseClass::submit_normal_derivative(grad_in, q_point);
7041}
7042
7043
7044template <typename Number, bool is_face, typename VectorizedArrayType>
7045inline DEAL_II_ALWAYS_INLINE void
7047 const Tensor<2, 1, VectorizedArrayType> hessian_in,
7048 const unsigned int q_point)
7049{
7051 hessian[0] = hessian_in;
7052 BaseClass::submit_hessian(hessian, q_point);
7053}
7054
7055
7056template <typename Number, bool is_face, typename VectorizedArrayType>
7057inline VectorizedArrayType
7059 integrate_value() const
7060{
7061 return BaseClass::integrate_value()[0];
7062}
7063
7064
7065
7066/*-------------------------- FEEvaluation -----------------------------------*/
7067
7068
7069template <int dim,
7070 int fe_degree,
7071 int n_q_points_1d,
7072 int n_components_,
7073 typename Number,
7074 typename VectorizedArrayType>
7075inline FEEvaluation<dim,
7076 fe_degree,
7077 n_q_points_1d,
7078 n_components_,
7079 Number,
7080 VectorizedArrayType>::
7081 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7082 const unsigned int fe_no,
7083 const unsigned int quad_no,
7084 const unsigned int first_selected_component,
7085 const unsigned int active_fe_index,
7086 const unsigned int active_quad_index)
7087 : BaseClass(matrix_free,
7088 fe_no,
7089 first_selected_component,
7090 quad_no,
7091 fe_degree,
7092 static_n_q_points,
7093 true /*note: this is not a face*/,
7094 active_fe_index,
7095 active_quad_index)
7096 , dofs_per_component(this->data->dofs_per_component_on_cell)
7097 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7098 , n_q_points(this->data->n_q_points)
7099{
7100 check_template_arguments(fe_no, 0);
7101}
7102
7103
7104
7105template <int dim,
7106 int fe_degree,
7107 int n_q_points_1d,
7108 int n_components_,
7109 typename Number,
7110 typename VectorizedArrayType>
7111inline FEEvaluation<dim,
7112 fe_degree,
7113 n_q_points_1d,
7114 n_components_,
7115 Number,
7116 VectorizedArrayType>::
7117 FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
7118 const std::pair<unsigned int, unsigned int> & range,
7119 const unsigned int dof_no,
7120 const unsigned int quad_no,
7121 const unsigned int first_selected_component)
7122 : FEEvaluation(matrix_free,
7123 dof_no,
7124 quad_no,
7125 first_selected_component,
7126 matrix_free.get_cell_active_fe_index(range))
7127{}
7128
7129
7130
7131template <int dim,
7132 int fe_degree,
7133 int n_q_points_1d,
7134 int n_components_,
7135 typename Number,
7136 typename VectorizedArrayType>
7137inline FEEvaluation<dim,
7138 fe_degree,
7139 n_q_points_1d,
7140 n_components_,
7141 Number,
7142 VectorizedArrayType>::
7143 FEEvaluation(const Mapping<dim> & mapping,
7144 const FiniteElement<dim> &fe,
7145 const Quadrature<1> & quadrature,
7146 const UpdateFlags update_flags,
7147 const unsigned int first_selected_component)
7148 : BaseClass(mapping,
7149 fe,
7150 quadrature,
7151 update_flags,
7152 first_selected_component,
7153 nullptr)
7154 , dofs_per_component(this->data->dofs_per_component_on_cell)
7155 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7156 , n_q_points(this->data->n_q_points)
7157{
7158 check_template_arguments(numbers::invalid_unsigned_int, 0);
7159}
7160
7161
7162
7163template <int dim,
7164 int fe_degree,
7165 int n_q_points_1d,
7166 int n_components_,
7167 typename Number,
7168 typename VectorizedArrayType>
7169inline FEEvaluation<dim,
7170 fe_degree,
7171 n_q_points_1d,
7172 n_components_,
7173 Number,
7174 VectorizedArrayType>::
7175 FEEvaluation(const FiniteElement<dim> &fe,
7176 const Quadrature<1> & quadrature,
7177 const UpdateFlags update_flags,
7178 const unsigned int first_selected_component)
7179 : BaseClass(StaticMappingQ1<dim>::mapping,
7180 fe,
7181 quadrature,
7182 update_flags,
7183 first_selected_component,
7184 nullptr)
7185 , dofs_per_component(this->data->dofs_per_component_on_cell)
7186 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7187 , n_q_points(this->data->n_q_points)
7188{
7189 check_template_arguments(numbers::invalid_unsigned_int, 0);
7190}
7191
7192
7193
7194template <int dim,
7195 int fe_degree,
7196 int n_q_points_1d,
7197 int n_components_,
7198 typename Number,
7199 typename VectorizedArrayType>
7200inline FEEvaluation<dim,
7201 fe_degree,
7202 n_q_points_1d,
7203 n_components_,
7204 Number,
7205 VectorizedArrayType>::
7206 FEEvaluation(const FiniteElement<dim> & fe,
7208 const unsigned int first_selected_component)
7209 : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
7210 fe,
7211 other.mapped_geometry->get_quadrature(),
7212 other.mapped_geometry->get_fe_values().get_update_flags(),
7213 first_selected_component,
7214 &other)
7215 , dofs_per_component(this->data->dofs_per_component_on_cell)
7216 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7217 , n_q_points(this->data->n_q_points)
7218{
7219 check_template_arguments(numbers::invalid_unsigned_int, 0);
7220}
7221
7222
7223
7224template <int dim,
7225 int fe_degree,
7226 int n_q_points_1d,
7227 int n_components_,
7228 typename Number,
7229 typename VectorizedArrayType>
7230inline FEEvaluation<dim,
7231 fe_degree,
7232 n_q_points_1d,
7233 n_components_,
7234 Number,
7235 VectorizedArrayType>::FEEvaluation(const FEEvaluation
7236 &other)
7237 : BaseClass(other)
7238 , dofs_per_component(this->data->dofs_per_component_on_cell)
7239 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
7240 , n_q_points(this->data->n_q_points)
7241{
7242 check_template_arguments(numbers::invalid_unsigned_int, 0);
7243}
7244
7245
7246
7247template <int dim,
7248 int fe_degree,
7249 int n_q_points_1d,
7250 int n_components_,
7251 typename Number,
7252 typename VectorizedArrayType>
7253inline FEEvaluation<dim,
7254 fe_degree,
7255 n_q_points_1d,
7256 n_components_,
7257 Number,
7258 VectorizedArrayType> &
7259FEEvaluation<dim,
7260 fe_degree,
7261 n_q_points_1d,
7262 n_components_,
7263 Number,
7264 VectorizedArrayType>::operator=(const FEEvaluation &other)
7265{
7266 BaseClass::operator=(other);
7267 check_template_arguments(numbers::invalid_unsigned_int, 0);
7268 return *this;
7269}
7270
7271
7272
7273template <int dim,
7274 int fe_degree,
7275 int n_q_points_1d,
7276 int n_components_,
7277 typename Number,
7278 typename VectorizedArrayType>
7279inline void
7280FEEvaluation<dim,
7281 fe_degree,
7282 n_q_points_1d,
7283 n_components_,
7284 Number,
7285 VectorizedArrayType>::
7286 check_template_arguments(const unsigned int dof_no,
7287 const unsigned int first_selected_component)
7288{
7289 (void)dof_no;
7290 (void)first_selected_component;
7291
7292 Assert(
7293 this->data->dofs_per_component_on_cell > 0,
7294 ExcMessage(
7295 "There is nothing useful you can do with an FEEvaluation object with "
7296 "FE_Nothing, i.e., without DoFs! If you have passed to "
7297 "MatrixFree::reinit() a collection of finite elements also containing "
7298 "FE_Nothing, please check - before creating FEEvaluation - the category "
7299 "of the current range by calling either "
7300 "MatrixFree::get_cell_range_category(range) or "
7301 "MatrixFree::get_face_range_category(range). The returned category "
7302 "is the index of the active FE, which you can use to exclude "
7303 "FE_Nothing."));
7304
7305# ifdef DEBUG
7306 // print error message when the dimensions do not match. Propose a possible
7307 // fix
7308 if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
7309 static_cast<unsigned int>(fe_degree) !=
7310 this->data->data.front().fe_degree) ||
7311 n_q_points != this->n_quadrature_points)
7312 {
7313 std::string message =
7314 "-------------------------------------------------------\n";
7315 message += "Illegal arguments in constructor/wrong template arguments!\n";
7316 message += " Called --> FEEvaluation<dim,";
7317 message += Utilities::int_to_string(fe_degree) + ",";
7318 message += Utilities::int_to_string(n_q_points_1d);
7319 message += "," + Utilities::int_to_string(n_components);
7320 message += ",Number>(data";
7321 if (first_selected_component != numbers::invalid_unsigned_int)
7322 {
7323 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7324 message += Utilities::int_to_string(this->quad_no) + ", ";
7325 message += Utilities::int_to_string(first_selected_component);
7326 }
7327 message += ")\n";
7328
7329 // check whether some other vector component has the correct number of
7330 // points
7331 unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
7332 proposed_fe_comp = numbers::invalid_unsigned_int,
7333 proposed_quad_comp = numbers::invalid_unsigned_int;
7334 if (dof_no != numbers::invalid_unsigned_int)
7335 {
7336 if (static_cast<unsigned int>(fe_degree) ==
7337 this->data->data.front().fe_degree)
7338 {
7339 proposed_dof_comp = dof_no;
7340 proposed_fe_comp = first_selected_component;
7341 }
7342 else
7343 for (unsigned int no = 0; no < this->matrix_free->n_components();
7344 ++no)
7345 for (unsigned int nf = 0;
7346 nf < this->matrix_free->n_base_elements(no);
7347 ++nf)
7348 if (this->matrix_free
7349 ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
7350 .data.front()
7351 .fe_degree == static_cast<unsigned int>(fe_degree))
7352 {
7353 proposed_dof_comp = no;
7354 proposed_fe_comp = nf;
7355 break;
7356 }
7357 if (n_q_points ==
7358 this->mapping_data->descriptor[this->active_quad_index]
7359 .n_q_points)
7360 proposed_quad_comp = this->quad_no;
7361 else
7362 for (unsigned int no = 0;
7363 no < this->matrix_free->get_mapping_info().cell_data.size();
7364 ++no)
7365 if (this->matrix_free->get_mapping_info()
7366 .cell_data[no]
7367 .descriptor[this->active_quad_index]
7368 .n_q_points == n_q_points)
7369 {
7370 proposed_quad_comp = no;
7371 break;
7372 }
7373 }
7374 if (proposed_dof_comp != numbers::invalid_unsigned_int &&
7375 proposed_quad_comp != numbers::invalid_unsigned_int)
7376 {
7377 if (proposed_dof_comp != first_selected_component)
7378 message += "Wrong vector component selection:\n";
7379 else
7380 message += "Wrong quadrature formula selection:\n";
7381 message += " Did you mean FEEvaluation<dim,";
7382 message += Utilities::int_to_string(fe_degree) + ",";
7383 message += Utilities::int_to_string(n_q_points_1d);
7384 message += "," + Utilities::int_to_string(n_components);
7385 message += ",Number>(data";
7386 if (dof_no != numbers::invalid_unsigned_int)
7387 {
7388 message +=
7389 ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
7390 message += Utilities::int_to_string(proposed_quad_comp) + ", ";
7391 message += Utilities::int_to_string(proposed_fe_comp);
7392 }
7393 message += ")?\n";
7394 std::string correct_pos;
7395 if (proposed_dof_comp != dof_no)
7396 correct_pos = " ^ ";
7397 else
7398 correct_pos = " ";
7399 if (proposed_quad_comp != this->quad_no)
7400 correct_pos += " ^ ";
7401 else
7402 correct_pos += " ";
7403 if (proposed_fe_comp != first_selected_component)
7404 correct_pos += " ^\n";
7405 else
7406 correct_pos += " \n";
7407 message += " " +
7408 correct_pos;
7409 }
7410 // ok, did not find the numbers specified by the template arguments in
7411 // the given list. Suggest correct template arguments
7412 const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(
7413 std::pow(1.001 * this->n_quadrature_points, 1. / dim));
7414 message += "Wrong template arguments:\n";
7415 message += " Did you mean FEEvaluation<dim,";
7416 message +=
7417 Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
7418 message += Utilities::int_to_string(proposed_n_q_points_1d);
7419 message += "," + Utilities::int_to_string(n_components);
7420 message += ",Number>(data";
7421 if (dof_no != numbers::invalid_unsigned_int)
7422 {
7423 message += ", " + Utilities::int_to_string(dof_no) + ", ";
7424 message += Utilities::int_to_string(this->quad_no);
7425 message += ", " + Utilities::int_to_string(first_selected_component);
7426 }
7427 message += ")?\n";
7428 std::string correct_pos;
7429 if (this->data->data.front().fe_degree !=
7430 static_cast<unsigned int>(fe_degree))
7431 correct_pos = " ^";
7432 else
7433 correct_pos = " ";
7434 if (proposed_n_q_points_1d != n_q_points_1d)
7435 correct_pos += " ^\n";
7436 else
7437 correct_pos += " \n";
7438 message += " " + correct_pos;
7439
7440 Assert(static_cast<unsigned int>(fe_degree) ==
7441 this->data->data.front().fe_degree &&
7442 n_q_points == this->n_quadrature_points,
7443 ExcMessage(message));
7444 }
7445 if (dof_no != numbers::invalid_unsigned_int)
7447 n_q_points,
7448 this->mapping_data->descriptor[this->active_quad_index].n_q_points);
7449# endif
7450}
7451
7452
7453
7454template <int dim,
7455 int fe_degree,
7456 int n_q_points_1d,
7457 int n_components_,
7458 typename Number,
7459 typename VectorizedArrayType>
7460inline void
7461FEEvaluation<dim,
7462 fe_degree,
7463 n_q_points_1d,
7464 n_components_,
7465 Number,
7466 VectorizedArrayType>::reinit(const unsigned int cell_index)
7467{
7468 Assert(this->matrix_free != nullptr,
7469 ExcMessage("FEEvaluation was initialized without a matrix-free object."
7470 " Integer indexing is not possible."));
7471
7472 Assert(this->dof_info != nullptr, ExcNotInitialized());
7473 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7474 this->cell = cell_index;
7475 this->cell_type =
7476 this->matrix_free->get_mapping_info().get_cell_type(cell_index);
7477
7478 const unsigned int offsets =
7479 this->mapping_data->data_index_offsets[cell_index];
7480 this->jacobian = &this->mapping_data->jacobians[0][offsets];
7481 this->J_value = &this->mapping_data->JxW_values[offsets];
7482 if (!this->mapping_data->jacobian_gradients[0].empty())
7483 {
7484 this->jacobian_gradients =
7485 this->mapping_data->jacobian_gradients[0].data() + offsets;
7486 this->jacobian_gradients_non_inverse =
7487 this->mapping_data->jacobian_gradients_non_inverse[0].data() + offsets;
7488 }
7489
7490 if (this->matrix_free->n_active_entries_per_cell_batch(this->cell) ==
7491 VectorizedArrayType::size())
7492 {
7494 for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
7495 this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7496 }
7497 else
7498 {
7499 unsigned int i = 0;
7500 for (; i < this->matrix_free->n_active_entries_per_cell_batch(this->cell);
7501 ++i)
7502 this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
7503 for (; i < VectorizedArrayType::size(); ++i)
7504 this->cell_ids[i] = numbers::invalid_unsigned_int;
7505 }
7506
7507 if (this->mapping_data->quadrature_points.empty() == false)
7508 this->quadrature_points =
7509 &this->mapping_data->quadrature_points
7510 [this->mapping_data->quadrature_point_offsets[this->cell]];
7511
7512# ifdef DEBUG
7513 this->is_reinitialized = true;
7514 this->dof_values_initialized = false;
7515 this->values_quad_initialized = false;
7516 this->gradients_quad_initialized = false;
7517 this->hessians_quad_initialized = false;
7518# endif
7519}
7520
7521
7522
7523template <int dim,
7524 int fe_degree,
7525 int n_q_points_1d,
7526 int n_components_,
7527 typename Number,
7528 typename VectorizedArrayType>
7529inline void
7530FEEvaluation<dim,
7531 fe_degree,
7532 n_q_points_1d,
7533 n_components_,
7534 Number,
7535 VectorizedArrayType>::
7536 reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
7537{
7538 Assert(this->dof_info != nullptr, ExcNotInitialized());
7539 Assert(this->mapping_data != nullptr, ExcNotInitialized());
7540
7541 this->cell = numbers::invalid_unsigned_int;
7542 this->cell_ids = cell_ids;
7543
7544 // determine type of cell batch
7546
7547 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7548 {
7549 const unsigned int cell_index = cell_ids[v];
7550
7552 continue;
7553
7554 this->cell_type =
7555 std::max(this->cell_type,
7556 this->matrix_free->get_mapping_info().get_cell_type(
7557 cell_index / VectorizedArrayType::size()));
7558 }
7559
7560 // allocate memory for internal data storage
7561 if (this->mapped_geometry == nullptr)
7562 this->mapped_geometry =
7563 std::make_shared<internal::MatrixFreeFunctions::
7564 MappingDataOnTheFly<dim, VectorizedArrayType>>();
7565
7566 auto &mapping_storage = this->mapped_geometry->get_data_storage();
7567
7568 auto &this_jacobian_data = mapping_storage.jacobians[0];
7569 auto &this_J_value_data = mapping_storage.JxW_values;
7570 auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
7571 auto &this_jacobian_gradients_non_inverse_data =
7572 mapping_storage.jacobian_gradients_non_inverse[0];
7573 auto &this_quadrature_points_data = mapping_storage.quadrature_points;
7574
7576 {
7577 if (this->mapping_data->jacobians[0].size() > 0)
7578 this_jacobian_data.resize_fast(2);
7579
7580 if (this->mapping_data->JxW_values.size() > 0)
7581 this_J_value_data.resize_fast(1);
7582
7583 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7584 this_jacobian_gradients_data.resize_fast(1);
7585
7586 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7587 this_jacobian_gradients_non_inverse_data.resize_fast(1);
7588
7589 if (this->mapping_data->quadrature_points.size() > 0)
7590 this_quadrature_points_data.resize_fast(1);
7591 }
7592 else
7593 {
7594 if (this->mapping_data->jacobians[0].size() > 0)
7595 this_jacobian_data.resize_fast(this->n_quadrature_points);
7596
7597 if (this->mapping_data->JxW_values.size() > 0)
7598 this_J_value_data.resize_fast(this->n_quadrature_points);
7599
7600 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7601 this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
7602
7603 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7604 this_jacobian_gradients_non_inverse_data.resize_fast(
7605 this->n_quadrature_points);
7606
7607 if (this->mapping_data->quadrature_points.size() > 0)
7608 this_quadrature_points_data.resize_fast(this->n_quadrature_points);
7609 }
7610
7611 // set pointers to internal data storage
7612 this->jacobian = this_jacobian_data.data();
7613 this->J_value = this_J_value_data.data();
7614 this->jacobian_gradients = this_jacobian_gradients_data.data();
7615 this->jacobian_gradients_non_inverse =
7616 this_jacobian_gradients_non_inverse_data.data();
7617 this->quadrature_points = this_quadrature_points_data.data();
7618
7619 // fill internal data storage lane by lane
7620 for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
7621 {
7622 const unsigned int cell_index = cell_ids[v];
7623
7625 continue;
7626
7627 const unsigned int cell_batch_index =
7628 cell_index / VectorizedArrayType::size();
7629 const unsigned int offsets =
7630 this->mapping_data->data_index_offsets[cell_batch_index];
7631 const unsigned int lane = cell_index % VectorizedArrayType::size();
7632
7633 if (this->cell_type <=
7635 {
7636 // case that all cells are Cartesian or affine
7637 const unsigned int q = 0;
7638
7639 if (this->mapping_data->JxW_values.size() > 0)
7640 this_J_value_data[q][v] =
7641 this->mapping_data->JxW_values[offsets + q][lane];
7642
7643 if (this->mapping_data->jacobians[0].size() > 0)
7644 for (unsigned int q = 0; q < 2; ++q)
7645 for (unsigned int i = 0; i < dim; ++i)
7646 for (unsigned int j = 0; j < dim; ++j)
7647 this_jacobian_data[q][i][j][v] =
7648 this->mapping_data->jacobians[0][offsets + q][i][j][lane];
7649
7650 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7651 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7652 for (unsigned int j = 0; j < dim; ++j)
7653 this_jacobian_gradients_data[q][i][j][v] =
7654 this->mapping_data
7655 ->jacobian_gradients[0][offsets + q][i][j][lane];
7656
7657 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() > 0)
7658 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7659 for (unsigned int j = 0; j < dim; ++j)
7660 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
7661 this->mapping_data
7662 ->jacobian_gradients_non_inverse[0][offsets + q][i][j]
7663 [lane];
7664
7665 if (this->mapping_data->quadrature_points.size() > 0)
7666 for (unsigned int i = 0; i < dim; ++i)
7667 this_quadrature_points_data[q][i][v] =
7668 this->mapping_data->quadrature_points
7669 [this->mapping_data
7670 ->quadrature_point_offsets[cell_batch_index] +
7671 q][i][lane];
7672 }
7673 else
7674 {
7675 // general case that at least one cell is not Cartesian or affine
7676 const auto cell_type =
7677 this->matrix_free->get_mapping_info().get_cell_type(
7678 cell_batch_index);
7679
7680 for (unsigned int q = 0; q < this->n_quadrature_points; ++q)
7681 {
7682 const unsigned int q_src =
7683 (cell_type <=
7685 0 :
7686 q;
7687
7688 if (this->mapping_data->JxW_values.size() > 0)
7689 this_J_value_data[q][v] =
7690 this->mapping_data->JxW_values[offsets + q_src][lane];
7691
7692 if (this->mapping_data->jacobians[0].size() > 0)
7693 for (unsigned int i = 0; i < dim; ++i)
7694 for (unsigned int j = 0; j < dim; ++j)
7695 this_jacobian_data[q][i][j][v] =
7696 this->mapping_data
7697 ->jacobians[0][offsets + q_src][i][j][lane];
7698
7699 if (this->mapping_data->jacobian_gradients[0].size() > 0)
7700 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7701 for (unsigned int j = 0; j < dim; ++j)
7702 this_jacobian_gradients_data[q][i][j][v] =
7703 this->mapping_data
7704 ->jacobian_gradients[0][offsets + q_src][i][j][lane];
7705
7706 if (this->mapping_data->jacobian_gradients_non_inverse[0].size() >
7707 0)
7708 for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
7709 for (unsigned int j = 0; j < dim; ++j)
7710 this_jacobian_gradients_non_inverse_data[q][i][j][v] =
7711 this->mapping_data
7712 ->jacobian_gradients_non_inverse[0][offsets + q_src][i]
7713 [j][lane];
7714
7715 if (this->mapping_data->quadrature_points.size() > 0)
7716 {
7717 if (cell_type <=
7719 {
7720 // affine case: quadrature points are not available but
7721 // have to be computed from the corner point and the
7722 // Jacobian
7724 this->mapping_data->quadrature_points
7725 [this->mapping_data
7726 ->quadrature_point_offsets[cell_batch_index] +
7727 0];
7728
7730 this->mapping_data->jacobians[0][offsets + 1];
7732 for (unsigned int d = 0; d < dim; ++d)
7733 point[d] +=
7734 jac[d][d] *
7735 static_cast<Number>(
7736 this->descriptor->quadrature.point(q)[d]);
7737 else
7738 for (unsigned int d = 0; d < dim; ++d)
7739 for (unsigned int e = 0; e < dim; ++e)
7740 point[d] +=
7741 jac[d][e] *
7742 static_cast<Number>(
7743 this->descriptor->quadrature.point(q)[e]);
7744
7745 for (unsigned int i = 0; i < dim; ++i)
7746 this_quadrature_points_data[q][i][v] = point[i][lane];
7747 }
7748 else
7749 {
7750 // general case: quadrature points are available
7751 for (unsigned int i = 0; i < dim; ++i)
7752 this_quadrature_points_data[q][i][v] =
7753 this->mapping_data->quadrature_points
7754 [this->mapping_data
7755 ->quadrature_point_offsets[cell_batch_index] +
7756 q][i][lane];
7757 }
7758 }
7759 }
7760 }
7761 }
7762
7763# ifdef DEBUG
7764 this->is_reinitialized = true;
7765 this->dof_values_initialized = false;
7766 this->values_quad_initialized = false;
7767 this->gradients_quad_initialized = false;
7768 this->hessians_quad_initialized = false;
7769# endif
7770}
7771
7772
7773
7774template <int dim,
7775 int fe_degree,
7776 int n_q_points_1d,
7777 int n_components_,
7778 typename Number,
7779 typename VectorizedArrayType>
7780template <bool level_dof_access>
7781inline void
7782FEEvaluation<dim,
7783 fe_degree,
7784 n_q_points_1d,
7785 n_components_,
7786 Number,
7787 VectorizedArrayType>::
7789{
7790 Assert(this->matrix_free == nullptr,
7791 ExcMessage("Cannot use initialization from cell iterator if "
7792 "initialized from MatrixFree object. Use variant for "
7793 "on the fly computation with arguments as for FEValues "
7794 "instead"));
7795 Assert(this->mapped_geometry.get() != nullptr, ExcNotInitialized());
7796 this->mapped_geometry->reinit(
7797 static_cast<typename Triangulation<dim>::cell_iterator>(cell));
7798 this->local_dof_indices.resize(cell->get_fe().n_dofs_per_cell());
7799 if (level_dof_access)
7800 cell->get_mg_dof_indices(this->local_dof_indices);
7801 else
7802 cell->get_dof_indices(this->local_dof_indices);
7803
7804# ifdef DEBUG
7805 this->is_reinitialized = true;
7806# endif
7807}
7808
7809
7810
7811template <int dim,
7812 int fe_degree,
7813 int n_q_points_1d,
7814 int n_components_,
7815 typename Number,
7816 typename VectorizedArrayType>
7817inline void
7818FEEvaluation<dim,
7819 fe_degree,
7820 n_q_points_1d,
7821 n_components_,
7822 Number,
7823 VectorizedArrayType>::
7824 reinit(const typename Triangulation<dim>::cell_iterator &cell)
7825{
7826 Assert(this->matrix_free == 0,
7827 ExcMessage("Cannot use initialization from cell iterator if "
7828 "initialized from MatrixFree object. Use variant for "
7829 "on the fly computation with arguments as for FEValues "
7830 "instead"));
7831 Assert(this->mapped_geometry.get() != 0, ExcNotInitialized());
7832 this->mapped_geometry->reinit(cell);
7833
7834# ifdef DEBUG
7835 this->is_reinitialized = true;
7836# endif
7837}
7838
7839
7840
7841template <int dim,
7842 int fe_degree,
7843 int n_q_points_1d,
7844 int n_components_,
7845 typename Number,
7846 typename VectorizedArrayType>
7847inline void
7848FEEvaluation<dim,
7849 fe_degree,
7850 n_q_points_1d,
7851 n_components_,
7852 Number,
7853 VectorizedArrayType>::evaluate(const bool evaluate_values,
7854 const bool evaluate_gradients,
7855 const bool evaluate_hessians)
7856{
7857# ifdef DEBUG
7858 Assert(this->dof_values_initialized == true,
7860# endif
7861 evaluate(this->values_dofs,
7862 evaluate_values,
7863 evaluate_gradients,
7864 evaluate_hessians);
7865}
7866
7867
7868template <int dim,
7869 int fe_degree,
7870 int n_q_points_1d,
7871 int n_components_,
7872 typename Number,
7873 typename VectorizedArrayType>
7874inline void
7875FEEvaluation<dim,
7876 fe_degree,
7877 n_q_points_1d,
7878 n_components_,
7879 Number,
7880 VectorizedArrayType>::
7881 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flags)
7882{
7883# ifdef DEBUG
7884 Assert(this->dof_values_initialized == true,
7886# endif
7887 evaluate(this->values_dofs, evaluation_flags);
7888}
7889
7890
7891
7892template <int dim,
7893 int fe_degree,
7894 int n_q_points_1d,
7895 int n_components_,
7896 typename Number,
7897 typename VectorizedArrayType>
7898inline void
7899FEEvaluation<dim,
7900 fe_degree,
7901 n_q_points_1d,
7902 n_components_,
7903 Number,
7904 VectorizedArrayType>::evaluate(const VectorizedArrayType
7905 * values_array,
7906 const bool evaluate_values,
7907 const bool evaluate_gradients,
7908 const bool evaluate_hessians)
7909{
7911 ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7912 ((evaluate_gradients) ? EvaluationFlags::gradients :
7914 ((evaluate_hessians) ? EvaluationFlags::hessians :
7916
7917 evaluate(values_array, flag);
7918}
7919
7920
7921
7922template <int dim,
7923 int fe_degree,
7924 int n_q_points_1d,
7925 int n_components_,
7926 typename Number,
7927 typename VectorizedArrayType>
7928inline void
7929FEEvaluation<dim,
7930 fe_degree,
7931 n_q_points_1d,
7932 n_components_,
7933 Number,
7934 VectorizedArrayType>::
7935 evaluate(const VectorizedArrayType * values_array,
7936 const EvaluationFlags::EvaluationFlags evaluation_flag)
7937{
7938 const bool hessians_on_general_cells =
7939 evaluation_flag & EvaluationFlags::hessians &&
7940 (this->cell_type > internal::MatrixFreeFunctions::affine);
7941 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
7942 if (hessians_on_general_cells)
7943 evaluation_flag_actual |= EvaluationFlags::gradients;
7944
7945 if (this->data->element_type ==
7947 evaluation_flag & EvaluationFlags::gradients &&
7948 (this->cell_type > internal::MatrixFreeFunctions::affine))
7949 evaluation_flag_actual |= EvaluationFlags::values;
7950
7951 if (fe_degree > -1)
7952 {
7954 evaluate(n_components, evaluation_flag_actual, values_array, *this);
7955 }
7956 else
7957 {
7959 n_components,
7960 evaluation_flag_actual,
7961 const_cast<VectorizedArrayType *>(values_array),
7962 *this);
7963 }
7964
7965# ifdef DEBUG
7966 if (evaluation_flag_actual & EvaluationFlags::values)
7967 this->values_quad_initialized = true;
7968 if (evaluation_flag_actual & EvaluationFlags::gradients)
7969 this->gradients_quad_initialized = true;
7970 if (evaluation_flag_actual & EvaluationFlags::hessians)
7971 this->hessians_quad_initialized = true;
7972# endif
7973}
7974
7975
7976
7977template <int dim,
7978 int fe_degree,
7979 int n_q_points_1d,
7980 int n_components_,
7981 typename Number,
7982 typename VectorizedArrayType>
7983template <typename VectorType>
7984inline void
7986 dim,
7987 fe_degree,
7988 n_q_points_1d,
7989 n_components_,
7990 Number,
7991 VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
7992 const bool evaluate_values,
7993 const bool evaluate_gradients,
7994 const bool evaluate_hessians)
7995{
7997 ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
7998 ((evaluate_gradients) ? EvaluationFlags::gradients :
8000 ((evaluate_hessians) ? EvaluationFlags::hessians :
8002
8003 gather_evaluate(input_vector, flag);
8004}
8005
8006
8007namespace internal
8008{
8012 template <typename Number,
8013 typename VectorizedArrayType,
8014 typename VectorType,
8015 typename EvaluatorType,
8016 std::enable_if_t<internal::has_begin<VectorType> &&
8018 VectorType> * = nullptr>
8019 VectorizedArrayType *
8020 check_vector_access_inplace(const EvaluatorType &fe_eval, VectorType &vector)
8021 {
8022 // for user-defined cell batches this functionality is not supported
8023 if (fe_eval.get_current_cell_index() == numbers::invalid_unsigned_int)
8024 return nullptr;
8025
8026 const unsigned int cell = fe_eval.get_cell_or_face_batch_id();
8027 const auto & dof_info = fe_eval.get_dof_info();
8028
8029 // If the index storage is interleaved and contiguous and the vector
8030 // storage has the correct alignment, we can directly pass the pointer
8031 // into the vector to the evaluate() and integrate() calls, without
8032 // reading the vector entries into a separate data field. This saves some
8033 // operations.
8034 if (std::is_same<typename VectorType::value_type, Number>::value &&
8035 dof_info.index_storage_variants
8038 interleaved_contiguous &&
8039 reinterpret_cast<std::size_t>(
8040 vector.begin() +
8041 dof_info.dof_indices_contiguous
8042 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
8043 [cell * VectorizedArrayType::size()]) %
8044 sizeof(VectorizedArrayType) ==
8045 0)
8046 {
8047 return reinterpret_cast<VectorizedArrayType *>(
8048 vector.begin() +
8049 dof_info.dof_indices_contiguous
8051 [cell * VectorizedArrayType::size()] +
8053 [fe_eval.get_active_fe_index()]
8054 [fe_eval.get_first_selected_component()] *
8055 VectorizedArrayType::size());
8056 }
8057 else
8058 return nullptr;
8059 }
8060
8064 template <typename Number,
8065 typename VectorizedArrayType,
8066 typename VectorType,
8067 typename EvaluatorType,
8068 std::enable_if_t<!internal::has_begin<VectorType> ||
8070 VectorType> * = nullptr>
8071 VectorizedArrayType *
8072 check_vector_access_inplace(const EvaluatorType &, VectorType &)
8073 {
8074 return nullptr;
8075 }
8076} // namespace internal
8077
8078
8079
8080template <int dim,
8081 int fe_degree,
8082 int n_q_points_1d,
8083 int n_components_,
8084 typename Number,
8085 typename VectorizedArrayType>
8086template <typename VectorType>
8087inline void
8088FEEvaluation<dim,
8089 fe_degree,
8090 n_q_points_1d,
8091 n_components_,
8092 Number,
8093 VectorizedArrayType>::
8094 gather_evaluate(const VectorType & input_vector,
8095 const EvaluationFlags::EvaluationFlags evaluation_flag)
8096{
8097 const VectorizedArrayType *src_ptr =
8098 internal::check_vector_access_inplace<Number, const VectorizedArrayType>(
8099 *this, input_vector);
8100 if (src_ptr != nullptr)
8101 evaluate(src_ptr, evaluation_flag);
8102 else
8103 {
8104 this->read_dof_values(input_vector);
8105 evaluate(this->begin_dof_values(), evaluation_flag);
8106 }
8107}
8108
8109
8110
8111template <int dim,
8112 int fe_degree,
8113 int n_q_points_1d,
8114 int n_components_,
8115 typename Number,
8116 typename VectorizedArrayType>
8117inline void
8118FEEvaluation<dim,
8119 fe_degree,
8120 n_q_points_1d,
8121 n_components_,
8122 Number,
8123 VectorizedArrayType>::integrate(const bool integrate_values,
8124 const bool integrate_gradients)
8125{
8126 integrate(integrate_values, integrate_gradients, this->values_dofs);
8127
8128# ifdef DEBUG
8129 this->dof_values_initialized = true;
8130# endif
8131}
8132
8133
8134
8135template <int dim,
8136 int fe_degree,
8137 int n_q_points_1d,
8138 int n_components_,
8139 typename Number,
8140 typename VectorizedArrayType>
8141inline void
8142FEEvaluation<dim,
8143 fe_degree,
8144 n_q_points_1d,
8145 n_components_,
8146 Number,
8147 VectorizedArrayType>::
8148 integrate(const EvaluationFlags::EvaluationFlags integration_flag)
8149{
8150 integrate(integration_flag, this->values_dofs);
8151
8152# ifdef DEBUG
8153 this->dof_values_initialized = true;
8154# endif
8155}
8156
8157
8158
8159template <int dim,
8160 int fe_degree,
8161 int n_q_points_1d,
8162 int n_components_,
8163 typename Number,
8164 typename VectorizedArrayType>
8165inline void
8166FEEvaluation<dim,
8167 fe_degree,
8168 n_q_points_1d,
8169 n_components_,
8170 Number,
8171 VectorizedArrayType>::integrate(const bool integrate_values,
8172 const bool integrate_gradients,
8173 VectorizedArrayType *values_array)
8174{
8176 (integrate_values ? EvaluationFlags::values : EvaluationFlags::nothing) |
8177 (integrate_gradients ? EvaluationFlags::gradients :
8179 integrate(flag, values_array);
8180}
8181
8182
8183
8184template <int dim,
8185 int fe_degree,
8186 int n_q_points_1d,
8187 int n_components_,
8188 typename Number,
8189 typename VectorizedArrayType>
8190inline void
8191FEEvaluation<dim,
8192 fe_degree,
8193 n_q_points_1d,
8194 n_components_,
8195 Number,
8196 VectorizedArrayType>::
8197 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
8198 VectorizedArrayType * values_array,
8199 const bool sum_into_values_array)
8200{
8201# ifdef DEBUG
8202 if (integration_flag & EvaluationFlags::values)
8203 Assert(this->values_quad_submitted == true,
8205 if (integration_flag & EvaluationFlags::gradients)
8206 Assert(this->gradients_quad_submitted == true,
8208 if ((integration_flag & EvaluationFlags::hessians) != 0u)
8209 Assert(this->hessians_quad_submitted == true,
8211# endif
8212 Assert(this->matrix_free != nullptr ||
8213 this->mapped_geometry->is_initialized(),
8215
8216 Assert(
8217 (integration_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients |
8219 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, and "
8220 "EvaluationFlags::hessians are supported."));
8221
8222 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
8223 if (integration_flag & EvaluationFlags::hessians &&
8224 (this->cell_type > internal::MatrixFreeFunctions::affine))
8225 {
8226 unsigned int size = n_components * dim * n_q_points;
8227 if ((integration_flag & EvaluationFlags::gradients) != 0u)
8228 {
8229 for (unsigned int i = 0; i < size; ++i)
8230 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
8231 }
8232 else
8233 {
8234 for (unsigned int i = 0; i < size; ++i)
8235 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
8236 integration_flag_actual |= EvaluationFlags::gradients;
8237 }
8238 }
8239
8240 if (this->data->element_type ==
8242 integration_flag & EvaluationFlags::gradients &&
8243 this->cell_type > internal::MatrixFreeFunctions::affine &&
8244 this->divergence_is_requested == false)
8245 {
8246 unsigned int size = n_components * n_q_points;
8247 if ((integration_flag & EvaluationFlags::values) != 0u)
8248 {
8249 for (unsigned int i = 0; i < size; ++i)
8250 this->values_quad[i] += this->values_from_gradients_quad[i];
8251 }
8252 else
8253 {
8254 for (unsigned int i = 0; i < size; ++i)
8255 this->values_quad[i] = this->values_from_gradients_quad[i];
8256 integration_flag_actual |= EvaluationFlags::values;
8257 }
8258 }
8259
8260 if (fe_degree > -1)
8261 {
8263 integrate(n_components,
8264 integration_flag_actual,
8265 values_array,
8266 *this,
8267 sum_into_values_array);
8268 }
8269 else
8270 {
8272 n_components,
8273 integration_flag_actual,
8274 values_array,
8275 *this,
8276 sum_into_values_array);
8277 }
8278
8279# ifdef DEBUG
8280 this->dof_values_initialized = true;
8281# endif
8282}
8283
8284
8285
8286template <int dim,
8287 int fe_degree,
8288 int n_q_points_1d,
8289 int n_components_,
8290 typename Number,
8291 typename VectorizedArrayType>
8292template <typename VectorType>
8293inline void
8295 dim,
8296 fe_degree,
8297 n_q_points_1d,
8298 n_components_,
8299 Number,
8300 VectorizedArrayType>::integrate_scatter(const bool integrate_values,
8301 const bool integrate_gradients,
8302 VectorType &destination)
8303{
8305 ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8306 ((integrate_gradients) ? EvaluationFlags::gradients :
8308
8309 integrate_scatter(flag, destination);
8310}
8311
8312
8313
8314template <int dim,
8315 int fe_degree,
8316 int n_q_points_1d,
8317 int n_components_,
8318 typename Number,
8319 typename VectorizedArrayType>
8320template <typename VectorType>
8321inline void
8322FEEvaluation<dim,
8323 fe_degree,
8324 n_q_points_1d,
8325 n_components_,
8326 Number,
8327 VectorizedArrayType>::
8328 integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
8329 VectorType & destination)
8330{
8331 VectorizedArrayType *dst_ptr =
8332 internal::check_vector_access_inplace<Number, VectorizedArrayType>(
8333 *this, destination);
8334 if (dst_ptr != nullptr)
8335 integrate(integration_flag, dst_ptr, true);
8336 else
8337 {
8338 integrate(integration_flag, this->begin_dof_values());
8339 this->distribute_local_to_global(destination);
8340 }
8341}
8342
8343
8344
8345template <int dim,
8346 int fe_degree,
8347 int n_q_points_1d,
8348 int n_components_,
8349 typename Number,
8350 typename VectorizedArrayType>
8352FEEvaluation<dim,
8353 fe_degree,
8354 n_q_points_1d,
8355 n_components_,
8356 Number,
8357 VectorizedArrayType>::dof_indices() const
8358{
8359 return {0U, dofs_per_cell};
8360}
8361
8362
8363
8364/*-------------------------- FEFaceEvaluation ---------------------------*/
8365
8366
8367
8368template <int dim,
8369 int fe_degree,
8370 int n_q_points_1d,
8371 int n_components_,
8372 typename Number,
8373 typename VectorizedArrayType>
8374inline FEFaceEvaluation<dim,
8375 fe_degree,
8376 n_q_points_1d,
8377 n_components_,
8378 Number,
8379 VectorizedArrayType>::
8380 FEFaceEvaluation(
8382 const bool is_interior_face,
8383 const unsigned int dof_no,
8384 const unsigned int quad_no,
8385 const unsigned int first_selected_component,
8386 const unsigned int active_fe_index,
8387 const unsigned int active_quad_index,
8388 const unsigned int face_type)
8389 : BaseClass(matrix_free,
8390 dof_no,
8391 first_selected_component,
8392 quad_no,
8393 fe_degree,
8394 static_n_q_points,
8395 is_interior_face,
8396 active_fe_index,
8397 active_quad_index,
8398 face_type)
8399 , dofs_per_component(this->data->dofs_per_component_on_cell)
8400 , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
8401 , n_q_points(this->n_quadrature_points)
8402{}
8403
8404
8405
8406template <int dim,
8407 int fe_degree,
8408 int n_q_points_1d,
8409 int n_components_,
8410 typename Number,
8411 typename VectorizedArrayType>
8412inline FEFaceEvaluation<dim,
8413 fe_degree,
8414 n_q_points_1d,
8415 n_components_,
8416 Number,
8417 VectorizedArrayType>::
8418 FEFaceEvaluation(
8420 const std::pair<unsigned int, unsigned int> & range,
8421 const bool is_interior_face,
8422 const unsigned int dof_no,
8423 const unsigned int quad_no,
8424 const unsigned int first_selected_component)
8425 : FEFaceEvaluation(matrix_free,
8426 is_interior_face,
8427 dof_no,
8428 quad_no,
8429 first_selected_component,
8430 matrix_free.get_face_active_fe_index(range,
8431 is_interior_face),
8432 numbers::invalid_unsigned_int,
8433 matrix_free.get_face_info(range.first).face_type)
8434{}
8435
8436
8437
8438template <int dim,
8439 int fe_degree,
8440 int n_q_points_1d,
8441 int n_components_,
8442 typename Number,
8443 typename VectorizedArrayType>
8444inline void
8446 fe_degree,
8447 n_q_points_1d,
8448 n_components_,
8449 Number,
8450 VectorizedArrayType>::reinit(const unsigned int face_index)
8451{
8452 Assert(this->mapped_geometry == nullptr,
8453 ExcMessage("FEEvaluation was initialized without a matrix-free object."
8454 " Integer indexing is not possible"));
8455 if (this->mapped_geometry != nullptr)
8456 return;
8457
8458 this->cell = face_index;
8459 this->dof_access_index =
8460 this->is_interior_face() ?
8463 Assert(this->mapping_data != nullptr, ExcNotInitialized());
8464
8465 if (face_index >=
8466 this->matrix_free->get_task_info().face_partition_data.back() &&
8467 face_index <
8468 this->matrix_free->get_task_info().boundary_partition_data.back())
8469 Assert(this->is_interior_face(),
8470 ExcMessage(
8471 "Boundary faces do not have a neighbor. When looping over "
8472 "boundary faces use FEFaceEvaluation with the parameter "
8473 "is_interior_face set to true. "));
8474
8475 this->reinit_face(this->matrix_free->get_face_info(face_index));
8476
8477 unsigned int i = 0;
8478 for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
8479 ++i)
8480 this->face_ids[i] = face_index * VectorizedArrayType::size() + i;
8481 for (; i < VectorizedArrayType::size(); ++i)
8482 this->face_ids[i] = numbers::invalid_unsigned_int;
8483
8484 this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
8485 const unsigned int offsets =
8486 this->mapping_data->data_index_offsets[face_index];
8487 this->J_value = &this->mapping_data->JxW_values[offsets];
8488 this->normal_vectors = &this->mapping_data->normal_vectors[offsets];
8489 this->jacobian =
8490 &this->mapping_data->jacobians[!this->is_interior_face()][offsets];
8491 this->normal_x_jacobian =
8492 &this->mapping_data
8493 ->normals_times_jacobians[!this->is_interior_face()][offsets];
8494 this->jacobian_gradients =
8495 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8496 offsets;
8497 this->jacobian_gradients_non_inverse =
8498 this->mapping_data
8499 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
8500 .data() +
8501 offsets;
8502
8503 if (this->mapping_data->quadrature_point_offsets.empty() == false)
8504 {
8505 AssertIndexRange(this->cell,
8506 this->mapping_data->quadrature_point_offsets.size());
8507 this->quadrature_points =
8508 this->mapping_data->quadrature_points.data() +
8509 this->mapping_data->quadrature_point_offsets[this->cell];
8510 }
8511
8512# ifdef DEBUG
8513 this->is_reinitialized = true;
8514 this->dof_values_initialized = false;
8515 this->values_quad_initialized = false;
8516 this->gradients_quad_initialized = false;
8517 this->hessians_quad_initialized = false;
8518# endif
8519}
8520
8521
8522
8523template <int dim,
8524 int fe_degree,
8525 int n_q_points_1d,
8526 int n_components_,
8527 typename Number,
8528 typename VectorizedArrayType>
8529inline void
8531 fe_degree,
8532 n_q_points_1d,
8533 n_components_,
8534 Number,
8535 VectorizedArrayType>::reinit(const unsigned int cell_index,
8536 const unsigned int face_number)
8537{
8538 Assert(
8539 this->quad_no <
8540 this->matrix_free->get_mapping_info().face_data_by_cells.size(),
8541 ExcMessage(
8542 "You must set MatrixFree::AdditionalData::mapping_update_flags_faces_by_cells to use the present reinit method."));
8545 this->matrix_free->get_mapping_info().cell_type.size());
8546 Assert(this->mapped_geometry == nullptr,
8547 ExcMessage("FEEvaluation was initialized without a matrix-free object."
8548 " Integer indexing is not possible"));
8549 if (this->mapped_geometry != nullptr)
8550 return;
8551 Assert(this->matrix_free != nullptr, ExcNotInitialized());
8552
8553 this->cell_type = this->matrix_free->get_mapping_info()
8554 .faces_by_cells_type[cell_index][face_number];
8555 this->cell = cell_index;
8556 this->subface_index = GeometryInfo<dim>::max_children_per_cell;
8557 this->dof_access_index =
8559
8560 constexpr unsigned int n_lanes = VectorizedArrayType::size();
8561
8562 if (this->is_interior_face() == false)
8563 {
8564 // for this case, we need to look into the FaceInfo field that collects
8565 // information from both sides of a face once for the global mesh, and
8566 // pick the face id that is not the local one (cell_this).
8567 for (unsigned int i = 0; i < n_lanes; ++i)
8568 {
8569 // compute actual (non vectorized) cell ID
8570 const unsigned int cell_this = cell_index * n_lanes + i;
8571 // compute face ID
8572 unsigned int face_index =
8573 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
8574 face_number,
8575 i);
8576
8577 this->face_ids[i] = face_index;
8578
8579 if (face_index == numbers::invalid_unsigned_int)
8580 {
8581 this->cell_ids[i] = numbers::invalid_unsigned_int;
8582 this->face_numbers[i] = static_cast<std::uint8_t>(-1);
8583 this->face_orientations[i] = static_cast<std::uint8_t>(-1);
8584 continue; // invalid face ID: no neighbor on boundary
8585 }
8586
8587 const auto &faces =
8588 this->matrix_free->get_face_info(face_index / n_lanes);
8589 // get cell ID on both sides of face
8590 auto cell_m = faces.cells_interior[face_index % n_lanes];
8591 auto cell_p = faces.cells_exterior[face_index % n_lanes];
8592
8593 const bool face_identifies_as_interior = cell_m != cell_this;
8594
8595 Assert(cell_m == cell_this || cell_p == cell_this,
8597
8598 // compare the IDs with the given cell ID
8599 if (face_identifies_as_interior)
8600 {
8601 this->cell_ids[i] = cell_m; // neighbor has the other ID
8602 this->face_numbers[i] = faces.interior_face_no;
8603 }
8604 else
8605 {
8606 this->cell_ids[i] = cell_p;
8607 this->face_numbers[i] = faces.exterior_face_no;
8608 }
8609
8610 const bool orientation_interior_face = faces.face_orientation >= 8;
8611 unsigned int face_orientation = faces.face_orientation % 8;
8612 if (face_identifies_as_interior != orientation_interior_face)
8613 {
8614 constexpr std::array<std::uint8_t, 8> table{
8615 {0, 1, 2, 3, 6, 5, 4, 7}};
8616 face_orientation = table[face_orientation];
8617 }
8618 this->face_orientations[i] = face_orientation;
8619 }
8620 }
8621 else
8622 {
8623 this->face_orientations[0] = 0;
8624 this->face_numbers[0] = face_number;
8625 for (unsigned int i = 0; i < n_lanes; ++i)
8626 this->cell_ids[i] = cell_index * n_lanes + i;
8627 for (unsigned int i = 0; i < n_lanes; ++i)
8628 this->face_ids[i] =
8629 this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
8630 face_number,
8631 i);
8632 }
8633
8634 const unsigned int offsets =
8635 this->matrix_free->get_mapping_info()
8636 .face_data_by_cells[this->quad_no]
8637 .data_index_offsets[cell_index * GeometryInfo<dim>::faces_per_cell +
8638 face_number];
8639 AssertIndexRange(offsets,
8640 this->matrix_free->get_mapping_info()
8641 .face_data_by_cells[this->quad_no]
8642 .JxW_values.size());
8643 this->J_value = &this->matrix_free->get_mapping_info()
8644 .face_data_by_cells[this->quad_no]
8645 .JxW_values[offsets];
8646 this->normal_vectors = &this->matrix_free->get_mapping_info()
8647 .face_data_by_cells[this->quad_no]
8648 .normal_vectors[offsets];
8649 this->jacobian = &this->matrix_free->get_mapping_info()
8650 .face_data_by_cells[this->quad_no]
8651 .jacobians[!this->is_interior_face()][offsets];
8652 this->normal_x_jacobian =
8653 &this->matrix_free->get_mapping_info()
8654 .face_data_by_cells[this->quad_no]
8655 .normals_times_jacobians[!this->is_interior_face()][offsets];
8656 this->jacobian_gradients =
8657 this->mapping_data->jacobian_gradients[!this->is_interior_face()].data() +
8658 offsets;
8659 this->jacobian_gradients_non_inverse =
8660 this->mapping_data
8661 ->jacobian_gradients_non_inverse[!this->is_interior_face()]
8662 .data() +
8663 offsets;
8664
8665 if (this->matrix_free->get_mapping_info()
8666 .face_data_by_cells[this->quad_no]
8667 .quadrature_point_offsets.empty() == false)
8668 {
8669 const unsigned int index =
8670 this->cell * GeometryInfo<dim>::faces_per_cell + this->face_numbers[0];
8672 this->matrix_free->get_mapping_info()
8673 .face_data_by_cells[this->quad_no]
8674 .quadrature_point_offsets.size());
8675 this->quadrature_points = this->matrix_free->get_mapping_info()
8676 .face_data_by_cells[this->quad_no]
8677 .quadrature_points.data() +
8678 this->matrix_free->get_mapping_info()
8679 .face_data_by_cells[this->quad_no]
8680 .quadrature_point_offsets[index];
8681 }
8682
8683# ifdef DEBUG
8684 this->is_reinitialized = true;
8685 this->dof_values_initialized = false;
8686 this->values_quad_initialized = false;
8687 this->gradients_quad_initialized = false;
8688 this->hessians_quad_initialized = false;
8689# endif
8690}
8691
8692
8693
8694template <int dim,
8695 int fe_degree,
8696 int n_q_points_1d,
8697 int n_components_,
8698 typename Number,
8699 typename VectorizedArrayType>
8700inline void
8702 fe_degree,
8703 n_q_points_1d,
8704 n_components_,
8705 Number,
8706 VectorizedArrayType>::evaluate(const bool evaluate_values,
8707 const bool evaluate_gradients)
8708{
8709# ifdef DEBUG
8710 Assert(this->dof_values_initialized, ExcNotInitialized());
8711# endif
8712
8713 evaluate(this->values_dofs, evaluate_values, evaluate_gradients);
8714}
8715
8716
8717
8718template <int dim,
8719 int fe_degree,
8720 int n_q_points_1d,
8721 int n_components_,
8722 typename Number,
8723 typename VectorizedArrayType>
8724inline void
8726 fe_degree,
8727 n_q_points_1d,
8728 n_components_,
8729 Number,
8730 VectorizedArrayType>::
8731 evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
8732{
8733# ifdef DEBUG
8734 Assert(this->dof_values_initialized, ExcNotInitialized());
8735# endif
8736
8737 evaluate(this->values_dofs, evaluation_flag);
8738}
8739
8740
8741
8742template <int dim,
8743 int fe_degree,
8744 int n_q_points_1d,
8745 int n_components_,
8746 typename Number,
8747 typename VectorizedArrayType>
8748inline void
8750 fe_degree,
8751 n_q_points_1d,
8752 n_components_,
8753 Number,
8754 VectorizedArrayType>::evaluate(const VectorizedArrayType
8755 * values_array,
8756 const bool evaluate_values,
8757 const bool evaluate_gradients)
8758{
8760 ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8761 ((evaluate_gradients) ? EvaluationFlags::gradients :
8763
8764 evaluate(values_array, flag);
8765}
8766
8767
8768
8769template <int dim,
8770 int fe_degree,
8771 int n_q_points_1d,
8772 int n_components_,
8773 typename Number,
8774 typename VectorizedArrayType>
8775inline void
8777 fe_degree,
8778 n_q_points_1d,
8779 n_components_,
8780 Number,
8781 VectorizedArrayType>::
8782 evaluate(const VectorizedArrayType * values_array,
8783 const EvaluationFlags::EvaluationFlags evaluation_flag)
8784{
8785 Assert((evaluation_flag &
8788 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8789 "and EvaluationFlags::hessians are supported."));
8790
8791 const bool hessians_on_general_cells =
8792 evaluation_flag & EvaluationFlags::hessians &&
8793 (this->cell_type > internal::MatrixFreeFunctions::affine);
8794 EvaluationFlags::EvaluationFlags evaluation_flag_actual = evaluation_flag;
8795 if (hessians_on_general_cells)
8796 evaluation_flag_actual |= EvaluationFlags::gradients;
8797
8798 if (this->data->element_type ==
8800 evaluation_flag & EvaluationFlags::gradients &&
8801 (this->cell_type > internal::MatrixFreeFunctions::affine))
8802 evaluation_flag_actual |= EvaluationFlags::values;
8803
8804 if (fe_degree > -1)
8806 template run<fe_degree, n_q_points_1d>(n_components,
8807 evaluation_flag_actual,
8808 values_array,
8809 *this);
8810 else
8812 n_components, evaluation_flag_actual, values_array, *this);
8813
8814# ifdef DEBUG
8815 if (evaluation_flag_actual & EvaluationFlags::values)
8816 this->values_quad_initialized = true;
8817 if (evaluation_flag_actual & EvaluationFlags::gradients)
8818 this->gradients_quad_initialized = true;
8819 if ((evaluation_flag_actual & EvaluationFlags::hessians) != 0u)
8820 this->hessians_quad_initialized = true;
8821# endif
8822}
8823
8824
8825
8826template <int dim,
8827 int fe_degree,
8828 int n_q_points_1d,
8829 int n_components_,
8830 typename Number,
8831 typename VectorizedArrayType>
8832inline void
8834 fe_degree,
8835 n_q_points_1d,
8836 n_components_,
8837 Number,
8838 VectorizedArrayType>::
8839 integrate(const EvaluationFlags::EvaluationFlags integration_flag)
8840{
8841 integrate(integration_flag, this->values_dofs);
8842
8843# ifdef DEBUG
8844 this->dof_values_initialized = true;
8845# endif
8846}
8847
8848
8849
8850template <int dim,
8851 int fe_degree,
8852 int n_q_points_1d,
8853 int n_components_,
8854 typename Number,
8855 typename VectorizedArrayType>
8856inline void
8858 fe_degree,
8859 n_q_points_1d,
8860 n_components_,
8861 Number,
8862 VectorizedArrayType>::integrate(const bool integrate_values,
8863 const bool integrate_gradients)
8864{
8865 integrate(integrate_values, integrate_gradients, this->values_dofs);
8866
8867# ifdef DEBUG
8868 this->dof_values_initialized = true;
8869# endif
8870}
8871
8872
8873
8874template <int dim,
8875 int fe_degree,
8876 int n_q_points_1d,
8877 int n_components_,
8878 typename Number,
8879 typename VectorizedArrayType>
8880inline void
8882 fe_degree,
8883 n_q_points_1d,
8884 n_components_,
8885 Number,
8886 VectorizedArrayType>::integrate(const bool integrate_values,
8887 const bool integrate_gradients,
8888 VectorizedArrayType
8889 *values_array)
8890{
8892 ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8893 ((integrate_gradients) ? EvaluationFlags::gradients :
8895
8896 integrate(flag, values_array);
8897}
8898
8899
8900
8901template <int dim,
8902 int fe_degree,
8903 int n_q_points_1d,
8904 int n_components_,
8905 typename Number,
8906 typename VectorizedArrayType>
8907inline void
8909 fe_degree,
8910 n_q_points_1d,
8911 n_components_,
8912 Number,
8913 VectorizedArrayType>::
8914 integrate(const EvaluationFlags::EvaluationFlags integration_flag,
8915 VectorizedArrayType * values_array)
8916{
8917 Assert((integration_flag &
8920 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
8921 "and EvaluationFlags::hessians are supported."));
8922
8923 EvaluationFlags::EvaluationFlags integration_flag_actual = integration_flag;
8924 if (integration_flag & EvaluationFlags::hessians &&
8925 (this->cell_type > internal::MatrixFreeFunctions::affine))
8926 {
8927 unsigned int size = n_components * dim * n_q_points;
8928 if ((integration_flag & EvaluationFlags::gradients) != 0u)
8929 {
8930 for (unsigned int i = 0; i < size; ++i)
8931 this->gradients_quad[i] += this->gradients_from_hessians_quad[i];
8932 }
8933 else
8934 {
8935 for (unsigned int i = 0; i < size; ++i)
8936 this->gradients_quad[i] = this->gradients_from_hessians_quad[i];
8937 integration_flag_actual |= EvaluationFlags::gradients;
8938 }
8939 }
8940
8941 if (this->data->element_type ==
8943 integration_flag & EvaluationFlags::gradients &&
8944 this->cell_type > internal::MatrixFreeFunctions::affine &&
8945 this->divergence_is_requested == false)
8946 {
8947 unsigned int size = n_components * n_q_points;
8948 if ((integration_flag & EvaluationFlags::values) != 0u)
8949 {
8950 for (unsigned int i = 0; i < size; ++i)
8951 this->values_quad[i] += this->values_from_gradients_quad[i];
8952 }
8953 else
8954 {
8955 for (unsigned int i = 0; i < size; ++i)
8956 this->values_quad[i] = this->values_from_gradients_quad[i];
8957 integration_flag_actual |= EvaluationFlags::values;
8958 }
8959 }
8960
8961 if (fe_degree > -1)
8963 template run<fe_degree, n_q_points_1d>(n_components,
8964 integration_flag_actual,
8965 values_array,
8966 *this);
8967 else
8969 n_components, integration_flag_actual, values_array, *this);
8970}
8971
8972
8973
8974template <int dim,
8975 int fe_degree,
8976 int n_q_points_1d,
8977 int n_components_,
8978 typename Number,
8979 typename VectorizedArrayType>
8980template <typename VectorType>
8981inline void
8983 dim,
8984 fe_degree,
8985 n_q_points_1d,
8986 n_components_,
8987 Number,
8988 VectorizedArrayType>::gather_evaluate(const VectorType &input_vector,
8989 const bool evaluate_values,
8990 const bool evaluate_gradients)
8991{
8993 ((evaluate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
8994 ((evaluate_gradients) ? EvaluationFlags::gradients :
8996
8997 gather_evaluate(input_vector, flag);
8998}
8999
9000
9001
9002template <int dim,
9003 int fe_degree,
9004 int n_q_points_1d,
9005 int n_components_,
9006 typename Number,
9007 typename VectorizedArrayType>
9008template <typename VectorType>
9009inline void
9011 fe_degree,
9012 n_q_points_1d,
9013 n_components_,
9014 Number,
9015 VectorizedArrayType>::
9016 gather_evaluate(const VectorType & input_vector,
9017 const EvaluationFlags::EvaluationFlags evaluation_flag)
9018{
9019 Assert((evaluation_flag &
9022 ExcMessage("Only EvaluationFlags::values, EvaluationFlags::gradients, "
9023 "and EvaluationFlags::hessians are supported."));
9024
9025 const auto shared_vector_data = internal::get_shared_vector_data(
9026 &input_vector,
9027 this->dof_access_index ==
9029 this->active_fe_index,
9030 this->dof_info);
9031
9032 if (this->data->data.front().fe_degree > 0 &&
9033 fast_evaluation_supported(this->data->data.front().fe_degree,
9034 this->data->data.front().n_q_points_1d) &&
9036 dim,
9037 typename VectorType::value_type,
9038 VectorizedArrayType>::
9039 supports(evaluation_flag,
9040 *this->data,
9041 internal::get_beginning<typename VectorType::value_type>(
9042 input_vector),
9043 this->dof_info->index_storage_variants[this->dof_access_index]
9044 [this->cell]))
9045 {
9046 if (fe_degree > -1)
9047 {
9049 dim,
9050 typename VectorType::value_type,
9051 VectorizedArrayType>::template run<fe_degree,
9052 n_q_points_1d>(
9053 n_components,
9054 evaluation_flag,
9055 internal::get_beginning<typename VectorType::value_type>(
9056 input_vector),
9057 shared_vector_data,
9058 *this);
9059 }
9060 else
9061 {
9063 dim,
9064 typename VectorType::value_type,
9065 VectorizedArrayType>::evaluate(n_components,
9066 evaluation_flag,
9067 internal::get_beginning<
9068 typename VectorType::value_type>(
9069 input_vector),
9070 shared_vector_data,
9071 *this);
9072 }
9073 }
9074 else
9075 {
9076 this->read_dof_values(input_vector);
9077 this->evaluate(evaluation_flag);
9078 }
9079
9080# ifdef DEBUG
9081 if (evaluation_flag & EvaluationFlags::values)
9082 this->values_quad_initialized = true;
9083 if (evaluation_flag & EvaluationFlags::gradients)
9084 this->gradients_quad_initialized = true;
9085 if (evaluation_flag & EvaluationFlags::hessians)
9086 this->hessians_quad_initialized = true;
9087# endif
9088}
9089
9090
9091
9092template <int dim,
9093 int fe_degree,
9094 int n_q_points_1d,
9095 int n_components_,
9096 typename Number,
9097 typename VectorizedArrayType>
9098template <typename VectorType>
9099inline void
9101 dim,
9102 fe_degree,
9103 n_q_points_1d,
9104 n_components_,
9105 Number,
9106 VectorizedArrayType>::integrate_scatter(const bool integrate_values,
9107 const bool integrate_gradients,
9108 VectorType &destination)
9109{
9111 ((integrate_values) ? EvaluationFlags::values : EvaluationFlags::nothing) |
9112 ((integrate_gradients) ? EvaluationFlags::gradients :
9114
9115 integrate_scatter(flag, destination);
9116}
9117
9118
9119
9120template <int dim,
9121 int fe_degree,
9122 int n_q_points_1d,
9123 int n_components_,
9124 typename Number,
9125 typename VectorizedArrayType>
9126template <typename VectorType>
9127inline void
9129 fe_degree,
9130 n_q_points_1d,
9131 n_components_,
9132 Number,
9133 VectorizedArrayType>::
9134 integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
9135 VectorType & destination)
9136{
9137 Assert((this->dof_access_index ==
9139 this->is_interior_face() == false) == false,
9141
9142 const auto shared_vector_data = internal::get_shared_vector_data(
9143 &destination,
9144 this->dof_access_index ==
9146 this->active_fe_index,
9147 this->dof_info);
9148
9149 if (this->data->data.front().fe_degree > 0 &&
9150 fast_evaluation_supported(this->data->data.front().fe_degree,
9151 this->data->data.front().n_q_points_1d) &&
9153 dim,
9154 typename VectorType::value_type,
9155 VectorizedArrayType>::
9156 supports(integration_flag,
9157 *this->data,
9158 internal::get_beginning<typename VectorType::value_type>(
9159 destination),
9160 this->dof_info->index_storage_variants[this->dof_access_index]
9161 [this->cell]))
9162 {
9163 if (fe_degree > -1)
9164 {
9166 dim,
9167 typename VectorType::value_type,
9168 VectorizedArrayType>::template run<fe_degree,
9169 n_q_points_1d>(
9170 n_components,
9171 integration_flag,
9172 internal::get_beginning<typename VectorType::value_type>(
9173 destination),
9174 shared_vector_data,
9175 *this);
9176 }
9177 else
9178 {
9180 dim,
9181 typename VectorType::value_type,
9182 VectorizedArrayType>::integrate(n_components,
9183 integration_flag,
9184 internal::get_beginning<
9185 typename VectorType::value_type>(
9186 destination),
9187 shared_vector_data,
9188 *this);
9189 }
9190 }
9191 else
9192 {
9193 integrate(integration_flag);
9194 this->distribute_local_to_global(destination);
9195 }
9196}
9197
9198
9199
9200template <int dim,
9201 int fe_degree,
9202 int n_q_points_1d,
9203 int n_components_,
9204 typename Number,
9205 typename VectorizedArrayType>
9208 fe_degree,
9209 n_q_points_1d,
9210 n_components_,
9211 Number,
9212 VectorizedArrayType>::dof_indices() const
9213{
9214 return {0U, dofs_per_cell};
9215}
9216
9217
9218
9219template <int dim,
9220 int fe_degree,
9221 int n_q_points_1d,
9222 int n_components_,
9223 typename Number,
9224 typename VectorizedArrayType>
9225bool
9226FEEvaluation<dim,
9227 fe_degree,
9228 n_q_points_1d,
9229 n_components_,
9230 Number,
9231 VectorizedArrayType>::
9232 fast_evaluation_supported(const unsigned int given_degree,
9233 const unsigned int given_n_q_points_1d)
9234{
9235 return fe_degree == -1 ?
9237 fast_evaluation_supported(given_degree, given_n_q_points_1d) :
9238 true;
9239}
9240
9241
9242
9243template <int dim,
9244 int fe_degree,
9245 int n_q_points_1d,
9246 int n_components_,
9247 typename Number,
9248 typename VectorizedArrayType>
9249bool
9251 fe_degree,
9252 n_q_points_1d,
9253 n_components_,
9254 Number,
9255 VectorizedArrayType>::
9256 fast_evaluation_supported(const unsigned int given_degree,
9257 const unsigned int given_n_q_points_1d)
9258{
9259 return fe_degree == -1 ?
9261 fast_evaluation_supported(given_degree, given_n_q_points_1d) :
9262 true;
9263}
9264
9265
9266
9267template <int dim,
9268 int fe_degree,
9269 int n_q_points_1d,
9270 int n_components_,
9271 typename Number,
9272 typename VectorizedArrayType>
9273bool
9275 fe_degree,
9276 n_q_points_1d,
9277 n_components_,
9278 Number,
9279 VectorizedArrayType>::at_boundary() const
9280{
9281 Assert(this->dof_access_index !=
9284
9285 if (this->is_interior_face() == false)
9286 return false;
9287 else if (this->cell < this->matrix_free->n_inner_face_batches())
9288 return false;
9289 else if (this->cell < (this->matrix_free->n_inner_face_batches() +
9290 this->matrix_free->n_boundary_face_batches()))
9291 return true;
9292 else
9293 return false;
9294}
9295
9296
9297
9298template <int dim,
9299 int fe_degree,
9300 int n_q_points_1d,
9301 int n_components_,
9302 typename Number,
9303 typename VectorizedArrayType>
9306 fe_degree,
9307 n_q_points_1d,
9308 n_components_,
9309 Number,
9310 VectorizedArrayType>::boundary_id() const
9311{
9312 Assert(this->dof_access_index !=
9315
9316 if (at_boundary())
9317 return this->matrix_free->get_boundary_id(this->cell);
9318 else
9320}
9321
9322
9323
9324/*------------------------- end FEFaceEvaluation ------------------------- */
9325
9326
9327#endif // ifndef DOXYGEN
9328
9329
9331
9332#endif
hessian_type get_hessian(unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
value_type get_laplacian(const unsigned int q_point) const
void submit_value(const gradient_type val_in, const unsigned int q_point)
value_type get_normal_derivative(const unsigned int q_point) const
void submit_normal_derivative(const gradient_type grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
void submit_gradient(const value_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
void submit_gradient(const Tensor< 2, 1, VectorizedArrayType > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< 1 > &mapping, const FiniteElement< 1 > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< 1, VectorizedArrayType, is_face > *other)
value_type get_divergence(const unsigned int q_point) const
value_type get_dof_value(const unsigned int dof) const
FEEvaluationAccess(const MatrixFree< 1, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
value_type get_dof_value(const unsigned int dof) const
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void submit_value(const Tensor< 1, 1, VectorizedArrayType > val_in, const unsigned int q_point)
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
value_type get_laplacian(const unsigned int q_point) const
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
value_type get_normal_derivative(const unsigned int q_point) const
void submit_value(const value_type val_in, const unsigned int q_point)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int dofs_per_cell, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
value_type get_value(const unsigned int q_point) const
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
Tensor< 3, dim, VectorizedArrayType > get_hessian(const unsigned int q_point) const
VectorizedArrayType get_divergence(const unsigned int q_point) const
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
gradient_type get_gradient(const unsigned int q_point) const
void submit_gradient(const Tensor< 1, dim, Tensor< 1, dim, VectorizedArrayType > > grad_in, const unsigned int q_point)
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
void submit_value(const Tensor< 1, dim, VectorizedArrayType > val_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
FEEvaluationAccess(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
static constexpr unsigned int dimension
Tensor< 1, n_components_, Tensor< 1, dim, VectorizedArrayType > > gradient_type
FEEvaluationAccess(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face=true, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
FEEvaluationAccess(const FEEvaluationAccess &other)
FEEvaluationAccess & operator=(const FEEvaluationAccess &other)
static constexpr unsigned int n_components
value_type get_dof_value(const unsigned int dof) const
Tensor< 1, n_components_, Tensor< 2, dim, VectorizedArrayType > > get_hessian(const unsigned int q_point) const
void read_write_operation_global(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors) const
void read_dof_values_plain(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
value_type get_laplacian(const unsigned int q_point) const
AlignedVector< VectorizedArrayType > * scratch_data_array
gradient_type get_gradient(const unsigned int q_point) const
void distribute_local_to_global(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
static constexpr unsigned int dimension
void submit_gradient(const gradient_type grad_in, const unsigned int q_point)
void submit_value(const value_type val_in, const unsigned int q_point)
Tensor< 1,(dim==2 ? 1 :dim), VectorizedArrayType > get_curl(const unsigned int q_point) const
void submit_dof_value(const value_type val_in, const unsigned int dof)
std::vector< types::global_dof_index > local_dof_indices
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
void submit_symmetric_gradient(const SymmetricTensor< 2, dim, VectorizedArrayType > grad_in, const unsigned int q_point)
VectorizedArrayType get_divergence(const unsigned int q_point) const
FEEvaluationBase(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component, const FEEvaluationData< dim, VectorizedArrayType, is_face > *other)
FEEvaluationBase & operator=(const FEEvaluationBase &other)
void read_dof_values(const VectorType &src, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip())
void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point)
static constexpr unsigned int n_components
gradient_type get_hessian_diagonal(const unsigned int q_point) const
void read_write_operation(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask, const bool apply_constraints=true) const
const MatrixFree< dim, Number, VectorizedArrayType > & get_matrix_free() const
void apply_hanging_node_constraints(const bool transpose) const
FEEvaluationBase(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no, const unsigned int first_selected_component, const unsigned int quad_no, const unsigned int fe_degree, const unsigned int n_q_points, const bool is_interior_face, const unsigned int active_fe_index, const unsigned int active_quad_index, const unsigned int face_type)
void submit_normal_derivative(const value_type grad_in, const unsigned int q_point)
void read_write_operation_contiguous(const VectorOperation &operation, const std::array< VectorType *, n_components_ > &vectors, const std::array< const std::vector< ArrayView< const typename VectorType::value_type > > *, n_components_ > &vectors_sm, const std::bitset< VectorizedArrayType::size()> &mask) const
value_type integrate_value() const
FEEvaluationBase(const FEEvaluationBase &other)
void submit_hessian(const hessian_type hessian_in, const unsigned int q_point)
SymmetricTensor< 2, dim, VectorizedArrayType > get_symmetric_gradient(const unsigned int q_point) const
const MatrixFree< dim, Number, VectorizedArrayType > * matrix_free
value_type get_normal_derivative(const unsigned int q_point) const
value_type get_value(const unsigned int q_point) const
void submit_curl(const Tensor< 1, dim==2 ? 1 :dim, VectorizedArrayType > curl_in, const unsigned int q_point)
void set_dof_values_plain(VectorType &dst, const unsigned int first_index=0, const std::bitset< VectorizedArrayType::size()> &mask=std::bitset< VectorizedArrayType::size()>().flip()) const
const MappingInfoStorageType::QuadratureDescriptor * descriptor
const MappingInfoStorageType * mapping_data
const ShapeInfoType * data
std::shared_ptr< internal::MatrixFreeFunctions::MappingDataOnTheFly< dim, Number > > mapped_geometry
const DoFInfo * dof_info
FEEvaluationData & operator=(const FEEvaluationData &other)
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
const unsigned int dofs_per_component
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
FEEvaluation(const Mapping< dim > &mapping, const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
const unsigned int n_q_points
void integrate(const bool integrate_values, const bool integrate_gradients)
void reinit(const std::array< unsigned int, VectorizedArrayType::size()> &cell_ids)
FEEvaluation(const FiniteElement< dim > &fe, const Quadrature< 1 > &quadrature, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void reinit(const typename Triangulation< dim >::cell_iterator &cell)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
void reinit(const unsigned int cell_batch_index)
FEEvaluation(const FiniteElement< dim > &fe, const FEEvaluationData< dim, VectorizedArrayType, false > &other, const unsigned int first_selected_component=0)
void reinit(const TriaIterator< DoFCellAccessor< dim, dim, level_dof_access > > &cell)
void evaluate(const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
static constexpr unsigned int tensor_dofs_per_cell
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
typename BaseClass::gradient_type gradient_type
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array, const bool sum_into_values=false)
static constexpr unsigned int dimension
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int static_dofs_per_cell
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
typename BaseClass::value_type value_type
FEEvaluation & operator=(const FEEvaluation &other)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int n_components
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients, const bool evaluate_hessians=false)
FEEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int)
static constexpr unsigned int static_n_q_points
FEEvaluation(const FEEvaluation &other)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void check_template_arguments(const unsigned int fe_no, const unsigned int first_selected_component)
static constexpr unsigned int static_dofs_per_component
typename BaseClass::value_type value_type
std_cxx20::ranges::iota_view< unsigned int, unsigned int > dof_indices() const
bool at_boundary() const
void integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag, VectorType &output_vector)
static constexpr unsigned int static_n_q_points_cell
static constexpr unsigned int tensor_dofs_per_cell
const unsigned int dofs_per_component
void reinit(const unsigned int face_batch_number)
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void reinit(const unsigned int cell_batch_number, const unsigned int face_number)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int given_n_q_points_1d)
void evaluate(const VectorizedArrayType *values_array, const EvaluationFlags::EvaluationFlags evaluation_flag)
const unsigned int dofs_per_cell
static constexpr unsigned int n_components
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
void integrate(const bool integrate_values, const bool integrate_gradients, VectorizedArrayType *values_array)
void evaluate(const VectorizedArrayType *values_array, const bool evaluate_values, const bool evaluate_gradients)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::pair< unsigned int, unsigned int > &range, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
void integrate(const EvaluationFlags::EvaluationFlags integration_flag, VectorizedArrayType *values_array)
static constexpr unsigned int static_dofs_per_component
static constexpr unsigned int static_n_q_points
FEFaceEvaluation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const bool is_interior_face=true, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0, const unsigned int active_fe_index=numbers::invalid_unsigned_int, const unsigned int active_quad_index=numbers::invalid_unsigned_int, const unsigned int face_type=numbers::invalid_unsigned_int)
static constexpr unsigned int dimension
typename BaseClass::gradient_type gradient_type
types::boundary_id boundary_id() const
void integrate_scatter(const bool integrate_values, const bool integrate_gradients, VectorType &output_vector)
static constexpr unsigned int static_dofs_per_cell
void integrate(const bool integrate_values, const bool integrate_gradients)
std::pair< unsigned int, unsigned int > component_to_base_index(const unsigned int component) const
unsigned int element_multiplicity(const unsigned int index) const
Abstract base class for mapping classes.
Definition mapping.h:317
const internal::MatrixFreeFunctions::TaskInfo & get_task_info() const
const internal::MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > & get_shape_info(const unsigned int dof_handler_index_component=0, const unsigned int quad_index=0, const unsigned int fe_base_element=0, const unsigned int hp_active_fe_index=0, const unsigned int hp_active_quad_index=0) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Table< 3, unsigned int > & get_cell_and_face_to_plain_faces() const
unsigned int n_inner_face_batches() const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
const internal::MatrixFreeFunctions::FaceToCellTopology< VectorizedArrayType::size()> & get_face_info(const unsigned int face_batch_index) const
const internal::MatrixFreeFunctions::DoFInfo & get_dof_info(const unsigned int dof_handler_index_component=0) const
AlignedVector< VectorizedArrayType > * acquire_scratch_data() const
const Number * constraint_pool_begin(const unsigned int pool_index) const
void release_scratch_data(const AlignedVector< VectorizedArrayType > *memory) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
bool indices_initialized() const
const Number * constraint_pool_end(const unsigned int pool_index) const
unsigned int n_components() const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
unsigned int n_base_elements(const unsigned int dof_handler_index) const
Definition point.h:112
DEAL_II_HOST constexpr const Number & access_raw_entry(const unsigned int unrolled_index) const
#define DEAL_II_ALWAYS_INLINE
Definition config.h:106
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:138
#define DEAL_II_DEPRECATED
Definition config.h:172
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcAccessToUninitializedField()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMatrixFreeAccessToUninitializedMappingField(std::string arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
constexpr compressed_constraint_kind unconstrained_compressed_constraint_kind
void check_vector_compatibility(const VectorType &vec, const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const internal::MatrixFreeFunctions::DoFInfo &dof_info)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const types::boundary_id internal_face_boundary_id
Definition types.h:277
static const unsigned int invalid_unsigned_int
Definition types.h:213
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:46
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &eval, const bool sum_into_values_array=false)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, false > &fe_eval, const bool sum_into_values_array)
static void apply(const unsigned int n_components, const unsigned int fe_degree, const MatrixFreeFunctions::ShapeInfo< VectorizedArrayType > &shape_info, const bool transpose, const std::array< MatrixFreeFunctions::compressed_constraint_kind, VectorizedArrayType::size()> &c_mask, VectorizedArrayType *values)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
static bool fast_evaluation_supported(const unsigned int given_degree, const unsigned int n_q_points_1d)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, Number *values_dofs, FEEvaluationData< dim, Number, true > &fe_eval)
std::vector< std::pair< unsigned short, unsigned short > > constraint_indicator
Definition dof_info.h:534
std::vector< std::pair< unsigned int, unsigned int > > row_starts
Definition dof_info.h:493
std::vector< std::vector< unsigned int > > component_dof_indices_offset
Definition dof_info.h:686
std::vector< std::vector< bool > > hanging_node_constraint_masks_comp
Definition dof_info.h:516
unsigned int fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const
std::vector< unsigned int > dof_indices
Definition dof_info.h:510
std::vector< compressed_constraint_kind > hanging_node_constraint_masks
Definition dof_info.h:522
std::array< std::vector< unsigned int >, 3 > dof_indices_interleave_strides
Definition dof_info.h:570
std::array< std::vector< std::pair< unsigned int, unsigned int > >, 3 > dof_indices_contiguous_sm
Definition dof_info.h:560
std::vector< unsigned int > row_starts_plain_indices
Definition dof_info.h:633
std::vector< unsigned int > component_to_base_index
Definition dof_info.h:673
std::array< std::vector< unsigned int >, 3 > dof_indices_contiguous
Definition dof_info.h:549
std::vector< unsigned int > plain_dof_indices
Definition dof_info.h:643
std::array< std::vector< unsigned char >, 3 > n_vectorization_lanes_filled
Definition dof_info.h:581
std::vector< unsigned int > dof_indices_interleaved
Definition dof_info.h:539
std::array< std::vector< IndexStorageVariants >, 3 > index_storage_variants
Definition dof_info.h:485
unsigned int quad_index_from_n_q_points(const unsigned int n_q_points) const
std::vector< unsigned int > face_partition_data
Definition task_info.h:497
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)