Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_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
97 static scalar_value_type
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 jacobian(const unsigned int point_index) const;
783
791 inverse_jacobian(const unsigned int point_index) const;
792
798 Number
799 JxW(const unsigned int point_index) const;
800
807 DEAL_II_DEPRECATED_EARLY Point<spacedim, Number>
808 real_point(const unsigned int point_index) const;
809
815 quadrature_point(const unsigned int point_index) const;
816
822 unit_point(const unsigned int point_index) const;
823
833
841
845 unsigned int
847
848protected:
849 static constexpr std::size_t n_lanes_user_interface =
851 static constexpr std::size_t n_lanes_internal =
853 static constexpr std::size_t stride =
855
864 void
865 setup(const unsigned int first_selected_component);
866
872 template <bool is_face, bool is_linear>
873 void
875
879 const unsigned int n_q_batches;
880
884 const unsigned int n_q_points;
885
889 const unsigned int n_q_points_scalar;
890
895
900
905 std::vector<Polynomials::Polynomial<double>> poly;
906
911
916 std::vector<unsigned int> renumber;
917
925 std::vector<scalar_value_type> solution_renumbered;
926
934
939
943 std::vector<value_type> values;
944
948 std::vector<gradient_type> gradients;
949
955
961
967
973
979
985
991 const Number *JxW_ptr;
992
997
1003
1010
1015
1020
1026 std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
1027
1032
1036 std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1037
1041 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1043
1049
1054
1059
1064
1072
1078
1084
1085 const bool is_interior;
1086};
1087
1088
1089
1120template <int n_components_,
1121 int dim,
1122 int spacedim = dim,
1123 typename Number = double>
1125 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1126{
1127public:
1128 static constexpr unsigned int dimension = dim;
1129 static constexpr unsigned int n_components = n_components_;
1130
1131 using number_type = Number;
1132
1135 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1137 using ETT = typename internal::FEPointEvaluation::
1138 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1139 using value_type = typename ETT::value_type;
1140 using scalar_value_type = typename ETT::scalar_value_type;
1141 using vectorized_value_type = typename ETT::vectorized_value_type;
1142 using unit_gradient_type = typename ETT::unit_gradient_type;
1143 using gradient_type = typename ETT::real_gradient_type;
1145 typename ETT::interface_vectorized_unit_gradient_type;
1146
1168 const unsigned int first_selected_component = 0);
1169
1189 const unsigned int first_selected_component = 0);
1190
1202 void
1204 const ArrayView<const Point<dim>> &unit_points);
1205
1210 void
1211 reinit();
1212
1217 void
1218 reinit(const unsigned int cell_index);
1219
1220
1232 template <std::size_t stride_view>
1233 void
1234 evaluate(
1236 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1237
1249 void
1250 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1251 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1252
1275 template <std::size_t stride_view>
1276 void
1278 const EvaluationFlags::EvaluationFlags &integration_flags,
1279 const bool sum_into_values = false);
1280
1303 void
1304 integrate(const ArrayView<ScalarNumber> &solution_values,
1305 const EvaluationFlags::EvaluationFlags &integration_flags,
1306 const bool sum_into_values = false);
1307
1334 template <std::size_t stride_view>
1335 void
1337 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1338 const EvaluationFlags::EvaluationFlags &integration_flags,
1339 const bool sum_into_values = false);
1340
1367 void
1368 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1369 const EvaluationFlags::EvaluationFlags &integration_flags,
1370 const bool sum_into_values = false);
1371
1378 normal_vector(const unsigned int point_index) const;
1379
1385 const value_type
1386 get_normal_derivative(const unsigned int point_index) const;
1387
1393 void
1394 submit_normal_derivative(const value_type &, const unsigned int point_index);
1395
1396private:
1397 static constexpr std::size_t n_lanes_user_interface =
1399 static constexpr std::size_t n_lanes_internal =
1401 static constexpr std::size_t stride =
1403
1408 template <bool is_linear, std::size_t stride_view>
1409 void
1412
1417 template <bool is_linear, std::size_t stride_view>
1418 void
1421 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1422 const unsigned int n_shapes,
1423 const unsigned int qb,
1424 vectorized_value_type &value,
1426
1430 template <bool is_linear, std::size_t stride_view>
1431 void
1434 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1435
1439 template <std::size_t stride_view>
1440 void
1443 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1444
1450 template <bool is_linear>
1451 void
1453 const EvaluationFlags::EvaluationFlags &integration_flags,
1454 const unsigned int n_shapes,
1455 const unsigned int qb,
1456 const vectorized_value_type value,
1458 vectorized_value_type *solution_values_vectorized_linear);
1459
1465 template <bool is_linear, std::size_t stride_view>
1466 void
1468 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1469 vectorized_value_type *solution_values_vectorized_linear,
1470 const bool sum_into_values);
1471
1475 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1476 void
1478 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1479 const EvaluationFlags::EvaluationFlags &integration_flags,
1480 const bool sum_into_values);
1481
1485 template <bool do_JxW, std::size_t stride_view>
1486 void
1488 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1489 const EvaluationFlags::EvaluationFlags &integration_flags,
1490 const bool sum_into_values);
1491
1495 template <bool do_JxW, std::size_t stride_view>
1496 void
1498 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1499 const EvaluationFlags::EvaluationFlags &integration_flags,
1500 const bool sum_into_values);
1501
1507 void
1509};
1510
1511
1512
1530template <int n_components_,
1531 int dim,
1532 int spacedim = dim,
1533 typename Number = double>
1535 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1536{
1537public:
1538 static constexpr unsigned int dimension = dim;
1539 static constexpr unsigned int n_components = n_components_;
1540
1541 using number_type = Number;
1542
1545 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1547 using ETT = typename internal::FEPointEvaluation::
1548 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1549 using value_type = typename ETT::value_type;
1550 using scalar_value_type = typename ETT::scalar_value_type;
1551 using vectorized_value_type = typename ETT::vectorized_value_type;
1552 using unit_gradient_type = typename ETT::unit_gradient_type;
1553 using gradient_type = typename ETT::real_gradient_type;
1555 typename ETT::interface_vectorized_unit_gradient_type;
1556
1563 const bool is_interior = true,
1564 const unsigned int first_selected_component = 0);
1565
1570 void
1571 reinit(const unsigned int cell_index, const unsigned int face_number);
1572
1577 void
1578 reinit(const unsigned int face_index);
1579
1591 template <std::size_t stride_view>
1592 void
1593 evaluate(
1595 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1596
1608 void
1609 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1610 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1611
1634 template <std::size_t stride_view>
1635 void
1637 const EvaluationFlags::EvaluationFlags &integration_flags,
1638 const bool sum_into_values = false);
1639
1662 void
1663 integrate(const ArrayView<ScalarNumber> &solution_values,
1664 const EvaluationFlags::EvaluationFlags &integration_flags,
1665 const bool sum_into_values = false);
1666
1689 template <std::size_t stride_view>
1690 void
1692 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1693 const EvaluationFlags::EvaluationFlags &integration_flags,
1694 const bool sum_into_values = false);
1695
1718 void
1719 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1720 const EvaluationFlags::EvaluationFlags &integration_flags,
1721 const bool sum_into_values = false);
1722
1729 template <int stride_face_dof = VectorizedArrayType::size()>
1730 void
1731 evaluate_in_face(const ScalarNumber *face_dof_values,
1732 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1733
1740 template <int stride_face_dof = VectorizedArrayType::size()>
1741 void
1742 integrate_in_face(ScalarNumber *face_dof_values,
1743 const EvaluationFlags::EvaluationFlags &integration_flags,
1744 const bool sum_into_values = false);
1745
1752 normal_vector(const unsigned int point_index) const;
1753
1759 const value_type
1760 get_normal_derivative(const unsigned int point_index) const;
1761
1767 void
1768 submit_normal_derivative(const value_type &, const unsigned int point_index);
1769
1770private:
1771 static constexpr std::size_t n_lanes_user_interface =
1773 static constexpr std::size_t n_lanes_internal =
1775 static constexpr std::size_t stride =
1777
1778 template <bool is_linear, std::size_t stride_view>
1779 void
1782 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1783
1784 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1785 void
1787 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1788 const EvaluationFlags::EvaluationFlags &integration_flags,
1789 const bool sum_into_values);
1790
1795 template <bool is_linear, int stride_face_dof>
1796 void
1797 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1798 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1799
1804 template <bool do_JxW, bool is_linear, int stride_face_dof>
1805 void
1807 ScalarNumber *face_dof_values,
1808 const EvaluationFlags::EvaluationFlags &integration_flags,
1809 const bool sum_into_values);
1810};
1811
1812
1813
1814// ----------------------- template and inline function ----------------------
1815
1816
1817template <int n_components_, int dim, int spacedim, typename Number>
1821 const UpdateFlags update_flags,
1822 const unsigned int first_selected_component)
1823 : n_q_batches(numbers::invalid_unsigned_int)
1824 , n_q_points(numbers::invalid_unsigned_int)
1825 , n_q_points_scalar(numbers::invalid_unsigned_int)
1826 , mapping(&mapping)
1827 , fe(&fe)
1828 , JxW_ptr(nullptr)
1829 , update_flags(update_flags)
1830 , mapping_info_on_the_fly(
1831 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1832 mapping,
1833 update_flags))
1834 , mapping_info(mapping_info_on_the_fly.get())
1835 , current_cell_index(numbers::invalid_unsigned_int)
1836 , current_face_number(numbers::invalid_unsigned_int)
1837 , must_reinitialize_pointers(false)
1838 , is_interior(true)
1839{
1840 setup(first_selected_component);
1841}
1842
1843
1844
1845template <int n_components_, int dim, int spacedim, typename Number>
1850 const unsigned int first_selected_component,
1851 const bool is_interior)
1852 : n_q_batches(numbers::invalid_unsigned_int)
1853 , n_q_points(numbers::invalid_unsigned_int)
1854 , n_q_points_scalar(numbers::invalid_unsigned_int)
1855 , mapping(&mapping_info.get_mapping())
1856 , fe(&fe)
1857 , JxW_ptr(nullptr)
1858 , update_flags(mapping_info.get_update_flags())
1859 , mapping_info(&mapping_info)
1860 , current_cell_index(numbers::invalid_unsigned_int)
1861 , current_face_number(numbers::invalid_unsigned_int)
1862 , must_reinitialize_pointers(true)
1863 , is_interior(is_interior)
1864{
1865 setup(first_selected_component);
1866}
1867
1868
1869
1870template <int n_components_, int dim, int spacedim, typename Number>
1874 : n_q_batches(other.n_q_batches)
1875 , n_q_points(other.n_q_points)
1876 , n_q_points_scalar(other.n_q_points_scalar)
1877 , mapping(other.mapping)
1878 , fe(other.fe)
1879 , poly(other.poly)
1880 , use_linear_path(other.use_linear_path)
1881 , renumber(other.renumber)
1882 , solution_renumbered(other.solution_renumbered)
1883 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1884 , values(other.values)
1885 , gradients(other.gradients)
1886 , dofs_per_component(other.dofs_per_component)
1887 , dofs_per_component_face(other.dofs_per_component_face)
1888 , component_in_base_element(other.component_in_base_element)
1889 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1890 , update_flags(other.update_flags)
1891 , fe_values(other.fe_values)
1892 , mapping_info_on_the_fly(
1893 other.mapping_info_on_the_fly ?
1895 *mapping,
1896 update_flags) :
1897 nullptr)
1898 , mapping_info(other.mapping_info)
1899 , current_cell_index(other.current_cell_index)
1900 , current_face_number(other.current_face_number)
1901 , fast_path(other.fast_path)
1902 , must_reinitialize_pointers(true)
1903 , is_interior(other.is_interior)
1904{}
1905
1906
1907
1908template <int n_components_, int dim, int spacedim, typename Number>
1911 : n_q_batches(other.n_q_batches)
1912 , n_q_points(other.n_q_points)
1913 , n_q_points_scalar(other.n_q_points_scalar)
1914 , mapping(other.mapping)
1915 , fe(other.fe)
1916 , poly(other.poly)
1917 , use_linear_path(other.use_linear_path)
1918 , renumber(other.renumber)
1919 , solution_renumbered(other.solution_renumbered)
1920 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1921 , values(other.values)
1922 , gradients(other.gradients)
1923 , dofs_per_component(other.dofs_per_component)
1924 , dofs_per_component_face(other.dofs_per_component_face)
1925 , component_in_base_element(other.component_in_base_element)
1926 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1927 , update_flags(other.update_flags)
1928 , fe_values(other.fe_values)
1929 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1930 , mapping_info(other.mapping_info)
1931 , current_cell_index(other.current_cell_index)
1932 , current_face_number(other.current_face_number)
1933 , fast_path(other.fast_path)
1934 , must_reinitialize_pointers(other.must_reinitialize_pointers)
1935 , is_interior(other.is_interior)
1936{}
1937
1938
1939
1940template <int n_components_, int dim, int spacedim, typename Number>
1941void
1943 const unsigned int first_selected_component)
1944{
1945 AssertIndexRange(first_selected_component + n_components,
1946 fe->n_components() + 1);
1947
1948 shapes.reserve(100);
1949
1950 bool same_base_element = true;
1951 unsigned int base_element_number = 0;
1952 component_in_base_element = 0;
1953 unsigned int component = 0;
1954 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1955 if (component + fe->element_multiplicity(base_element_number) >
1956 first_selected_component)
1957 {
1958 if (first_selected_component + n_components >
1959 component + fe->element_multiplicity(base_element_number))
1960 same_base_element = false;
1961 component_in_base_element = first_selected_component - component;
1962 break;
1963 }
1964 else
1965 component += fe->element_multiplicity(base_element_number);
1966
1969 *fe, base_element_number) &&
1970 same_base_element)
1971 {
1972 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1973 renumber = shape_info.lexicographic_numbering;
1974 dofs_per_component = shape_info.dofs_per_component_on_cell;
1975 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1977 fe->base_element(base_element_number));
1978
1979 bool is_lexicographic = true;
1980 for (unsigned int i = 0; i < renumber.size(); ++i)
1981 if (i != renumber[i])
1982 is_lexicographic = false;
1983
1984 if (is_lexicographic)
1985 renumber.clear();
1986
1987 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1988 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1989 poly[1].value(1.) == 1.) &&
1990 (fe->n_components() == n_components);
1991
1992 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1993 const unsigned int size_cell = dofs_per_component * n_components;
1994 scratch_data_scalar.resize(size_face + size_cell);
1995
1996 solution_renumbered.resize(dofs_per_component);
1997 solution_renumbered_vectorized.resize(dofs_per_component);
1998
1999 fast_path = true;
2000 }
2001 else
2002 {
2003 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
2004 for (unsigned int d = 0; d < n_components; ++d)
2005 {
2006 const unsigned int component = first_selected_component + d;
2007 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
2008 {
2009 const bool is_primitive =
2010 fe->is_primitive() || fe->is_primitive(i);
2011 if (is_primitive)
2012 nonzero_shape_function_component[i][d] =
2013 (component == fe->system_to_component_index(i).first);
2014 else
2015 nonzero_shape_function_component[i][d] =
2016 (fe->get_nonzero_components(i)[component] == true);
2017 }
2018 }
2019
2020 fast_path = false;
2021 }
2022}
2023
2024
2025
2026template <int n_components_, int dim, int spacedim, typename Number>
2027template <bool is_face, bool is_linear>
2028inline void
2030{
2031 const unsigned int geometry_index =
2032 mapping_info->template compute_geometry_index_offset<is_face>(
2033 current_cell_index, current_face_number);
2034
2035 cell_type = mapping_info->get_cell_type(geometry_index);
2036
2037 const_cast<unsigned int &>(n_q_points_scalar) =
2038 mapping_info->get_n_q_points_unvectorized(geometry_index);
2039
2040 // round up n_q_points_scalar / n_lanes_internal
2041 const_cast<unsigned int &>(n_q_batches) =
2042 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
2043
2044 const unsigned int n_q_points_before = n_q_points;
2045
2046 const_cast<unsigned int &>(n_q_points) =
2047 (stride == 1) ? n_q_batches : n_q_points_scalar;
2048
2049 if (n_q_points != n_q_points_before)
2050 {
2051 if (update_flags & update_values)
2052 values.resize(n_q_points);
2053 if (update_flags & update_gradients)
2054 gradients.resize(n_q_points);
2055 }
2056
2057 if (n_q_points == 0)
2058 return;
2059
2060 // set unit point pointer
2061 const unsigned int unit_point_offset =
2062 mapping_info->compute_unit_point_index_offset(geometry_index);
2063
2064 if (is_face)
2065 unit_point_faces_ptr =
2066 mapping_info->get_unit_point_faces(unit_point_offset);
2067 else
2068 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2069
2070 // set data pointers
2071 const unsigned int data_offset =
2072 mapping_info->compute_data_index_offset(geometry_index);
2073 const unsigned int compressed_data_offset =
2074 mapping_info->compute_compressed_data_index_offset(geometry_index);
2075#ifdef DEBUG
2076 const UpdateFlags update_flags_mapping =
2077 mapping_info->get_update_flags_mapping();
2078 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2079 real_point_ptr = mapping_info->get_real_point(data_offset);
2080 if (update_flags_mapping & UpdateFlags::update_jacobians)
2081 jacobian_ptr =
2082 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2083 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2084 inverse_jacobian_ptr =
2085 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2086 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2087 normal_ptr = mapping_info->get_normal_vector(data_offset);
2088 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2089 JxW_ptr = mapping_info->get_JxW(data_offset);
2090#else
2091 real_point_ptr = mapping_info->get_real_point(data_offset);
2092 jacobian_ptr =
2093 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2094 inverse_jacobian_ptr =
2095 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2096 normal_ptr = mapping_info->get_normal_vector(data_offset);
2097 JxW_ptr = mapping_info->get_JxW(data_offset);
2098#endif
2099
2100 if (!is_linear && fast_path)
2101 {
2102 const std::size_t n_shapes = poly.size();
2103 if (is_face)
2104 shapes_faces.resize_fast(n_q_batches * n_shapes);
2105 else
2106 shapes.resize_fast(n_q_batches * n_shapes);
2107
2108 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2109 if (is_face)
2110 {
2111 if (dim > 1)
2112 {
2114 shapes_faces.data() + qb * n_shapes,
2115 poly,
2116 unit_point_faces_ptr[qb],
2117 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2118 }
2119 }
2120 else
2121 {
2122 if (update_flags & UpdateFlags::update_gradients)
2123 {
2124 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2125 poly,
2126 unit_point_ptr[qb],
2127 1);
2128 }
2129 else if (qb + 1 < n_q_batches)
2130 {
2131 // Use function with reduced overhead to compute for two
2132 // points at once
2134 shapes.data() + qb * n_shapes,
2135 poly,
2136 unit_point_ptr[qb],
2137 unit_point_ptr[qb + 1]);
2138 ++qb;
2139 }
2140 else
2141 {
2142 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2143 poly,
2144 unit_point_ptr[qb],
2145 0);
2146 }
2147 }
2148 }
2149}
2150
2151
2152
2153template <int n_components_, int dim, int spacedim, typename Number>
2154inline const typename FEPointEvaluationBase<n_components_,
2155 dim,
2156 spacedim,
2157 Number>::value_type &
2159 const unsigned int point_index) const
2160{
2161 AssertIndexRange(point_index, values.size());
2162 return values[point_index];
2163}
2164
2165
2166
2167template <int n_components_, int dim, int spacedim, typename Number>
2168inline const typename FEPointEvaluationBase<n_components_,
2169 dim,
2170 spacedim,
2171 Number>::gradient_type &
2173 const unsigned int point_index) const
2174{
2175 AssertIndexRange(point_index, gradients.size());
2176 return gradients[point_index];
2177}
2178
2179
2180
2181template <int n_components_, int dim, int spacedim, typename Number>
2182inline Number
2184 const unsigned int point_index) const
2185{
2186 static_assert(n_components == dim,
2187 "Only makes sense for a vector field with dim components");
2188
2189 AssertIndexRange(point_index, values.size());
2190 return trace(gradients[point_index]);
2191}
2192
2193
2194
2195template <int n_components_, int dim, int spacedim, typename Number>
2196inline void
2198 const value_type &value,
2199 const unsigned int point_index)
2200{
2201 AssertIndexRange(point_index, n_q_points);
2202 values[point_index] = value;
2203}
2204
2205
2206
2207template <int n_components_, int dim, int spacedim, typename Number>
2208inline void
2210 const gradient_type &gradient,
2211 const unsigned int point_index)
2212{
2213 AssertIndexRange(point_index, n_q_points);
2214 gradients[point_index] = gradient;
2215}
2216
2217
2218
2219template <int n_components_, int dim, int spacedim, typename Number>
2220inline void
2222 const Number &value,
2223 const unsigned int point_index)
2224{
2225 static_assert(n_components == dim,
2226 "Only makes sense for a vector field with dim components");
2227
2228 AssertIndexRange(point_index, n_q_points);
2229 gradients[point_index] = gradient_type();
2230 for (unsigned int d = 0; d < dim; ++d)
2231 gradients[point_index][d][d] = value;
2232}
2233
2234
2235
2236template <int n_components_, int dim, int spacedim, typename Number>
2239 const unsigned int point_index) const
2240{
2241 AssertIndexRange(point_index, n_q_points);
2242 Assert(jacobian_ptr != nullptr,
2244 ExcFEPointEvaluationAccessToUninitializedMappingField(
2245 "update_jacobians"));
2246 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2247 GeometryType::affine ?
2248 0 :
2249 point_index];
2250}
2251
2252
2253
2254template <int n_components_, int dim, int spacedim, typename Number>
2257 const unsigned int point_index) const
2258{
2259 AssertIndexRange(point_index, n_q_points);
2260 Assert(inverse_jacobian_ptr != nullptr,
2262 ExcFEPointEvaluationAccessToUninitializedMappingField(
2263 "update_inverse_jacobians"));
2264 return inverse_jacobian_ptr
2265 [cell_type <=
2267 0 :
2268 point_index];
2269}
2270
2271
2272
2273template <int n_components_, int dim, int spacedim, typename Number>
2274inline Number
2276 const unsigned int point_index) const
2277{
2278 AssertIndexRange(point_index, n_q_points);
2279 Assert(JxW_ptr != nullptr,
2281 ExcFEPointEvaluationAccessToUninitializedMappingField(
2282 "update_JxW_values"));
2283 return JxW_ptr[point_index];
2284}
2285
2286
2287
2288template <int n_components_, int dim, int spacedim, typename Number>
2291 const unsigned int point_index) const
2292{
2293 return quadrature_point(point_index);
2294 ;
2295}
2296
2297
2298
2299template <int n_components_, int dim, int spacedim, typename Number>
2302 const unsigned int point_index) const
2303{
2304 AssertIndexRange(point_index, n_q_points);
2305 Assert(real_point_ptr != nullptr,
2307 ExcFEPointEvaluationAccessToUninitializedMappingField(
2308 "update_quadrature_points"));
2309 return real_point_ptr[point_index];
2310}
2311
2312
2313
2314template <int n_components_, int dim, int spacedim, typename Number>
2315inline Point<dim, Number>
2317 const unsigned int point_index) const
2318{
2319 AssertIndexRange(point_index, n_q_points);
2320 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2321 Point<dim, Number> unit_point;
2322 for (unsigned int d = 0; d < dim; ++d)
2324 unit_point_ptr[point_index / stride][d], point_index % stride);
2325 return unit_point;
2326}
2327
2328
2329
2330template <int n_components_, int dim, int spacedim, typename Number>
2337
2338
2339
2340template <int n_components_, int dim, int spacedim, typename Number>
2344 const unsigned int first_selected_component)
2345 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2346 mapping_info,
2347 fe,
2348 first_selected_component)
2349{}
2350
2351
2352
2353template <int n_components_, int dim, int spacedim, typename Number>
2355 const Mapping<dim, spacedim> &mapping,
2357 const UpdateFlags update_flags,
2358 const unsigned int first_selected_component)
2359 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2360 mapping,
2361 fe,
2362 update_flags,
2363 first_selected_component)
2364{}
2365
2366
2367
2368template <int n_components_, int dim, int spacedim, typename Number>
2369inline void
2372{
2373 this->current_cell_index = numbers::invalid_unsigned_int;
2374 this->current_face_number = numbers::invalid_unsigned_int;
2375
2376 if (this->use_linear_path)
2377 this->template do_reinit<false, true>();
2378 else
2379 this->template do_reinit<false, false>();
2380}
2381
2382
2383
2384template <int n_components_, int dim, int spacedim, typename Number>
2385inline void
2387{
2388 internal_reinit_single_cell_state_mapping_info();
2389 this->must_reinitialize_pointers = false;
2390}
2391
2392
2393
2394template <int n_components_, int dim, int spacedim, typename Number>
2395inline void
2398 const ArrayView<const Point<dim>> &unit_points)
2399{
2400 // reinit is only allowed for mapping computation on the fly
2401 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2403
2404 this->mapping_info->reinit(cell, unit_points);
2405 this->must_reinitialize_pointers = false;
2406
2407 if (!this->fast_path)
2408 {
2409 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2410 *this->mapping,
2411 *this->fe,
2413 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2414 this->update_flags);
2415 this->fe_values->reinit(cell);
2416 }
2417
2418 if (this->use_linear_path)
2419 this->template do_reinit<false, true>();
2420 else
2421 this->template do_reinit<false, false>();
2422}
2423
2424
2425
2426template <int n_components_, int dim, int spacedim, typename Number>
2427inline void
2429 const unsigned int cell_index)
2430{
2431 this->current_cell_index = cell_index;
2432 this->current_face_number = numbers::invalid_unsigned_int;
2433 this->must_reinitialize_pointers = false;
2434
2435 if (this->use_linear_path)
2436 this->template do_reinit<false, true>();
2437 else
2438 this->template do_reinit<false, false>();
2439
2440 if (!this->fast_path)
2441 {
2442 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2443
2444 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2445 for (unsigned int d = 0; d < dim; ++d)
2446 unit_points[v][d] =
2447 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2448
2449 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2450 *this->mapping,
2451 *this->fe,
2453 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2454 this->update_flags);
2455
2456 this->fe_values->reinit(
2457 this->mapping_info->get_cell_iterator(this->current_cell_index));
2458 }
2459}
2460
2461
2462
2463template <int n_components_, int dim, int spacedim, typename Number>
2464template <std::size_t stride_view>
2465void
2468 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2469{
2470 if (this->must_reinitialize_pointers)
2471 internal_reinit_single_cell_state_mapping_info();
2472
2473 if (this->n_q_points == 0)
2474 return;
2475
2476 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2477
2478 if (!((evaluation_flags & EvaluationFlags::values) ||
2479 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2480 return;
2481
2482 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2483 if (this->fast_path)
2484 {
2485 if (this->use_linear_path)
2486 evaluate_fast<true>(solution_values, evaluation_flags);
2487 else
2488 evaluate_fast<false>(solution_values, evaluation_flags);
2489 }
2490 else
2491 evaluate_slow(solution_values, evaluation_flags);
2492}
2493
2494
2495
2496template <int n_components_, int dim, int spacedim, typename Number>
2497void
2499 const ArrayView<const ScalarNumber> &solution_values,
2500 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2501{
2502 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2503 solution_values.size()),
2504 evaluation_flags);
2505}
2506
2507
2508
2509template <int n_components_, int dim, int spacedim, typename Number>
2510template <std::size_t stride_view>
2511void
2513 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2514 const EvaluationFlags::EvaluationFlags &integration_flags,
2515 const bool sum_into_values)
2516{
2517 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2518}
2519
2520
2521
2522template <int n_components_, int dim, int spacedim, typename Number>
2523void
2525 const ArrayView<ScalarNumber> &solution_values,
2526 const EvaluationFlags::EvaluationFlags &integration_flags,
2527 const bool sum_into_values)
2528{
2529 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2530 solution_values.size()),
2531 integration_flags,
2532 sum_into_values);
2533}
2534
2535
2536
2537template <int n_components_, int dim, int spacedim, typename Number>
2539 scalar_value_type
2541 const
2542{
2543 value_type return_value = {};
2544
2545 for (const auto point_index : this->quadrature_point_indices())
2546 return_value += values[point_index] * this->JxW(point_index);
2547
2548 return ETT::sum_value(return_value);
2549}
2550
2551
2552
2553template <int n_components_, int dim, int spacedim, typename Number>
2554unsigned int
2557{
2558 Assert(stride == 1,
2559 ExcMessage(
2560 "Calling this function only makes sense in fully vectorized mode."));
2561 if (q == n_q_batches - 1)
2562 {
2563 const unsigned int n_filled_lanes =
2564 n_q_points_scalar & (n_lanes_user_interface - 1);
2565 if (n_filled_lanes == 0)
2566 return n_lanes_user_interface;
2567 else
2568 return n_filled_lanes;
2569 }
2570 else
2571 return n_lanes_user_interface;
2572}
2573
2574
2575
2576template <int n_components_, int dim, int spacedim, typename Number>
2577template <std::size_t stride_view>
2578void
2580 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2581 const EvaluationFlags::EvaluationFlags &integration_flags,
2582 const bool sum_into_values)
2583{
2584 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2585}
2586
2587
2588
2589template <int n_components_, int dim, int spacedim, typename Number>
2590void
2592 const ArrayView<ScalarNumber> &solution_values,
2593 const EvaluationFlags::EvaluationFlags &integration_flags,
2594 const bool sum_into_values)
2595{
2596 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2597 solution_values.size()),
2598 integration_flags,
2599 sum_into_values);
2600}
2601
2602
2603
2604template <int n_components_, int dim, int spacedim, typename Number>
2605template <bool is_linear, std::size_t stride_view>
2606inline void
2609{
2610 const unsigned int dofs_per_comp =
2611 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2612
2613 for (unsigned int comp = 0; comp < n_components; ++comp)
2614 {
2615 const std::size_t offset =
2616 (this->component_in_base_element + comp) * dofs_per_comp;
2617
2618 if ((is_linear && n_components == 1) || this->renumber.empty())
2619 {
2620 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2621 ETT::read_value(solution_values[i + offset],
2622 comp,
2623 this->solution_renumbered[i]);
2624 }
2625 else
2626 {
2627 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2628 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2629 ETT::read_value(solution_values[renumber_ptr[i]],
2630 comp,
2631 this->solution_renumbered[i]);
2632 }
2633 }
2634}
2635
2636
2637
2638template <int n_components_, int dim, int spacedim, typename Number>
2639template <bool is_linear, std::size_t stride_view>
2640inline void
2643 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2644 const unsigned int n_shapes,
2645 const unsigned int qb,
2646 vectorized_value_type &value,
2648{
2649 if (evaluation_flags & EvaluationFlags::gradients)
2650 {
2651 std::array<vectorized_value_type, dim + 1> result;
2652 if constexpr (is_linear)
2653 {
2654 if constexpr (n_components == 1)
2655 result =
2657 dim,
2660 1,
2661 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2662 else
2663 result =
2665 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2666 }
2667 else
2669 dim,
2672 1,
2673 false>(this->shapes.data() + qb * n_shapes,
2674 n_shapes,
2675 this->solution_renumbered.data());
2676 gradient[0] = result[0];
2677 if (dim > 1)
2678 gradient[1] = result[1];
2679 if (dim > 2)
2680 gradient[2] = result[2];
2681 value = result[dim];
2682 }
2683 else
2684 {
2685 if constexpr (is_linear)
2686 {
2687 if constexpr (n_components == 1)
2689 dim,
2692 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2693 else
2695 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2696 }
2697 else
2698 value =
2702 false>(
2703 this->shapes.data() + qb * n_shapes,
2704 n_shapes,
2705 this->solution_renumbered.data());
2706 }
2707}
2708
2709
2710
2711template <int n_components_, int dim, int spacedim, typename Number>
2712template <bool is_linear, std::size_t stride_view>
2713inline void
2716 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2717{
2718 if (!(is_linear && n_components == 1))
2719 prepare_evaluate_fast<is_linear>(solution_values);
2720
2721 // loop over quadrature batches qb
2722 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2723
2724 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2725 {
2728
2729 compute_evaluate_fast<is_linear>(
2730 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2731
2732 if (evaluation_flags & EvaluationFlags::values)
2733 {
2734 for (unsigned int v = 0, offset = qb * stride;
2735 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2736 ++v, ++offset)
2737 ETT::set_value(value, v, this->values[offset]);
2738 }
2739 if (evaluation_flags & EvaluationFlags::gradients)
2740 {
2741 Assert(this->update_flags & update_gradients ||
2742 this->update_flags & update_inverse_jacobians,
2744
2745 for (unsigned int v = 0, offset = qb * stride;
2746 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2747 ++v, ++offset)
2748 {
2749 unit_gradient_type unit_gradient;
2750 ETT::set_gradient(gradient, v, unit_gradient);
2751 this->gradients[offset] =
2752 this->cell_type <=
2755 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2757 this
2758 ->inverse_jacobian_ptr[this->cell_type <=
2760 GeometryType::affine ?
2761 0 :
2762 offset]
2763 .transpose(),
2764 unit_gradient);
2765 }
2766 }
2767 }
2768}
2769
2770
2771
2772template <int n_components_, int dim, int spacedim, typename Number>
2773template <std::size_t stride_view>
2774inline void
2777 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2778{
2779 // slow path with FEValues
2780 Assert(this->fe_values.get() != nullptr,
2781 ExcMessage(
2782 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2783
2784 const std::size_t n_points = this->fe_values->get_quadrature().size();
2785
2786 if (evaluation_flags & EvaluationFlags::values)
2787 {
2788 this->values.resize(this->n_q_points);
2789 std::fill(this->values.begin(), this->values.end(), value_type());
2790 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2791 {
2792 const ScalarNumber value = solution_values[i];
2793 for (unsigned int d = 0; d < n_components; ++d)
2794 if (this->nonzero_shape_function_component[i][d] &&
2795 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2796 for (unsigned int qb = 0, q = 0; q < n_points;
2797 ++qb, q += n_lanes_user_interface)
2798 for (unsigned int v = 0;
2799 v < n_lanes_user_interface && q + v < n_points;
2800 ++v)
2801 ETT::access(this->values[qb],
2802 v,
2803 d,
2804 this->fe_values->shape_value(i, q + v) * value);
2805 else if (this->nonzero_shape_function_component[i][d])
2806 for (unsigned int qb = 0, q = 0; q < n_points;
2807 ++qb, q += n_lanes_user_interface)
2808 for (unsigned int v = 0;
2809 v < n_lanes_user_interface && q + v < n_points;
2810 ++v)
2811 ETT::access(this->values[qb],
2812 v,
2813 d,
2814 this->fe_values->shape_value_component(i,
2815 q + v,
2816 d) *
2817 value);
2818 }
2819 }
2820
2821 if (evaluation_flags & EvaluationFlags::gradients)
2822 {
2823 this->gradients.resize(this->n_q_points);
2824 std::fill(this->gradients.begin(),
2825 this->gradients.end(),
2826 gradient_type());
2827 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2828 {
2829 const ScalarNumber value = solution_values[i];
2830 for (unsigned int d = 0; d < n_components; ++d)
2831 if (this->nonzero_shape_function_component[i][d] &&
2832 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2833 for (unsigned int qb = 0, q = 0; q < n_points;
2834 ++qb, q += n_lanes_user_interface)
2835 for (unsigned int v = 0;
2836 v < n_lanes_user_interface && q + v < n_points;
2837 ++v)
2838 ETT::access(this->gradients[qb],
2839 v,
2840 d,
2841 this->fe_values->shape_grad(i, q + v) * value);
2842 else if (this->nonzero_shape_function_component[i][d])
2843 for (unsigned int qb = 0, q = 0; q < n_points;
2844 ++qb, q += n_lanes_user_interface)
2845 for (unsigned int v = 0;
2846 v < n_lanes_user_interface && q + v < n_points;
2847 ++v)
2848 ETT::access(
2849 this->gradients[qb],
2850 v,
2851 d,
2852 this->fe_values->shape_grad_component(i, q + v, d) * value);
2853 }
2854 }
2855}
2856
2857
2858
2859template <int n_components_, int dim, int spacedim, typename Number>
2860template <bool is_linear>
2861inline void
2863 const EvaluationFlags::EvaluationFlags &integration_flags,
2864 const unsigned int n_shapes,
2865 const unsigned int qb,
2866 const vectorized_value_type value,
2868 vectorized_value_type *solution_values_vectorized_linear)
2869{
2870 if (integration_flags & EvaluationFlags::gradients)
2872 is_linear,
2873 dim,
2875 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2876 n_shapes,
2877 &value,
2878 gradient,
2879 is_linear ?
2880 solution_values_vectorized_linear :
2881 this->solution_renumbered_vectorized.data(),
2882 this->unit_point_ptr[qb],
2883 qb != 0);
2884 else
2886 dim,
2889 this->shapes.data() + qb * n_shapes,
2890 n_shapes,
2891 value,
2892 is_linear ? solution_values_vectorized_linear :
2893 this->solution_renumbered_vectorized.data(),
2894 this->unit_point_ptr[qb],
2895 qb != 0);
2896}
2897
2898
2899
2900template <int n_components_, int dim, int spacedim, typename Number>
2901template <bool is_linear, std::size_t stride_view>
2902inline void
2904 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2905 vectorized_value_type *solution_values_vectorized_linear,
2906 const bool sum_into_values)
2907{
2908 if (!sum_into_values && this->fe->n_components() > n_components)
2909 for (unsigned int i = 0; i < solution_values.size(); ++i)
2910 solution_values[i] = 0;
2911
2912 const unsigned int dofs_per_comp =
2913 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2914
2915 for (unsigned int comp = 0; comp < n_components; ++comp)
2916 {
2917 const std::size_t offset =
2918 (this->component_in_base_element + comp) * dofs_per_comp;
2919
2920 if (is_linear || this->renumber.empty())
2921 {
2922 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2923 if (sum_into_values)
2924 solution_values[i + offset] +=
2925 ETT::sum_value(comp,
2926 is_linear ?
2927 *(solution_values_vectorized_linear + i) :
2928 this->solution_renumbered_vectorized[i]);
2929 else
2930 solution_values[i + offset] =
2931 ETT::sum_value(comp,
2932 is_linear ?
2933 *(solution_values_vectorized_linear + i) :
2934 this->solution_renumbered_vectorized[i]);
2935 }
2936 else
2937 {
2938 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2939 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2940 if (sum_into_values)
2941 solution_values[renumber_ptr[i]] +=
2942 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2943 else
2944 solution_values[renumber_ptr[i]] =
2945 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2946 }
2947 }
2948}
2949
2950
2951
2952template <int n_components_, int dim, int spacedim, typename Number>
2953template <bool do_JxW, bool is_linear, std::size_t stride_view>
2954inline void
2956 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2957 const EvaluationFlags::EvaluationFlags &integration_flags,
2958 const bool sum_into_values)
2959{
2960 // zero out lanes of incomplete last quadrature point batch
2961 if constexpr (stride == 1)
2962 if (const unsigned int n_filled_lanes =
2963 this->n_q_points_scalar & (n_lanes_internal - 1);
2964 n_filled_lanes > 0)
2965 {
2966 if (integration_flags & EvaluationFlags::values)
2967 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2968 ETT::set_zero_value(this->values.back(), v);
2969 if (integration_flags & EvaluationFlags::gradients)
2970 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2971 ETT::set_zero_gradient(this->gradients.back(), v);
2972 }
2973
2974 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2975 solution_values_vectorized_linear = {};
2976
2977 // loop over quadrature batches qb
2978 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2979
2980 const bool cartesian_cell =
2982 const bool affine_cell =
2984 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2985 {
2986 vectorized_value_type value = {};
2988
2989 if (integration_flags & EvaluationFlags::values)
2990 for (unsigned int v = 0, offset = qb * stride;
2991 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2992 ++v, ++offset)
2993 ETT::get_value(value,
2994 v,
2995 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2996 this->values[offset]);
2997
2998 if (integration_flags & EvaluationFlags::gradients)
2999 for (unsigned int v = 0, offset = qb * stride;
3000 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3001 ++v, ++offset)
3002 {
3003 const gradient_type grad_w =
3004 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3005 this->gradients[offset];
3006 ETT::get_gradient(
3007 gradient,
3008 v,
3009 cartesian_cell ?
3010 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3011 grad_w) :
3013 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3014 grad_w));
3015 }
3016
3017 compute_integrate_fast<is_linear>(
3018 integration_flags,
3019 n_shapes,
3020 qb,
3021 value,
3022 gradient,
3023 solution_values_vectorized_linear.data());
3024 }
3025
3026 // add between the lanes and write into the result
3027 finish_integrate_fast<is_linear>(solution_values,
3028 solution_values_vectorized_linear.data(),
3029 sum_into_values);
3030}
3031
3032
3033
3034template <int n_components_, int dim, int spacedim, typename Number>
3035template <bool do_JxW, std::size_t stride_view>
3036inline void
3038 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3039 const EvaluationFlags::EvaluationFlags &integration_flags,
3040 const bool sum_into_values)
3041{
3042 // slow path with FEValues
3043 Assert(this->fe_values.get() != nullptr,
3044 ExcMessage(
3045 "Not initialized. Please call FEPointEvaluation::reinit()!"));
3046 if (!sum_into_values)
3047 for (unsigned int i = 0; i < solution_values.size(); ++i)
3048 solution_values[i] = 0;
3049
3050 const std::size_t n_points = this->fe_values->get_quadrature().size();
3051
3052 if (integration_flags & EvaluationFlags::values)
3053 {
3054 AssertIndexRange(this->n_q_points, this->values.size() + 1);
3055 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3056 {
3057 for (unsigned int d = 0; d < n_components; ++d)
3058 if (this->nonzero_shape_function_component[i][d] &&
3059 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3060 for (unsigned int qb = 0, q = 0; q < n_points;
3061 ++qb, q += n_lanes_user_interface)
3062 for (unsigned int v = 0;
3063 v < n_lanes_user_interface && q + v < n_points;
3064 ++v)
3065 solution_values[i] +=
3066 this->fe_values->shape_value(i, q + v) *
3067 ETT::access(this->values[qb], v, d) *
3068 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3069 else if (this->nonzero_shape_function_component[i][d])
3070 for (unsigned int qb = 0, q = 0; q < n_points;
3071 ++qb, q += n_lanes_user_interface)
3072 for (unsigned int v = 0;
3073 v < n_lanes_user_interface && q + v < n_points;
3074 ++v)
3075 solution_values[i] +=
3076 this->fe_values->shape_value_component(i, q + v, d) *
3077 ETT::access(this->values[qb], v, d) *
3078 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3079 }
3080 }
3081
3082 if (integration_flags & EvaluationFlags::gradients)
3083 {
3084 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
3085 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
3086 {
3087 for (unsigned int d = 0; d < n_components; ++d)
3088 if (this->nonzero_shape_function_component[i][d] &&
3089 (this->fe->is_primitive(i) || this->fe->is_primitive()))
3090 for (unsigned int qb = 0, q = 0; q < n_points;
3091 ++qb, q += n_lanes_user_interface)
3092 for (unsigned int v = 0;
3093 v < n_lanes_user_interface && q + v < n_points;
3094 ++v)
3095 solution_values[i] +=
3096 this->fe_values->shape_grad(i, q + v) *
3097 ETT::access(this->gradients[qb], v, d) *
3098 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3099 else if (this->nonzero_shape_function_component[i][d])
3100 for (unsigned int qb = 0, q = 0; q < n_points;
3101 ++qb, q += n_lanes_user_interface)
3102 for (unsigned int v = 0;
3103 v < n_lanes_user_interface && q + v < n_points;
3104 ++v)
3105 solution_values[i] +=
3106 this->fe_values->shape_grad_component(i, q + v, d) *
3107 ETT::access(this->gradients[qb], v, d) *
3108 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
3109 }
3110 }
3111}
3112
3113
3114
3115template <int n_components_, int dim, int spacedim, typename Number>
3116template <bool do_JxW, std::size_t stride_view>
3117void
3119 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3120 const EvaluationFlags::EvaluationFlags &integration_flags,
3121 const bool sum_into_values)
3122{
3123 if (this->must_reinitialize_pointers)
3124 internal_reinit_single_cell_state_mapping_info();
3125
3126 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3127
3128 if (this->n_q_points == 0 || // no evaluation points provided
3129 !((integration_flags & EvaluationFlags::values) ||
3130 (integration_flags &
3131 EvaluationFlags::gradients))) // no integration flags
3132 {
3133 if (!sum_into_values)
3134 for (unsigned int i = 0; i < solution_values.size(); ++i)
3135 solution_values[i] = 0;
3136 return;
3137 }
3138
3139 Assert(
3140 !do_JxW || this->JxW_ptr != nullptr,
3141 ExcMessage(
3142 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3143
3144 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3145 if (this->fast_path)
3146 {
3147 if (this->use_linear_path)
3148 integrate_fast<do_JxW, true>(solution_values,
3149 integration_flags,
3150 sum_into_values);
3151 else
3152 integrate_fast<do_JxW, false>(solution_values,
3153 integration_flags,
3154 sum_into_values);
3155 }
3156 else
3157 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3158}
3159
3160
3161
3162template <int n_components_, int dim, int spacedim, typename Number>
3165 const unsigned int point_index) const
3166{
3167 AssertIndexRange(point_index, this->n_q_points);
3168 Assert(this->normal_ptr != nullptr,
3170 ExcFEPointEvaluationAccessToUninitializedMappingField(
3171 "update_normal_vectors"));
3172 return this->normal_ptr[point_index];
3173}
3174
3175
3176
3177template <int n_components_, int dim, int spacedim, typename Number>
3178inline const typename FEPointEvaluation<n_components_,
3179 dim,
3180 spacedim,
3181 Number>::value_type
3183 const unsigned int point_index) const
3184{
3185 AssertIndexRange(point_index, this->gradients.size());
3186 return this->gradients[point_index] * normal_vector(point_index);
3187}
3188
3189
3190
3191template <int n_components_, int dim, int spacedim, typename Number>
3192inline void
3195 const unsigned int point_index)
3196{
3197 AssertIndexRange(point_index, this->gradients.size());
3198 if constexpr (n_components == 1)
3199 this->gradients[point_index] = value * normal_vector(point_index);
3200 else
3201 this->gradients[point_index] =
3202 outer_product(value, normal_vector(point_index));
3203}
3204
3205
3206
3207template <int n_components_, int dim, int spacedim, typename Number>
3212 const bool is_interior,
3213 const unsigned int first_selected_component)
3214 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3215 mapping_info,
3216 fe,
3217 first_selected_component,
3218 is_interior)
3219{}
3220
3221
3222
3223template <int n_components_, int dim, int spacedim, typename Number>
3224inline void
3226 const unsigned int cell_index,
3227 const unsigned int face_number)
3228{
3229 this->current_cell_index = cell_index;
3230 this->current_face_number = face_number;
3231 this->must_reinitialize_pointers = false;
3232
3233 if (this->use_linear_path)
3234 this->template do_reinit<true, true>();
3235 else
3236 this->template do_reinit<true, false>();
3237}
3238
3239
3240
3241template <int n_components_, int dim, int spacedim, typename Number>
3242inline void
3244 const unsigned int face_index)
3245{
3246 this->current_cell_index = face_index;
3247 this->current_face_number =
3248 this->mapping_info->get_face_number(face_index, this->is_interior);
3249 this->must_reinitialize_pointers = false;
3250
3251 if (this->use_linear_path)
3252 this->template do_reinit<true, true>();
3253 else
3254 this->template do_reinit<true, false>();
3255}
3256
3257
3258
3259template <int n_components_, int dim, int spacedim, typename Number>
3260template <std::size_t stride_view>
3261void
3264 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3265{
3266 Assert(!this->must_reinitialize_pointers,
3267 ExcMessage("Object has not been reinitialized!"));
3268
3269 if (this->n_q_points == 0)
3270 return;
3271
3272 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3273
3274 if (!((evaluation_flags & EvaluationFlags::values) ||
3275 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3276 return;
3277
3278 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3279
3280 if (this->use_linear_path)
3281 do_evaluate<true>(solution_values, evaluation_flags);
3282 else
3283 do_evaluate<false>(solution_values, evaluation_flags);
3284}
3285
3286
3287
3288template <int n_components_, int dim, int spacedim, typename Number>
3289void
3291 const ArrayView<const ScalarNumber> &solution_values,
3292 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3293{
3294 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3295 solution_values.size()),
3296 evaluation_flags);
3297}
3298
3299
3300
3301template <int n_components_, int dim, int spacedim, typename Number>
3302template <bool is_linear, std::size_t stride_view>
3303void
3306 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3307{
3308 const unsigned int dofs_per_comp =
3309 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3310
3311 const ScalarNumber *input;
3312 if (stride_view == 1 && this->component_in_base_element == 0 &&
3313 (is_linear || this->renumber.empty()))
3314 input = solution_values.data();
3315 else
3316 {
3317 for (unsigned int comp = 0; comp < n_components; ++comp)
3318 {
3319 const std::size_t offset =
3320 (this->component_in_base_element + comp) * dofs_per_comp;
3321
3322 if (is_linear || this->renumber.empty())
3323 {
3324 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3325 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3326 solution_values[i + offset];
3327 }
3328 else
3329 {
3330 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3331 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3332 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3333 solution_values[renumber_ptr[i]];
3334 }
3335 }
3336 input = this->scratch_data_scalar.data();
3337 }
3338
3339 ScalarNumber *output =
3340 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3341
3342 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3343 template interpolate<true, false>(n_components,
3344 evaluation_flags,
3345 this->shape_info,
3346 input,
3347 output,
3348 this->current_face_number);
3349
3350 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3351}
3352
3353
3354
3355template <int n_components_, int dim, int spacedim, typename Number>
3356template <std::size_t stride_view>
3357void
3359 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3360 const EvaluationFlags::EvaluationFlags &integration_flags,
3361 const bool sum_into_values)
3362{
3363 Assert(!this->must_reinitialize_pointers,
3364 ExcMessage("Object has not been reinitialized!"));
3365
3366 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3367
3368 if (this->n_q_points == 0 || // no evaluation points provided
3369 !((integration_flags & EvaluationFlags::values) ||
3370 (integration_flags &
3371 EvaluationFlags::gradients))) // no integration flags
3372 {
3373 if (!sum_into_values)
3374 for (unsigned int i = 0; i < solution_values.size(); ++i)
3375 solution_values[i] = 0;
3376 return;
3377 }
3378
3379 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3380
3381 if (this->use_linear_path)
3382 do_integrate<true, true>(solution_values,
3383 integration_flags,
3384 sum_into_values);
3385 else
3386 do_integrate<true, false>(solution_values,
3387 integration_flags,
3388 sum_into_values);
3389}
3390
3391
3392
3393template <int n_components_, int dim, int spacedim, typename Number>
3394void
3396 const ArrayView<ScalarNumber> &solution_values,
3397 const EvaluationFlags::EvaluationFlags &integration_flags,
3398 const bool sum_into_values)
3399{
3400 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3401 solution_values.size()),
3402 integration_flags,
3403 sum_into_values);
3404}
3405
3406
3407
3408template <int n_components_, int dim, int spacedim, typename Number>
3409template <std::size_t stride_view>
3410void
3412 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3413 const EvaluationFlags::EvaluationFlags &integration_flags,
3414 const bool sum_into_values)
3415{
3416 Assert(!this->must_reinitialize_pointers,
3417 ExcMessage("Object has not been reinitialized!"));
3418
3419 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3420
3421 if (this->n_q_points == 0 || // no evaluation points provided
3422 !((integration_flags & EvaluationFlags::values) ||
3423 (integration_flags &
3424 EvaluationFlags::gradients))) // no integration flags
3425 {
3426 if (!sum_into_values)
3427 for (unsigned int i = 0; i < solution_values.size(); ++i)
3428 solution_values[i] = 0;
3429 return;
3430 }
3431
3432 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3433
3434 if (this->use_linear_path)
3435 do_integrate<false, true>(solution_values,
3436 integration_flags,
3437 sum_into_values);
3438 else
3439 do_integrate<false, false>(solution_values,
3440 integration_flags,
3441 sum_into_values);
3442}
3443
3444
3445
3446template <int n_components_, int dim, int spacedim, typename Number>
3447void
3449 const ArrayView<ScalarNumber> &solution_values,
3450 const EvaluationFlags::EvaluationFlags &integration_flags,
3451 const bool sum_into_values)
3452{
3453 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3454 solution_values.size()),
3455 integration_flags,
3456 sum_into_values);
3457}
3458
3459
3460
3461template <int n_components_, int dim, int spacedim, typename Number>
3462template <bool do_JxW, bool is_linear, std::size_t stride_view>
3463void
3465 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3466 const EvaluationFlags::EvaluationFlags &integration_flags,
3467 const bool sum_into_values)
3468{
3469 if (!sum_into_values && this->fe->n_components() > n_components)
3470 for (unsigned int i = 0; i < solution_values.size(); ++i)
3471 solution_values[i] = 0;
3472
3473 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3474 integration_flags,
3475 false);
3476
3477 ScalarNumber *input = this->scratch_data_scalar.begin();
3478
3479 if (stride_view == 1 && this->component_in_base_element == 0 &&
3480 (is_linear || this->renumber.empty()))
3481 {
3482 if (sum_into_values)
3483 internal::
3484 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3485 template interpolate<false, true>(n_components,
3486 integration_flags,
3487 this->shape_info,
3488 input,
3489 solution_values.data(),
3490 this->current_face_number);
3491 else
3492 internal::
3493 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3494 template interpolate<false, false>(n_components,
3495 integration_flags,
3496 this->shape_info,
3497 input,
3498 solution_values.data(),
3499 this->current_face_number);
3500 }
3501 else
3502 {
3503 const unsigned int dofs_per_comp_face =
3504 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3505
3506 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3507 ScalarNumber *output = input + size_input;
3508
3509 internal::
3510 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3511 template interpolate<false, false>(n_components,
3512 integration_flags,
3513 this->shape_info,
3514 input,
3515 output,
3516 this->current_face_number);
3517
3518 const unsigned int dofs_per_comp =
3519 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3520
3521 for (unsigned int comp = 0; comp < n_components; ++comp)
3522 {
3523 const std::size_t offset =
3524 (this->component_in_base_element + comp) * dofs_per_comp;
3525
3526 if (is_linear || this->renumber.empty())
3527 {
3528 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3529 if (sum_into_values)
3530 solution_values[i + offset] +=
3531 output[i + comp * dofs_per_comp];
3532 else
3533 solution_values[i + offset] =
3534 output[i + comp * dofs_per_comp];
3535 }
3536 else
3537 {
3538 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3539 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3540 if (sum_into_values)
3541 solution_values[renumber_ptr[i]] +=
3542 output[i + comp * dofs_per_comp];
3543 else
3544 solution_values[renumber_ptr[i]] =
3545 output[i + comp * dofs_per_comp];
3546 }
3547 }
3548 }
3549}
3550
3551
3552
3553template <int n_components_, int dim, int spacedim, typename Number>
3554template <int stride_face_dof>
3555void
3557 const ScalarNumber *face_dof_values,
3558 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3559{
3560 if (this->use_linear_path)
3561 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3562 evaluation_flags);
3563 else
3564 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3565 evaluation_flags);
3566}
3567
3568
3569
3570template <int n_components_, int dim, int spacedim, typename Number>
3571template <bool is_linear, int stride_face_dof>
3572inline void
3574 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3575 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3576{
3577 const scalar_value_type *face_dof_values_ptr;
3578 if constexpr (n_components == 1)
3579 face_dof_values_ptr = face_dof_values;
3580 else
3581 {
3582 const unsigned int dofs_per_comp_face =
3583 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3584 for (unsigned int comp = 0; comp < n_components; ++comp)
3585 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3586 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3587 stride_face_dof],
3588 comp,
3589 this->solution_renumbered[i]);
3590
3591 face_dof_values_ptr = this->solution_renumbered.data();
3592 }
3593
3594 constexpr int stride_face_dof_actual =
3595 n_components == 1 ? stride_face_dof : 1;
3596
3597 // loop over quadrature batches qb
3598 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3599
3600 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3601 {
3604
3605 if (evaluation_flags & EvaluationFlags::gradients)
3606 {
3607 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3608 is_linear ?
3610 dim - 1,
3613 2,
3614 stride_face_dof_actual>(face_dof_values_ptr,
3615 this->unit_point_faces_ptr[qb]) :
3617 dim - 1,
3620 2,
3621 false,
3622 stride_face_dof_actual>(this->shapes_faces.data() +
3623 qb * n_shapes,
3624 n_shapes,
3625 face_dof_values_ptr);
3626
3627 value = interpolated_value[dim - 1];
3628 // reorder derivative from tangential/normal derivatives into tensor
3629 // in physical coordinates
3630 if (this->current_face_number / 2 == 0)
3631 {
3632 gradient[0] = interpolated_value[dim];
3633 if (dim > 1)
3634 gradient[1] = interpolated_value[0];
3635 if (dim > 2)
3636 gradient[2] = interpolated_value[1];
3637 }
3638 else if (this->current_face_number / 2 == 1)
3639 {
3640 if (dim > 1)
3641 gradient[1] = interpolated_value[dim];
3642 if (dim == 3)
3643 {
3644 gradient[0] = interpolated_value[1];
3645 gradient[2] = interpolated_value[0];
3646 }
3647 else if (dim == 2)
3648 gradient[0] = interpolated_value[0];
3649 else
3651 }
3652 else if (this->current_face_number / 2 == 2)
3653 {
3654 if (dim > 2)
3655 {
3656 gradient[0] = interpolated_value[0];
3657 gradient[1] = interpolated_value[1];
3658 gradient[2] = interpolated_value[dim];
3659 }
3660 else
3662 }
3663 else
3665 }
3666 else
3667 {
3668 value = is_linear ?
3670 dim - 1,
3673 stride_face_dof_actual>(face_dof_values_ptr,
3674 this->unit_point_faces_ptr[qb]) :
3676 dim - 1,
3679 false,
3680 stride_face_dof_actual>(this->shapes_faces.data() +
3681 qb * n_shapes,
3682 n_shapes,
3683 face_dof_values_ptr);
3684 }
3685
3686 if (evaluation_flags & EvaluationFlags::values)
3687 {
3688 for (unsigned int v = 0, offset = qb * stride;
3689 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3690 ++v, ++offset)
3691 ETT::set_value(value, v, this->values[offset]);
3692 }
3693 if (evaluation_flags & EvaluationFlags::gradients)
3694 {
3695 Assert(this->update_flags & update_gradients ||
3696 this->update_flags & update_inverse_jacobians,
3698
3699 for (unsigned int v = 0, offset = qb * stride;
3700 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3701 ++v, ++offset)
3702 {
3703 unit_gradient_type unit_gradient;
3704 ETT::set_gradient(gradient, v, unit_gradient);
3705 this->gradients[offset] =
3706 this->cell_type <=
3709 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3711 this
3712 ->inverse_jacobian_ptr[this->cell_type <=
3714 GeometryType::affine ?
3715 0 :
3716 offset]
3717 .transpose(),
3718 unit_gradient);
3719 }
3720 }
3721 }
3722}
3723
3724
3725
3726template <int n_components_, int dim, int spacedim, typename Number>
3727template <int stride_face_dof>
3728void
3730 ScalarNumber *face_dof_values,
3731 const EvaluationFlags::EvaluationFlags &integration_flags,
3732 const bool sum_into_values)
3733{
3734 if (this->use_linear_path)
3735 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3736 integration_flags,
3737 sum_into_values);
3738 else
3739 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3740 integration_flags,
3741 sum_into_values);
3742}
3743
3744
3745
3746template <int n_components_, int dim, int spacedim, typename Number>
3747template <bool do_JxW, bool is_linear, int stride_face_dof>
3748inline void
3751 ScalarNumber *face_dof_values,
3752 const EvaluationFlags::EvaluationFlags &integration_flags,
3753 const bool sum_into_values)
3754{
3755 // zero out lanes of incomplete last quadrature point batch
3756 if constexpr (stride == 1)
3757 if (const unsigned int n_filled_lanes =
3758 this->n_q_points_scalar & (n_lanes_internal - 1);
3759 n_filled_lanes > 0)
3760 {
3761 if (integration_flags & EvaluationFlags::values)
3762 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3763 ETT::set_zero_value(this->values.back(), v);
3764 if (integration_flags & EvaluationFlags::gradients)
3765 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3766 ETT::set_zero_gradient(this->gradients.back(), v);
3767 }
3768
3769 std::array<vectorized_value_type,
3770 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3771 solution_values_vectorized_linear = {};
3772
3773 // loop over quadrature batches qb
3774 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3775
3776 const bool cartesian_cell =
3778 const bool affine_cell =
3780 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3781 {
3782 vectorized_value_type value = {};
3784
3785 if (integration_flags & EvaluationFlags::values)
3786 for (unsigned int v = 0, offset = qb * stride;
3787 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3788 ++v, ++offset)
3789 ETT::get_value(value,
3790 v,
3791 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3792 this->values[offset]);
3793
3794 if (integration_flags & EvaluationFlags::gradients)
3795 for (unsigned int v = 0, offset = qb * stride;
3796 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3797 ++v, ++offset)
3798 {
3799 const auto grad_w =
3800 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3801 this->gradients[offset];
3802 ETT::get_gradient(
3803 gradient,
3804 v,
3805 cartesian_cell ?
3806 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3807 grad_w) :
3809 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3810 grad_w));
3811 }
3812
3813 if (integration_flags & EvaluationFlags::gradients)
3814 {
3815 std::array<vectorized_value_type, 2> value_face = {};
3816 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3817
3818 value_face[0] = value;
3819 // fill derivative in physical coordinates into tangential/normal
3820 // derivatives
3821 if (this->current_face_number / 2 == 0)
3822 {
3823 value_face[1] = gradient[0];
3824 if (dim > 1)
3825 gradient_in_face[0] = gradient[1];
3826 if (dim > 2)
3827 gradient_in_face[1] = gradient[2];
3828 }
3829 else if (this->current_face_number / 2 == 1)
3830 {
3831 if (dim > 1)
3832 value_face[1] = gradient[1];
3833 if (dim == 3)
3834 {
3835 gradient_in_face[0] = gradient[2];
3836 gradient_in_face[1] = gradient[0];
3837 }
3838 else if (dim == 2)
3839 gradient_in_face[0] = gradient[0];
3840 else
3842 }
3843 else if (this->current_face_number / 2 == 2)
3844 {
3845 if (dim > 2)
3846 {
3847 value_face[1] = gradient[2];
3848 gradient_in_face[0] = gradient[0];
3849 gradient_in_face[1] = gradient[1];
3850 }
3851 else
3853 }
3854 else
3856
3858 is_linear,
3859 dim - 1,
3862 2>(this->shapes_faces.data() + qb * n_shapes,
3863 n_shapes,
3864 value_face.data(),
3865 gradient_in_face,
3866 is_linear ? solution_values_vectorized_linear.data() :
3867 this->solution_renumbered_vectorized.data(),
3868 this->unit_point_faces_ptr[qb],
3869 qb != 0);
3870 }
3871 else
3873 dim - 1,
3876 this->shapes_faces.data() + qb * n_shapes,
3877 n_shapes,
3878 value,
3879 is_linear ? solution_values_vectorized_linear.data() :
3880 this->solution_renumbered_vectorized.data(),
3881 this->unit_point_faces_ptr[qb],
3882 qb != 0);
3883 }
3884
3885 const unsigned int dofs_per_comp_face =
3886 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3887
3888 for (unsigned int comp = 0; comp < n_components; ++comp)
3889 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3890 if (sum_into_values)
3891 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3892 stride_face_dof] +=
3893 ETT::sum_value(comp,
3894 is_linear ?
3895 *(solution_values_vectorized_linear.data() + i) :
3896 this->solution_renumbered_vectorized[i]);
3897 else
3898 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3899 ETT::sum_value(comp,
3900 is_linear ?
3901 *(solution_values_vectorized_linear.data() + i) :
3902 this->solution_renumbered_vectorized[i]);
3903}
3904
3905
3906
3907template <int n_components_, int dim, int spacedim, typename Number>
3910 const unsigned int point_index) const
3911{
3912 AssertIndexRange(point_index, this->n_q_points);
3913 Assert(this->normal_ptr != nullptr,
3915 ExcFEPointEvaluationAccessToUninitializedMappingField(
3916 "update_normal_vectors"));
3917 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3918 {
3920 for (unsigned int d = 0; d < dim; ++d)
3921 normal[d] =
3922 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3923 0);
3924
3925 return normal;
3926 }
3927 else
3928 {
3929 return this->normal_ptr[point_index];
3930 }
3931}
3932
3933
3934
3935template <int n_components_, int dim, int spacedim, typename Number>
3936inline const typename FEFacePointEvaluation<n_components_,
3937 dim,
3938 spacedim,
3939 Number>::value_type
3941 get_normal_derivative(const unsigned int point_index) const
3942{
3943 AssertIndexRange(point_index, this->gradients.size());
3944 return this->gradients[point_index] * normal_vector(point_index);
3945}
3946
3947
3948
3949template <int n_components_, int dim, int spacedim, typename Number>
3950inline void
3953 const unsigned int point_index)
3954{
3955 AssertIndexRange(point_index, this->gradients.size());
3956 if constexpr (n_components == 1)
3957 this->gradients[point_index] = value * normal_vector(point_index);
3958 else
3959 this->gradients[point_index] =
3960 outer_product(value, normal_vector(point_index));
3961}
3962
3964
3965#endif
value_type * data() const noexcept
Definition array_view.h:661
std::size_t size() const
Definition array_view.h:684
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
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)
static constexpr std::size_t stride
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
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 internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
FEFacePointEvaluation(NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const bool is_interior=true, const unsigned int first_selected_component=0)
typename ETT::unit_gradient_type unit_gradient_type
typename internal::FEPointEvaluation:: EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
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
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
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
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
const unsigned int n_q_points_scalar
const Point< dim, VectorizedArrayType > * unit_point_ptr
typename internal::FEPointEvaluation:: EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
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
SmartPointer< NonMatching::MappingInfo< dim, spacedim, Number > > mapping_info
const unsigned int n_q_batches
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)
SmartPointer< const FiniteElement< dim, spacedim > > fe
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
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)
FEPointEvaluationBase(NonMatching::MappingInfo< dim, spacedim, Number > &mapping_info, const FiniteElement< dim, spacedim > &fe, const unsigned int first_selected_component=0, const bool is_interior=true)
typename ETT::scalar_value_type scalar_value_type
static constexpr unsigned int dimension
typename ETT::vectorized_value_type vectorized_value_type
static constexpr unsigned int n_components
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
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
SmartPointer< const Mapping< dim, spacedim > > mapping
internal::MatrixFreeFunctions::ShapeInfo< ScalarNumber > shape_info
void submit_value(const value_type &value, const unsigned int point_index)
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
static constexpr std::size_t stride
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
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 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 ETT::real_gradient_type gradient_type
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
typename internal::FEPointEvaluation:: EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
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)
void internal_reinit_single_cell_state_mapping_info()
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
value_type * data() const noexcept
Definition array_view.h:846
std::size_t size() const
Definition array_view.h:867
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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)
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)
Definition exceptions.h:516
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.
#define DEAL_II_ASSERT_UNREACHABLE()
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
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)
static const unsigned int invalid_unsigned_int
Definition types.h:220
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)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
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)
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)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
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)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
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)
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
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)
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, Tensor< 1, dim, VectorizedArrayType > > vectorized_unit_gradient_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)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
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)
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
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 SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)