Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
fe_point_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_fe_point_evaluation_h
16#define dealii_fe_point_evaluation_h
17
18#include <deal.II/base/config.h>
19
24#include <deal.II/base/tensor.h>
26
28#include <deal.II/fe/mapping.h>
29
35
37
39
40namespace internal
41{
43 {
46 std::string,
47 << "You are requesting information from an FEPointEvaluationBase "
48 << "object for which this kind of information has not been computed. "
49 << "What information these objects compute is determined by the update_* "
50 << "flags you pass to MappingInfo() in the Constructor. Here, "
51 << "the operation you are attempting requires the <" << arg1
52 << "> flag to be set, but it was apparently not specified "
53 << "upon initialization.");
54
59 template <int dim,
60 int spacedim,
61 int n_components,
62 typename Number,
63 typename Enable = void>
65 {
69 typename ::internal::VectorizedArrayTrait<
77 using real_gradient_type = std::conditional_t<
78 n_components == spacedim,
87
88 static void
89 read_value(const ScalarNumber vector_entry,
90 const unsigned int component,
91 scalar_value_type &result)
92 {
93 AssertIndexRange(component, n_components);
94 result[component] = vector_entry;
95 }
96
99 {
100 return result;
101 }
102
103 static scalar_value_type
105 {
106 scalar_value_type result_scalar = {};
107
108 for (unsigned int c = 0; c < n_components; ++c)
109 result_scalar[c] = result[c].sum();
110
111 return result_scalar;
112 }
113
114 static ScalarNumber
115 sum_value(const unsigned int component,
116 const vectorized_value_type &result)
117 {
118 AssertIndexRange(component, n_components);
119 return result[component].sum();
120 }
121
122 static void
124 const unsigned int vector_lane,
125 unit_gradient_type &result)
126 {
127 for (unsigned int i = 0; i < n_components; ++i)
128 for (unsigned int d = 0; d < dim; ++d)
129 result[i][d] =
131 value[d][i], vector_lane);
132 }
133
134 static void
136 const unsigned int vector_lane,
137 const unit_gradient_type &result)
138 {
139 for (unsigned int i = 0; i < n_components; ++i)
140 for (unsigned int d = 0; d < dim; ++d)
142 value[d][i], vector_lane) = result[i][d];
143 }
144
145 static void
147 const unsigned int vector_lane,
149 {
150 for (unsigned int i = 0; i < n_components; ++i)
151 for (unsigned int d = 0; d < dim; ++d)
153 value[d][i], vector_lane) = result[i][d];
154 }
155
156 static void
158 const unsigned int vector_lane)
159 {
160 for (unsigned int i = 0; i < n_components; ++i)
161 for (unsigned int d = 0; d < spacedim; ++d)
163 vector_lane) = 0.;
164 }
165
166 static void
168 const unsigned int vector_lane,
169 scalar_value_type &result)
170 {
171 for (unsigned int i = 0; i < n_components; ++i)
172 result[i] = value[i][vector_lane];
173 }
174
175 static void
177 const unsigned int,
178 vectorized_value_type &result)
179 {
180 result = value;
181 }
182
183 static void
185 const unsigned int vector_lane,
186 const scalar_value_type &result)
187 {
188 for (unsigned int i = 0; i < n_components; ++i)
189 value[i][vector_lane] = result[i];
190 }
191
192 static void
194 const unsigned int,
195 const vectorized_value_type &result)
196 {
197 value = result;
198 }
199
200 static void
201 set_zero_value(value_type &value, const unsigned int vector_lane)
202 {
203 for (unsigned int i = 0; i < n_components; ++i)
205 0.;
206 }
207
208 static void
210 const unsigned int vector_lane,
211 const unsigned int component,
212 const ScalarNumber &shape_value)
213 {
215 vector_lane) += shape_value;
216 }
217
218 static ScalarNumber
220 const unsigned int vector_lane,
221 const unsigned int component)
222 {
224 vector_lane);
225 }
226
227 static void
229 const unsigned int vector_lane,
230 const unsigned int component,
231 const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
232 {
233 for (unsigned int d = 0; d < spacedim; ++d)
235 vector_lane) +=
236 shape_gradient[d];
237 }
238
241 const unsigned int vector_lane,
242 const unsigned int component)
243 {
245 for (unsigned int d = 0; d < spacedim; ++d)
246 result[d] =
248 vector_lane);
249 return result;
250 }
251 };
252
253 template <int dim, int spacedim, typename Number>
254 struct EvaluatorTypeTraits<dim, spacedim, 1, Number>
255 {
259 typename ::internal::VectorizedArrayTrait<
261 using value_type = Number;
270
271 static void
272 read_value(const ScalarNumber vector_entry,
273 const unsigned int,
274 scalar_value_type &result)
275 {
276 result = vector_entry;
277 }
278
279 static scalar_value_type
281 {
282 return result;
283 }
284
285 static scalar_value_type
287 {
288 return result.sum();
289 }
290
291 static ScalarNumber
292 sum_value(const unsigned int, const vectorized_value_type &result)
293 {
294 return result.sum();
295 }
296
297 static void
299 const unsigned int vector_lane,
301 {
302 for (unsigned int d = 0; d < dim; ++d)
303 result[d] = value[d][vector_lane];
304 }
305
306 static void
308 const unsigned int,
310 {
311 result = value;
312 }
313
314 static void
316 const unsigned int vector_lane,
317 const scalar_unit_gradient_type &result)
318 {
319 for (unsigned int d = 0; d < dim; ++d)
320 value[d][vector_lane] = result[d];
321 }
322
323 static void
325 const unsigned int,
326 const vectorized_unit_gradient_type &result)
327 {
328 value = result;
329 }
330
331 static void
333 const unsigned int vector_lane)
334 {
335 for (unsigned int d = 0; d < spacedim; ++d)
337 0.;
338 }
339
340 static void
342 const unsigned int vector_lane,
343 scalar_value_type &result)
344 {
345 result = value[vector_lane];
346 }
347
348 static void
350 const unsigned int,
351 vectorized_value_type &result)
352 {
353 result = value;
354 }
355
356 static void
358 const unsigned int vector_lane,
359 const scalar_value_type &result)
360 {
361 value[vector_lane] = result;
362 }
363
364 static void
366 const unsigned int,
367 const vectorized_value_type &result)
368 {
369 value = result;
370 }
371
372 static void
373 set_zero_value(value_type &value, const unsigned int vector_lane)
374 {
376 }
377
378 static void
380 const unsigned int vector_lane,
381 const unsigned int,
382 const ScalarNumber &shape_value)
383 {
385 shape_value;
386 }
387
388 static ScalarNumber
390 const unsigned int vector_lane,
391 const unsigned int)
392 {
394 }
395
396 static void
398 const unsigned int vector_lane,
399 const unsigned int,
400 const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
401 {
402 for (unsigned int d = 0; d < spacedim; ++d)
404 shape_gradient[d];
405 }
406
409 const unsigned int vector_lane,
410 const unsigned int)
411 {
413 for (unsigned int d = 0; d < spacedim; ++d)
414 result[d] =
416 return result;
417 }
418 };
419
420 template <int dim, typename Number>
422 dim,
423 dim,
424 Number,
425 std::enable_if_t<dim != 1>>
426 {
430 typename ::internal::VectorizedArrayTrait<
441
442 static void
443 read_value(const ScalarNumber vector_entry,
444 const unsigned int component,
445 scalar_value_type &result)
446 {
447 AssertIndexRange(component, dim);
448 result[component] = vector_entry;
449 }
450
451 static scalar_value_type
453 {
454 return result;
455 }
456
457 static scalar_value_type
459 {
460 scalar_value_type result_scalar = {};
461
462 for (unsigned int c = 0; c < dim; ++c)
463 result_scalar[c] = result[c].sum();
464
465 return result_scalar;
466 }
467
468 static ScalarNumber
469 sum_value(const unsigned int component,
470 const vectorized_value_type &result)
471 {
472 AssertIndexRange(component, dim);
473 return result[component].sum();
474 }
475
476 static void
478 const unsigned int vector_lane,
479 unit_gradient_type &result)
480 {
481 for (unsigned int i = 0; i < dim; ++i)
482 for (unsigned int d = 0; d < dim; ++d)
483 result[i][d] =
485 value[d][i], vector_lane);
486 }
487
488 static void
490 const unsigned int vector_lane,
491 const unit_gradient_type &result)
492 {
493 for (unsigned int i = 0; i < dim; ++i)
494 for (unsigned int d = 0; d < dim; ++d)
496 value[d][i], vector_lane) = result[i][d];
497 }
498
499 static void
501 const unsigned int vector_lane)
502 {
503 for (unsigned int i = 0; i < dim; ++i)
504 for (unsigned int d = 0; d < dim; ++d)
506 vector_lane) = 0.;
507 }
508
509 static void
511 const unsigned int vector_lane,
512 scalar_value_type &result)
513 {
514 for (unsigned int i = 0; i < dim; ++i)
515 result[i] = value[i][vector_lane];
516 }
517
518 static void
520 const unsigned int,
521 vectorized_value_type &result)
522 {
523 result = value;
524 }
525
526 static void
528 const unsigned int vector_lane,
529 const scalar_value_type &result)
530 {
531 for (unsigned int i = 0; i < dim; ++i)
532 value[i][vector_lane] = result[i];
533 }
534
535 static void
537 const unsigned int,
538 const vectorized_value_type &result)
539 {
540 value = result;
541 }
542
543 static void
544 set_zero_value(value_type &value, const unsigned int vector_lane)
545 {
546 for (unsigned int i = 0; i < dim; ++i)
548 0.;
549 }
550
551 static void
553 const unsigned int vector_lane,
554 const unsigned int component,
555 const ScalarNumber &shape_value)
556 {
558 vector_lane) += shape_value;
559 }
560
561 static ScalarNumber
563 const unsigned int vector_lane,
564 const unsigned int component)
565 {
567 vector_lane);
568 }
569
570 static void
572 const unsigned int vector_lane,
573 const unsigned int component,
574 const Tensor<1, dim, ScalarNumber> &shape_gradient)
575 {
576 for (unsigned int d = 0; d < dim; ++d)
578 vector_lane) +=
579 shape_gradient[d];
580 }
581
584 const unsigned int vector_lane,
585 const unsigned int component)
586 {
588 for (unsigned int d = 0; d < dim; ++d)
589 result[d] =
591 vector_lane);
592 return result;
593 }
594 };
595
596 template <int dim, int spacedim>
597 bool
599 const unsigned int base_element_number);
600
601 template <int dim, int spacedim>
602 bool
604
605 template <int dim, int spacedim>
606 std::vector<Polynomials::Polynomial<double>>
608 } // namespace FEPointEvaluation
609} // namespace internal
610
611
612
619template <int n_components_,
620 int dim,
621 int spacedim = dim,
622 typename Number = double>
624{
625public:
626 static constexpr unsigned int dimension = dim;
627 static constexpr unsigned int n_components = n_components_;
628
629 using number_type = Number;
630
633 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
635 using ETT = typename internal::FEPointEvaluation::
636 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
637 using value_type = typename ETT::value_type;
638 using scalar_value_type = typename ETT::scalar_value_type;
639 using vectorized_value_type = typename ETT::vectorized_value_type;
640 using gradient_type = typename ETT::real_gradient_type;
642 typename ETT::interface_vectorized_unit_gradient_type;
643
644protected:
666 const unsigned int first_selected_component = 0);
667
690 const unsigned int first_selected_component = 0,
691 const bool is_interior = true);
692
697
702
703public:
711 const value_type &
712 get_value(const unsigned int point_index) const;
713
722 void
723 submit_value(const value_type &value, const unsigned int point_index);
724
734 const gradient_type &
735 get_gradient(const unsigned int point_index) const;
736
745 void
746 submit_gradient(const gradient_type &, const unsigned int point_index);
747
755 Number
756 get_divergence(const unsigned int point_index) const;
757
773 void
774 submit_divergence(const Number &value, const unsigned int point_index);
775
782 Tensor<1, (dim == 2 ? 1 : dim), Number>
783 get_curl(const unsigned int point_index) const;
784
791 jacobian(const unsigned int point_index) const;
792
800 inverse_jacobian(const unsigned int point_index) const;
801
807 Number
808 JxW(const unsigned int point_index) const;
809
817 real_point(const unsigned int point_index) const;
818
824 quadrature_point(const unsigned int point_index) const;
825
831 unit_point(const unsigned int point_index) const;
832
842
850
854 unsigned int
856
857protected:
858 static constexpr std::size_t n_lanes_user_interface =
860 static constexpr std::size_t n_lanes_internal =
862 static constexpr std::size_t stride =
864
873 void
874 setup(const unsigned int first_selected_component);
875
881 template <bool is_face, bool is_linear>
882 void
884
888 const unsigned int n_q_batches;
889
893 const unsigned int n_q_points;
894
898 const unsigned int n_q_points_scalar;
899
904
909
914 std::vector<Polynomials::Polynomial<double>> poly;
915
920
925 std::vector<unsigned int> renumber;
926
934 std::vector<scalar_value_type> solution_renumbered;
935
943
948
952 std::vector<value_type> values;
953
957 std::vector<gradient_type> gradients;
958
964
970
976
982
988
994
1000 const Number *JxW_ptr;
1001
1006
1012
1019
1024
1029
1035 std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
1036
1041
1045 std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1046
1050 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1052
1059
1064
1069
1074
1082
1088
1094
1095 const bool is_interior;
1096};
1097
1098
1099
1130template <int n_components_,
1131 int dim,
1132 int spacedim = dim,
1133 typename Number = double>
1135 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1136{
1137public:
1138 static constexpr unsigned int dimension = dim;
1139 static constexpr unsigned int n_components = n_components_;
1140
1141 using number_type = Number;
1142
1145 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1147 using ETT = typename internal::FEPointEvaluation::
1148 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1149 using value_type = typename ETT::value_type;
1150 using scalar_value_type = typename ETT::scalar_value_type;
1151 using vectorized_value_type = typename ETT::vectorized_value_type;
1152 using unit_gradient_type = typename ETT::unit_gradient_type;
1153 using gradient_type = typename ETT::real_gradient_type;
1155 typename ETT::interface_vectorized_unit_gradient_type;
1156
1178 const unsigned int first_selected_component = 0);
1179
1199 const unsigned int first_selected_component = 0);
1200
1212 void
1214 const ArrayView<const Point<dim>> &unit_points);
1215
1220 void
1221 reinit();
1222
1227 void
1228 reinit(const unsigned int cell_index);
1229
1230
1242 template <std::size_t stride_view>
1243 void
1244 evaluate(
1246 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1247
1259 void
1260 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1261 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1262
1285 template <std::size_t stride_view>
1286 void
1288 const EvaluationFlags::EvaluationFlags &integration_flags,
1289 const bool sum_into_values = false);
1290
1313 void
1314 integrate(const ArrayView<ScalarNumber> &solution_values,
1315 const EvaluationFlags::EvaluationFlags &integration_flags,
1316 const bool sum_into_values = false);
1317
1344 template <std::size_t stride_view>
1345 void
1347 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1348 const EvaluationFlags::EvaluationFlags &integration_flags,
1349 const bool sum_into_values = false);
1350
1377 void
1378 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1379 const EvaluationFlags::EvaluationFlags &integration_flags,
1380 const bool sum_into_values = false);
1381
1388 normal_vector(const unsigned int point_index) const;
1389
1395 const value_type
1396 get_normal_derivative(const unsigned int point_index) const;
1397
1403 void
1404 submit_normal_derivative(const value_type &, const unsigned int point_index);
1405
1406private:
1407 static constexpr std::size_t n_lanes_user_interface =
1409 static constexpr std::size_t n_lanes_internal =
1411 static constexpr std::size_t stride =
1413
1418 template <bool is_linear, std::size_t stride_view>
1419 void
1422
1427 template <bool is_linear, std::size_t stride_view>
1428 void
1431 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1432 const unsigned int n_shapes,
1433 const unsigned int qb,
1434 vectorized_value_type &value,
1436
1440 template <bool is_linear, std::size_t stride_view>
1441 void
1444 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1445
1449 template <std::size_t stride_view>
1450 void
1453 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1454
1460 template <bool is_linear>
1461 void
1463 const EvaluationFlags::EvaluationFlags &integration_flags,
1464 const unsigned int n_shapes,
1465 const unsigned int qb,
1466 const vectorized_value_type value,
1468 vectorized_value_type *solution_values_vectorized_linear);
1469
1475 template <bool is_linear, std::size_t stride_view>
1476 void
1478 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1479 vectorized_value_type *solution_values_vectorized_linear,
1480 const bool sum_into_values);
1481
1485 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1486 void
1488 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1489 const EvaluationFlags::EvaluationFlags &integration_flags,
1490 const bool sum_into_values);
1491
1495 template <bool do_JxW, std::size_t stride_view>
1496 void
1498 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1499 const EvaluationFlags::EvaluationFlags &integration_flags,
1500 const bool sum_into_values);
1501
1505 template <bool do_JxW, std::size_t stride_view>
1506 void
1508 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1509 const EvaluationFlags::EvaluationFlags &integration_flags,
1510 const bool sum_into_values);
1511
1517 void
1519};
1520
1521
1522
1540template <int n_components_,
1541 int dim,
1542 int spacedim = dim,
1543 typename Number = double>
1545 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1546{
1547public:
1548 static constexpr unsigned int dimension = dim;
1549 static constexpr unsigned int n_components = n_components_;
1550
1551 using number_type = Number;
1552
1555 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1557 using ETT = typename internal::FEPointEvaluation::
1558 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1559 using value_type = typename ETT::value_type;
1560 using scalar_value_type = typename ETT::scalar_value_type;
1561 using vectorized_value_type = typename ETT::vectorized_value_type;
1562 using unit_gradient_type = typename ETT::unit_gradient_type;
1563 using gradient_type = typename ETT::real_gradient_type;
1565 typename ETT::interface_vectorized_unit_gradient_type;
1566
1573 const bool is_interior = true,
1574 const unsigned int first_selected_component = 0);
1575
1580 void
1581 reinit(const unsigned int cell_index, const unsigned int face_number);
1582
1587 void
1588 reinit(const unsigned int face_index);
1589
1601 template <std::size_t stride_view>
1602 void
1603 evaluate(
1605 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1606
1618 void
1619 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1620 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1621
1644 template <std::size_t stride_view>
1645 void
1647 const EvaluationFlags::EvaluationFlags &integration_flags,
1648 const bool sum_into_values = false);
1649
1672 void
1673 integrate(const ArrayView<ScalarNumber> &solution_values,
1674 const EvaluationFlags::EvaluationFlags &integration_flags,
1675 const bool sum_into_values = false);
1676
1699 template <std::size_t stride_view>
1700 void
1702 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1703 const EvaluationFlags::EvaluationFlags &integration_flags,
1704 const bool sum_into_values = false);
1705
1728 void
1729 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1730 const EvaluationFlags::EvaluationFlags &integration_flags,
1731 const bool sum_into_values = false);
1732
1739 template <int stride_face_dof = VectorizedArrayType::size()>
1740 void
1741 evaluate_in_face(const ScalarNumber *face_dof_values,
1742 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1743
1750 template <int stride_face_dof = VectorizedArrayType::size()>
1751 void
1752 integrate_in_face(ScalarNumber *face_dof_values,
1753 const EvaluationFlags::EvaluationFlags &integration_flags,
1754 const bool sum_into_values = false);
1755
1762 normal_vector(const unsigned int point_index) const;
1763
1769 const value_type
1770 get_normal_derivative(const unsigned int point_index) const;
1771
1777 void
1778 submit_normal_derivative(const value_type &, const unsigned int point_index);
1779
1780private:
1781 static constexpr std::size_t n_lanes_user_interface =
1783 static constexpr std::size_t n_lanes_internal =
1785 static constexpr std::size_t stride =
1787
1788 template <bool is_linear, std::size_t stride_view>
1789 void
1792 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1793
1794 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1795 void
1797 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1798 const EvaluationFlags::EvaluationFlags &integration_flags,
1799 const bool sum_into_values);
1800
1805 template <bool is_linear, int stride_face_dof>
1806 void
1807 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1808 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1809
1814 template <bool do_JxW, bool is_linear, int stride_face_dof>
1815 void
1817 ScalarNumber *face_dof_values,
1818 const EvaluationFlags::EvaluationFlags &integration_flags,
1819 const bool sum_into_values);
1820};
1821
1822
1823
1824// ----------------------- template and inline function ----------------------
1825
1826
1827template <int n_components_, int dim, int spacedim, typename Number>
1831 const UpdateFlags update_flags,
1832 const unsigned int first_selected_component)
1833 : n_q_batches(numbers::invalid_unsigned_int)
1834 , n_q_points(numbers::invalid_unsigned_int)
1835 , n_q_points_scalar(numbers::invalid_unsigned_int)
1836 , mapping(&mapping)
1837 , fe(&fe)
1838 , JxW_ptr(nullptr)
1839 , update_flags(update_flags)
1840 , mapping_info_on_the_fly(
1841 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1842 mapping,
1843 update_flags))
1844 , mapping_info(mapping_info_on_the_fly.get())
1845 , current_cell_index(numbers::invalid_unsigned_int)
1846 , current_face_number(numbers::invalid_unsigned_int)
1847 , must_reinitialize_pointers(false)
1848 , is_interior(true)
1849{
1850 setup(first_selected_component);
1851}
1852
1853
1854
1855template <int n_components_, int dim, int spacedim, typename Number>
1860 const unsigned int first_selected_component,
1861 const bool is_interior)
1862 : n_q_batches(numbers::invalid_unsigned_int)
1863 , n_q_points(numbers::invalid_unsigned_int)
1864 , n_q_points_scalar(numbers::invalid_unsigned_int)
1865 , mapping(&mapping_info.get_mapping())
1866 , fe(&fe)
1867 , JxW_ptr(nullptr)
1868 , update_flags(mapping_info.get_update_flags())
1869 , mapping_info(&mapping_info)
1870 , current_cell_index(numbers::invalid_unsigned_int)
1871 , current_face_number(numbers::invalid_unsigned_int)
1872 , must_reinitialize_pointers(true)
1873 , is_interior(is_interior)
1874{
1875 setup(first_selected_component);
1876}
1877
1878
1879
1880template <int n_components_, int dim, int spacedim, typename Number>
1884 : n_q_batches(other.n_q_batches)
1885 , n_q_points(other.n_q_points)
1886 , n_q_points_scalar(other.n_q_points_scalar)
1887 , mapping(other.mapping)
1888 , fe(other.fe)
1889 , poly(other.poly)
1890 , use_linear_path(other.use_linear_path)
1891 , renumber(other.renumber)
1892 , solution_renumbered(other.solution_renumbered)
1893 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1894 , values(other.values)
1895 , gradients(other.gradients)
1896 , dofs_per_component(other.dofs_per_component)
1897 , dofs_per_component_face(other.dofs_per_component_face)
1898 , component_in_base_element(other.component_in_base_element)
1899 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1900 , update_flags(other.update_flags)
1901 , fe_values(other.fe_values)
1902 , mapping_info_on_the_fly(
1903 other.mapping_info_on_the_fly ?
1905 *mapping,
1906 update_flags) :
1907 nullptr)
1908 , mapping_info(other.mapping_info)
1909 , current_cell_index(other.current_cell_index)
1910 , current_face_number(other.current_face_number)
1911 , fast_path(other.fast_path)
1912 , must_reinitialize_pointers(true)
1913 , is_interior(other.is_interior)
1914{}
1915
1916
1917
1918template <int n_components_, int dim, int spacedim, typename Number>
1920 FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept = default;
1921
1922
1923
1924template <int n_components_, int dim, int spacedim, typename Number>
1925void
1927 const unsigned int first_selected_component)
1928{
1929 AssertIndexRange(first_selected_component + n_components,
1930 fe->n_components() + 1);
1931
1932 shapes.reserve(100);
1933
1934 bool same_base_element = true;
1935 unsigned int base_element_number = 0;
1936 component_in_base_element = 0;
1937 unsigned int component = 0;
1938 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1939 if (component + fe->element_multiplicity(base_element_number) >
1940 first_selected_component)
1941 {
1942 if (first_selected_component + n_components >
1943 component + fe->element_multiplicity(base_element_number))
1944 same_base_element = false;
1945 component_in_base_element = first_selected_component - component;
1946 break;
1947 }
1948 else
1949 component += fe->element_multiplicity(base_element_number);
1950
1953 *fe, base_element_number) &&
1954 same_base_element)
1955 {
1956 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1957 renumber = shape_info.lexicographic_numbering;
1958 dofs_per_component = shape_info.dofs_per_component_on_cell;
1959 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1961 fe->base_element(base_element_number));
1962
1963 bool is_lexicographic = true;
1964 for (unsigned int i = 0; i < renumber.size(); ++i)
1965 if (i != renumber[i])
1966 is_lexicographic = false;
1967
1968 if (is_lexicographic)
1969 renumber.clear();
1970
1971 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1972 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1973 poly[1].value(1.) == 1.) &&
1974 (fe->n_components() == n_components);
1975
1976 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1977 const unsigned int size_cell = dofs_per_component * n_components;
1978 scratch_data_scalar.resize(size_face + size_cell);
1979
1980 solution_renumbered.resize(dofs_per_component);
1981 solution_renumbered_vectorized.resize(dofs_per_component);
1982
1983 fast_path = true;
1984 }
1985 else
1986 {
1987 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1988 for (unsigned int d = 0; d < n_components; ++d)
1989 {
1990 const unsigned int component = first_selected_component + d;
1991 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1992 {
1993 const bool is_primitive =
1994 fe->is_primitive() || fe->is_primitive(i);
1995 if (is_primitive)
1996 nonzero_shape_function_component[i][d] =
1997 (component == fe->system_to_component_index(i).first);
1998 else
1999 nonzero_shape_function_component[i][d] =
2000 (fe->get_nonzero_components(i)[component] == true);
2001 }
2002 }
2003
2004 fast_path = false;
2005 }
2006}
2007
2008
2009
2010template <int n_components_, int dim, int spacedim, typename Number>
2011template <bool is_face, bool is_linear>
2012inline void
2014{
2015 const unsigned int geometry_index =
2016 mapping_info->template compute_geometry_index_offset<is_face>(
2017 current_cell_index, current_face_number);
2018
2019 cell_type = mapping_info->get_cell_type(geometry_index);
2020
2021 const_cast<unsigned int &>(n_q_points_scalar) =
2022 mapping_info->get_n_q_points_unvectorized(geometry_index);
2023
2024 // round up n_q_points_scalar / n_lanes_internal
2025 const_cast<unsigned int &>(n_q_batches) =
2026 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
2027
2028 const unsigned int n_q_points_before = n_q_points;
2029
2030 const_cast<unsigned int &>(n_q_points) =
2031 (stride == 1) ? n_q_batches : n_q_points_scalar;
2032
2033 if (n_q_points != n_q_points_before)
2034 {
2035 if (update_flags & update_values)
2036 values.resize(n_q_points);
2037 if (update_flags & update_gradients)
2038 gradients.resize(n_q_points);
2039 }
2040
2041 if (n_q_points == 0)
2042 return;
2043
2044 // set unit point pointer
2045 const unsigned int unit_point_offset =
2046 mapping_info->compute_unit_point_index_offset(geometry_index);
2047
2048 if (is_face)
2049 unit_point_faces_ptr =
2050 mapping_info->get_unit_point_faces(unit_point_offset);
2051 else
2052 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2053
2054 // set data pointers
2055 const unsigned int data_offset =
2056 mapping_info->compute_data_index_offset(geometry_index);
2057 const unsigned int compressed_data_offset =
2058 mapping_info->compute_compressed_data_index_offset(geometry_index);
2059 if constexpr (running_in_debug_mode())
2060 {
2061 const UpdateFlags update_flags_mapping =
2062 mapping_info->get_update_flags_mapping();
2063 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2064 real_point_ptr = mapping_info->get_real_point(data_offset);
2065 if (update_flags_mapping & UpdateFlags::update_jacobians)
2066 jacobian_ptr =
2067 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2068 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2069 inverse_jacobian_ptr =
2070 mapping_info->get_inverse_jacobian(compressed_data_offset,
2071 is_interior);
2072 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2073 normal_ptr = mapping_info->get_normal_vector(data_offset);
2074 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2075 JxW_ptr = mapping_info->get_JxW(data_offset);
2076 }
2077 else
2078 {
2079 real_point_ptr = mapping_info->get_real_point(data_offset);
2080 jacobian_ptr =
2081 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2082 inverse_jacobian_ptr =
2083 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2084 normal_ptr = mapping_info->get_normal_vector(data_offset);
2085 JxW_ptr = mapping_info->get_JxW(data_offset);
2086 }
2087
2088 if (!is_linear && fast_path)
2089 {
2090 const std::size_t n_shapes = poly.size();
2091 if (is_face)
2092 shapes_faces.resize_fast(n_q_batches * n_shapes);
2093 else
2094 shapes.resize_fast(n_q_batches * n_shapes);
2095
2096 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2097 if (is_face)
2098 {
2099 if (dim > 1)
2100 {
2102 shapes_faces.data() + qb * n_shapes,
2103 poly,
2104 unit_point_faces_ptr[qb],
2105 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2106 }
2107 }
2108 else
2109 {
2110 if (update_flags & UpdateFlags::update_gradients)
2111 {
2112 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2113 poly,
2114 unit_point_ptr[qb],
2115 1);
2116 }
2117 else if (qb + 1 < n_q_batches)
2118 {
2119 // Use function with reduced overhead to compute for two
2120 // points at once
2122 shapes.data() + qb * n_shapes,
2123 poly,
2124 unit_point_ptr[qb],
2125 unit_point_ptr[qb + 1]);
2126 ++qb;
2127 }
2128 else
2129 {
2130 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2131 poly,
2132 unit_point_ptr[qb],
2133 0);
2134 }
2135 }
2136 }
2137}
2138
2139
2140
2141template <int n_components_, int dim, int spacedim, typename Number>
2142inline const typename FEPointEvaluationBase<n_components_,
2143 dim,
2144 spacedim,
2145 Number>::value_type &
2147 const unsigned int point_index) const
2148{
2149 AssertIndexRange(point_index, values.size());
2150 return values[point_index];
2151}
2152
2153
2154
2155template <int n_components_, int dim, int spacedim, typename Number>
2156inline const typename FEPointEvaluationBase<n_components_,
2157 dim,
2158 spacedim,
2159 Number>::gradient_type &
2161 const unsigned int point_index) const
2162{
2163 AssertIndexRange(point_index, gradients.size());
2164 return gradients[point_index];
2165}
2166
2167
2168
2169template <int n_components_, int dim, int spacedim, typename Number>
2170inline Number
2172 const unsigned int point_index) const
2173{
2174 static_assert(n_components == dim,
2175 "Only makes sense for a vector field with dim components");
2176
2177 AssertIndexRange(point_index, values.size());
2178 return trace(gradients[point_index]);
2179}
2180
2181
2182
2183template <int n_components_, int dim, int spacedim, typename Number>
2184inline void
2186 const value_type &value,
2187 const unsigned int point_index)
2188{
2189 AssertIndexRange(point_index, n_q_points);
2190 values[point_index] = value;
2191}
2192
2193
2194
2195template <int n_components_, int dim, int spacedim, typename Number>
2196inline void
2198 const gradient_type &gradient,
2199 const unsigned int point_index)
2200{
2201 AssertIndexRange(point_index, n_q_points);
2202 gradients[point_index] = gradient;
2203}
2204
2205
2206
2207template <int n_components_, int dim, int spacedim, typename Number>
2208inline void
2210 const Number &value,
2211 const unsigned int point_index)
2212{
2213 static_assert(n_components == dim,
2214 "Only makes sense for a vector field with dim components");
2215
2216 AssertIndexRange(point_index, n_q_points);
2217 gradients[point_index] = gradient_type();
2218 for (unsigned int d = 0; d < dim; ++d)
2219 gradients[point_index][d][d] = value;
2220}
2221
2222
2223
2224template <int n_components_, int dim, int spacedim, typename Number>
2225Tensor<1, (dim == 2 ? 1 : dim), Number>
2227 const unsigned int point_index) const
2228{
2229 static_assert(
2230 dim > 1 && n_components == dim,
2231 "Only makes sense for a vector field with dim components and dim > 1");
2232
2233 const Tensor<2, dim, Number> grad = get_gradient(point_index);
2234 Tensor<1, (dim == 2 ? 1 : dim), Number> curl;
2235 switch (dim)
2236 {
2237 case 2:
2238 curl[0] = grad[1][0] - grad[0][1];
2239 break;
2240 case 3:
2241 curl[0] = grad[2][1] - grad[1][2];
2242 curl[1] = grad[0][2] - grad[2][0];
2243 curl[2] = grad[1][0] - grad[0][1];
2244 break;
2245 default:
2247 }
2248 return curl;
2249}
2250
2251
2252
2253template <int n_components_, int dim, int spacedim, typename Number>
2256 const unsigned int point_index) const
2257{
2258 AssertIndexRange(point_index, n_q_points);
2259 Assert(jacobian_ptr != nullptr,
2261 ExcFEPointEvaluationAccessToUninitializedMappingField(
2262 "update_jacobians"));
2263 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2264 GeometryType::affine ?
2265 0 :
2266 point_index];
2267}
2268
2269
2270
2271template <int n_components_, int dim, int spacedim, typename Number>
2274 const unsigned int point_index) const
2275{
2276 AssertIndexRange(point_index, n_q_points);
2277 Assert(inverse_jacobian_ptr != nullptr,
2279 ExcFEPointEvaluationAccessToUninitializedMappingField(
2280 "update_inverse_jacobians"));
2281 return inverse_jacobian_ptr
2282 [cell_type <=
2284 0 :
2285 point_index];
2286}
2287
2288
2289
2290template <int n_components_, int dim, int spacedim, typename Number>
2291inline Number
2293 const unsigned int point_index) const
2294{
2295 AssertIndexRange(point_index, n_q_points);
2296 Assert(JxW_ptr != nullptr,
2298 ExcFEPointEvaluationAccessToUninitializedMappingField(
2299 "update_JxW_values"));
2300 return JxW_ptr[point_index];
2301}
2302
2303
2304
2305template <int n_components_, int dim, int spacedim, typename Number>
2308 const unsigned int point_index) const
2309{
2310 return quadrature_point(point_index);
2311}
2312
2313
2314
2315template <int n_components_, int dim, int spacedim, typename Number>
2318 const unsigned int point_index) const
2319{
2320 AssertIndexRange(point_index, n_q_points);
2321 Assert(real_point_ptr != nullptr,
2323 ExcFEPointEvaluationAccessToUninitializedMappingField(
2324 "update_quadrature_points"));
2325 return real_point_ptr[point_index];
2326}
2327
2328
2329
2330template <int n_components_, int dim, int spacedim, typename Number>
2331inline Point<dim, Number>
2333 const unsigned int point_index) const
2334{
2335 AssertIndexRange(point_index, n_q_points);
2336 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2337 Point<dim, Number> unit_point;
2338 for (unsigned int d = 0; d < dim; ++d)
2340 unit_point_ptr[point_index / stride][d], point_index % stride);
2341 return unit_point;
2342}
2343
2344
2345
2346template <int n_components_, int dim, int spacedim, typename Number>
2354
2355
2356
2357template <int n_components_, int dim, int spacedim, typename Number>
2361 const unsigned int first_selected_component)
2362 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2363 mapping_info,
2364 fe,
2365 first_selected_component)
2366{}
2367
2368
2369
2370template <int n_components_, int dim, int spacedim, typename Number>
2372 const Mapping<dim, spacedim> &mapping,
2374 const UpdateFlags update_flags,
2375 const unsigned int first_selected_component)
2376 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2377 mapping,
2378 fe,
2379 update_flags,
2380 first_selected_component)
2381{}
2382
2383
2384
2385template <int n_components_, int dim, int spacedim, typename Number>
2386inline void
2389{
2390 this->current_cell_index = numbers::invalid_unsigned_int;
2391 this->current_face_number = numbers::invalid_unsigned_int;
2392
2393 if (this->use_linear_path)
2394 this->template do_reinit<false, true>();
2395 else
2396 this->template do_reinit<false, false>();
2397}
2398
2399
2400
2401template <int n_components_, int dim, int spacedim, typename Number>
2402inline void
2404{
2405 internal_reinit_single_cell_state_mapping_info();
2406 this->must_reinitialize_pointers = false;
2407}
2408
2409
2410
2411template <int n_components_, int dim, int spacedim, typename Number>
2412inline void
2415 const ArrayView<const Point<dim>> &unit_points)
2416{
2417 // reinit is only allowed for mapping computation on the fly
2418 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2420
2421 this->mapping_info_on_the_fly->reinit(cell, unit_points);
2422 this->must_reinitialize_pointers = false;
2423
2424 if (!this->fast_path)
2425 {
2426 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2427 *this->mapping,
2428 *this->fe,
2430 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2431 this->update_flags);
2432 this->fe_values->reinit(cell);
2433 }
2434
2435 if (this->use_linear_path)
2436 this->template do_reinit<false, true>();
2437 else
2438 this->template do_reinit<false, false>();
2439}
2440
2441
2442
2443template <int n_components_, int dim, int spacedim, typename Number>
2444inline void
2446 const unsigned int cell_index)
2447{
2448 this->current_cell_index = cell_index;
2449 this->current_face_number = numbers::invalid_unsigned_int;
2450 this->must_reinitialize_pointers = false;
2451
2452 if (this->use_linear_path)
2453 this->template do_reinit<false, true>();
2454 else
2455 this->template do_reinit<false, false>();
2456
2457 if (!this->fast_path)
2458 {
2459 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2460
2461 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2462 for (unsigned int d = 0; d < dim; ++d)
2463 unit_points[v][d] =
2464 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2465
2466 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2467 *this->mapping,
2468 *this->fe,
2470 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2471 this->update_flags);
2472
2473 this->fe_values->reinit(
2474 this->mapping_info->get_cell_iterator(this->current_cell_index));
2475 }
2476}
2477
2478
2479
2480template <int n_components_, int dim, int spacedim, typename Number>
2481template <std::size_t stride_view>
2482void
2485 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2486{
2487 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2488
2489 if (!((evaluation_flags & EvaluationFlags::values) ||
2490 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2491 return;
2492
2493 if (this->must_reinitialize_pointers)
2494 internal_reinit_single_cell_state_mapping_info();
2495
2496 if (this->n_q_points == 0)
2497 return;
2498
2499 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2500 if (this->fast_path)
2501 {
2502 if (this->use_linear_path)
2503 evaluate_fast<true>(solution_values, evaluation_flags);
2504 else
2505 evaluate_fast<false>(solution_values, evaluation_flags);
2506 }
2507 else
2508 evaluate_slow(solution_values, evaluation_flags);
2509}
2510
2511
2512
2513template <int n_components_, int dim, int spacedim, typename Number>
2514void
2516 const ArrayView<const ScalarNumber> &solution_values,
2517 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2518{
2519 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2520 solution_values.size()),
2521 evaluation_flags);
2522}
2523
2524
2525
2526template <int n_components_, int dim, int spacedim, typename Number>
2527template <std::size_t stride_view>
2528void
2530 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2531 const EvaluationFlags::EvaluationFlags &integration_flags,
2532 const bool sum_into_values)
2533{
2534 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2535}
2536
2537
2538
2539template <int n_components_, int dim, int spacedim, typename Number>
2540void
2542 const ArrayView<ScalarNumber> &solution_values,
2543 const EvaluationFlags::EvaluationFlags &integration_flags,
2544 const bool sum_into_values)
2545{
2546 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2547 solution_values.size()),
2548 integration_flags,
2549 sum_into_values);
2550}
2551
2552
2553
2554template <int n_components_, int dim, int spacedim, typename Number>
2556 scalar_value_type
2558 const
2559{
2560 value_type return_value = {};
2561
2562 for (const auto point_index : this->quadrature_point_indices())
2563 return_value += values[point_index] * this->JxW(point_index);
2564
2565 return ETT::sum_value(return_value);
2566}
2567
2568
2569
2570template <int n_components_, int dim, int spacedim, typename Number>
2571unsigned int
2574{
2575 Assert(stride == 1,
2576 ExcMessage(
2577 "Calling this function only makes sense in fully vectorized mode."));
2578 if (q == n_q_batches - 1)
2579 {
2580 const unsigned int n_filled_lanes =
2581 n_q_points_scalar & (n_lanes_user_interface - 1);
2582 if (n_filled_lanes == 0)
2583 return n_lanes_user_interface;
2584 else
2585 return n_filled_lanes;
2586 }
2587 else
2588 return n_lanes_user_interface;
2589}
2590
2591
2592
2593template <int n_components_, int dim, int spacedim, typename Number>
2594template <std::size_t stride_view>
2595void
2597 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2598 const EvaluationFlags::EvaluationFlags &integration_flags,
2599 const bool sum_into_values)
2600{
2601 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2602}
2603
2604
2605
2606template <int n_components_, int dim, int spacedim, typename Number>
2607void
2609 const ArrayView<ScalarNumber> &solution_values,
2610 const EvaluationFlags::EvaluationFlags &integration_flags,
2611 const bool sum_into_values)
2612{
2613 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2614 solution_values.size()),
2615 integration_flags,
2616 sum_into_values);
2617}
2618
2619
2620
2621template <int n_components_, int dim, int spacedim, typename Number>
2622template <bool is_linear, std::size_t stride_view>
2623inline void
2626{
2627 const unsigned int dofs_per_comp =
2628 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2629
2630 for (unsigned int comp = 0; comp < n_components; ++comp)
2631 {
2632 const std::size_t offset =
2633 (this->component_in_base_element + comp) * dofs_per_comp;
2634
2635 if ((is_linear && n_components == 1) || this->renumber.empty())
2636 {
2637 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2638 ETT::read_value(solution_values[i + offset],
2639 comp,
2640 this->solution_renumbered[i]);
2641 }
2642 else
2643 {
2644 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2645 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2646 ETT::read_value(solution_values[renumber_ptr[i]],
2647 comp,
2648 this->solution_renumbered[i]);
2649 }
2650 }
2651}
2652
2653
2654
2655template <int n_components_, int dim, int spacedim, typename Number>
2656template <bool is_linear, std::size_t stride_view>
2657inline void
2660 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2661 const unsigned int n_shapes,
2662 const unsigned int qb,
2663 vectorized_value_type &value,
2665{
2666 if (evaluation_flags & EvaluationFlags::gradients)
2667 {
2668 std::array<vectorized_value_type, dim + 1> result;
2669 if constexpr (is_linear)
2670 {
2671 if constexpr (n_components == 1)
2672 result =
2674 dim,
2677 1,
2678 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2679 else
2680 result =
2682 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2683 }
2684 else
2686 dim,
2689 1,
2690 false>(this->shapes.data() + qb * n_shapes,
2691 n_shapes,
2692 this->solution_renumbered.data());
2693 gradient[0] = result[0];
2694 if (dim > 1)
2695 gradient[1] = result[1];
2696 if (dim > 2)
2697 gradient[2] = result[2];
2698 value = result[dim];
2699 }
2700 else
2701 {
2702 if constexpr (is_linear)
2703 {
2704 if constexpr (n_components == 1)
2706 dim,
2709 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2710 else
2712 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2713 }
2714 else
2715 value =
2719 false>(
2720 this->shapes.data() + qb * n_shapes,
2721 n_shapes,
2722 this->solution_renumbered.data());
2723 }
2724}
2725
2726
2727
2728template <int n_components_, int dim, int spacedim, typename Number>
2729template <bool is_linear, std::size_t stride_view>
2730inline void
2733 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2734{
2735 if (!(is_linear && n_components == 1))
2736 prepare_evaluate_fast<is_linear>(solution_values);
2737
2738 // loop over quadrature batches qb
2739 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2740
2741 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2742 {
2745
2746 compute_evaluate_fast<is_linear>(
2747 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2748
2749 if (evaluation_flags & EvaluationFlags::values)
2750 {
2751 for (unsigned int v = 0, offset = qb * stride;
2752 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2753 ++v, ++offset)
2754 ETT::set_value(value, v, this->values[offset]);
2755 }
2756 if (evaluation_flags & EvaluationFlags::gradients)
2757 {
2758 Assert(this->update_flags & update_gradients ||
2759 this->update_flags & update_inverse_jacobians,
2761
2762 for (unsigned int v = 0, offset = qb * stride;
2763 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2764 ++v, ++offset)
2765 {
2766 unit_gradient_type unit_gradient;
2767 ETT::set_gradient(gradient, v, unit_gradient);
2768 this->gradients[offset] =
2769 this->cell_type <=
2772 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2774 this
2775 ->inverse_jacobian_ptr[this->cell_type <=
2777 GeometryType::affine ?
2778 0 :
2779 offset]
2780 .transpose(),
2781 unit_gradient);
2782 }
2783 }
2784 }
2785}
2786
2787
2788
2789template <int n_components_, int dim, int spacedim, typename Number>
2790template <std::size_t stride_view>
2791inline void
2794 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2795{
2796 // slow path with FEValues
2797 Assert(this->fe_values.get() != nullptr,
2798 ExcMessage(
2799 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2800
2801 const std::size_t n_points = this->fe_values->get_quadrature().size();
2802
2803 if (evaluation_flags & EvaluationFlags::values)
2804 {
2805 this->values.resize(this->n_q_points);
2806 std::fill(this->values.begin(), this->values.end(), value_type());
2807 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2808 {
2809 const ScalarNumber value = solution_values[i];
2810 for (unsigned int d = 0; d < n_components; ++d)
2811 if (this->nonzero_shape_function_component[i][d] &&
2812 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2813 for (unsigned int qb = 0, q = 0; q < n_points;
2814 ++qb, q += n_lanes_user_interface)
2815 for (unsigned int v = 0;
2816 v < n_lanes_user_interface && q + v < n_points;
2817 ++v)
2818 ETT::access(this->values[qb],
2819 v,
2820 d,
2821 this->fe_values->shape_value(i, q + v) * value);
2822 else if (this->nonzero_shape_function_component[i][d])
2823 for (unsigned int qb = 0, q = 0; q < n_points;
2824 ++qb, q += n_lanes_user_interface)
2825 for (unsigned int v = 0;
2826 v < n_lanes_user_interface && q + v < n_points;
2827 ++v)
2828 ETT::access(this->values[qb],
2829 v,
2830 d,
2831 this->fe_values->shape_value_component(i,
2832 q + v,
2833 d) *
2834 value);
2835 }
2836 }
2837
2838 if (evaluation_flags & EvaluationFlags::gradients)
2839 {
2840 this->gradients.resize(this->n_q_points);
2841 std::fill(this->gradients.begin(),
2842 this->gradients.end(),
2843 gradient_type());
2844 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2845 {
2846 const ScalarNumber value = solution_values[i];
2847 for (unsigned int d = 0; d < n_components; ++d)
2848 if (this->nonzero_shape_function_component[i][d] &&
2849 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2850 for (unsigned int qb = 0, q = 0; q < n_points;
2851 ++qb, q += n_lanes_user_interface)
2852 for (unsigned int v = 0;
2853 v < n_lanes_user_interface && q + v < n_points;
2854 ++v)
2855 ETT::access(this->gradients[qb],
2856 v,
2857 d,
2858 this->fe_values->shape_grad(i, q + v) * value);
2859 else if (this->nonzero_shape_function_component[i][d])
2860 for (unsigned int qb = 0, q = 0; q < n_points;
2861 ++qb, q += n_lanes_user_interface)
2862 for (unsigned int v = 0;
2863 v < n_lanes_user_interface && q + v < n_points;
2864 ++v)
2865 ETT::access(
2866 this->gradients[qb],
2867 v,
2868 d,
2869 this->fe_values->shape_grad_component(i, q + v, d) * value);
2870 }
2871 }
2872}
2873
2874
2875
2876template <int n_components_, int dim, int spacedim, typename Number>
2877template <bool is_linear>
2878inline void
2880 const EvaluationFlags::EvaluationFlags &integration_flags,
2881 const unsigned int n_shapes,
2882 const unsigned int qb,
2883 const vectorized_value_type value,
2885 vectorized_value_type *solution_values_vectorized_linear)
2886{
2887 if (integration_flags & EvaluationFlags::gradients)
2889 is_linear,
2890 dim,
2892 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2893 n_shapes,
2894 &value,
2895 gradient,
2896 is_linear ?
2897 solution_values_vectorized_linear :
2898 this->solution_renumbered_vectorized.data(),
2899 this->unit_point_ptr[qb],
2900 qb != 0);
2901 else
2903 dim,
2906 this->shapes.data() + qb * n_shapes,
2907 n_shapes,
2908 value,
2909 is_linear ? solution_values_vectorized_linear :
2910 this->solution_renumbered_vectorized.data(),
2911 this->unit_point_ptr[qb],
2912 qb != 0);
2913}
2914
2915
2916
2917template <int n_components_, int dim, int spacedim, typename Number>
2918template <bool is_linear, std::size_t stride_view>
2919inline void
2921 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2922 vectorized_value_type *solution_values_vectorized_linear,
2923 const bool sum_into_values)
2924{
2925 if (!sum_into_values && this->fe->n_components() > n_components)
2926 for (unsigned int i = 0; i < solution_values.size(); ++i)
2927 solution_values[i] = 0;
2928
2929 const unsigned int dofs_per_comp =
2930 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2931
2932 for (unsigned int comp = 0; comp < n_components; ++comp)
2933 {
2934 const std::size_t offset =
2935 (this->component_in_base_element + comp) * dofs_per_comp;
2936
2937 if (is_linear || this->renumber.empty())
2938 {
2939 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2940 if (sum_into_values)
2941 solution_values[i + offset] +=
2942 ETT::sum_value(comp,
2943 is_linear ?
2944 *(solution_values_vectorized_linear + i) :
2945 this->solution_renumbered_vectorized[i]);
2946 else
2947 solution_values[i + offset] =
2948 ETT::sum_value(comp,
2949 is_linear ?
2950 *(solution_values_vectorized_linear + i) :
2951 this->solution_renumbered_vectorized[i]);
2952 }
2953 else
2954 {
2955 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2956 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2957 if (sum_into_values)
2958 solution_values[renumber_ptr[i]] +=
2959 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2960 else
2961 solution_values[renumber_ptr[i]] =
2962 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2963 }
2964 }
2965}
2966
2967
2968
2969template <int n_components_, int dim, int spacedim, typename Number>
2970template <bool do_JxW, bool is_linear, std::size_t stride_view>
2971inline void
2973 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2974 const EvaluationFlags::EvaluationFlags &integration_flags,
2975 const bool sum_into_values)
2976{
2977 // zero out lanes of incomplete last quadrature point batch
2978 if constexpr (stride == 1)
2979 if (const unsigned int n_filled_lanes =
2980 this->n_q_points_scalar & (n_lanes_internal - 1);
2981 n_filled_lanes > 0)
2982 {
2983 if (integration_flags & EvaluationFlags::values)
2984 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2985 ETT::set_zero_value(this->values.back(), v);
2986 if (integration_flags & EvaluationFlags::gradients)
2987 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2988 ETT::set_zero_gradient(this->gradients.back(), v);
2989 }
2990
2991 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2992 solution_values_vectorized_linear = {};
2993
2994 // loop over quadrature batches qb
2995 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2996
2997 const bool cartesian_cell =
2999 const bool affine_cell =
3001 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3002 {
3003 vectorized_value_type value = {};
3005
3006 if (integration_flags & EvaluationFlags::values)
3007 for (unsigned int v = 0, offset = qb * stride;
3008 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3009 ++v, ++offset)
3010 ETT::get_value(value,
3011 v,
3012 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3013 this->values[offset]);
3014
3015 if (integration_flags & EvaluationFlags::gradients)
3016 for (unsigned int v = 0, offset = qb * stride;
3017 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3018 ++v, ++offset)
3019 {
3020 const gradient_type grad_w =
3021 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3022 this->gradients[offset];
3023 ETT::get_gradient(
3024 gradient,
3025 v,
3026 cartesian_cell ?
3027 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3028 grad_w) :
3030 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3031 grad_w));
3032 }
3033
3034 compute_integrate_fast<is_linear>(
3035 integration_flags,
3036 n_shapes,
3037 qb,
3038 value,
3039 gradient,
3040 solution_values_vectorized_linear.data());
3041 }
3042
3043 // add between the lanes and write into the result
3044 finish_integrate_fast<is_linear>(solution_values,
3045 solution_values_vectorized_linear.data(),
3046 sum_into_values);
3047}
3048
3049
3050
3051template <int n_components_, int dim, int spacedim, typename Number>
3052template <bool do_JxW, std::size_t stride_view>
3053inline void
3055 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3056 const EvaluationFlags::EvaluationFlags &integration_flags,
3057 const bool sum_into_values)
3058{
3059 // slow path with FEValues
3060 Assert(this->fe_values.get() != nullptr,
3061 ExcMessage(
3062 "Not initialized. Please call FEPointEvaluation::reinit()!"));
3063 if (!sum_into_values)
3064 for (unsigned int i = 0; i < solution_values.size(); ++i)
3065 solution_values[i] = 0;
3066
3067 const std::size_t n_points = this->fe_values->get_quadrature().size();
3068
3069 if (integration_flags & EvaluationFlags::values)
3070 {
3071 AssertIndexRange(this->n_q_points, this->values.size() + 1);
3072 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3073 {
3074 for (unsigned int d = 0; d < n_components; ++d)
3075 if (this->nonzero_shape_function_component[i][d] &&
3076 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3077 for (unsigned int qb = 0, q = 0; q < n_points;
3078 ++qb, q += n_lanes_user_interface)
3079 for (unsigned int v = 0;
3080 v < n_lanes_user_interface && q + v < n_points;
3081 ++v)
3082 solution_values[i] +=
3083 this->fe_values->shape_value(i, q + v) *
3084 ETT::access(this->values[qb], v, d) *
3085 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3086 else if (this->nonzero_shape_function_component[i][d])
3087 for (unsigned int qb = 0, q = 0; q < n_points;
3088 ++qb, q += n_lanes_user_interface)
3089 for (unsigned int v = 0;
3090 v < n_lanes_user_interface && q + v < n_points;
3091 ++v)
3092 solution_values[i] +=
3093 this->fe_values->shape_value_component(i, q + v, d) *
3094 ETT::access(this->values[qb], v, d) *
3095 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3096 }
3097 }
3098
3099 if (integration_flags & EvaluationFlags::gradients)
3100 {
3101 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
3102 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3103 {
3104 for (unsigned int d = 0; d < n_components; ++d)
3105 if (this->nonzero_shape_function_component[i][d] &&
3106 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3107 for (unsigned int qb = 0, q = 0; q < n_points;
3108 ++qb, q += n_lanes_user_interface)
3109 for (unsigned int v = 0;
3110 v < n_lanes_user_interface && q + v < n_points;
3111 ++v)
3112 solution_values[i] +=
3113 this->fe_values->shape_grad(i, q + v) *
3114 ETT::access(this->gradients[qb], v, d) *
3115 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3116 else if (this->nonzero_shape_function_component[i][d])
3117 for (unsigned int qb = 0, q = 0; q < n_points;
3118 ++qb, q += n_lanes_user_interface)
3119 for (unsigned int v = 0;
3120 v < n_lanes_user_interface && q + v < n_points;
3121 ++v)
3122 solution_values[i] +=
3123 this->fe_values->shape_grad_component(i, q + v, d) *
3124 ETT::access(this->gradients[qb], v, d) *
3125 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3126 }
3127 }
3128}
3129
3130
3131
3132template <int n_components_, int dim, int spacedim, typename Number>
3133template <bool do_JxW, std::size_t stride_view>
3134void
3136 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3137 const EvaluationFlags::EvaluationFlags &integration_flags,
3138 const bool sum_into_values)
3139{
3140 if (this->must_reinitialize_pointers)
3141 internal_reinit_single_cell_state_mapping_info();
3142
3143 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3144
3145 if (this->n_q_points == 0 || // no evaluation points provided
3146 !((integration_flags & EvaluationFlags::values) ||
3147 (integration_flags &
3148 EvaluationFlags::gradients))) // no integration flags
3149 {
3150 if (!sum_into_values)
3151 for (unsigned int i = 0; i < solution_values.size(); ++i)
3152 solution_values[i] = 0;
3153 return;
3154 }
3155
3156 Assert(
3157 !do_JxW || this->JxW_ptr != nullptr,
3158 ExcMessage(
3159 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3160
3161 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3162 if (this->fast_path)
3163 {
3164 if (this->use_linear_path)
3165 integrate_fast<do_JxW, true>(solution_values,
3166 integration_flags,
3167 sum_into_values);
3168 else
3169 integrate_fast<do_JxW, false>(solution_values,
3170 integration_flags,
3171 sum_into_values);
3172 }
3173 else
3174 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3175}
3176
3177
3178
3179template <int n_components_, int dim, int spacedim, typename Number>
3182 const unsigned int point_index) const
3183{
3184 AssertIndexRange(point_index, this->n_q_points);
3185 Assert(this->normal_ptr != nullptr,
3187 ExcFEPointEvaluationAccessToUninitializedMappingField(
3188 "update_normal_vectors"));
3189 return this->normal_ptr[point_index];
3190}
3191
3192
3193
3194template <int n_components_, int dim, int spacedim, typename Number>
3195inline const typename FEPointEvaluation<n_components_,
3196 dim,
3197 spacedim,
3198 Number>::value_type
3200 const unsigned int point_index) const
3201{
3202 AssertIndexRange(point_index, this->gradients.size());
3203
3204 value_type normal_derivative;
3205 if constexpr (n_components == 1)
3206 normal_derivative =
3207 this->gradients[point_index] * normal_vector(point_index);
3208 else
3209 for (unsigned int comp = 0; comp < n_components; ++comp)
3210 normal_derivative[comp] =
3211 this->gradients[point_index][comp] * normal_vector(point_index);
3212
3213 return normal_derivative;
3214}
3215
3216
3217
3218template <int n_components_, int dim, int spacedim, typename Number>
3219inline void
3222 const unsigned int point_index)
3223{
3224 AssertIndexRange(point_index, this->gradients.size());
3225 if constexpr (n_components == 1)
3226 this->gradients[point_index] = value * normal_vector(point_index);
3227 else
3228 for (unsigned int comp = 0; comp < n_components; ++comp)
3229 this->gradients[point_index][comp] =
3230 value[comp] * normal_vector(point_index);
3231}
3232
3233
3234
3235template <int n_components_, int dim, int spacedim, typename Number>
3240 const bool is_interior,
3241 const unsigned int first_selected_component)
3242 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3243 mapping_info,
3244 fe,
3245 first_selected_component,
3246 is_interior)
3247{}
3248
3249
3250
3251template <int n_components_, int dim, int spacedim, typename Number>
3252inline void
3254 const unsigned int cell_index,
3255 const unsigned int face_number)
3256{
3257 this->current_cell_index = cell_index;
3258 this->current_face_number = face_number;
3259 this->must_reinitialize_pointers = false;
3260
3261 if (this->use_linear_path)
3262 this->template do_reinit<true, true>();
3263 else
3264 this->template do_reinit<true, false>();
3265}
3266
3267
3268
3269template <int n_components_, int dim, int spacedim, typename Number>
3270inline void
3272 const unsigned int face_index)
3273{
3274 this->current_cell_index = face_index;
3275 this->current_face_number =
3276 this->mapping_info->get_face_number(face_index, this->is_interior);
3277 this->must_reinitialize_pointers = false;
3278
3279 if (this->use_linear_path)
3280 this->template do_reinit<true, true>();
3281 else
3282 this->template do_reinit<true, false>();
3283}
3284
3285
3286
3287template <int n_components_, int dim, int spacedim, typename Number>
3288template <std::size_t stride_view>
3289void
3292 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3293{
3294 Assert(!this->must_reinitialize_pointers,
3295 ExcMessage("Object has not been reinitialized!"));
3296
3297 if (this->n_q_points == 0)
3298 return;
3299
3300 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3301
3302 if (!((evaluation_flags & EvaluationFlags::values) ||
3303 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3304 return;
3305
3306 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3307
3308 if (this->use_linear_path)
3309 do_evaluate<true>(solution_values, evaluation_flags);
3310 else
3311 do_evaluate<false>(solution_values, evaluation_flags);
3312}
3313
3314
3315
3316template <int n_components_, int dim, int spacedim, typename Number>
3317void
3319 const ArrayView<const ScalarNumber> &solution_values,
3320 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3321{
3322 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3323 solution_values.size()),
3324 evaluation_flags);
3325}
3326
3327
3328
3329template <int n_components_, int dim, int spacedim, typename Number>
3330template <bool is_linear, std::size_t stride_view>
3331void
3334 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3335{
3336 const unsigned int dofs_per_comp =
3337 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3338
3339 const ScalarNumber *input;
3340 if (stride_view == 1 && this->component_in_base_element == 0 &&
3341 (is_linear || this->renumber.empty()))
3342 input = solution_values.data();
3343 else
3344 {
3345 for (unsigned int comp = 0; comp < n_components; ++comp)
3346 {
3347 const std::size_t offset =
3348 (this->component_in_base_element + comp) * dofs_per_comp;
3349
3350 if (is_linear || this->renumber.empty())
3351 {
3352 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3353 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3354 solution_values[i + offset];
3355 }
3356 else
3357 {
3358 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3359 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3360 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3361 solution_values[renumber_ptr[i]];
3362 }
3363 }
3364 input = this->scratch_data_scalar.data();
3365 }
3366
3367 ScalarNumber *output =
3368 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3369
3370 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3371 template interpolate<true, false>(n_components,
3372 evaluation_flags,
3373 this->shape_info,
3374 input,
3375 output,
3376 this->current_face_number);
3377
3378 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3379}
3380
3381
3382
3383template <int n_components_, int dim, int spacedim, typename Number>
3384template <std::size_t stride_view>
3385void
3387 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3388 const EvaluationFlags::EvaluationFlags &integration_flags,
3389 const bool sum_into_values)
3390{
3391 Assert(!this->must_reinitialize_pointers,
3392 ExcMessage("Object has not been reinitialized!"));
3393
3394 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3395
3396 if (this->n_q_points == 0 || // no evaluation points provided
3397 !((integration_flags & EvaluationFlags::values) ||
3398 (integration_flags &
3399 EvaluationFlags::gradients))) // no integration flags
3400 {
3401 if (!sum_into_values)
3402 for (unsigned int i = 0; i < solution_values.size(); ++i)
3403 solution_values[i] = 0;
3404 return;
3405 }
3406
3407 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3408
3409 if (this->use_linear_path)
3410 do_integrate<true, true>(solution_values,
3411 integration_flags,
3412 sum_into_values);
3413 else
3414 do_integrate<true, false>(solution_values,
3415 integration_flags,
3416 sum_into_values);
3417}
3418
3419
3420
3421template <int n_components_, int dim, int spacedim, typename Number>
3422void
3424 const ArrayView<ScalarNumber> &solution_values,
3425 const EvaluationFlags::EvaluationFlags &integration_flags,
3426 const bool sum_into_values)
3427{
3428 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3429 solution_values.size()),
3430 integration_flags,
3431 sum_into_values);
3432}
3433
3434
3435
3436template <int n_components_, int dim, int spacedim, typename Number>
3437template <std::size_t stride_view>
3438void
3440 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3441 const EvaluationFlags::EvaluationFlags &integration_flags,
3442 const bool sum_into_values)
3443{
3444 Assert(!this->must_reinitialize_pointers,
3445 ExcMessage("Object has not been reinitialized!"));
3446
3447 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3448
3449 if (this->n_q_points == 0 || // no evaluation points provided
3450 !((integration_flags & EvaluationFlags::values) ||
3451 (integration_flags &
3452 EvaluationFlags::gradients))) // no integration flags
3453 {
3454 if (!sum_into_values)
3455 for (unsigned int i = 0; i < solution_values.size(); ++i)
3456 solution_values[i] = 0;
3457 return;
3458 }
3459
3460 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3461
3462 if (this->use_linear_path)
3463 do_integrate<false, true>(solution_values,
3464 integration_flags,
3465 sum_into_values);
3466 else
3467 do_integrate<false, false>(solution_values,
3468 integration_flags,
3469 sum_into_values);
3470}
3471
3472
3473
3474template <int n_components_, int dim, int spacedim, typename Number>
3475void
3477 const ArrayView<ScalarNumber> &solution_values,
3478 const EvaluationFlags::EvaluationFlags &integration_flags,
3479 const bool sum_into_values)
3480{
3481 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3482 solution_values.size()),
3483 integration_flags,
3484 sum_into_values);
3485}
3486
3487
3488
3489template <int n_components_, int dim, int spacedim, typename Number>
3490template <bool do_JxW, bool is_linear, std::size_t stride_view>
3491void
3493 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3494 const EvaluationFlags::EvaluationFlags &integration_flags,
3495 const bool sum_into_values)
3496{
3497 if (!sum_into_values && this->fe->n_components() > n_components)
3498 for (unsigned int i = 0; i < solution_values.size(); ++i)
3499 solution_values[i] = 0;
3500
3501 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3502 integration_flags,
3503 false);
3504
3505 ScalarNumber *input = this->scratch_data_scalar.begin();
3506
3507 if (stride_view == 1 && this->component_in_base_element == 0 &&
3508 (is_linear || this->renumber.empty()))
3509 {
3510 if (sum_into_values)
3511 internal::
3512 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3513 template interpolate<false, true>(n_components,
3514 integration_flags,
3515 this->shape_info,
3516 input,
3517 solution_values.data(),
3518 this->current_face_number);
3519 else
3520 internal::
3521 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3522 template interpolate<false, false>(n_components,
3523 integration_flags,
3524 this->shape_info,
3525 input,
3526 solution_values.data(),
3527 this->current_face_number);
3528 }
3529 else
3530 {
3531 const unsigned int dofs_per_comp_face =
3532 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3533
3534 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3535 ScalarNumber *output = input + size_input;
3536
3537 internal::
3538 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3539 template interpolate<false, false>(n_components,
3540 integration_flags,
3541 this->shape_info,
3542 input,
3543 output,
3544 this->current_face_number);
3545
3546 const unsigned int dofs_per_comp =
3547 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3548
3549 for (unsigned int comp = 0; comp < n_components; ++comp)
3550 {
3551 const std::size_t offset =
3552 (this->component_in_base_element + comp) * dofs_per_comp;
3553
3554 if (is_linear || this->renumber.empty())
3555 {
3556 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3557 if (sum_into_values)
3558 solution_values[i + offset] +=
3559 output[i + comp * dofs_per_comp];
3560 else
3561 solution_values[i + offset] =
3562 output[i + comp * dofs_per_comp];
3563 }
3564 else
3565 {
3566 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3567 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3568 if (sum_into_values)
3569 solution_values[renumber_ptr[i]] +=
3570 output[i + comp * dofs_per_comp];
3571 else
3572 solution_values[renumber_ptr[i]] =
3573 output[i + comp * dofs_per_comp];
3574 }
3575 }
3576 }
3577}
3578
3579
3580
3581template <int n_components_, int dim, int spacedim, typename Number>
3582template <int stride_face_dof>
3583void
3585 const ScalarNumber *face_dof_values,
3586 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3587{
3588 if (this->use_linear_path)
3589 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3590 evaluation_flags);
3591 else
3592 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3593 evaluation_flags);
3594}
3595
3596
3597
3598template <int n_components_, int dim, int spacedim, typename Number>
3599template <bool is_linear, int stride_face_dof>
3600inline void
3602 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3603 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3604{
3605 const scalar_value_type *face_dof_values_ptr;
3606 if constexpr (n_components == 1)
3607 face_dof_values_ptr = face_dof_values;
3608 else
3609 {
3610 const unsigned int dofs_per_comp_face =
3611 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3612 for (unsigned int comp = 0; comp < n_components; ++comp)
3613 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3614 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3615 stride_face_dof],
3616 comp,
3617 this->solution_renumbered[i]);
3618
3619 face_dof_values_ptr = this->solution_renumbered.data();
3620 }
3621
3622 constexpr int stride_face_dof_actual =
3623 n_components == 1 ? stride_face_dof : 1;
3624
3625 // loop over quadrature batches qb
3626 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3627
3628 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3629 {
3632
3633 if (evaluation_flags & EvaluationFlags::gradients)
3634 {
3635 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3636 is_linear ?
3638 dim - 1,
3641 2,
3642 stride_face_dof_actual>(face_dof_values_ptr,
3643 this->unit_point_faces_ptr[qb]) :
3645 dim - 1,
3648 2,
3649 false,
3650 stride_face_dof_actual>(this->shapes_faces.data() +
3651 qb * n_shapes,
3652 n_shapes,
3653 face_dof_values_ptr);
3654
3655 value = interpolated_value[dim - 1];
3656 // reorder derivative from tangential/normal derivatives into tensor
3657 // in physical coordinates
3658 if (this->current_face_number / 2 == 0)
3659 {
3660 gradient[0] = interpolated_value[dim];
3661 if (dim > 1)
3662 gradient[1] = interpolated_value[0];
3663 if (dim > 2)
3664 gradient[2] = interpolated_value[1];
3665 }
3666 else if (this->current_face_number / 2 == 1)
3667 {
3668 if (dim > 1)
3669 gradient[1] = interpolated_value[dim];
3670 if (dim == 3)
3671 {
3672 gradient[0] = interpolated_value[1];
3673 gradient[2] = interpolated_value[0];
3674 }
3675 else if (dim == 2)
3676 gradient[0] = interpolated_value[0];
3677 else
3679 }
3680 else if (this->current_face_number / 2 == 2)
3681 {
3682 if (dim > 2)
3683 {
3684 gradient[0] = interpolated_value[0];
3685 gradient[1] = interpolated_value[1];
3686 gradient[2] = interpolated_value[dim];
3687 }
3688 else
3690 }
3691 else
3693 }
3694 else
3695 {
3696 value = is_linear ?
3698 dim - 1,
3701 stride_face_dof_actual>(face_dof_values_ptr,
3702 this->unit_point_faces_ptr[qb]) :
3704 dim - 1,
3707 false,
3708 stride_face_dof_actual>(this->shapes_faces.data() +
3709 qb * n_shapes,
3710 n_shapes,
3711 face_dof_values_ptr);
3712 }
3713
3714 if (evaluation_flags & EvaluationFlags::values)
3715 {
3716 for (unsigned int v = 0, offset = qb * stride;
3717 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3718 ++v, ++offset)
3719 ETT::set_value(value, v, this->values[offset]);
3720 }
3721 if (evaluation_flags & EvaluationFlags::gradients)
3722 {
3723 Assert(this->update_flags & update_gradients ||
3724 this->update_flags & update_inverse_jacobians,
3726
3727 for (unsigned int v = 0, offset = qb * stride;
3728 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3729 ++v, ++offset)
3730 {
3731 unit_gradient_type unit_gradient;
3732 ETT::set_gradient(gradient, v, unit_gradient);
3733 this->gradients[offset] =
3734 this->cell_type <=
3737 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3739 this
3740 ->inverse_jacobian_ptr[this->cell_type <=
3742 GeometryType::affine ?
3743 0 :
3744 offset]
3745 .transpose(),
3746 unit_gradient);
3747 }
3748 }
3749 }
3750}
3751
3752
3753
3754template <int n_components_, int dim, int spacedim, typename Number>
3755template <int stride_face_dof>
3756void
3758 ScalarNumber *face_dof_values,
3759 const EvaluationFlags::EvaluationFlags &integration_flags,
3760 const bool sum_into_values)
3761{
3762 if (this->use_linear_path)
3763 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3764 integration_flags,
3765 sum_into_values);
3766 else
3767 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3768 integration_flags,
3769 sum_into_values);
3770}
3771
3772
3773
3774template <int n_components_, int dim, int spacedim, typename Number>
3775template <bool do_JxW, bool is_linear, int stride_face_dof>
3776inline void
3779 ScalarNumber *face_dof_values,
3780 const EvaluationFlags::EvaluationFlags &integration_flags,
3781 const bool sum_into_values)
3782{
3783 // zero out lanes of incomplete last quadrature point batch
3784 if constexpr (stride == 1)
3785 if (const unsigned int n_filled_lanes =
3786 this->n_q_points_scalar & (n_lanes_internal - 1);
3787 n_filled_lanes > 0)
3788 {
3789 if (integration_flags & EvaluationFlags::values)
3790 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3791 ETT::set_zero_value(this->values.back(), v);
3792 if (integration_flags & EvaluationFlags::gradients)
3793 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3794 ETT::set_zero_gradient(this->gradients.back(), v);
3795 }
3796
3797 std::array<vectorized_value_type,
3798 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3799 solution_values_vectorized_linear = {};
3800
3801 // loop over quadrature batches qb
3802 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3803
3804 const bool cartesian_cell =
3806 const bool affine_cell =
3808 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3809 {
3810 vectorized_value_type value = {};
3812
3813 if (integration_flags & EvaluationFlags::values)
3814 for (unsigned int v = 0, offset = qb * stride;
3815 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3816 ++v, ++offset)
3817 ETT::get_value(value,
3818 v,
3819 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3820 this->values[offset]);
3821
3822 if (integration_flags & EvaluationFlags::gradients)
3823 for (unsigned int v = 0, offset = qb * stride;
3824 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3825 ++v, ++offset)
3826 {
3827 const auto grad_w =
3828 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3829 this->gradients[offset];
3830 ETT::get_gradient(
3831 gradient,
3832 v,
3833 cartesian_cell ?
3834 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3835 grad_w) :
3837 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3838 grad_w));
3839 }
3840
3841 if (integration_flags & EvaluationFlags::gradients)
3842 {
3843 std::array<vectorized_value_type, 2> value_face = {};
3844 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3845
3846 value_face[0] = value;
3847 // fill derivative in physical coordinates into tangential/normal
3848 // derivatives
3849 if (this->current_face_number / 2 == 0)
3850 {
3851 value_face[1] = gradient[0];
3852 if (dim > 1)
3853 gradient_in_face[0] = gradient[1];
3854 if (dim > 2)
3855 gradient_in_face[1] = gradient[2];
3856 }
3857 else if (this->current_face_number / 2 == 1)
3858 {
3859 if (dim > 1)
3860 value_face[1] = gradient[1];
3861 if (dim == 3)
3862 {
3863 gradient_in_face[0] = gradient[2];
3864 gradient_in_face[1] = gradient[0];
3865 }
3866 else if (dim == 2)
3867 gradient_in_face[0] = gradient[0];
3868 else
3870 }
3871 else if (this->current_face_number / 2 == 2)
3872 {
3873 if (dim > 2)
3874 {
3875 value_face[1] = gradient[2];
3876 gradient_in_face[0] = gradient[0];
3877 gradient_in_face[1] = gradient[1];
3878 }
3879 else
3881 }
3882 else
3884
3886 is_linear,
3887 dim - 1,
3890 2>(this->shapes_faces.data() + qb * n_shapes,
3891 n_shapes,
3892 value_face.data(),
3893 gradient_in_face,
3894 is_linear ? solution_values_vectorized_linear.data() :
3895 this->solution_renumbered_vectorized.data(),
3896 this->unit_point_faces_ptr[qb],
3897 qb != 0);
3898 }
3899 else
3901 dim - 1,
3904 this->shapes_faces.data() + qb * n_shapes,
3905 n_shapes,
3906 value,
3907 is_linear ? solution_values_vectorized_linear.data() :
3908 this->solution_renumbered_vectorized.data(),
3909 this->unit_point_faces_ptr[qb],
3910 qb != 0);
3911 }
3912
3913 const unsigned int dofs_per_comp_face =
3914 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3915
3916 for (unsigned int comp = 0; comp < n_components; ++comp)
3917 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3918 if (sum_into_values)
3919 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3920 stride_face_dof] +=
3921 ETT::sum_value(comp,
3922 is_linear ?
3923 *(solution_values_vectorized_linear.data() + i) :
3924 this->solution_renumbered_vectorized[i]);
3925 else
3926 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3927 ETT::sum_value(comp,
3928 is_linear ?
3929 *(solution_values_vectorized_linear.data() + i) :
3930 this->solution_renumbered_vectorized[i]);
3931}
3932
3933
3934
3935template <int n_components_, int dim, int spacedim, typename Number>
3938 const unsigned int point_index) const
3939{
3940 AssertIndexRange(point_index, this->n_q_points);
3941 Assert(this->normal_ptr != nullptr,
3943 ExcFEPointEvaluationAccessToUninitializedMappingField(
3944 "update_normal_vectors"));
3945 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3946 {
3948 for (unsigned int d = 0; d < dim; ++d)
3949 normal[d] =
3950 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3951 0);
3952
3953 return normal;
3954 }
3955 else
3956 {
3957 return this->normal_ptr[point_index];
3958 }
3959}
3960
3961
3962
3963template <int n_components_, int dim, int spacedim, typename Number>
3964inline const typename FEFacePointEvaluation<n_components_,
3965 dim,
3966 spacedim,
3967 Number>::value_type
3969 get_normal_derivative(const unsigned int point_index) const
3970{
3971 AssertIndexRange(point_index, this->gradients.size());
3972
3973 value_type normal_derivative;
3974 if constexpr (n_components == 1)
3975 normal_derivative =
3976 this->gradients[point_index] * normal_vector(point_index);
3977 else
3978 for (unsigned int comp = 0; comp < n_components; ++comp)
3979 normal_derivative[comp] =
3980 this->gradients[point_index][comp] * normal_vector(point_index);
3981
3982 return normal_derivative;
3983}
3984
3985
3986
3987template <int n_components_, int dim, int spacedim, typename Number>
3988inline void
3991 const unsigned int point_index)
3992{
3993 AssertIndexRange(point_index, this->gradients.size());
3994 if constexpr (n_components == 1)
3995 this->gradients[point_index] = value * normal_vector(point_index);
3996 else
3997 for (unsigned int comp = 0; comp < n_components; ++comp)
3998 this->gradients[point_index][comp] =
3999 value[comp] * normal_vector(point_index);
4000}
4001
4003
4004#endif
value_type * data() const noexcept
Definition array_view.h:666
std::size_t size() const
Definition array_view.h:689
typename ETT::vectorized_value_type vectorized_value_type
static constexpr std::size_t n_lanes_internal
void do_integrate_in_face(ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void integrate_in_face(ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename ETT::scalar_value_type scalar_value_type
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void do_evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void submit_normal_derivative(const value_type &, const unsigned int point_index)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
static constexpr std::size_t stride
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
void do_evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr unsigned int n_components
static constexpr std::size_t n_lanes_user_interface
typename ETT::real_gradient_type gradient_type
FEFacePointEvaluation(const NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const bool is_interior=true, const unsigned int first_selected_component=0)
void evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void reinit(const unsigned int cell_index, const unsigned int face_number)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename ETT::unit_gradient_type unit_gradient_type
const value_type get_normal_derivative(const unsigned int point_index) const
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
static constexpr unsigned int dimension
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
typename ETT::value_type value_type
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
ObserverPointer< const Mapping< dim, spacedim > > mapping
const DerivativeForm< 1, dim, spacedim, Number > * jacobian_ptr
std::unique_ptr< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info_on_the_fly
std::vector< gradient_type > gradients
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Number get_divergence(const unsigned int point_index) const
Number JxW(const unsigned int point_index) const
const UpdateFlags update_flags
static constexpr std::size_t n_lanes_user_interface
static constexpr std::size_t stride
internal::MatrixFreeFunctions::GeometryType cell_type
std::vector< Polynomials::Polynomial< double > > poly
Point< spacedim, Number > real_point(const unsigned int point_index) const
Tensor< 1,(dim==2 ? 1 :dim), Number > get_curl(const unsigned int point_index) const
const unsigned int n_q_points_scalar
const Point< dim, VectorizedArrayType > * unit_point_ptr
FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept
AlignedVector< ScalarNumber > scratch_data_scalar
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
const value_type & get_value(const unsigned int point_index) const
const Point< dim - 1, VectorizedArrayType > * unit_point_faces_ptr
std::vector< scalar_value_type > solution_renumbered
const unsigned int n_q_batches
ObserverPointer< const FiniteElement< dim, spacedim > > fe
std::vector< std::array< bool, n_components > > nonzero_shape_function_component
Point< spacedim, Number > quadrature_point(const unsigned int point_index) const
unsigned int n_active_entries_per_quadrature_batch(unsigned int q)
const gradient_type & get_gradient(const unsigned int point_index) const
std::shared_ptr< FEValues< dim, spacedim > > fe_values
FEPointEvaluationBase(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
std::vector< unsigned int > renumber
std::vector< value_type > values
AlignedVector<::ndarray< VectorizedArrayType, 2, dim - 1 > > shapes_faces
Point< dim, Number > unit_point(const unsigned int point_index) const
const unsigned int n_q_points
void submit_divergence(const Number &value, const unsigned int point_index)
const Tensor< 1, spacedim, Number > * normal_ptr
DerivativeForm< 1, spacedim, dim, Number > inverse_jacobian(const unsigned int point_index) const
AlignedVector< vectorized_value_type > solution_renumbered_vectorized
static constexpr std::size_t n_lanes_internal
ObserverPointer< const NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info
const DerivativeForm< 1, spacedim, dim, Number > * inverse_jacobian_ptr
AlignedVector<::ndarray< VectorizedArrayType, 2, dim > > shapes
typename ETT::value_type value_type
scalar_value_type integrate_value() const
void submit_gradient(const gradient_type &, const unsigned int point_index)
void setup(const unsigned int first_selected_component)
typename ETT::scalar_value_type scalar_value_type
static constexpr unsigned int dimension
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::vectorized_value_type vectorized_value_type
static constexpr unsigned int n_components
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept
DerivativeForm< 1, dim, spacedim, Number > jacobian(const unsigned int point_index) const
typename ETT::real_gradient_type gradient_type
const Point< spacedim, Number > * real_point_ptr
FEPointEvaluationBase(const NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const unsigned int first_selected_component=0, const bool is_interior=true)
internal::MatrixFreeFunctions::ShapeInfo< ScalarNumber > shape_info
void submit_value(const value_type &value, const unsigned int point_index)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static constexpr std::size_t stride
void integrate_slow(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr unsigned int dimension
typename ETT::value_type value_type
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename ETT::scalar_value_type scalar_value_type
typename ETT::unit_gradient_type unit_gradient_type
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
const value_type get_normal_derivative(const unsigned int point_index) const
void prepare_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
typename ETT::vectorized_value_type vectorized_value_type
void compute_integrate_fast(const EvaluationFlags::EvaluationFlags &integration_flags, const unsigned int n_shapes, const unsigned int qb, const vectorized_value_type value, const interface_vectorized_unit_gradient_type gradient, vectorized_value_type *solution_values_vectorized_linear)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::real_gradient_type gradient_type
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void submit_normal_derivative(const value_type &, const unsigned int point_index)
void evaluate_slow(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr std::size_t n_lanes_internal
static constexpr unsigned int n_components
void finish_integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, vectorized_value_type *solution_values_vectorized_linear, const bool sum_into_values)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
void internal_reinit_single_cell_state_mapping_info()
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
Abstract base class for mapping classes.
Definition mapping.h:320
Definition point.h:113
value_type * data() const noexcept
Definition array_view.h:851
std::size_t size() const
Definition array_view.h:872
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_DEPRECATED
Definition config.h:280
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_diagonal_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int cell_index
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcFEPointEvaluationAccessToUninitializedMappingField(std::string arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_jacobians
Volume element.
@ update_inverse_jacobians
Volume element.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
std::vector< Polynomials::Polynomial< double > > get_polynomial_space(const FiniteElement< dim, spacedim > &fe)
bool is_fast_path_supported(const FiniteElement< dim, spacedim > &fe, const unsigned int base_element_number)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void integrate_tensor_product_value_and_gradient(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 *value, const Tensor< 1, dim, Number2 > &gradient, Number2 *values, const Point< dim, Number > &p, const bool do_add)
void compute_values_of_array_in_pairs(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p0, const Point< dim, Number > &p1)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_shapes(const ::ndarray< Number2, 2, dim > *shapes, const int n_shapes, const Number *values, const std::vector< unsigned int > &renumber={})
void compute_values_of_array(::ndarray< Number, 2, dim > *shapes, const std::vector< Polynomials::Polynomial< double > > &poly, const Point< dim, Number > &p, const unsigned int derivative=1)
ProductTypeNoPoint< Number, Number2 >::type evaluate_tensor_product_value_linear(const Number *values, const Point< dim, Number2 > &p)
std::array< typename ProductTypeNoPoint< Number, Number2 >::type, dim+n_values > evaluate_tensor_product_value_and_gradient_linear(const Number *values, const Point< dim, Number2 > &p)
void integrate_tensor_product_value(const ::ndarray< Number, 2, dim > *shapes, const unsigned int n_shapes, const Number2 &value, Number2 *values, const Point< dim, Number > &p, const bool do_add)
constexpr unsigned int invalid_unsigned_int
Definition types.h:238
boost::integer_range< IncrementableType > iota_view
Definition iota_view.h:45
STL namespace.
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
Definition ndarray.h:107
static void set_zero_gradient(unit_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static Tensor< 1, dim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int, const vectorized_unit_gradient_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int vector_lane, scalar_unit_gradient_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static scalar_value_type sum_value(const vectorized_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int, vectorized_unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int vector_lane, const scalar_unit_gradient_type &result)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const DerivativeForm< 1, dim, n_components, Number > &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static scalar_value_type sum_value(const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
Tensor< 1, n_components, ScalarNumber > scalar_value_type
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, VectorizedArrayType > > vectorized_unit_gradient_type
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
std::conditional_t< n_components==spacedim, Tensor< 2, spacedim, Number >, Tensor< 1, n_components, Tensor< 1, spacedim, Number > > > real_gradient_type
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
Tensor< 1, n_components, Tensor< 1, dim, Number > > unit_gradient_type
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static constexpr std::size_t width()
static constexpr std::size_t stride()
static value_type & get(value_type &value, unsigned int c)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)