deal.II version GIT relicensing-3315-g061cc17459 2025-05-16 03: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\}}\)
Loading...
Searching...
No Matches
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
1186 const unsigned int first_selected_component = 0,
1187 const bool force_lexicographic_numbering = false);
1188
1216 const unsigned int first_selected_component = 0,
1217 const bool force_lexicographic_numbering = false);
1218
1230 void
1232 const ArrayView<const Point<dim>> &unit_points);
1233
1238 void
1239 reinit();
1240
1245 void
1246 reinit(const unsigned int cell_index);
1247
1248
1260 template <std::size_t stride_view>
1261 void
1262 evaluate(
1264 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1265
1277 void
1278 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1279 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1280
1303 template <std::size_t stride_view>
1304 void
1306 const EvaluationFlags::EvaluationFlags &integration_flags,
1307 const bool sum_into_values = false);
1308
1331 void
1332 integrate(const ArrayView<ScalarNumber> &solution_values,
1333 const EvaluationFlags::EvaluationFlags &integration_flags,
1334 const bool sum_into_values = false);
1335
1362 template <std::size_t stride_view>
1363 void
1365 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1366 const EvaluationFlags::EvaluationFlags &integration_flags,
1367 const bool sum_into_values = false);
1368
1395 void
1396 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1397 const EvaluationFlags::EvaluationFlags &integration_flags,
1398 const bool sum_into_values = false);
1399
1406 normal_vector(const unsigned int point_index) const;
1407
1413 const value_type
1414 get_normal_derivative(const unsigned int point_index) const;
1415
1421 void
1422 submit_normal_derivative(const value_type &, const unsigned int point_index);
1423
1424private:
1425 static constexpr std::size_t n_lanes_user_interface =
1427 static constexpr std::size_t n_lanes_internal =
1429 static constexpr std::size_t stride =
1431
1433
1438 template <bool is_linear, std::size_t stride_view>
1439 void
1442
1447 template <bool is_linear, std::size_t stride_view>
1448 void
1451 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1452 const unsigned int n_shapes,
1453 const unsigned int qb,
1454 vectorized_value_type &value,
1456
1460 template <bool is_linear, std::size_t stride_view>
1461 void
1464 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1465
1469 template <std::size_t stride_view>
1470 void
1473 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1474
1480 template <bool is_linear>
1481 void
1483 const EvaluationFlags::EvaluationFlags &integration_flags,
1484 const unsigned int n_shapes,
1485 const unsigned int qb,
1486 const vectorized_value_type value,
1488 vectorized_value_type *solution_values_vectorized_linear);
1489
1495 template <bool is_linear, std::size_t stride_view>
1496 void
1498 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1499 vectorized_value_type *solution_values_vectorized_linear,
1500 const bool sum_into_values);
1501
1505 template <bool do_JxW, bool is_linear, 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
1515 template <bool do_JxW, std::size_t stride_view>
1516 void
1518 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1519 const EvaluationFlags::EvaluationFlags &integration_flags,
1520 const bool sum_into_values);
1521
1525 template <bool do_JxW, std::size_t stride_view>
1526 void
1528 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1529 const EvaluationFlags::EvaluationFlags &integration_flags,
1530 const bool sum_into_values);
1531
1537 void
1539};
1540
1541
1542
1560template <int n_components_,
1561 int dim,
1562 int spacedim = dim,
1563 typename Number = double>
1565 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1566{
1567public:
1568 static constexpr unsigned int dimension = dim;
1569 static constexpr unsigned int n_components = n_components_;
1570
1571 using number_type = Number;
1572
1575 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1577 using ETT = typename internal::FEPointEvaluation::
1578 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1579 using value_type = typename ETT::value_type;
1580 using scalar_value_type = typename ETT::scalar_value_type;
1581 using vectorized_value_type = typename ETT::vectorized_value_type;
1582 using unit_gradient_type = typename ETT::unit_gradient_type;
1583 using gradient_type = typename ETT::real_gradient_type;
1585 typename ETT::interface_vectorized_unit_gradient_type;
1586
1593 const bool is_interior = true,
1594 const unsigned int first_selected_component = 0);
1595
1600 void
1601 reinit(const unsigned int cell_index, const unsigned int face_number);
1602
1607 void
1608 reinit(const unsigned int face_index);
1609
1621 template <std::size_t stride_view>
1622 void
1623 evaluate(
1625 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1626
1638 void
1639 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1640 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1641
1664 template <std::size_t stride_view>
1665 void
1667 const EvaluationFlags::EvaluationFlags &integration_flags,
1668 const bool sum_into_values = false);
1669
1692 void
1693 integrate(const ArrayView<ScalarNumber> &solution_values,
1694 const EvaluationFlags::EvaluationFlags &integration_flags,
1695 const bool sum_into_values = false);
1696
1719 template <std::size_t stride_view>
1720 void
1722 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1723 const EvaluationFlags::EvaluationFlags &integration_flags,
1724 const bool sum_into_values = false);
1725
1748 void
1749 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1750 const EvaluationFlags::EvaluationFlags &integration_flags,
1751 const bool sum_into_values = false);
1752
1759 template <int stride_face_dof = VectorizedArrayType::size()>
1760 void
1761 evaluate_in_face(const ScalarNumber *face_dof_values,
1762 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1763
1770 template <int stride_face_dof = VectorizedArrayType::size()>
1771 void
1772 integrate_in_face(ScalarNumber *face_dof_values,
1773 const EvaluationFlags::EvaluationFlags &integration_flags,
1774 const bool sum_into_values = false);
1775
1782 normal_vector(const unsigned int point_index) const;
1783
1789 const value_type
1790 get_normal_derivative(const unsigned int point_index) const;
1791
1797 void
1798 submit_normal_derivative(const value_type &, const unsigned int point_index);
1799
1800private:
1801 static constexpr std::size_t n_lanes_user_interface =
1803 static constexpr std::size_t n_lanes_internal =
1805 static constexpr std::size_t stride =
1807
1808 template <bool is_linear, std::size_t stride_view>
1809 void
1812 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1813
1814 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1815 void
1817 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1818 const EvaluationFlags::EvaluationFlags &integration_flags,
1819 const bool sum_into_values);
1820
1825 template <bool is_linear, int stride_face_dof>
1826 void
1827 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1828 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1829
1834 template <bool do_JxW, bool is_linear, int stride_face_dof>
1835 void
1837 ScalarNumber *face_dof_values,
1838 const EvaluationFlags::EvaluationFlags &integration_flags,
1839 const bool sum_into_values);
1840};
1841
1842
1843
1844// ----------------------- template and inline function ----------------------
1845
1846
1847template <int n_components_, int dim, int spacedim, typename Number>
1851 const UpdateFlags update_flags,
1852 const unsigned int first_selected_component)
1853 : n_q_batches(numbers::invalid_unsigned_int)
1854 , n_q_points(numbers::invalid_unsigned_int)
1855 , n_q_points_scalar(numbers::invalid_unsigned_int)
1856 , mapping(&mapping)
1857 , fe(&fe)
1858 , JxW_ptr(nullptr)
1859 , update_flags(update_flags)
1860 , mapping_info_on_the_fly(
1861 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1862 mapping,
1863 update_flags))
1864 , mapping_info(mapping_info_on_the_fly.get())
1865 , current_cell_index(numbers::invalid_unsigned_int)
1866 , current_face_number(numbers::invalid_unsigned_int)
1867 , must_reinitialize_pointers(false)
1868 , is_interior(true)
1869{
1870 setup(first_selected_component);
1871}
1872
1873
1874
1875template <int n_components_, int dim, int spacedim, typename Number>
1880 const unsigned int first_selected_component,
1881 const bool is_interior)
1882 : n_q_batches(numbers::invalid_unsigned_int)
1883 , n_q_points(numbers::invalid_unsigned_int)
1884 , n_q_points_scalar(numbers::invalid_unsigned_int)
1885 , mapping(&mapping_info.get_mapping())
1886 , fe(&fe)
1887 , JxW_ptr(nullptr)
1888 , update_flags(mapping_info.get_update_flags())
1889 , mapping_info(&mapping_info)
1890 , current_cell_index(numbers::invalid_unsigned_int)
1891 , current_face_number(numbers::invalid_unsigned_int)
1892 , must_reinitialize_pointers(true)
1893 , is_interior(is_interior)
1894{
1895 setup(first_selected_component);
1896}
1897
1898
1899
1900template <int n_components_, int dim, int spacedim, typename Number>
1904 : n_q_batches(other.n_q_batches)
1905 , n_q_points(other.n_q_points)
1906 , n_q_points_scalar(other.n_q_points_scalar)
1907 , mapping(other.mapping)
1908 , fe(other.fe)
1909 , poly(other.poly)
1910 , use_linear_path(other.use_linear_path)
1911 , renumber(other.renumber)
1912 , solution_renumbered(other.solution_renumbered)
1913 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1914 , values(other.values)
1915 , gradients(other.gradients)
1916 , dofs_per_component(other.dofs_per_component)
1917 , dofs_per_component_face(other.dofs_per_component_face)
1918 , component_in_base_element(other.component_in_base_element)
1919 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1920 , update_flags(other.update_flags)
1921 , fe_values(other.fe_values)
1922 , mapping_info_on_the_fly(
1923 other.mapping_info_on_the_fly ?
1925 *mapping,
1926 update_flags) :
1927 nullptr)
1928 , mapping_info(other.mapping_info)
1929 , current_cell_index(other.current_cell_index)
1930 , current_face_number(other.current_face_number)
1931 , fast_path(other.fast_path)
1932 , must_reinitialize_pointers(true)
1933 , is_interior(other.is_interior)
1934{}
1935
1936
1937
1938template <int n_components_, int dim, int spacedim, typename Number>
1940 FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept = default;
1941
1942
1943
1944template <int n_components_, int dim, int spacedim, typename Number>
1945void
1947 const unsigned int first_selected_component)
1948{
1949 AssertIndexRange(first_selected_component + n_components,
1950 fe->n_components() + 1);
1951
1952 shapes.reserve(100);
1953
1954 bool same_base_element = true;
1955 unsigned int base_element_number = 0;
1956 component_in_base_element = 0;
1957 unsigned int component = 0;
1958 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1959 if (component + fe->element_multiplicity(base_element_number) >
1960 first_selected_component)
1961 {
1962 if (first_selected_component + n_components >
1963 component + fe->element_multiplicity(base_element_number))
1964 same_base_element = false;
1965 component_in_base_element = first_selected_component - component;
1966 break;
1967 }
1968 else
1969 component += fe->element_multiplicity(base_element_number);
1970
1973 *fe, base_element_number) &&
1974 same_base_element)
1975 {
1976 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1977 renumber = shape_info.lexicographic_numbering;
1978 dofs_per_component = shape_info.dofs_per_component_on_cell;
1979 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1981 fe->base_element(base_element_number));
1982
1983 bool is_lexicographic = true;
1984 for (unsigned int i = 0; i < renumber.size(); ++i)
1985 if (i != renumber[i])
1986 is_lexicographic = false;
1987
1988 if (is_lexicographic)
1989 renumber.clear();
1990
1991 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1992 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1993 poly[1].value(1.) == 1.) &&
1994 (fe->n_components() == n_components);
1995
1996 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1997 const unsigned int size_cell = dofs_per_component * n_components;
1998 scratch_data_scalar.resize(size_face + size_cell);
1999
2000 solution_renumbered.resize(dofs_per_component);
2001 solution_renumbered_vectorized.resize(dofs_per_component);
2002
2003 fast_path = true;
2004 }
2005 else
2006 {
2007 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
2008 for (unsigned int d = 0; d < n_components; ++d)
2009 {
2010 const unsigned int component = first_selected_component + d;
2011 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2012 {
2013 const bool is_primitive =
2014 fe->is_primitive() || fe->is_primitive(i);
2015 if (is_primitive)
2016 nonzero_shape_function_component[i][d] =
2017 (component == fe->system_to_component_index(i).first);
2018 else
2019 nonzero_shape_function_component[i][d] =
2020 (fe->get_nonzero_components(i)[component] == true);
2021 }
2022 }
2023
2024 fast_path = false;
2025 }
2026}
2027
2028
2029
2030template <int n_components_, int dim, int spacedim, typename Number>
2031template <bool is_face, bool is_linear>
2032inline void
2034{
2035 const unsigned int geometry_index =
2036 mapping_info->template compute_geometry_index_offset<is_face>(
2037 current_cell_index, current_face_number);
2038
2039 cell_type = mapping_info->get_cell_type(geometry_index);
2040
2041 const_cast<unsigned int &>(n_q_points_scalar) =
2042 mapping_info->get_n_q_points_unvectorized(geometry_index);
2043
2044 // round up n_q_points_scalar / n_lanes_internal
2045 const_cast<unsigned int &>(n_q_batches) =
2046 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
2047
2048 const unsigned int n_q_points_before = n_q_points;
2049
2050 const_cast<unsigned int &>(n_q_points) =
2051 (stride == 1) ? n_q_batches : n_q_points_scalar;
2052
2053 if (n_q_points != n_q_points_before)
2054 {
2055 if (update_flags & update_values)
2056 values.resize(n_q_points);
2057 if (update_flags & update_gradients)
2058 gradients.resize(n_q_points);
2059 }
2060
2061 if (n_q_points == 0)
2062 return;
2063
2064 // set unit point pointer
2065 const unsigned int unit_point_offset =
2066 mapping_info->compute_unit_point_index_offset(geometry_index);
2067
2068 if (is_face)
2069 unit_point_faces_ptr =
2070 mapping_info->get_unit_point_faces(unit_point_offset);
2071 else
2072 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2073
2074 // set data pointers
2075 const unsigned int data_offset =
2076 mapping_info->compute_data_index_offset(geometry_index);
2077 const unsigned int compressed_data_offset =
2078 mapping_info->compute_compressed_data_index_offset(geometry_index);
2079 if constexpr (running_in_debug_mode())
2080 {
2081 const UpdateFlags update_flags_mapping =
2082 mapping_info->get_update_flags_mapping();
2083 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2084 real_point_ptr = mapping_info->get_real_point(data_offset);
2085 if (update_flags_mapping & UpdateFlags::update_jacobians)
2086 jacobian_ptr =
2087 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2088 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2089 inverse_jacobian_ptr =
2090 mapping_info->get_inverse_jacobian(compressed_data_offset,
2091 is_interior);
2092 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2093 normal_ptr = mapping_info->get_normal_vector(data_offset);
2094 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2095 JxW_ptr = mapping_info->get_JxW(data_offset);
2096 }
2097 else
2098 {
2099 real_point_ptr = mapping_info->get_real_point(data_offset);
2100 jacobian_ptr =
2101 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2102 inverse_jacobian_ptr =
2103 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2104 normal_ptr = mapping_info->get_normal_vector(data_offset);
2105 JxW_ptr = mapping_info->get_JxW(data_offset);
2106 }
2107
2108 if (!is_linear && fast_path)
2109 {
2110 const std::size_t n_shapes = poly.size();
2111 if (is_face)
2112 shapes_faces.resize_fast(n_q_batches * n_shapes);
2113 else
2114 shapes.resize_fast(n_q_batches * n_shapes);
2115
2116 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2117 if (is_face)
2118 {
2119 if (dim > 1)
2120 {
2122 shapes_faces.data() + qb * n_shapes,
2123 poly,
2124 unit_point_faces_ptr[qb],
2125 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2126 }
2127 }
2128 else
2129 {
2130 if (update_flags & UpdateFlags::update_gradients)
2131 {
2132 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2133 poly,
2134 unit_point_ptr[qb],
2135 1);
2136 }
2137 else if (qb + 1 < n_q_batches)
2138 {
2139 // Use function with reduced overhead to compute for two
2140 // points at once
2142 shapes.data() + qb * n_shapes,
2143 poly,
2144 unit_point_ptr[qb],
2145 unit_point_ptr[qb + 1]);
2146 ++qb;
2147 }
2148 else
2149 {
2150 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2151 poly,
2152 unit_point_ptr[qb],
2153 0);
2154 }
2155 }
2156 }
2157}
2158
2159
2160
2161template <int n_components_, int dim, int spacedim, typename Number>
2162inline const typename FEPointEvaluationBase<n_components_,
2163 dim,
2164 spacedim,
2165 Number>::value_type &
2167 const unsigned int point_index) const
2168{
2169 AssertIndexRange(point_index, values.size());
2170 return values[point_index];
2171}
2172
2173
2174
2175template <int n_components_, int dim, int spacedim, typename Number>
2176inline const typename FEPointEvaluationBase<n_components_,
2177 dim,
2178 spacedim,
2179 Number>::gradient_type &
2181 const unsigned int point_index) const
2182{
2183 AssertIndexRange(point_index, gradients.size());
2184 return gradients[point_index];
2185}
2186
2187
2188
2189template <int n_components_, int dim, int spacedim, typename Number>
2190inline Number
2192 const unsigned int point_index) const
2193{
2194 static_assert(n_components == dim,
2195 "Only makes sense for a vector field with dim components");
2196
2197 AssertIndexRange(point_index, values.size());
2198 return trace(gradients[point_index]);
2199}
2200
2201
2202
2203template <int n_components_, int dim, int spacedim, typename Number>
2204inline void
2206 const value_type &value,
2207 const unsigned int point_index)
2208{
2209 AssertIndexRange(point_index, n_q_points);
2210 values[point_index] = value;
2211}
2212
2213
2214
2215template <int n_components_, int dim, int spacedim, typename Number>
2216inline void
2218 const gradient_type &gradient,
2219 const unsigned int point_index)
2220{
2221 AssertIndexRange(point_index, n_q_points);
2222 gradients[point_index] = gradient;
2223}
2224
2225
2226
2227template <int n_components_, int dim, int spacedim, typename Number>
2228inline void
2230 const Number &value,
2231 const unsigned int point_index)
2232{
2233 static_assert(n_components == dim,
2234 "Only makes sense for a vector field with dim components");
2235
2236 AssertIndexRange(point_index, n_q_points);
2237 gradients[point_index] = gradient_type();
2238 for (unsigned int d = 0; d < dim; ++d)
2239 gradients[point_index][d][d] = value;
2240}
2241
2242
2243
2244template <int n_components_, int dim, int spacedim, typename Number>
2245Tensor<1, (dim == 2 ? 1 : dim), Number>
2247 const unsigned int point_index) const
2248{
2249 static_assert(
2250 dim > 1 && n_components == dim,
2251 "Only makes sense for a vector field with dim components and dim > 1");
2252
2253 const Tensor<2, dim, Number> grad = get_gradient(point_index);
2254 Tensor<1, (dim == 2 ? 1 : dim), Number> curl;
2255 switch (dim)
2256 {
2257 case 2:
2258 curl[0] = grad[1][0] - grad[0][1];
2259 break;
2260 case 3:
2261 curl[0] = grad[2][1] - grad[1][2];
2262 curl[1] = grad[0][2] - grad[2][0];
2263 curl[2] = grad[1][0] - grad[0][1];
2264 break;
2265 default:
2267 }
2268 return curl;
2269}
2270
2271
2272
2273template <int n_components_, int dim, int spacedim, typename Number>
2276 const unsigned int point_index) const
2277{
2278 AssertIndexRange(point_index, n_q_points);
2279 Assert(jacobian_ptr != nullptr,
2281 ExcFEPointEvaluationAccessToUninitializedMappingField(
2282 "update_jacobians"));
2283 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2284 GeometryType::affine ?
2285 0 :
2286 point_index];
2287}
2288
2289
2290
2291template <int n_components_, int dim, int spacedim, typename Number>
2294 const unsigned int point_index) const
2295{
2296 AssertIndexRange(point_index, n_q_points);
2297 Assert(inverse_jacobian_ptr != nullptr,
2299 ExcFEPointEvaluationAccessToUninitializedMappingField(
2300 "update_inverse_jacobians"));
2301 return inverse_jacobian_ptr
2302 [cell_type <=
2304 0 :
2305 point_index];
2306}
2307
2308
2309
2310template <int n_components_, int dim, int spacedim, typename Number>
2311inline Number
2313 const unsigned int point_index) const
2314{
2315 AssertIndexRange(point_index, n_q_points);
2316 Assert(JxW_ptr != nullptr,
2318 ExcFEPointEvaluationAccessToUninitializedMappingField(
2319 "update_JxW_values"));
2320 return JxW_ptr[point_index];
2321}
2322
2323
2324
2325template <int n_components_, int dim, int spacedim, typename Number>
2328 const unsigned int point_index) const
2329{
2330 return quadrature_point(point_index);
2331}
2332
2333
2334
2335template <int n_components_, int dim, int spacedim, typename Number>
2338 const unsigned int point_index) const
2339{
2340 AssertIndexRange(point_index, n_q_points);
2341 Assert(real_point_ptr != nullptr,
2343 ExcFEPointEvaluationAccessToUninitializedMappingField(
2344 "update_quadrature_points"));
2345 return real_point_ptr[point_index];
2346}
2347
2348
2349
2350template <int n_components_, int dim, int spacedim, typename Number>
2351inline Point<dim, Number>
2353 const unsigned int point_index) const
2354{
2355 AssertIndexRange(point_index, n_q_points);
2356 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2357 Point<dim, Number> unit_point;
2358 for (unsigned int d = 0; d < dim; ++d)
2360 unit_point_ptr[point_index / stride][d], point_index % stride);
2361 return unit_point;
2362}
2363
2364
2365
2366template <int n_components_, int dim, int spacedim, typename Number>
2374
2375
2376
2377template <int n_components_, int dim, int spacedim, typename Number>
2381 const unsigned int first_selected_component,
2382 const bool force_lexicographic_numbering)
2383 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2384 mapping_info,
2385 fe,
2386 first_selected_component)
2387 , lexicographic_numbering(force_lexicographic_numbering ||
2388 this->renumber.empty())
2389{}
2390
2391
2392
2393template <int n_components_, int dim, int spacedim, typename Number>
2395 const Mapping<dim, spacedim> &mapping,
2397 const UpdateFlags update_flags,
2398 const unsigned int first_selected_component,
2399 const bool force_lexicographic_numbering)
2400 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2401 mapping,
2402 fe,
2403 update_flags,
2404 first_selected_component)
2405 , lexicographic_numbering(force_lexicographic_numbering ||
2406 this->renumber.empty())
2407{}
2408
2409
2410
2411template <int n_components_, int dim, int spacedim, typename Number>
2412inline void
2415{
2416 this->current_cell_index = numbers::invalid_unsigned_int;
2417 this->current_face_number = numbers::invalid_unsigned_int;
2418
2419 if (this->use_linear_path)
2420 this->template do_reinit<false, true>();
2421 else
2422 this->template do_reinit<false, false>();
2423}
2424
2425
2426
2427template <int n_components_, int dim, int spacedim, typename Number>
2428inline void
2430{
2431 internal_reinit_single_cell_state_mapping_info();
2432 this->must_reinitialize_pointers = false;
2433}
2434
2435
2436
2437template <int n_components_, int dim, int spacedim, typename Number>
2438inline void
2441 const ArrayView<const Point<dim>> &unit_points)
2442{
2443 // reinit is only allowed for mapping computation on the fly
2444 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2446
2447 this->mapping_info_on_the_fly->reinit(cell, unit_points);
2448 this->must_reinitialize_pointers = false;
2449
2450 if (!this->fast_path)
2451 {
2452 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2453 *this->mapping,
2454 *this->fe,
2456 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2457 this->update_flags);
2458 this->fe_values->reinit(cell);
2459 }
2460
2461 if (this->use_linear_path)
2462 this->template do_reinit<false, true>();
2463 else
2464 this->template do_reinit<false, false>();
2465}
2466
2467
2468
2469template <int n_components_, int dim, int spacedim, typename Number>
2470inline void
2472 const unsigned int cell_index)
2473{
2474 this->current_cell_index = cell_index;
2475 this->current_face_number = numbers::invalid_unsigned_int;
2476 this->must_reinitialize_pointers = false;
2477
2478 if (this->use_linear_path)
2479 this->template do_reinit<false, true>();
2480 else
2481 this->template do_reinit<false, false>();
2482
2483 if (!this->fast_path)
2484 {
2485 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2486
2487 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2488 for (unsigned int d = 0; d < dim; ++d)
2489 unit_points[v][d] =
2490 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2491
2492 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2493 *this->mapping,
2494 *this->fe,
2496 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2497 this->update_flags);
2498
2499 this->fe_values->reinit(
2500 this->mapping_info->get_cell_iterator(this->current_cell_index));
2501 }
2502}
2503
2504
2505
2506template <int n_components_, int dim, int spacedim, typename Number>
2507template <std::size_t stride_view>
2508void
2511 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2512{
2513 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2514
2515 if (!((evaluation_flags & EvaluationFlags::values) ||
2516 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2517 return;
2518
2519 if (this->must_reinitialize_pointers)
2520 internal_reinit_single_cell_state_mapping_info();
2521
2522 if (this->n_q_points == 0)
2523 return;
2524
2525 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2526 if (this->fast_path)
2527 {
2528 if (this->use_linear_path)
2529 evaluate_fast<true>(solution_values, evaluation_flags);
2530 else
2531 evaluate_fast<false>(solution_values, evaluation_flags);
2532 }
2533 else
2534 evaluate_slow(solution_values, evaluation_flags);
2535}
2536
2537
2538
2539template <int n_components_, int dim, int spacedim, typename Number>
2540void
2542 const ArrayView<const ScalarNumber> &solution_values,
2543 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2544{
2545 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2546 solution_values.size()),
2547 evaluation_flags);
2548}
2549
2550
2551
2552template <int n_components_, int dim, int spacedim, typename Number>
2553template <std::size_t stride_view>
2554void
2556 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2557 const EvaluationFlags::EvaluationFlags &integration_flags,
2558 const bool sum_into_values)
2559{
2560 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2561}
2562
2563
2564
2565template <int n_components_, int dim, int spacedim, typename Number>
2566void
2568 const ArrayView<ScalarNumber> &solution_values,
2569 const EvaluationFlags::EvaluationFlags &integration_flags,
2570 const bool sum_into_values)
2571{
2572 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2573 solution_values.size()),
2574 integration_flags,
2575 sum_into_values);
2576}
2577
2578
2579
2580template <int n_components_, int dim, int spacedim, typename Number>
2582 scalar_value_type
2584 const
2585{
2586 value_type return_value = {};
2587
2588 for (const auto point_index : this->quadrature_point_indices())
2589 return_value += values[point_index] * this->JxW(point_index);
2590
2591 return ETT::sum_value(return_value);
2592}
2593
2594
2595
2596template <int n_components_, int dim, int spacedim, typename Number>
2597unsigned int
2600{
2601 Assert(stride == 1,
2602 ExcMessage(
2603 "Calling this function only makes sense in fully vectorized mode."));
2604 if (q == n_q_batches - 1)
2605 {
2606 const unsigned int n_filled_lanes =
2607 n_q_points_scalar & (n_lanes_user_interface - 1);
2608 if (n_filled_lanes == 0)
2609 return n_lanes_user_interface;
2610 else
2611 return n_filled_lanes;
2612 }
2613 else
2614 return n_lanes_user_interface;
2615}
2616
2617
2618
2619template <int n_components_, int dim, int spacedim, typename Number>
2620template <std::size_t stride_view>
2621void
2623 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2624 const EvaluationFlags::EvaluationFlags &integration_flags,
2625 const bool sum_into_values)
2626{
2627 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2628}
2629
2630
2631
2632template <int n_components_, int dim, int spacedim, typename Number>
2633void
2635 const ArrayView<ScalarNumber> &solution_values,
2636 const EvaluationFlags::EvaluationFlags &integration_flags,
2637 const bool sum_into_values)
2638{
2639 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2640 solution_values.size()),
2641 integration_flags,
2642 sum_into_values);
2643}
2644
2645
2646
2647template <int n_components_, int dim, int spacedim, typename Number>
2648template <bool is_linear, std::size_t stride_view>
2649inline void
2652{
2653 const unsigned int dofs_per_comp =
2654 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2655
2656 for (unsigned int comp = 0; comp < n_components; ++comp)
2657 {
2658 const std::size_t offset =
2659 (this->component_in_base_element + comp) * dofs_per_comp;
2660
2661 if ((is_linear && n_components == 1) || lexicographic_numbering)
2662 {
2663 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2664 ETT::read_value(solution_values[i + offset],
2665 comp,
2666 this->solution_renumbered[i]);
2667 }
2668 else
2669 {
2670 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2671 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2672 ETT::read_value(solution_values[renumber_ptr[i]],
2673 comp,
2674 this->solution_renumbered[i]);
2675 }
2676 }
2677}
2678
2679
2680
2681template <int n_components_, int dim, int spacedim, typename Number>
2682template <bool is_linear, std::size_t stride_view>
2683inline void
2686 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2687 const unsigned int n_shapes,
2688 const unsigned int qb,
2689 vectorized_value_type &value,
2691{
2692 if (evaluation_flags & EvaluationFlags::gradients)
2693 {
2694 std::array<vectorized_value_type, dim + 1> result;
2695 if constexpr (is_linear)
2696 {
2697 if constexpr (n_components == 1)
2698 result =
2700 dim,
2703 1,
2704 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2705 else
2706 result =
2708 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2709 }
2710 else
2712 dim,
2715 1,
2716 false>(this->shapes.data() + qb * n_shapes,
2717 n_shapes,
2718 this->solution_renumbered.data());
2719 gradient[0] = result[0];
2720 if (dim > 1)
2721 gradient[1] = result[1];
2722 if (dim > 2)
2723 gradient[2] = result[2];
2724 value = result[dim];
2725 }
2726 else
2727 {
2728 if constexpr (is_linear)
2729 {
2730 if constexpr (n_components == 1)
2732 dim,
2735 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2736 else
2738 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2739 }
2740 else
2741 value =
2745 false>(
2746 this->shapes.data() + qb * n_shapes,
2747 n_shapes,
2748 this->solution_renumbered.data());
2749 }
2750}
2751
2752
2753
2754template <int n_components_, int dim, int spacedim, typename Number>
2755template <bool is_linear, std::size_t stride_view>
2756inline void
2759 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2760{
2761 if (!(is_linear && n_components == 1))
2762 prepare_evaluate_fast<is_linear>(solution_values);
2763
2764 // loop over quadrature batches qb
2765 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2766
2767 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2768 {
2771
2772 compute_evaluate_fast<is_linear>(
2773 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2774
2775 if (evaluation_flags & EvaluationFlags::values)
2776 {
2777 for (unsigned int v = 0, offset = qb * stride;
2778 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2779 ++v, ++offset)
2780 ETT::set_value(value, v, this->values[offset]);
2781 }
2782 if (evaluation_flags & EvaluationFlags::gradients)
2783 {
2784 Assert(this->update_flags & update_gradients ||
2785 this->update_flags & update_inverse_jacobians,
2787
2788 for (unsigned int v = 0, offset = qb * stride;
2789 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2790 ++v, ++offset)
2791 {
2792 unit_gradient_type unit_gradient;
2793 ETT::set_gradient(gradient, v, unit_gradient);
2794 this->gradients[offset] =
2795 this->cell_type <=
2798 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2800 this
2801 ->inverse_jacobian_ptr[this->cell_type <=
2803 GeometryType::affine ?
2804 0 :
2805 offset]
2806 .transpose(),
2807 unit_gradient);
2808 }
2809 }
2810 }
2811}
2812
2813
2814
2815template <int n_components_, int dim, int spacedim, typename Number>
2816template <std::size_t stride_view>
2817inline void
2820 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2821{
2822 // slow path with FEValues
2823 Assert(this->fe_values.get() != nullptr,
2824 ExcMessage(
2825 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2826
2827 const std::size_t n_points = this->fe_values->get_quadrature().size();
2828
2829 if (evaluation_flags & EvaluationFlags::values)
2830 {
2831 this->values.resize(this->n_q_points);
2832 std::fill(this->values.begin(), this->values.end(), value_type());
2833 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2834 {
2835 const ScalarNumber value = solution_values[i];
2836 for (unsigned int d = 0; d < n_components; ++d)
2837 if (this->nonzero_shape_function_component[i][d] &&
2838 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2839 for (unsigned int qb = 0, q = 0; q < n_points;
2840 ++qb, q += n_lanes_user_interface)
2841 for (unsigned int v = 0;
2842 v < n_lanes_user_interface && q + v < n_points;
2843 ++v)
2844 ETT::access(this->values[qb],
2845 v,
2846 d,
2847 this->fe_values->shape_value(i, q + v) * value);
2848 else if (this->nonzero_shape_function_component[i][d])
2849 for (unsigned int qb = 0, q = 0; q < n_points;
2850 ++qb, q += n_lanes_user_interface)
2851 for (unsigned int v = 0;
2852 v < n_lanes_user_interface && q + v < n_points;
2853 ++v)
2854 ETT::access(this->values[qb],
2855 v,
2856 d,
2857 this->fe_values->shape_value_component(i,
2858 q + v,
2859 d) *
2860 value);
2861 }
2862 }
2863
2864 if (evaluation_flags & EvaluationFlags::gradients)
2865 {
2866 this->gradients.resize(this->n_q_points);
2867 std::fill(this->gradients.begin(),
2868 this->gradients.end(),
2869 gradient_type());
2870 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2871 {
2872 const ScalarNumber value = solution_values[i];
2873 for (unsigned int d = 0; d < n_components; ++d)
2874 if (this->nonzero_shape_function_component[i][d] &&
2875 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2876 for (unsigned int qb = 0, q = 0; q < n_points;
2877 ++qb, q += n_lanes_user_interface)
2878 for (unsigned int v = 0;
2879 v < n_lanes_user_interface && q + v < n_points;
2880 ++v)
2881 ETT::access(this->gradients[qb],
2882 v,
2883 d,
2884 this->fe_values->shape_grad(i, q + v) * value);
2885 else if (this->nonzero_shape_function_component[i][d])
2886 for (unsigned int qb = 0, q = 0; q < n_points;
2887 ++qb, q += n_lanes_user_interface)
2888 for (unsigned int v = 0;
2889 v < n_lanes_user_interface && q + v < n_points;
2890 ++v)
2891 ETT::access(
2892 this->gradients[qb],
2893 v,
2894 d,
2895 this->fe_values->shape_grad_component(i, q + v, d) * value);
2896 }
2897 }
2898}
2899
2900
2901
2902template <int n_components_, int dim, int spacedim, typename Number>
2903template <bool is_linear>
2904inline void
2906 const EvaluationFlags::EvaluationFlags &integration_flags,
2907 const unsigned int n_shapes,
2908 const unsigned int qb,
2909 const vectorized_value_type value,
2911 vectorized_value_type *solution_values_vectorized_linear)
2912{
2913 if (integration_flags & EvaluationFlags::gradients)
2915 is_linear,
2916 dim,
2918 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2919 n_shapes,
2920 &value,
2921 gradient,
2922 is_linear ?
2923 solution_values_vectorized_linear :
2924 this->solution_renumbered_vectorized.data(),
2925 this->unit_point_ptr[qb],
2926 qb != 0);
2927 else
2929 dim,
2932 this->shapes.data() + qb * n_shapes,
2933 n_shapes,
2934 value,
2935 is_linear ? solution_values_vectorized_linear :
2936 this->solution_renumbered_vectorized.data(),
2937 this->unit_point_ptr[qb],
2938 qb != 0);
2939}
2940
2941
2942
2943template <int n_components_, int dim, int spacedim, typename Number>
2944template <bool is_linear, std::size_t stride_view>
2945inline void
2947 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2948 vectorized_value_type *solution_values_vectorized_linear,
2949 const bool sum_into_values)
2950{
2951 if (!sum_into_values && this->fe->n_components() > n_components)
2952 for (unsigned int i = 0; i < solution_values.size(); ++i)
2953 solution_values[i] = 0;
2954
2955 const unsigned int dofs_per_comp =
2956 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2957
2958 for (unsigned int comp = 0; comp < n_components; ++comp)
2959 {
2960 const std::size_t offset =
2961 (this->component_in_base_element + comp) * dofs_per_comp;
2962
2963 if (is_linear || lexicographic_numbering)
2964 {
2965 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2966 if (sum_into_values)
2967 solution_values[i + offset] +=
2968 ETT::sum_value(comp,
2969 is_linear ?
2970 *(solution_values_vectorized_linear + i) :
2971 this->solution_renumbered_vectorized[i]);
2972 else
2973 solution_values[i + offset] =
2974 ETT::sum_value(comp,
2975 is_linear ?
2976 *(solution_values_vectorized_linear + i) :
2977 this->solution_renumbered_vectorized[i]);
2978 }
2979 else
2980 {
2981 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2982 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2983 if (sum_into_values)
2984 solution_values[renumber_ptr[i]] +=
2985 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2986 else
2987 solution_values[renumber_ptr[i]] =
2988 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2989 }
2990 }
2991}
2992
2993
2994
2995template <int n_components_, int dim, int spacedim, typename Number>
2996template <bool do_JxW, bool is_linear, std::size_t stride_view>
2997inline void
2999 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3000 const EvaluationFlags::EvaluationFlags &integration_flags,
3001 const bool sum_into_values)
3002{
3003 // zero out lanes of incomplete last quadrature point batch
3004 if constexpr (stride == 1)
3005 if (const unsigned int n_filled_lanes =
3006 this->n_q_points_scalar & (n_lanes_internal - 1);
3007 n_filled_lanes > 0)
3008 {
3009 if (integration_flags & EvaluationFlags::values)
3010 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3011 ETT::set_zero_value(this->values.back(), v);
3012 if (integration_flags & EvaluationFlags::gradients)
3013 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3014 ETT::set_zero_gradient(this->gradients.back(), v);
3015 }
3016
3017 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
3018 solution_values_vectorized_linear = {};
3019
3020 // loop over quadrature batches qb
3021 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3022
3023 const bool cartesian_cell =
3025 const bool affine_cell =
3027 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3028 {
3029 vectorized_value_type value = {};
3031
3032 if (integration_flags & EvaluationFlags::values)
3033 for (unsigned int v = 0, offset = qb * stride;
3034 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3035 ++v, ++offset)
3036 ETT::get_value(value,
3037 v,
3038 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3039 this->values[offset]);
3040
3041 if (integration_flags & EvaluationFlags::gradients)
3042 for (unsigned int v = 0, offset = qb * stride;
3043 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3044 ++v, ++offset)
3045 {
3046 const gradient_type grad_w =
3047 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3048 this->gradients[offset];
3049 ETT::get_gradient(
3050 gradient,
3051 v,
3052 cartesian_cell ?
3053 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3054 grad_w) :
3056 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3057 grad_w));
3058 }
3059
3060 compute_integrate_fast<is_linear>(
3061 integration_flags,
3062 n_shapes,
3063 qb,
3064 value,
3065 gradient,
3066 solution_values_vectorized_linear.data());
3067 }
3068
3069 // add between the lanes and write into the result
3070 finish_integrate_fast<is_linear>(solution_values,
3071 solution_values_vectorized_linear.data(),
3072 sum_into_values);
3073}
3074
3075
3076
3077template <int n_components_, int dim, int spacedim, typename Number>
3078template <bool do_JxW, std::size_t stride_view>
3079inline void
3081 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3082 const EvaluationFlags::EvaluationFlags &integration_flags,
3083 const bool sum_into_values)
3084{
3085 // slow path with FEValues
3086 Assert(this->fe_values.get() != nullptr,
3087 ExcMessage(
3088 "Not initialized. Please call FEPointEvaluation::reinit()!"));
3089 if (!sum_into_values)
3090 for (unsigned int i = 0; i < solution_values.size(); ++i)
3091 solution_values[i] = 0;
3092
3093 const std::size_t n_points = this->fe_values->get_quadrature().size();
3094
3095 if (integration_flags & EvaluationFlags::values)
3096 {
3097 AssertIndexRange(this->n_q_points, this->values.size() + 1);
3098 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3099 {
3100 for (unsigned int d = 0; d < n_components; ++d)
3101 if (this->nonzero_shape_function_component[i][d] &&
3102 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3103 for (unsigned int qb = 0, q = 0; q < n_points;
3104 ++qb, q += n_lanes_user_interface)
3105 for (unsigned int v = 0;
3106 v < n_lanes_user_interface && q + v < n_points;
3107 ++v)
3108 solution_values[i] +=
3109 this->fe_values->shape_value(i, q + v) *
3110 ETT::access(this->values[qb], v, d) *
3111 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3112 else if (this->nonzero_shape_function_component[i][d])
3113 for (unsigned int qb = 0, q = 0; q < n_points;
3114 ++qb, q += n_lanes_user_interface)
3115 for (unsigned int v = 0;
3116 v < n_lanes_user_interface && q + v < n_points;
3117 ++v)
3118 solution_values[i] +=
3119 this->fe_values->shape_value_component(i, q + v, d) *
3120 ETT::access(this->values[qb], v, d) *
3121 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3122 }
3123 }
3124
3125 if (integration_flags & EvaluationFlags::gradients)
3126 {
3127 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
3128 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3129 {
3130 for (unsigned int d = 0; d < n_components; ++d)
3131 if (this->nonzero_shape_function_component[i][d] &&
3132 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3133 for (unsigned int qb = 0, q = 0; q < n_points;
3134 ++qb, q += n_lanes_user_interface)
3135 for (unsigned int v = 0;
3136 v < n_lanes_user_interface && q + v < n_points;
3137 ++v)
3138 solution_values[i] +=
3139 this->fe_values->shape_grad(i, q + v) *
3140 ETT::access(this->gradients[qb], v, d) *
3141 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3142 else if (this->nonzero_shape_function_component[i][d])
3143 for (unsigned int qb = 0, q = 0; q < n_points;
3144 ++qb, q += n_lanes_user_interface)
3145 for (unsigned int v = 0;
3146 v < n_lanes_user_interface && q + v < n_points;
3147 ++v)
3148 solution_values[i] +=
3149 this->fe_values->shape_grad_component(i, q + v, d) *
3150 ETT::access(this->gradients[qb], v, d) *
3151 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3152 }
3153 }
3154}
3155
3156
3157
3158template <int n_components_, int dim, int spacedim, typename Number>
3159template <bool do_JxW, std::size_t stride_view>
3160void
3162 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3163 const EvaluationFlags::EvaluationFlags &integration_flags,
3164 const bool sum_into_values)
3165{
3166 if (this->must_reinitialize_pointers)
3167 internal_reinit_single_cell_state_mapping_info();
3168
3169 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3170
3171 if (this->n_q_points == 0 || // no evaluation points provided
3172 !((integration_flags & EvaluationFlags::values) ||
3173 (integration_flags &
3174 EvaluationFlags::gradients))) // no integration flags
3175 {
3176 if (!sum_into_values)
3177 for (unsigned int i = 0; i < solution_values.size(); ++i)
3178 solution_values[i] = 0;
3179 return;
3180 }
3181
3182 Assert(
3183 !do_JxW || this->JxW_ptr != nullptr,
3184 ExcMessage(
3185 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3186
3187 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3188 if (this->fast_path)
3189 {
3190 if (this->use_linear_path)
3191 integrate_fast<do_JxW, true>(solution_values,
3192 integration_flags,
3193 sum_into_values);
3194 else
3195 integrate_fast<do_JxW, false>(solution_values,
3196 integration_flags,
3197 sum_into_values);
3198 }
3199 else
3200 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3201}
3202
3203
3204
3205template <int n_components_, int dim, int spacedim, typename Number>
3208 const unsigned int point_index) const
3209{
3210 AssertIndexRange(point_index, this->n_q_points);
3211 Assert(this->normal_ptr != nullptr,
3213 ExcFEPointEvaluationAccessToUninitializedMappingField(
3214 "update_normal_vectors"));
3215 return this->normal_ptr[point_index];
3216}
3217
3218
3219
3220template <int n_components_, int dim, int spacedim, typename Number>
3221inline const typename FEPointEvaluation<n_components_,
3222 dim,
3223 spacedim,
3224 Number>::value_type
3226 const unsigned int point_index) const
3227{
3228 AssertIndexRange(point_index, this->gradients.size());
3229
3230 value_type normal_derivative;
3231 if constexpr (n_components == 1)
3232 normal_derivative =
3233 this->gradients[point_index] * normal_vector(point_index);
3234 else
3235 for (unsigned int comp = 0; comp < n_components; ++comp)
3236 normal_derivative[comp] =
3237 this->gradients[point_index][comp] * normal_vector(point_index);
3238
3239 return normal_derivative;
3240}
3241
3242
3243
3244template <int n_components_, int dim, int spacedim, typename Number>
3245inline void
3248 const unsigned int point_index)
3249{
3250 AssertIndexRange(point_index, this->gradients.size());
3251 if constexpr (n_components == 1)
3252 this->gradients[point_index] = value * normal_vector(point_index);
3253 else
3254 for (unsigned int comp = 0; comp < n_components; ++comp)
3255 this->gradients[point_index][comp] =
3256 value[comp] * normal_vector(point_index);
3257}
3258
3259
3260
3261template <int n_components_, int dim, int spacedim, typename Number>
3266 const bool is_interior,
3267 const unsigned int first_selected_component)
3268 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3269 mapping_info,
3270 fe,
3271 first_selected_component,
3272 is_interior)
3273{}
3274
3275
3276
3277template <int n_components_, int dim, int spacedim, typename Number>
3278inline void
3280 const unsigned int cell_index,
3281 const unsigned int face_number)
3282{
3283 this->current_cell_index = cell_index;
3284 this->current_face_number = face_number;
3285 this->must_reinitialize_pointers = false;
3286
3287 if (this->use_linear_path)
3288 this->template do_reinit<true, true>();
3289 else
3290 this->template do_reinit<true, false>();
3291}
3292
3293
3294
3295template <int n_components_, int dim, int spacedim, typename Number>
3296inline void
3298 const unsigned int face_index)
3299{
3300 this->current_cell_index = face_index;
3301 this->current_face_number =
3302 this->mapping_info->get_face_number(face_index, this->is_interior);
3303 this->must_reinitialize_pointers = false;
3304
3305 if (this->use_linear_path)
3306 this->template do_reinit<true, true>();
3307 else
3308 this->template do_reinit<true, false>();
3309}
3310
3311
3312
3313template <int n_components_, int dim, int spacedim, typename Number>
3314template <std::size_t stride_view>
3315void
3318 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3319{
3320 Assert(!this->must_reinitialize_pointers,
3321 ExcMessage("Object has not been reinitialized!"));
3322
3323 if (this->n_q_points == 0)
3324 return;
3325
3326 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3327
3328 if (!((evaluation_flags & EvaluationFlags::values) ||
3329 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3330 return;
3331
3332 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3333
3334 if (this->use_linear_path)
3335 do_evaluate<true>(solution_values, evaluation_flags);
3336 else
3337 do_evaluate<false>(solution_values, evaluation_flags);
3338}
3339
3340
3341
3342template <int n_components_, int dim, int spacedim, typename Number>
3343void
3345 const ArrayView<const ScalarNumber> &solution_values,
3346 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3347{
3348 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3349 solution_values.size()),
3350 evaluation_flags);
3351}
3352
3353
3354
3355template <int n_components_, int dim, int spacedim, typename Number>
3356template <bool is_linear, std::size_t stride_view>
3357void
3360 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3361{
3362 const unsigned int dofs_per_comp =
3363 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3364
3365 const ScalarNumber *input;
3366 if (stride_view == 1 && this->component_in_base_element == 0 &&
3367 (is_linear || this->renumber.empty()))
3368 input = solution_values.data();
3369 else
3370 {
3371 for (unsigned int comp = 0; comp < n_components; ++comp)
3372 {
3373 const std::size_t offset =
3374 (this->component_in_base_element + comp) * dofs_per_comp;
3375
3376 if (is_linear || this->renumber.empty())
3377 {
3378 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3379 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3380 solution_values[i + offset];
3381 }
3382 else
3383 {
3384 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3385 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3386 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3387 solution_values[renumber_ptr[i]];
3388 }
3389 }
3390 input = this->scratch_data_scalar.data();
3391 }
3392
3393 ScalarNumber *output =
3394 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3395
3396 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3397 template interpolate<true, false>(n_components,
3398 evaluation_flags,
3399 this->shape_info,
3400 input,
3401 output,
3402 this->current_face_number);
3403
3404 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3405}
3406
3407
3408
3409template <int n_components_, int dim, int spacedim, typename Number>
3410template <std::size_t stride_view>
3411void
3413 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3414 const EvaluationFlags::EvaluationFlags &integration_flags,
3415 const bool sum_into_values)
3416{
3417 Assert(!this->must_reinitialize_pointers,
3418 ExcMessage("Object has not been reinitialized!"));
3419
3420 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3421
3422 if (this->n_q_points == 0 || // no evaluation points provided
3423 !((integration_flags & EvaluationFlags::values) ||
3424 (integration_flags &
3425 EvaluationFlags::gradients))) // no integration flags
3426 {
3427 if (!sum_into_values)
3428 for (unsigned int i = 0; i < solution_values.size(); ++i)
3429 solution_values[i] = 0;
3430 return;
3431 }
3432
3433 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3434
3435 if (this->use_linear_path)
3436 do_integrate<true, true>(solution_values,
3437 integration_flags,
3438 sum_into_values);
3439 else
3440 do_integrate<true, false>(solution_values,
3441 integration_flags,
3442 sum_into_values);
3443}
3444
3445
3446
3447template <int n_components_, int dim, int spacedim, typename Number>
3448void
3450 const ArrayView<ScalarNumber> &solution_values,
3451 const EvaluationFlags::EvaluationFlags &integration_flags,
3452 const bool sum_into_values)
3453{
3454 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3455 solution_values.size()),
3456 integration_flags,
3457 sum_into_values);
3458}
3459
3460
3461
3462template <int n_components_, int dim, int spacedim, typename Number>
3463template <std::size_t stride_view>
3464void
3466 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3467 const EvaluationFlags::EvaluationFlags &integration_flags,
3468 const bool sum_into_values)
3469{
3470 Assert(!this->must_reinitialize_pointers,
3471 ExcMessage("Object has not been reinitialized!"));
3472
3473 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3474
3475 if (this->n_q_points == 0 || // no evaluation points provided
3476 !((integration_flags & EvaluationFlags::values) ||
3477 (integration_flags &
3478 EvaluationFlags::gradients))) // no integration flags
3479 {
3480 if (!sum_into_values)
3481 for (unsigned int i = 0; i < solution_values.size(); ++i)
3482 solution_values[i] = 0;
3483 return;
3484 }
3485
3486 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3487
3488 if (this->use_linear_path)
3489 do_integrate<false, true>(solution_values,
3490 integration_flags,
3491 sum_into_values);
3492 else
3493 do_integrate<false, false>(solution_values,
3494 integration_flags,
3495 sum_into_values);
3496}
3497
3498
3499
3500template <int n_components_, int dim, int spacedim, typename Number>
3501void
3503 const ArrayView<ScalarNumber> &solution_values,
3504 const EvaluationFlags::EvaluationFlags &integration_flags,
3505 const bool sum_into_values)
3506{
3507 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3508 solution_values.size()),
3509 integration_flags,
3510 sum_into_values);
3511}
3512
3513
3514
3515template <int n_components_, int dim, int spacedim, typename Number>
3516template <bool do_JxW, bool is_linear, std::size_t stride_view>
3517void
3519 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3520 const EvaluationFlags::EvaluationFlags &integration_flags,
3521 const bool sum_into_values)
3522{
3523 if (!sum_into_values && this->fe->n_components() > n_components)
3524 for (unsigned int i = 0; i < solution_values.size(); ++i)
3525 solution_values[i] = 0;
3526
3527 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3528 integration_flags,
3529 false);
3530
3531 ScalarNumber *input = this->scratch_data_scalar.begin();
3532
3533 if (stride_view == 1 && this->component_in_base_element == 0 &&
3534 (is_linear || this->renumber.empty()))
3535 {
3536 if (sum_into_values)
3537 internal::
3538 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3539 template interpolate<false, true>(n_components,
3540 integration_flags,
3541 this->shape_info,
3542 input,
3543 solution_values.data(),
3544 this->current_face_number);
3545 else
3546 internal::
3547 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3548 template interpolate<false, false>(n_components,
3549 integration_flags,
3550 this->shape_info,
3551 input,
3552 solution_values.data(),
3553 this->current_face_number);
3554 }
3555 else
3556 {
3557 const unsigned int dofs_per_comp_face =
3558 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3559
3560 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3561 ScalarNumber *output = input + size_input;
3562
3563 internal::
3564 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3565 template interpolate<false, false>(n_components,
3566 integration_flags,
3567 this->shape_info,
3568 input,
3569 output,
3570 this->current_face_number);
3571
3572 const unsigned int dofs_per_comp =
3573 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3574
3575 for (unsigned int comp = 0; comp < n_components; ++comp)
3576 {
3577 const std::size_t offset =
3578 (this->component_in_base_element + comp) * dofs_per_comp;
3579
3580 if (is_linear || this->renumber.empty())
3581 {
3582 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3583 if (sum_into_values)
3584 solution_values[i + offset] +=
3585 output[i + comp * dofs_per_comp];
3586 else
3587 solution_values[i + offset] =
3588 output[i + comp * dofs_per_comp];
3589 }
3590 else
3591 {
3592 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3593 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3594 if (sum_into_values)
3595 solution_values[renumber_ptr[i]] +=
3596 output[i + comp * dofs_per_comp];
3597 else
3598 solution_values[renumber_ptr[i]] =
3599 output[i + comp * dofs_per_comp];
3600 }
3601 }
3602 }
3603}
3604
3605
3606
3607template <int n_components_, int dim, int spacedim, typename Number>
3608template <int stride_face_dof>
3609void
3611 const ScalarNumber *face_dof_values,
3612 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3613{
3614 if (this->use_linear_path)
3615 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3616 evaluation_flags);
3617 else
3618 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3619 evaluation_flags);
3620}
3621
3622
3623
3624template <int n_components_, int dim, int spacedim, typename Number>
3625template <bool is_linear, int stride_face_dof>
3626inline void
3628 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3629 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3630{
3631 const scalar_value_type *face_dof_values_ptr;
3632 if constexpr (n_components == 1)
3633 face_dof_values_ptr = face_dof_values;
3634 else
3635 {
3636 const unsigned int dofs_per_comp_face =
3637 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3638 for (unsigned int comp = 0; comp < n_components; ++comp)
3639 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3640 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3641 stride_face_dof],
3642 comp,
3643 this->solution_renumbered[i]);
3644
3645 face_dof_values_ptr = this->solution_renumbered.data();
3646 }
3647
3648 constexpr int stride_face_dof_actual =
3649 n_components == 1 ? stride_face_dof : 1;
3650
3651 // loop over quadrature batches qb
3652 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3653
3654 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3655 {
3658
3659 if (evaluation_flags & EvaluationFlags::gradients)
3660 {
3661 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3662 is_linear ?
3664 dim - 1,
3667 2,
3668 stride_face_dof_actual>(face_dof_values_ptr,
3669 this->unit_point_faces_ptr[qb]) :
3671 dim - 1,
3674 2,
3675 false,
3676 stride_face_dof_actual>(this->shapes_faces.data() +
3677 qb * n_shapes,
3678 n_shapes,
3679 face_dof_values_ptr);
3680
3681 value = interpolated_value[dim - 1];
3682 // reorder derivative from tangential/normal derivatives into tensor
3683 // in physical coordinates
3684 if (this->current_face_number / 2 == 0)
3685 {
3686 gradient[0] = interpolated_value[dim];
3687 if (dim > 1)
3688 gradient[1] = interpolated_value[0];
3689 if (dim > 2)
3690 gradient[2] = interpolated_value[1];
3691 }
3692 else if (this->current_face_number / 2 == 1)
3693 {
3694 if (dim > 1)
3695 gradient[1] = interpolated_value[dim];
3696 if (dim == 3)
3697 {
3698 gradient[0] = interpolated_value[1];
3699 gradient[2] = interpolated_value[0];
3700 }
3701 else if (dim == 2)
3702 gradient[0] = interpolated_value[0];
3703 else
3705 }
3706 else if (this->current_face_number / 2 == 2)
3707 {
3708 if (dim > 2)
3709 {
3710 gradient[0] = interpolated_value[0];
3711 gradient[1] = interpolated_value[1];
3712 gradient[2] = interpolated_value[dim];
3713 }
3714 else
3716 }
3717 else
3719 }
3720 else
3721 {
3722 value = is_linear ?
3724 dim - 1,
3727 stride_face_dof_actual>(face_dof_values_ptr,
3728 this->unit_point_faces_ptr[qb]) :
3730 dim - 1,
3733 false,
3734 stride_face_dof_actual>(this->shapes_faces.data() +
3735 qb * n_shapes,
3736 n_shapes,
3737 face_dof_values_ptr);
3738 }
3739
3740 if (evaluation_flags & EvaluationFlags::values)
3741 {
3742 for (unsigned int v = 0, offset = qb * stride;
3743 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3744 ++v, ++offset)
3745 ETT::set_value(value, v, this->values[offset]);
3746 }
3747 if (evaluation_flags & EvaluationFlags::gradients)
3748 {
3749 Assert(this->update_flags & update_gradients ||
3750 this->update_flags & update_inverse_jacobians,
3752
3753 for (unsigned int v = 0, offset = qb * stride;
3754 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3755 ++v, ++offset)
3756 {
3757 unit_gradient_type unit_gradient;
3758 ETT::set_gradient(gradient, v, unit_gradient);
3759 this->gradients[offset] =
3760 this->cell_type <=
3763 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3765 this
3766 ->inverse_jacobian_ptr[this->cell_type <=
3768 GeometryType::affine ?
3769 0 :
3770 offset]
3771 .transpose(),
3772 unit_gradient);
3773 }
3774 }
3775 }
3776}
3777
3778
3779
3780template <int n_components_, int dim, int spacedim, typename Number>
3781template <int stride_face_dof>
3782void
3784 ScalarNumber *face_dof_values,
3785 const EvaluationFlags::EvaluationFlags &integration_flags,
3786 const bool sum_into_values)
3787{
3788 if (this->use_linear_path)
3789 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3790 integration_flags,
3791 sum_into_values);
3792 else
3793 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3794 integration_flags,
3795 sum_into_values);
3796}
3797
3798
3799
3800template <int n_components_, int dim, int spacedim, typename Number>
3801template <bool do_JxW, bool is_linear, int stride_face_dof>
3802inline void
3805 ScalarNumber *face_dof_values,
3806 const EvaluationFlags::EvaluationFlags &integration_flags,
3807 const bool sum_into_values)
3808{
3809 // zero out lanes of incomplete last quadrature point batch
3810 if constexpr (stride == 1)
3811 if (const unsigned int n_filled_lanes =
3812 this->n_q_points_scalar & (n_lanes_internal - 1);
3813 n_filled_lanes > 0)
3814 {
3815 if (integration_flags & EvaluationFlags::values)
3816 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3817 ETT::set_zero_value(this->values.back(), v);
3818 if (integration_flags & EvaluationFlags::gradients)
3819 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3820 ETT::set_zero_gradient(this->gradients.back(), v);
3821 }
3822
3823 std::array<vectorized_value_type,
3824 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3825 solution_values_vectorized_linear = {};
3826
3827 // loop over quadrature batches qb
3828 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3829
3830 const bool cartesian_cell =
3832 const bool affine_cell =
3834 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3835 {
3836 vectorized_value_type value = {};
3838
3839 if (integration_flags & EvaluationFlags::values)
3840 for (unsigned int v = 0, offset = qb * stride;
3841 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3842 ++v, ++offset)
3843 ETT::get_value(value,
3844 v,
3845 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3846 this->values[offset]);
3847
3848 if (integration_flags & EvaluationFlags::gradients)
3849 for (unsigned int v = 0, offset = qb * stride;
3850 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3851 ++v, ++offset)
3852 {
3853 const auto grad_w =
3854 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3855 this->gradients[offset];
3856 ETT::get_gradient(
3857 gradient,
3858 v,
3859 cartesian_cell ?
3860 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3861 grad_w) :
3863 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3864 grad_w));
3865 }
3866
3867 if (integration_flags & EvaluationFlags::gradients)
3868 {
3869 std::array<vectorized_value_type, 2> value_face = {};
3870 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3871
3872 value_face[0] = value;
3873 // fill derivative in physical coordinates into tangential/normal
3874 // derivatives
3875 if (this->current_face_number / 2 == 0)
3876 {
3877 value_face[1] = gradient[0];
3878 if (dim > 1)
3879 gradient_in_face[0] = gradient[1];
3880 if (dim > 2)
3881 gradient_in_face[1] = gradient[2];
3882 }
3883 else if (this->current_face_number / 2 == 1)
3884 {
3885 if (dim > 1)
3886 value_face[1] = gradient[1];
3887 if (dim == 3)
3888 {
3889 gradient_in_face[0] = gradient[2];
3890 gradient_in_face[1] = gradient[0];
3891 }
3892 else if (dim == 2)
3893 gradient_in_face[0] = gradient[0];
3894 else
3896 }
3897 else if (this->current_face_number / 2 == 2)
3898 {
3899 if (dim > 2)
3900 {
3901 value_face[1] = gradient[2];
3902 gradient_in_face[0] = gradient[0];
3903 gradient_in_face[1] = gradient[1];
3904 }
3905 else
3907 }
3908 else
3910
3912 is_linear,
3913 dim - 1,
3916 2>(this->shapes_faces.data() + qb * n_shapes,
3917 n_shapes,
3918 value_face.data(),
3919 gradient_in_face,
3920 is_linear ? solution_values_vectorized_linear.data() :
3921 this->solution_renumbered_vectorized.data(),
3922 this->unit_point_faces_ptr[qb],
3923 qb != 0);
3924 }
3925 else
3927 dim - 1,
3930 this->shapes_faces.data() + qb * n_shapes,
3931 n_shapes,
3932 value,
3933 is_linear ? solution_values_vectorized_linear.data() :
3934 this->solution_renumbered_vectorized.data(),
3935 this->unit_point_faces_ptr[qb],
3936 qb != 0);
3937 }
3938
3939 const unsigned int dofs_per_comp_face =
3940 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3941
3942 for (unsigned int comp = 0; comp < n_components; ++comp)
3943 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3944 if (sum_into_values)
3945 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3946 stride_face_dof] +=
3947 ETT::sum_value(comp,
3948 is_linear ?
3949 *(solution_values_vectorized_linear.data() + i) :
3950 this->solution_renumbered_vectorized[i]);
3951 else
3952 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3953 ETT::sum_value(comp,
3954 is_linear ?
3955 *(solution_values_vectorized_linear.data() + i) :
3956 this->solution_renumbered_vectorized[i]);
3957}
3958
3959
3960
3961template <int n_components_, int dim, int spacedim, typename Number>
3964 const unsigned int point_index) const
3965{
3966 AssertIndexRange(point_index, this->n_q_points);
3967 Assert(this->normal_ptr != nullptr,
3969 ExcFEPointEvaluationAccessToUninitializedMappingField(
3970 "update_normal_vectors"));
3971 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3972 {
3974 for (unsigned int d = 0; d < dim; ++d)
3975 normal[d] =
3976 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3977 0);
3978
3979 return normal;
3980 }
3981 else
3982 {
3983 return this->normal_ptr[point_index];
3984 }
3985}
3986
3987
3988
3989template <int n_components_, int dim, int spacedim, typename Number>
3990inline const typename FEFacePointEvaluation<n_components_,
3991 dim,
3992 spacedim,
3993 Number>::value_type
3995 get_normal_derivative(const unsigned int point_index) const
3996{
3997 AssertIndexRange(point_index, this->gradients.size());
3998
3999 value_type normal_derivative;
4000 if constexpr (n_components == 1)
4001 normal_derivative =
4002 this->gradients[point_index] * normal_vector(point_index);
4003 else
4004 for (unsigned int comp = 0; comp < n_components; ++comp)
4005 normal_derivative[comp] =
4006 this->gradients[point_index][comp] * normal_vector(point_index);
4007
4008 return normal_derivative;
4009}
4010
4011
4012
4013template <int n_components_, int dim, int spacedim, typename Number>
4014inline void
4017 const unsigned int point_index)
4018{
4019 AssertIndexRange(point_index, this->gradients.size());
4020 if constexpr (n_components == 1)
4021 this->gradients[point_index] = value * normal_vector(point_index);
4022 else
4023 for (unsigned int comp = 0; comp < n_components; ++comp)
4024 this->gradients[point_index][comp] =
4025 value[comp] * normal_vector(point_index);
4026}
4027
4029
4030#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
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
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0, const bool force_lexicographic_numbering=false)
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:281
#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 > &)