Reference documentation for deal.II version GIT relicensing-559-gd86056192f 2024-05-09 06:30:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
fe_point_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_fe_point_evaluation_h
16#define dealii_fe_point_evaluation_h
17
18#include <deal.II/base/config.h>
19
24#include <deal.II/base/tensor.h>
26
28#include <deal.II/fe/mapping.h>
29
35
37
39
40namespace internal
41{
43 {
46 std::string,
47 << "You are requesting information from an FEPointEvaluationBase "
48 << "object for which this kind of information has not been computed. "
49 << "What information these objects compute is determined by the update_* "
50 << "flags you pass to MappingInfo() in the Constructor. Here, "
51 << "the operation you are attempting requires the <" << arg1
52 << "> flag to be set, but it was apparently not specified "
53 << "upon initialization.");
54
59 template <int dim,
60 int spacedim,
61 int n_components,
62 typename Number,
63 typename Enable = void>
65 {
69 typename ::internal::VectorizedArrayTrait<
77 using real_gradient_type = std::conditional_t<
78 n_components == spacedim,
87
88 static void
89 read_value(const ScalarNumber vector_entry,
90 const unsigned int component,
91 scalar_value_type &result)
92 {
93 AssertIndexRange(component, n_components);
94 result[component] = vector_entry;
95 }
96
99 {
100 return result;
101 }
102
103 static scalar_value_type
105 {
106 scalar_value_type result_scalar = {};
107
108 for (unsigned int c = 0; c < n_components; ++c)
109 result_scalar[c] = result[c].sum();
110
111 return result_scalar;
112 }
113
114 static ScalarNumber
115 sum_value(const unsigned int component,
116 const vectorized_value_type &result)
117 {
118 AssertIndexRange(component, n_components);
119 return result[component].sum();
120 }
121
122 static void
124 const unsigned int vector_lane,
125 unit_gradient_type &result)
126 {
127 for (unsigned int i = 0; i < n_components; ++i)
128 for (unsigned int d = 0; d < dim; ++d)
129 result[i][d] =
131 value[d][i], vector_lane);
132 }
133
134 static void
136 const unsigned int vector_lane,
137 const unit_gradient_type &result)
138 {
139 for (unsigned int i = 0; i < n_components; ++i)
140 for (unsigned int d = 0; d < dim; ++d)
142 value[d][i], vector_lane) = result[i][d];
143 }
144
145 static void
147 const unsigned int vector_lane,
149 {
150 for (unsigned int i = 0; i < n_components; ++i)
151 for (unsigned int d = 0; d < dim; ++d)
153 value[d][i], vector_lane) = result[i][d];
154 }
155
156 static void
158 const unsigned int vector_lane)
159 {
160 for (unsigned int i = 0; i < n_components; ++i)
161 for (unsigned int d = 0; d < spacedim; ++d)
163 vector_lane) = 0.;
164 }
165
166 static void
168 const unsigned int vector_lane,
169 scalar_value_type &result)
170 {
171 for (unsigned int i = 0; i < n_components; ++i)
172 result[i] = value[i][vector_lane];
173 }
174
175 static void
177 const unsigned int,
178 vectorized_value_type &result)
179 {
180 result = value;
181 }
182
183 static void
185 const unsigned int vector_lane,
186 const scalar_value_type &result)
187 {
188 for (unsigned int i = 0; i < n_components; ++i)
189 value[i][vector_lane] = result[i];
190 }
191
192 static void
194 const unsigned int,
195 const vectorized_value_type &result)
196 {
197 value = result;
198 }
199
200 static void
201 set_zero_value(value_type &value, const unsigned int vector_lane)
202 {
203 for (unsigned int i = 0; i < n_components; ++i)
205 0.;
206 }
207
208 static void
210 const unsigned int vector_lane,
211 const unsigned int component,
212 const ScalarNumber &shape_value)
213 {
215 vector_lane) += shape_value;
216 }
217
218 static ScalarNumber
220 const unsigned int vector_lane,
221 const unsigned int component)
222 {
224 vector_lane);
225 }
226
227 static void
229 const unsigned int vector_lane,
230 const unsigned int component,
231 const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
232 {
233 for (unsigned int d = 0; d < spacedim; ++d)
235 vector_lane) +=
236 shape_gradient[d];
237 }
238
241 const unsigned int vector_lane,
242 const unsigned int component)
243 {
245 for (unsigned int d = 0; d < spacedim; ++d)
246 result[d] =
248 vector_lane);
249 return result;
250 }
251 };
252
253 template <int dim, int spacedim, typename Number>
254 struct EvaluatorTypeTraits<dim, spacedim, 1, Number>
255 {
259 typename ::internal::VectorizedArrayTrait<
261 using value_type = Number;
270
271 static void
272 read_value(const ScalarNumber vector_entry,
273 const unsigned int,
274 scalar_value_type &result)
275 {
276 result = vector_entry;
277 }
278
279 static scalar_value_type
281 {
282 return result;
283 }
284
285 static scalar_value_type
287 {
288 return result.sum();
289 }
290
291 static ScalarNumber
292 sum_value(const unsigned int, const vectorized_value_type &result)
293 {
294 return result.sum();
295 }
296
297 static void
299 const unsigned int vector_lane,
301 {
302 for (unsigned int d = 0; d < dim; ++d)
303 result[d] = value[d][vector_lane];
304 }
305
306 static void
308 const unsigned int,
310 {
311 result = value;
312 }
313
314 static void
316 const unsigned int vector_lane,
317 const scalar_unit_gradient_type &result)
318 {
319 for (unsigned int d = 0; d < dim; ++d)
320 value[d][vector_lane] = result[d];
321 }
322
323 static void
325 const unsigned int,
326 const vectorized_unit_gradient_type &result)
327 {
328 value = result;
329 }
330
331 static void
333 const unsigned int vector_lane)
334 {
335 for (unsigned int d = 0; d < spacedim; ++d)
337 0.;
338 }
339
340 static void
342 const unsigned int vector_lane,
343 scalar_value_type &result)
344 {
345 result = value[vector_lane];
346 }
347
348 static void
350 const unsigned int,
351 vectorized_value_type &result)
352 {
353 result = value;
354 }
355
356 static void
358 const unsigned int vector_lane,
359 const scalar_value_type &result)
360 {
361 value[vector_lane] = result;
362 }
363
364 static void
366 const unsigned int,
367 const vectorized_value_type &result)
368 {
369 value = result;
370 }
371
372 static void
373 set_zero_value(value_type &value, const unsigned int vector_lane)
374 {
376 }
377
378 static void
380 const unsigned int vector_lane,
381 const unsigned int,
382 const ScalarNumber &shape_value)
383 {
385 shape_value;
386 }
387
388 static ScalarNumber
390 const unsigned int vector_lane,
391 const unsigned int)
392 {
394 }
395
396 static void
398 const unsigned int vector_lane,
399 const unsigned int,
400 const Tensor<1, spacedim, ScalarNumber> &shape_gradient)
401 {
402 for (unsigned int d = 0; d < spacedim; ++d)
404 shape_gradient[d];
405 }
406
409 const unsigned int vector_lane,
410 const unsigned int)
411 {
413 for (unsigned int d = 0; d < spacedim; ++d)
414 result[d] =
416 return result;
417 }
418 };
419
420 template <int dim, typename Number>
422 dim,
423 dim,
424 Number,
425 std::enable_if_t<dim != 1>>
426 {
430 typename ::internal::VectorizedArrayTrait<
441
442 static void
443 read_value(const ScalarNumber vector_entry,
444 const unsigned int component,
445 scalar_value_type &result)
446 {
447 AssertIndexRange(component, dim);
448 result[component] = vector_entry;
449 }
450
451 static scalar_value_type
453 {
454 return result;
455 }
456
457 static scalar_value_type
459 {
460 scalar_value_type result_scalar = {};
461
462 for (unsigned int c = 0; c < dim; ++c)
463 result_scalar[c] = result[c].sum();
464
465 return result_scalar;
466 }
467
468 static ScalarNumber
469 sum_value(const unsigned int component,
470 const vectorized_value_type &result)
471 {
472 AssertIndexRange(component, dim);
473 return result[component].sum();
474 }
475
476 static void
478 const unsigned int vector_lane,
479 unit_gradient_type &result)
480 {
481 for (unsigned int i = 0; i < dim; ++i)
482 for (unsigned int d = 0; d < dim; ++d)
483 result[i][d] =
485 value[d][i], vector_lane);
486 }
487
488 static void
490 const unsigned int vector_lane,
491 const unit_gradient_type &result)
492 {
493 for (unsigned int i = 0; i < dim; ++i)
494 for (unsigned int d = 0; d < dim; ++d)
496 value[d][i], vector_lane) = result[i][d];
497 }
498
499 static void
501 const unsigned int vector_lane)
502 {
503 for (unsigned int i = 0; i < dim; ++i)
504 for (unsigned int d = 0; d < dim; ++d)
506 vector_lane) = 0.;
507 }
508
509 static void
511 const unsigned int vector_lane,
512 scalar_value_type &result)
513 {
514 for (unsigned int i = 0; i < dim; ++i)
515 result[i] = value[i][vector_lane];
516 }
517
518 static void
520 const unsigned int,
521 vectorized_value_type &result)
522 {
523 result = value;
524 }
525
526 static void
528 const unsigned int vector_lane,
529 const scalar_value_type &result)
530 {
531 for (unsigned int i = 0; i < dim; ++i)
532 value[i][vector_lane] = result[i];
533 }
534
535 static void
537 const unsigned int,
538 const vectorized_value_type &result)
539 {
540 value = result;
541 }
542
543 static void
544 set_zero_value(value_type &value, const unsigned int vector_lane)
545 {
546 for (unsigned int i = 0; i < dim; ++i)
548 0.;
549 }
550
551 static void
553 const unsigned int vector_lane,
554 const unsigned int component,
555 const ScalarNumber &shape_value)
556 {
558 vector_lane) += shape_value;
559 }
560
561 static ScalarNumber
563 const unsigned int vector_lane,
564 const unsigned int component)
565 {
567 vector_lane);
568 }
569
570 static void
572 const unsigned int vector_lane,
573 const unsigned int component,
574 const Tensor<1, dim, ScalarNumber> &shape_gradient)
575 {
576 for (unsigned int d = 0; d < dim; ++d)
578 vector_lane) +=
579 shape_gradient[d];
580 }
581
584 const unsigned int vector_lane,
585 const unsigned int component)
586 {
588 for (unsigned int d = 0; d < dim; ++d)
589 result[d] =
591 vector_lane);
592 return result;
593 }
594 };
595
596 template <int dim, int spacedim>
597 bool
599 const unsigned int base_element_number);
600
601 template <int dim, int spacedim>
602 bool
604
605 template <int dim, int spacedim>
606 std::vector<Polynomials::Polynomial<double>>
608 } // namespace FEPointEvaluation
609} // namespace internal
610
611
612
619template <int n_components_,
620 int dim,
621 int spacedim = dim,
622 typename Number = double>
624{
625public:
626 static constexpr unsigned int dimension = dim;
627 static constexpr unsigned int n_components = n_components_;
628
629 using number_type = Number;
630
633 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
635 using ETT = typename internal::FEPointEvaluation::
636 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
637 using value_type = typename ETT::value_type;
638 using scalar_value_type = typename ETT::scalar_value_type;
639 using vectorized_value_type = typename ETT::vectorized_value_type;
640 using gradient_type = typename ETT::real_gradient_type;
642 typename ETT::interface_vectorized_unit_gradient_type;
643
644protected:
666 const unsigned int first_selected_component = 0);
667
690 const unsigned int first_selected_component = 0,
691 const bool is_interior = true);
692
697
702
703public:
711 const value_type &
712 get_value(const unsigned int point_index) const;
713
722 void
723 submit_value(const value_type &value, const unsigned int point_index);
724
734 const gradient_type &
735 get_gradient(const unsigned int point_index) const;
736
745 void
746 submit_gradient(const gradient_type &, const unsigned int point_index);
747
754 jacobian(const unsigned int point_index) const;
755
763 inverse_jacobian(const unsigned int point_index) const;
764
770 Number
771 JxW(const unsigned int point_index) const;
772
779 DEAL_II_DEPRECATED_EARLY Point<spacedim, Number>
780 real_point(const unsigned int point_index) const;
781
787 quadrature_point(const unsigned int point_index) const;
788
794 unit_point(const unsigned int point_index) const;
795
805
813
817 unsigned int
819
820protected:
821 static constexpr std::size_t n_lanes_user_interface =
823 static constexpr std::size_t n_lanes_internal =
825 static constexpr std::size_t stride =
827
836 void
837 setup(const unsigned int first_selected_component);
838
844 template <bool is_face, bool is_linear>
845 void
847
851 const unsigned int n_q_batches;
852
856 const unsigned int n_q_points;
857
861 const unsigned int n_q_points_scalar;
862
867
872
877 std::vector<Polynomials::Polynomial<double>> poly;
878
883
888 std::vector<unsigned int> renumber;
889
897 std::vector<scalar_value_type> solution_renumbered;
898
906
911
915 std::vector<value_type> values;
916
920 std::vector<gradient_type> gradients;
921
927
933
939
945
951
957
963 const Number *JxW_ptr;
964
969
974 unsigned int dofs_per_component;
975
982
987
992
998 std::vector<std::array<bool, n_components>> nonzero_shape_function_component;
999
1004
1008 std::shared_ptr<FEValues<dim, spacedim>> fe_values;
1009
1013 std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
1015
1021
1026
1031
1036
1044
1050
1056
1057 const bool is_interior;
1058};
1059
1060
1061
1092template <int n_components_,
1093 int dim,
1094 int spacedim = dim,
1095 typename Number = double>
1097 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1098{
1099public:
1100 static constexpr unsigned int dimension = dim;
1101 static constexpr unsigned int n_components = n_components_;
1102
1103 using number_type = Number;
1104
1107 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1109 using ETT = typename internal::FEPointEvaluation::
1110 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1111 using value_type = typename ETT::value_type;
1112 using scalar_value_type = typename ETT::scalar_value_type;
1113 using vectorized_value_type = typename ETT::vectorized_value_type;
1114 using unit_gradient_type = typename ETT::unit_gradient_type;
1115 using gradient_type = typename ETT::real_gradient_type;
1117 typename ETT::interface_vectorized_unit_gradient_type;
1118
1140 const unsigned int first_selected_component = 0);
1141
1161 const unsigned int first_selected_component = 0);
1162
1174 void
1176 const ArrayView<const Point<dim>> &unit_points);
1177
1182 void
1183 reinit();
1184
1189 void
1190 reinit(const unsigned int cell_index);
1191
1192
1204 template <std::size_t stride_view>
1205 void
1206 evaluate(
1208 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1209
1221 void
1222 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1223 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1224
1247 template <std::size_t stride_view>
1248 void
1250 const EvaluationFlags::EvaluationFlags &integration_flags,
1251 const bool sum_into_values = false);
1252
1275 void
1276 integrate(const ArrayView<ScalarNumber> &solution_values,
1277 const EvaluationFlags::EvaluationFlags &integration_flags,
1278 const bool sum_into_values = false);
1279
1306 template <std::size_t stride_view>
1307 void
1309 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1310 const EvaluationFlags::EvaluationFlags &integration_flags,
1311 const bool sum_into_values = false);
1312
1339 void
1340 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1341 const EvaluationFlags::EvaluationFlags &integration_flags,
1342 const bool sum_into_values = false);
1343
1350 normal_vector(const unsigned int point_index) const;
1351
1352private:
1353 static constexpr std::size_t n_lanes_user_interface =
1355 static constexpr std::size_t n_lanes_internal =
1357 static constexpr std::size_t stride =
1359
1364 template <bool is_linear, std::size_t stride_view>
1365 void
1368
1373 template <bool is_linear, std::size_t stride_view>
1374 void
1377 const EvaluationFlags::EvaluationFlags &evaluation_flags,
1378 const unsigned int n_shapes,
1379 const unsigned int qb,
1380 vectorized_value_type &value,
1382
1386 template <bool is_linear, std::size_t stride_view>
1387 void
1390 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1391
1395 template <std::size_t stride_view>
1396 void
1399 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1400
1406 template <bool is_linear>
1407 void
1409 const EvaluationFlags::EvaluationFlags &integration_flags,
1410 const unsigned int n_shapes,
1411 const unsigned int qb,
1412 const vectorized_value_type value,
1414 vectorized_value_type *solution_values_vectorized_linear);
1415
1421 template <bool is_linear, std::size_t stride_view>
1422 void
1424 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1425 vectorized_value_type *solution_values_vectorized_linear,
1426 const bool sum_into_values);
1427
1431 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1432 void
1434 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1435 const EvaluationFlags::EvaluationFlags &integration_flags,
1436 const bool sum_into_values);
1437
1441 template <bool do_JxW, std::size_t stride_view>
1442 void
1444 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1445 const EvaluationFlags::EvaluationFlags &integration_flags,
1446 const bool sum_into_values);
1447
1451 template <bool do_JxW, std::size_t stride_view>
1452 void
1454 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1455 const EvaluationFlags::EvaluationFlags &integration_flags,
1456 const bool sum_into_values);
1457};
1458
1459
1460
1478template <int n_components_,
1479 int dim,
1480 int spacedim = dim,
1481 typename Number = double>
1483 : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
1484{
1485public:
1486 static constexpr unsigned int dimension = dim;
1487 static constexpr unsigned int n_components = n_components_;
1488
1489 using number_type = Number;
1490
1493 using VectorizedArrayType = typename ::internal::VectorizedArrayTrait<
1495 using ETT = typename internal::FEPointEvaluation::
1496 EvaluatorTypeTraits<dim, spacedim, n_components, Number>;
1497 using value_type = typename ETT::value_type;
1498 using scalar_value_type = typename ETT::scalar_value_type;
1499 using vectorized_value_type = typename ETT::vectorized_value_type;
1500 using unit_gradient_type = typename ETT::unit_gradient_type;
1501 using gradient_type = typename ETT::real_gradient_type;
1503 typename ETT::interface_vectorized_unit_gradient_type;
1504
1511 const bool is_interior = true,
1512 const unsigned int first_selected_component = 0);
1513
1518 void
1519 reinit(const unsigned int cell_index, const unsigned int face_number);
1520
1525 void
1526 reinit(const unsigned int face_index);
1527
1539 template <std::size_t stride_view>
1540 void
1541 evaluate(
1543 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1544
1556 void
1557 evaluate(const ArrayView<const ScalarNumber> &solution_values,
1558 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1559
1582 template <std::size_t stride_view>
1583 void
1585 const EvaluationFlags::EvaluationFlags &integration_flags,
1586 const bool sum_into_values = false);
1587
1610 void
1611 integrate(const ArrayView<ScalarNumber> &solution_values,
1612 const EvaluationFlags::EvaluationFlags &integration_flags,
1613 const bool sum_into_values = false);
1614
1637 template <std::size_t stride_view>
1638 void
1640 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1641 const EvaluationFlags::EvaluationFlags &integration_flags,
1642 const bool sum_into_values = false);
1643
1666 void
1667 test_and_sum(const ArrayView<ScalarNumber> &solution_values,
1668 const EvaluationFlags::EvaluationFlags &integration_flags,
1669 const bool sum_into_values = false);
1670
1677 template <int stride_face_dof = VectorizedArrayType::size()>
1678 void
1679 evaluate_in_face(const ScalarNumber *face_dof_values,
1680 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1681
1688 template <int stride_face_dof = VectorizedArrayType::size()>
1689 void
1690 integrate_in_face(ScalarNumber *face_dof_values,
1691 const EvaluationFlags::EvaluationFlags &integration_flags,
1692 const bool sum_into_values = false);
1693
1700 normal_vector(const unsigned int point_index) const;
1701
1702private:
1703 static constexpr std::size_t n_lanes_user_interface =
1705 static constexpr std::size_t n_lanes_internal =
1707 static constexpr std::size_t stride =
1709
1710 template <bool is_linear, std::size_t stride_view>
1711 void
1714 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1715
1716 template <bool do_JxW, bool is_linear, std::size_t stride_view>
1717 void
1719 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
1720 const EvaluationFlags::EvaluationFlags &integration_flags,
1721 const bool sum_into_values);
1722
1727 template <bool is_linear, int stride_face_dof>
1728 void
1729 do_evaluate_in_face(const ScalarNumber *face_dof_values,
1730 const EvaluationFlags::EvaluationFlags &evaluation_flags);
1731
1736 template <bool do_JxW, bool is_linear, int stride_face_dof>
1737 void
1739 ScalarNumber *face_dof_values,
1740 const EvaluationFlags::EvaluationFlags &integration_flags,
1741 const bool sum_into_values);
1742};
1743
1744
1745
1746// ----------------------- template and inline function ----------------------
1747
1748
1749template <int n_components_, int dim, int spacedim, typename Number>
1753 const UpdateFlags update_flags,
1754 const unsigned int first_selected_component)
1755 : n_q_batches(numbers::invalid_unsigned_int)
1756 , n_q_points(numbers::invalid_unsigned_int)
1757 , n_q_points_scalar(numbers::invalid_unsigned_int)
1758 , mapping(&mapping)
1759 , fe(&fe)
1760 , JxW_ptr(nullptr)
1761 , update_flags(update_flags)
1762 , mapping_info_on_the_fly(
1763 std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
1764 mapping,
1765 update_flags))
1766 , mapping_info(mapping_info_on_the_fly.get())
1767 , current_cell_index(numbers::invalid_unsigned_int)
1768 , current_face_number(numbers::invalid_unsigned_int)
1769 , must_reinitialize_pointers(false)
1770 , is_interior(true)
1771{
1772 setup(first_selected_component);
1773}
1774
1775
1776
1777template <int n_components_, int dim, int spacedim, typename Number>
1782 const unsigned int first_selected_component,
1783 const bool is_interior)
1784 : n_q_batches(numbers::invalid_unsigned_int)
1785 , n_q_points(numbers::invalid_unsigned_int)
1786 , n_q_points_scalar(numbers::invalid_unsigned_int)
1787 , mapping(&mapping_info.get_mapping())
1788 , fe(&fe)
1789 , JxW_ptr(nullptr)
1790 , update_flags(mapping_info.get_update_flags())
1791 , mapping_info(&mapping_info)
1792 , current_cell_index(numbers::invalid_unsigned_int)
1793 , current_face_number(numbers::invalid_unsigned_int)
1794 , must_reinitialize_pointers(true)
1795 , is_interior(is_interior)
1796{
1797 setup(first_selected_component);
1798}
1799
1800
1801
1802template <int n_components_, int dim, int spacedim, typename Number>
1806 : n_q_batches(other.n_q_batches)
1807 , n_q_points(other.n_q_points)
1808 , n_q_points_scalar(other.n_q_points_scalar)
1809 , mapping(other.mapping)
1810 , fe(other.fe)
1811 , poly(other.poly)
1812 , use_linear_path(other.use_linear_path)
1813 , renumber(other.renumber)
1814 , solution_renumbered(other.solution_renumbered)
1815 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1816 , values(other.values)
1817 , gradients(other.gradients)
1818 , dofs_per_component(other.dofs_per_component)
1819 , dofs_per_component_face(other.dofs_per_component_face)
1820 , component_in_base_element(other.component_in_base_element)
1821 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1822 , update_flags(other.update_flags)
1823 , fe_values(other.fe_values)
1824 , mapping_info_on_the_fly(
1825 other.mapping_info_on_the_fly ?
1827 *mapping,
1828 update_flags) :
1829 nullptr)
1830 , mapping_info(other.mapping_info)
1831 , current_cell_index(other.current_cell_index)
1832 , current_face_number(other.current_face_number)
1833 , fast_path(other.fast_path)
1834 , must_reinitialize_pointers(true)
1835 , is_interior(other.is_interior)
1836{}
1837
1838
1839
1840template <int n_components_, int dim, int spacedim, typename Number>
1843 : n_q_batches(other.n_q_batches)
1844 , n_q_points(other.n_q_points)
1845 , n_q_points_scalar(other.n_q_points_scalar)
1846 , mapping(other.mapping)
1847 , fe(other.fe)
1848 , poly(other.poly)
1849 , use_linear_path(other.use_linear_path)
1850 , renumber(other.renumber)
1851 , solution_renumbered(other.solution_renumbered)
1852 , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
1853 , values(other.values)
1854 , gradients(other.gradients)
1855 , dofs_per_component(other.dofs_per_component)
1856 , dofs_per_component_face(other.dofs_per_component_face)
1857 , component_in_base_element(other.component_in_base_element)
1858 , nonzero_shape_function_component(other.nonzero_shape_function_component)
1859 , update_flags(other.update_flags)
1860 , fe_values(other.fe_values)
1861 , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
1862 , mapping_info(other.mapping_info)
1863 , current_cell_index(other.current_cell_index)
1864 , current_face_number(other.current_face_number)
1865 , fast_path(other.fast_path)
1866 , must_reinitialize_pointers(other.must_reinitialize_pointers)
1867 , is_interior(other.is_interior)
1868{}
1869
1870
1871
1872template <int n_components_, int dim, int spacedim, typename Number>
1873void
1875 const unsigned int first_selected_component)
1876{
1877 AssertIndexRange(first_selected_component + n_components,
1878 fe->n_components() + 1);
1879
1880 shapes.reserve(100);
1881
1882 bool same_base_element = true;
1883 unsigned int base_element_number = 0;
1884 component_in_base_element = 0;
1885 unsigned int component = 0;
1886 for (; base_element_number < fe->n_base_elements(); ++base_element_number)
1887 if (component + fe->element_multiplicity(base_element_number) >
1888 first_selected_component)
1889 {
1890 if (first_selected_component + n_components >
1891 component + fe->element_multiplicity(base_element_number))
1892 same_base_element = false;
1893 component_in_base_element = first_selected_component - component;
1894 break;
1895 }
1896 else
1897 component += fe->element_multiplicity(base_element_number);
1898
1901 *fe, base_element_number) &&
1902 same_base_element)
1903 {
1904 shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
1905 renumber = shape_info.lexicographic_numbering;
1906 dofs_per_component = shape_info.dofs_per_component_on_cell;
1907 dofs_per_component_face = shape_info.dofs_per_component_on_face;
1909 fe->base_element(base_element_number));
1910
1911 bool is_lexicographic = true;
1912 for (unsigned int i = 0; i < renumber.size(); ++i)
1913 if (i != renumber[i])
1914 is_lexicographic = false;
1915
1916 if (is_lexicographic)
1917 renumber.clear();
1918
1919 use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
1920 poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
1921 poly[1].value(1.) == 1.) &&
1922 (fe->n_components() == n_components);
1923
1924 const unsigned int size_face = 3 * dofs_per_component_face * n_components;
1925 const unsigned int size_cell = dofs_per_component * n_components;
1926 scratch_data_scalar.resize(size_face + size_cell);
1927
1928 solution_renumbered.resize(dofs_per_component);
1929 solution_renumbered_vectorized.resize(dofs_per_component);
1930
1931 fast_path = true;
1932 }
1933 else
1934 {
1935 nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
1936 for (unsigned int d = 0; d < n_components; ++d)
1937 {
1938 const unsigned int component = first_selected_component + d;
1939 for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
1940 {
1941 const bool is_primitive =
1942 fe->is_primitive() || fe->is_primitive(i);
1943 if (is_primitive)
1944 nonzero_shape_function_component[i][d] =
1945 (component == fe->system_to_component_index(i).first);
1946 else
1947 nonzero_shape_function_component[i][d] =
1948 (fe->get_nonzero_components(i)[component] == true);
1949 }
1950 }
1951
1952 fast_path = false;
1953 }
1954}
1955
1956
1957
1958template <int n_components_, int dim, int spacedim, typename Number>
1959template <bool is_face, bool is_linear>
1960inline void
1962{
1963 const unsigned int geometry_index =
1964 mapping_info->template compute_geometry_index_offset<is_face>(
1965 current_cell_index, current_face_number);
1966
1967 cell_type = mapping_info->get_cell_type(geometry_index);
1968
1969 const_cast<unsigned int &>(n_q_points_scalar) =
1970 mapping_info->get_n_q_points_unvectorized(geometry_index);
1971
1972 // round up n_q_points_scalar / n_lanes_internal
1973 const_cast<unsigned int &>(n_q_batches) =
1974 (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
1975
1976 const unsigned int n_q_points_before = n_q_points;
1977
1978 const_cast<unsigned int &>(n_q_points) =
1979 (stride == 1) ? n_q_batches : n_q_points_scalar;
1980
1981 if (n_q_points != n_q_points_before)
1982 {
1983 if (update_flags & update_values)
1984 values.resize(n_q_points);
1985 if (update_flags & update_gradients)
1986 gradients.resize(n_q_points);
1987 }
1988
1989 if (n_q_points == 0)
1990 return;
1991
1992 // set unit point pointer
1993 const unsigned int unit_point_offset =
1994 mapping_info->compute_unit_point_index_offset(geometry_index);
1995
1996 if (is_face)
1997 unit_point_faces_ptr =
1998 mapping_info->get_unit_point_faces(unit_point_offset);
1999 else
2000 unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
2001
2002 // set data pointers
2003 const unsigned int data_offset =
2004 mapping_info->compute_data_index_offset(geometry_index);
2005 const unsigned int compressed_data_offset =
2006 mapping_info->compute_compressed_data_index_offset(geometry_index);
2007#ifdef DEBUG
2008 const UpdateFlags update_flags_mapping =
2009 mapping_info->get_update_flags_mapping();
2010 if (update_flags_mapping & UpdateFlags::update_quadrature_points)
2011 real_point_ptr = mapping_info->get_real_point(data_offset);
2012 if (update_flags_mapping & UpdateFlags::update_jacobians)
2013 jacobian_ptr =
2014 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2015 if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
2016 inverse_jacobian_ptr =
2017 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2018 if (update_flags_mapping & UpdateFlags::update_normal_vectors)
2019 normal_ptr = mapping_info->get_normal_vector(data_offset);
2020 if (update_flags_mapping & UpdateFlags::update_JxW_values)
2021 JxW_ptr = mapping_info->get_JxW(data_offset);
2022#else
2023 real_point_ptr = mapping_info->get_real_point(data_offset);
2024 jacobian_ptr =
2025 mapping_info->get_jacobian(compressed_data_offset, is_interior);
2026 inverse_jacobian_ptr =
2027 mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
2028 normal_ptr = mapping_info->get_normal_vector(data_offset);
2029 JxW_ptr = mapping_info->get_JxW(data_offset);
2030#endif
2031
2032 if (!is_linear && fast_path)
2033 {
2034 const std::size_t n_shapes = poly.size();
2035 if (is_face)
2036 shapes_faces.resize_fast(n_q_batches * n_shapes);
2037 else
2038 shapes.resize_fast(n_q_batches * n_shapes);
2039
2040 for (unsigned int qb = 0; qb < n_q_batches; ++qb)
2041 if (is_face)
2042 {
2043 if (dim > 1)
2044 {
2046 shapes_faces.data() + qb * n_shapes,
2047 poly,
2048 unit_point_faces_ptr[qb],
2049 update_flags & UpdateFlags::update_gradients ? 1 : 0);
2050 }
2051 }
2052 else
2053 {
2054 if (update_flags & UpdateFlags::update_gradients)
2055 {
2056 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2057 poly,
2058 unit_point_ptr[qb],
2059 1);
2060 }
2061 else if (qb + 1 < n_q_batches)
2062 {
2063 // Use function with reduced overhead to compute for two
2064 // points at once
2066 shapes.data() + qb * n_shapes,
2067 poly,
2068 unit_point_ptr[qb],
2069 unit_point_ptr[qb + 1]);
2070 ++qb;
2071 }
2072 else
2073 {
2074 internal::compute_values_of_array(shapes.data() + qb * n_shapes,
2075 poly,
2076 unit_point_ptr[qb],
2077 0);
2078 }
2079 }
2080 }
2081}
2082
2083
2084
2085template <int n_components_, int dim, int spacedim, typename Number>
2086inline const typename FEPointEvaluationBase<n_components_,
2087 dim,
2088 spacedim,
2089 Number>::value_type &
2091 const unsigned int point_index) const
2092{
2093 AssertIndexRange(point_index, values.size());
2094 return values[point_index];
2095}
2096
2097
2098
2099template <int n_components_, int dim, int spacedim, typename Number>
2100inline const typename FEPointEvaluationBase<n_components_,
2101 dim,
2102 spacedim,
2103 Number>::gradient_type &
2105 const unsigned int point_index) const
2106{
2107 AssertIndexRange(point_index, gradients.size());
2108 return gradients[point_index];
2109}
2110
2111
2112
2113template <int n_components_, int dim, int spacedim, typename Number>
2114inline void
2116 const value_type &value,
2117 const unsigned int point_index)
2118{
2119 AssertIndexRange(point_index, n_q_points);
2120 values[point_index] = value;
2121}
2122
2123
2124
2125template <int n_components_, int dim, int spacedim, typename Number>
2126inline void
2128 const gradient_type &gradient,
2129 const unsigned int point_index)
2130{
2131 AssertIndexRange(point_index, n_q_points);
2132 gradients[point_index] = gradient;
2133}
2134
2135
2136
2137template <int n_components_, int dim, int spacedim, typename Number>
2140 const unsigned int point_index) const
2141{
2142 AssertIndexRange(point_index, n_q_points);
2143 Assert(jacobian_ptr != nullptr,
2145 ExcFEPointEvaluationAccessToUninitializedMappingField(
2146 "update_jacobians"));
2147 return jacobian_ptr[cell_type <= ::internal::MatrixFreeFunctions::
2148 GeometryType::affine ?
2149 0 :
2150 point_index];
2151}
2152
2153
2154
2155template <int n_components_, int dim, int spacedim, typename Number>
2158 const unsigned int point_index) const
2159{
2160 AssertIndexRange(point_index, n_q_points);
2161 Assert(inverse_jacobian_ptr != nullptr,
2163 ExcFEPointEvaluationAccessToUninitializedMappingField(
2164 "update_inverse_jacobians"));
2165 return inverse_jacobian_ptr
2166 [cell_type <=
2168 0 :
2169 point_index];
2170}
2171
2172
2173
2174template <int n_components_, int dim, int spacedim, typename Number>
2175inline Number
2177 const unsigned int point_index) const
2178{
2179 AssertIndexRange(point_index, n_q_points);
2180 Assert(JxW_ptr != nullptr,
2182 ExcFEPointEvaluationAccessToUninitializedMappingField(
2183 "update_JxW_values"));
2184 return JxW_ptr[point_index];
2185}
2186
2187
2188
2189template <int n_components_, int dim, int spacedim, typename Number>
2192 const unsigned int point_index) const
2193{
2194 return quadrature_point(point_index);
2195 ;
2196}
2197
2198
2199
2200template <int n_components_, int dim, int spacedim, typename Number>
2203 const unsigned int point_index) const
2204{
2205 AssertIndexRange(point_index, n_q_points);
2206 Assert(real_point_ptr != nullptr,
2208 ExcFEPointEvaluationAccessToUninitializedMappingField(
2209 "update_quadrature_points"));
2210 return real_point_ptr[point_index];
2211}
2212
2213
2214
2215template <int n_components_, int dim, int spacedim, typename Number>
2216inline Point<dim, Number>
2218 const unsigned int point_index) const
2219{
2220 AssertIndexRange(point_index, n_q_points);
2221 Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
2222 Point<dim, Number> unit_point;
2223 for (unsigned int d = 0; d < dim; ++d)
2225 unit_point_ptr[point_index / stride][d], point_index % stride);
2226 return unit_point;
2227}
2228
2229
2230
2231template <int n_components_, int dim, int spacedim, typename Number>
2238
2239
2240
2241template <int n_components_, int dim, int spacedim, typename Number>
2245 const unsigned int first_selected_component)
2246 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2247 mapping_info,
2248 fe,
2249 first_selected_component)
2250{}
2251
2252
2253
2254template <int n_components_, int dim, int spacedim, typename Number>
2256 const Mapping<dim, spacedim> &mapping,
2258 const UpdateFlags update_flags,
2259 const unsigned int first_selected_component)
2260 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
2261 mapping,
2262 fe,
2263 update_flags,
2264 first_selected_component)
2265{}
2266
2267
2268
2269template <int n_components_, int dim, int spacedim, typename Number>
2270inline void
2272{
2273 this->current_cell_index = numbers::invalid_unsigned_int;
2274 this->current_face_number = numbers::invalid_unsigned_int;
2275
2276 if (this->use_linear_path)
2277 this->template do_reinit<false, true>();
2278 else
2279 this->template do_reinit<false, false>();
2280}
2281
2282
2283
2284template <int n_components_, int dim, int spacedim, typename Number>
2285inline void
2288 const ArrayView<const Point<dim>> &unit_points)
2289{
2290 // reinit is only allowed for mapping computation on the fly
2291 AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
2293
2294 this->mapping_info->reinit(cell, unit_points);
2295 this->must_reinitialize_pointers = false;
2296
2297 if (!this->fast_path)
2298 {
2299 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2300 *this->mapping,
2301 *this->fe,
2303 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2304 this->update_flags);
2305 this->fe_values->reinit(cell);
2306 }
2307
2308 if (this->use_linear_path)
2309 this->template do_reinit<false, true>();
2310 else
2311 this->template do_reinit<false, false>();
2312}
2313
2314
2315
2316template <int n_components_, int dim, int spacedim, typename Number>
2317inline void
2319 const unsigned int cell_index)
2320{
2321 this->current_cell_index = cell_index;
2322 this->current_face_number = numbers::invalid_unsigned_int;
2323 this->must_reinitialize_pointers = false;
2324
2325 if (this->use_linear_path)
2326 this->template do_reinit<false, true>();
2327 else
2328 this->template do_reinit<false, false>();
2329
2330 if (!this->fast_path)
2331 {
2332 std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
2333
2334 for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
2335 for (unsigned int d = 0; d < dim; ++d)
2336 unit_points[v][d] =
2337 this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
2338
2339 this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
2340 *this->mapping,
2341 *this->fe,
2343 std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
2344 this->update_flags);
2345
2346 this->fe_values->reinit(
2347 this->mapping_info->get_cell_iterator(this->current_cell_index));
2348 }
2349}
2350
2351
2352
2353template <int n_components_, int dim, int spacedim, typename Number>
2354template <std::size_t stride_view>
2355void
2358 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2359{
2360 if (this->must_reinitialize_pointers)
2361 reinit();
2362
2363 if (this->n_q_points == 0)
2364 return;
2365
2366 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
2367
2368 if (!((evaluation_flags & EvaluationFlags::values) ||
2369 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
2370 return;
2371
2372 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
2373 if (this->fast_path)
2374 {
2375 if (this->use_linear_path)
2376 evaluate_fast<true>(solution_values, evaluation_flags);
2377 else
2378 evaluate_fast<false>(solution_values, evaluation_flags);
2379 }
2380 else
2381 evaluate_slow(solution_values, evaluation_flags);
2382}
2383
2384
2385
2386template <int n_components_, int dim, int spacedim, typename Number>
2387void
2389 const ArrayView<const ScalarNumber> &solution_values,
2390 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2391{
2392 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
2393 solution_values.size()),
2394 evaluation_flags);
2395}
2396
2397
2398
2399template <int n_components_, int dim, int spacedim, typename Number>
2400template <std::size_t stride_view>
2401void
2403 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2404 const EvaluationFlags::EvaluationFlags &integration_flags,
2405 const bool sum_into_values)
2406{
2407 do_integrate<true>(solution_values, integration_flags, sum_into_values);
2408}
2409
2410
2411
2412template <int n_components_, int dim, int spacedim, typename Number>
2413void
2415 const ArrayView<ScalarNumber> &solution_values,
2416 const EvaluationFlags::EvaluationFlags &integration_flags,
2417 const bool sum_into_values)
2418{
2419 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2420 solution_values.size()),
2421 integration_flags,
2422 sum_into_values);
2423}
2424
2425
2426
2427template <int n_components_, int dim, int spacedim, typename Number>
2429 scalar_value_type
2431 const
2432{
2433 value_type return_value = {};
2434
2435 for (const auto point_index : this->quadrature_point_indices())
2436 return_value += values[point_index] * this->JxW(point_index);
2437
2438 return ETT::sum_value(return_value);
2439}
2440
2441
2442
2443template <int n_components_, int dim, int spacedim, typename Number>
2444unsigned int
2447{
2448 Assert(stride == 1,
2449 ExcMessage(
2450 "Calling this function only makes sense in fully vectorized mode."));
2451 if (q == n_q_batches - 1)
2452 {
2453 const unsigned int n_filled_lanes =
2454 n_q_points_scalar & (n_lanes_user_interface - 1);
2455 if (n_filled_lanes == 0)
2456 return n_lanes_user_interface;
2457 else
2458 return n_filled_lanes;
2459 }
2460 else
2461 return n_lanes_user_interface;
2462}
2463
2464
2465
2466template <int n_components_, int dim, int spacedim, typename Number>
2467template <std::size_t stride_view>
2468void
2470 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2471 const EvaluationFlags::EvaluationFlags &integration_flags,
2472 const bool sum_into_values)
2473{
2474 do_integrate<false>(solution_values, integration_flags, sum_into_values);
2475}
2476
2477
2478
2479template <int n_components_, int dim, int spacedim, typename Number>
2480void
2482 const ArrayView<ScalarNumber> &solution_values,
2483 const EvaluationFlags::EvaluationFlags &integration_flags,
2484 const bool sum_into_values)
2485{
2486 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
2487 solution_values.size()),
2488 integration_flags,
2489 sum_into_values);
2490}
2491
2492
2493
2494template <int n_components_, int dim, int spacedim, typename Number>
2495template <bool is_linear, std::size_t stride_view>
2496inline void
2499{
2500 const unsigned int dofs_per_comp =
2501 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2502
2503 for (unsigned int comp = 0; comp < n_components; ++comp)
2504 {
2505 const std::size_t offset =
2506 (this->component_in_base_element + comp) * dofs_per_comp;
2507
2508 if ((is_linear && n_components == 1) || this->renumber.empty())
2509 {
2510 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2511 ETT::read_value(solution_values[i + offset],
2512 comp,
2513 this->solution_renumbered[i]);
2514 }
2515 else
2516 {
2517 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2518 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2519 ETT::read_value(solution_values[renumber_ptr[i]],
2520 comp,
2521 this->solution_renumbered[i]);
2522 }
2523 }
2524}
2525
2526
2527
2528template <int n_components_, int dim, int spacedim, typename Number>
2529template <bool is_linear, std::size_t stride_view>
2530inline void
2533 const EvaluationFlags::EvaluationFlags &evaluation_flags,
2534 const unsigned int n_shapes,
2535 const unsigned int qb,
2536 vectorized_value_type &value,
2538{
2539 if (evaluation_flags & EvaluationFlags::gradients)
2540 {
2541 std::array<vectorized_value_type, dim + 1> result;
2542 if constexpr (is_linear)
2543 {
2544 if constexpr (n_components == 1)
2545 result =
2547 dim,
2550 1,
2551 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2552 else
2553 result =
2555 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2556 }
2557 else
2559 dim,
2562 1,
2563 false>(this->shapes.data() + qb * n_shapes,
2564 n_shapes,
2565 this->solution_renumbered.data());
2566 gradient[0] = result[0];
2567 if (dim > 1)
2568 gradient[1] = result[1];
2569 if (dim > 2)
2570 gradient[2] = result[2];
2571 value = result[dim];
2572 }
2573 else
2574 {
2575 if constexpr (is_linear)
2576 {
2577 if constexpr (n_components == 1)
2579 dim,
2582 stride_view>(solution_values.data(), this->unit_point_ptr[qb]);
2583 else
2585 this->solution_renumbered.data(), this->unit_point_ptr[qb]);
2586 }
2587 else
2588 value =
2592 false>(
2593 this->shapes.data() + qb * n_shapes,
2594 n_shapes,
2595 this->solution_renumbered.data());
2596 }
2597}
2598
2599
2600
2601template <int n_components_, int dim, int spacedim, typename Number>
2602template <bool is_linear, std::size_t stride_view>
2603inline void
2606 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2607{
2608 if (!(is_linear && n_components == 1))
2609 prepare_evaluate_fast<is_linear>(solution_values);
2610
2611 // loop over quadrature batches qb
2612 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2613
2614 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2615 {
2618
2619 compute_evaluate_fast<is_linear>(
2620 solution_values, evaluation_flags, n_shapes, qb, value, gradient);
2621
2622 if (evaluation_flags & EvaluationFlags::values)
2623 {
2624 for (unsigned int v = 0, offset = qb * stride;
2625 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2626 ++v, ++offset)
2627 ETT::set_value(value, v, this->values[offset]);
2628 }
2629 if (evaluation_flags & EvaluationFlags::gradients)
2630 {
2631 Assert(this->update_flags & update_gradients ||
2632 this->update_flags & update_inverse_jacobians,
2634
2635 for (unsigned int v = 0, offset = qb * stride;
2636 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2637 ++v, ++offset)
2638 {
2639 unit_gradient_type unit_gradient;
2640 ETT::set_gradient(gradient, v, unit_gradient);
2641 this->gradients[offset] =
2642 this->cell_type <=
2645 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
2647 this
2648 ->inverse_jacobian_ptr[this->cell_type <=
2650 GeometryType::affine ?
2651 0 :
2652 offset]
2653 .transpose(),
2654 unit_gradient);
2655 }
2656 }
2657 }
2658}
2659
2660
2661
2662template <int n_components_, int dim, int spacedim, typename Number>
2663template <std::size_t stride_view>
2664inline void
2667 const EvaluationFlags::EvaluationFlags &evaluation_flags)
2668{
2669 // slow path with FEValues
2670 Assert(this->fe_values.get() != nullptr,
2671 ExcMessage(
2672 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2673
2674 const std::size_t n_points = this->fe_values->get_quadrature().size();
2675
2676 if (evaluation_flags & EvaluationFlags::values)
2677 {
2678 this->values.resize(this->n_q_points);
2679 std::fill(this->values.begin(), this->values.end(), value_type());
2680 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2681 {
2682 const ScalarNumber value = solution_values[i];
2683 for (unsigned int d = 0; d < n_components; ++d)
2684 if (this->nonzero_shape_function_component[i][d] &&
2685 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2686 for (unsigned int qb = 0, q = 0; q < n_points;
2687 ++qb, q += n_lanes_user_interface)
2688 for (unsigned int v = 0;
2689 v < n_lanes_user_interface && q + v < n_points;
2690 ++v)
2691 ETT::access(this->values[qb],
2692 v,
2693 d,
2694 this->fe_values->shape_value(i, q + v) * value);
2695 else if (this->nonzero_shape_function_component[i][d])
2696 for (unsigned int qb = 0, q = 0; q < n_points;
2697 ++qb, q += n_lanes_user_interface)
2698 for (unsigned int v = 0;
2699 v < n_lanes_user_interface && q + v < n_points;
2700 ++v)
2701 ETT::access(this->values[qb],
2702 v,
2703 d,
2704 this->fe_values->shape_value_component(i,
2705 q + v,
2706 d) *
2707 value);
2708 }
2709 }
2710
2711 if (evaluation_flags & EvaluationFlags::gradients)
2712 {
2713 this->gradients.resize(this->n_q_points);
2714 std::fill(this->gradients.begin(),
2715 this->gradients.end(),
2716 gradient_type());
2717 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2718 {
2719 const ScalarNumber value = solution_values[i];
2720 for (unsigned int d = 0; d < n_components; ++d)
2721 if (this->nonzero_shape_function_component[i][d] &&
2722 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2723 for (unsigned int qb = 0, q = 0; q < n_points;
2724 ++qb, q += n_lanes_user_interface)
2725 for (unsigned int v = 0;
2726 v < n_lanes_user_interface && q + v < n_points;
2727 ++v)
2728 ETT::access(this->gradients[qb],
2729 v,
2730 d,
2731 this->fe_values->shape_grad(i, q + v) * value);
2732 else if (this->nonzero_shape_function_component[i][d])
2733 for (unsigned int qb = 0, q = 0; q < n_points;
2734 ++qb, q += n_lanes_user_interface)
2735 for (unsigned int v = 0;
2736 v < n_lanes_user_interface && q + v < n_points;
2737 ++v)
2738 ETT::access(
2739 this->gradients[qb],
2740 v,
2741 d,
2742 this->fe_values->shape_grad_component(i, q + v, d) * value);
2743 }
2744 }
2745}
2746
2747
2748
2749template <int n_components_, int dim, int spacedim, typename Number>
2750template <bool is_linear>
2751inline void
2753 const EvaluationFlags::EvaluationFlags &integration_flags,
2754 const unsigned int n_shapes,
2755 const unsigned int qb,
2756 const vectorized_value_type value,
2758 vectorized_value_type *solution_values_vectorized_linear)
2759{
2760 if (integration_flags & EvaluationFlags::gradients)
2762 is_linear,
2763 dim,
2765 vectorized_value_type>(this->shapes.data() + qb * n_shapes,
2766 n_shapes,
2767 &value,
2768 gradient,
2769 is_linear ?
2770 solution_values_vectorized_linear :
2771 this->solution_renumbered_vectorized.data(),
2772 this->unit_point_ptr[qb],
2773 qb != 0);
2774 else
2776 dim,
2779 this->shapes.data() + qb * n_shapes,
2780 n_shapes,
2781 value,
2782 is_linear ? solution_values_vectorized_linear :
2783 this->solution_renumbered_vectorized.data(),
2784 this->unit_point_ptr[qb],
2785 qb != 0);
2786}
2787
2788
2789
2790template <int n_components_, int dim, int spacedim, typename Number>
2791template <bool is_linear, std::size_t stride_view>
2792inline void
2794 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2795 vectorized_value_type *solution_values_vectorized_linear,
2796 const bool sum_into_values)
2797{
2798 if (!sum_into_values && this->fe->n_components() > n_components)
2799 for (unsigned int i = 0; i < solution_values.size(); ++i)
2800 solution_values[i] = 0;
2801
2802 const unsigned int dofs_per_comp =
2803 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
2804
2805 for (unsigned int comp = 0; comp < n_components; ++comp)
2806 {
2807 const std::size_t offset =
2808 (this->component_in_base_element + comp) * dofs_per_comp;
2809
2810 if (is_linear || this->renumber.empty())
2811 {
2812 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2813 if (sum_into_values)
2814 solution_values[i + offset] +=
2815 ETT::sum_value(comp,
2816 is_linear ?
2817 *(solution_values_vectorized_linear + i) :
2818 this->solution_renumbered_vectorized[i]);
2819 else
2820 solution_values[i + offset] =
2821 ETT::sum_value(comp,
2822 is_linear ?
2823 *(solution_values_vectorized_linear + i) :
2824 this->solution_renumbered_vectorized[i]);
2825 }
2826 else
2827 {
2828 const unsigned int *renumber_ptr = this->renumber.data() + offset;
2829 for (unsigned int i = 0; i < dofs_per_comp; ++i)
2830 if (sum_into_values)
2831 solution_values[renumber_ptr[i]] +=
2832 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2833 else
2834 solution_values[renumber_ptr[i]] =
2835 ETT::sum_value(comp, this->solution_renumbered_vectorized[i]);
2836 }
2837 }
2838}
2839
2840
2841
2842template <int n_components_, int dim, int spacedim, typename Number>
2843template <bool do_JxW, bool is_linear, std::size_t stride_view>
2844inline void
2846 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2847 const EvaluationFlags::EvaluationFlags &integration_flags,
2848 const bool sum_into_values)
2849{
2850 // zero out lanes of incomplete last quadrature point batch
2851 if constexpr (stride == 1)
2852 if (const unsigned int n_filled_lanes =
2853 this->n_q_points_scalar & (n_lanes_internal - 1);
2854 n_filled_lanes > 0)
2855 {
2856 if (integration_flags & EvaluationFlags::values)
2857 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2858 ETT::set_zero_value(this->values.back(), v);
2859 if (integration_flags & EvaluationFlags::gradients)
2860 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
2861 ETT::set_zero_gradient(this->gradients.back(), v);
2862 }
2863
2864 std::array<vectorized_value_type, is_linear ? Utilities::pow(2, dim) : 0>
2865 solution_values_vectorized_linear = {};
2866
2867 // loop over quadrature batches qb
2868 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
2869
2870 const bool cartesian_cell =
2872 const bool affine_cell =
2874 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
2875 {
2876 vectorized_value_type value = {};
2878
2879 if (integration_flags & EvaluationFlags::values)
2880 for (unsigned int v = 0, offset = qb * stride;
2881 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2882 ++v, ++offset)
2883 ETT::get_value(value,
2884 v,
2885 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
2886 this->values[offset]);
2887
2888 if (integration_flags & EvaluationFlags::gradients)
2889 for (unsigned int v = 0, offset = qb * stride;
2890 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
2891 ++v, ++offset)
2892 {
2893 const gradient_type grad_w =
2894 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
2895 this->gradients[offset];
2896 ETT::get_gradient(
2897 gradient,
2898 v,
2899 cartesian_cell ?
2900 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
2901 grad_w) :
2903 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
2904 grad_w));
2905 }
2906
2907 compute_integrate_fast<is_linear>(
2908 integration_flags,
2909 n_shapes,
2910 qb,
2911 value,
2912 gradient,
2913 solution_values_vectorized_linear.data());
2914 }
2915
2916 // add between the lanes and write into the result
2917 finish_integrate_fast<is_linear>(solution_values,
2918 solution_values_vectorized_linear.data(),
2919 sum_into_values);
2920}
2921
2922
2923
2924template <int n_components_, int dim, int spacedim, typename Number>
2925template <bool do_JxW, std::size_t stride_view>
2926inline void
2928 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
2929 const EvaluationFlags::EvaluationFlags &integration_flags,
2930 const bool sum_into_values)
2931{
2932 // slow path with FEValues
2933 Assert(this->fe_values.get() != nullptr,
2934 ExcMessage(
2935 "Not initialized. Please call FEPointEvaluation::reinit()!"));
2936 if (!sum_into_values)
2937 for (unsigned int i = 0; i < solution_values.size(); ++i)
2938 solution_values[i] = 0;
2939
2940 const std::size_t n_points = this->fe_values->get_quadrature().size();
2941
2942 if (integration_flags & EvaluationFlags::values)
2943 {
2944 AssertIndexRange(this->n_q_points, this->values.size() + 1);
2945 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2946 {
2947 for (unsigned int d = 0; d < n_components; ++d)
2948 if (this->nonzero_shape_function_component[i][d] &&
2949 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2950 for (unsigned int qb = 0, q = 0; q < n_points;
2951 ++qb, q += n_lanes_user_interface)
2952 for (unsigned int v = 0;
2953 v < n_lanes_user_interface && q + v < n_points;
2954 ++v)
2955 solution_values[i] +=
2956 this->fe_values->shape_value(i, q + v) *
2957 ETT::access(this->values[qb], v, d) *
2958 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2959 else if (this->nonzero_shape_function_component[i][d])
2960 for (unsigned int qb = 0, q = 0; q < n_points;
2961 ++qb, q += n_lanes_user_interface)
2962 for (unsigned int v = 0;
2963 v < n_lanes_user_interface && q + v < n_points;
2964 ++v)
2965 solution_values[i] +=
2966 this->fe_values->shape_value_component(i, q + v, d) *
2967 ETT::access(this->values[qb], v, d) *
2968 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2969 }
2970 }
2971
2972 if (integration_flags & EvaluationFlags::gradients)
2973 {
2974 AssertIndexRange(this->n_q_points, this->gradients.size() + 1);
2975 for (unsigned int i = 0; i < this->fe->n_dofs_per_cell(); ++i)
2976 {
2977 for (unsigned int d = 0; d < n_components; ++d)
2978 if (this->nonzero_shape_function_component[i][d] &&
2979 (this->fe->is_primitive(i) || this->fe->is_primitive()))
2980 for (unsigned int qb = 0, q = 0; q < n_points;
2981 ++qb, q += n_lanes_user_interface)
2982 for (unsigned int v = 0;
2983 v < n_lanes_user_interface && q + v < n_points;
2984 ++v)
2985 solution_values[i] +=
2986 this->fe_values->shape_grad(i, q + v) *
2987 ETT::access(this->gradients[qb], v, d) *
2988 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2989 else if (this->nonzero_shape_function_component[i][d])
2990 for (unsigned int qb = 0, q = 0; q < n_points;
2991 ++qb, q += n_lanes_user_interface)
2992 for (unsigned int v = 0;
2993 v < n_lanes_user_interface && q + v < n_points;
2994 ++v)
2995 solution_values[i] +=
2996 this->fe_values->shape_grad_component(i, q + v, d) *
2997 ETT::access(this->gradients[qb], v, d) *
2998 (do_JxW ? this->fe_values->JxW(q + v) : 1.);
2999 }
3000 }
3001}
3002
3003
3004
3005template <int n_components_, int dim, int spacedim, typename Number>
3006template <bool do_JxW, std::size_t stride_view>
3007void
3009 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3010 const EvaluationFlags::EvaluationFlags &integration_flags,
3011 const bool sum_into_values)
3012{
3013 if (this->must_reinitialize_pointers)
3014 reinit();
3015
3016 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3017
3018 if (this->n_q_points == 0 || // no evaluation points provided
3019 !((integration_flags & EvaluationFlags::values) ||
3020 (integration_flags &
3021 EvaluationFlags::gradients))) // no integration flags
3022 {
3023 if (!sum_into_values)
3024 for (unsigned int i = 0; i < solution_values.size(); ++i)
3025 solution_values[i] = 0;
3026 return;
3027 }
3028
3029 Assert(
3030 !do_JxW || this->JxW_ptr != nullptr,
3031 ExcMessage(
3032 "JxW pointer is not set! If you do not want to integrate() use test_and_sum()"));
3033
3034 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3035 if (this->fast_path)
3036 {
3037 if (this->use_linear_path)
3038 integrate_fast<do_JxW, true>(solution_values,
3039 integration_flags,
3040 sum_into_values);
3041 else
3042 integrate_fast<do_JxW, false>(solution_values,
3043 integration_flags,
3044 sum_into_values);
3045 }
3046 else
3047 integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
3048}
3049
3050
3051
3052template <int n_components_, int dim, int spacedim, typename Number>
3055 const unsigned int point_index) const
3056{
3057 AssertIndexRange(point_index, this->n_q_points);
3058 Assert(this->normal_ptr != nullptr,
3060 ExcFEPointEvaluationAccessToUninitializedMappingField(
3061 "update_normal_vectors"));
3062 if (this->is_interior)
3063 return this->normal_ptr[point_index];
3064 else
3065 return -this->normal_ptr[point_index];
3066}
3067
3068
3069
3070template <int n_components_, int dim, int spacedim, typename Number>
3075 const bool is_interior,
3076 const unsigned int first_selected_component)
3077 : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
3078 mapping_info,
3079 fe,
3080 first_selected_component,
3081 is_interior)
3082{}
3083
3084
3085
3086template <int n_components_, int dim, int spacedim, typename Number>
3087inline void
3089 const unsigned int cell_index,
3090 const unsigned int face_number)
3091{
3092 this->current_cell_index = cell_index;
3093 this->current_face_number = face_number;
3094 this->must_reinitialize_pointers = false;
3095
3096 if (this->use_linear_path)
3097 this->template do_reinit<true, true>();
3098 else
3099 this->template do_reinit<true, false>();
3100}
3101
3102
3103
3104template <int n_components_, int dim, int spacedim, typename Number>
3105inline void
3107 const unsigned int face_index)
3108{
3109 this->current_cell_index = face_index;
3110 this->current_face_number =
3111 this->mapping_info->get_face_number(face_index, this->is_interior);
3112 this->must_reinitialize_pointers = false;
3113
3114 if (this->use_linear_path)
3115 this->template do_reinit<true, true>();
3116 else
3117 this->template do_reinit<true, false>();
3118}
3119
3120
3121
3122template <int n_components_, int dim, int spacedim, typename Number>
3123template <std::size_t stride_view>
3124void
3127 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3128{
3129 Assert(!this->must_reinitialize_pointers,
3130 ExcMessage("Object has not been reinitialized!"));
3131
3132 if (this->n_q_points == 0)
3133 return;
3134
3135 Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
3136
3137 if (!((evaluation_flags & EvaluationFlags::values) ||
3138 (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
3139 return;
3140
3141 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3142
3143 if (this->use_linear_path)
3144 do_evaluate<true>(solution_values, evaluation_flags);
3145 else
3146 do_evaluate<false>(solution_values, evaluation_flags);
3147}
3148
3149
3150
3151template <int n_components_, int dim, int spacedim, typename Number>
3152void
3154 const ArrayView<const ScalarNumber> &solution_values,
3155 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3156{
3157 evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
3158 solution_values.size()),
3159 evaluation_flags);
3160}
3161
3162
3163
3164template <int n_components_, int dim, int spacedim, typename Number>
3165template <bool is_linear, std::size_t stride_view>
3166void
3169 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3170{
3171 const unsigned int dofs_per_comp =
3172 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3173
3174 const ScalarNumber *input;
3175 if (stride_view == 1 && this->component_in_base_element == 0 &&
3176 (is_linear || this->renumber.empty()))
3177 input = solution_values.data();
3178 else
3179 {
3180 for (unsigned int comp = 0; comp < n_components; ++comp)
3181 {
3182 const std::size_t offset =
3183 (this->component_in_base_element + comp) * dofs_per_comp;
3184
3185 if (is_linear || this->renumber.empty())
3186 {
3187 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3188 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3189 solution_values[i + offset];
3190 }
3191 else
3192 {
3193 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3194 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3195 this->scratch_data_scalar[i + comp * dofs_per_comp] =
3196 solution_values[renumber_ptr[i]];
3197 }
3198 }
3199 input = this->scratch_data_scalar.data();
3200 }
3201
3202 ScalarNumber *output =
3203 this->scratch_data_scalar.begin() + dofs_per_comp * n_components;
3204
3205 internal::FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3206 template interpolate<true, false>(n_components,
3207 evaluation_flags,
3208 this->shape_info,
3209 input,
3210 output,
3211 this->current_face_number);
3212
3213 do_evaluate_in_face<is_linear, 1>(output, evaluation_flags);
3214}
3215
3216
3217
3218template <int n_components_, int dim, int spacedim, typename Number>
3219template <std::size_t stride_view>
3220void
3222 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3223 const EvaluationFlags::EvaluationFlags &integration_flags,
3224 const bool sum_into_values)
3225{
3226 Assert(!this->must_reinitialize_pointers,
3227 ExcMessage("Object has not been reinitialized!"));
3228
3229 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3230
3231 if (this->n_q_points == 0 || // no evaluation points provided
3232 !((integration_flags & EvaluationFlags::values) ||
3233 (integration_flags &
3234 EvaluationFlags::gradients))) // no integration flags
3235 {
3236 if (!sum_into_values)
3237 for (unsigned int i = 0; i < solution_values.size(); ++i)
3238 solution_values[i] = 0;
3239 return;
3240 }
3241
3242 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3243
3244 if (this->use_linear_path)
3245 do_integrate<true, true>(solution_values,
3246 integration_flags,
3247 sum_into_values);
3248 else
3249 do_integrate<true, false>(solution_values,
3250 integration_flags,
3251 sum_into_values);
3252}
3253
3254
3255
3256template <int n_components_, int dim, int spacedim, typename Number>
3257void
3259 const ArrayView<ScalarNumber> &solution_values,
3260 const EvaluationFlags::EvaluationFlags &integration_flags,
3261 const bool sum_into_values)
3262{
3263 integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3264 solution_values.size()),
3265 integration_flags,
3266 sum_into_values);
3267}
3268
3269
3270
3271template <int n_components_, int dim, int spacedim, typename Number>
3272template <std::size_t stride_view>
3273void
3275 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3276 const EvaluationFlags::EvaluationFlags &integration_flags,
3277 const bool sum_into_values)
3278{
3279 Assert(!this->must_reinitialize_pointers,
3280 ExcMessage("Object has not been reinitialized!"));
3281
3282 Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
3283
3284 if (this->n_q_points == 0 || // no evaluation points provided
3285 !((integration_flags & EvaluationFlags::values) ||
3286 (integration_flags &
3287 EvaluationFlags::gradients))) // no integration flags
3288 {
3289 if (!sum_into_values)
3290 for (unsigned int i = 0; i < solution_values.size(); ++i)
3291 solution_values[i] = 0;
3292 return;
3293 }
3294
3295 AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
3296
3297 if (this->use_linear_path)
3298 do_integrate<false, true>(solution_values,
3299 integration_flags,
3300 sum_into_values);
3301 else
3302 do_integrate<false, false>(solution_values,
3303 integration_flags,
3304 sum_into_values);
3305}
3306
3307
3308
3309template <int n_components_, int dim, int spacedim, typename Number>
3310void
3312 const ArrayView<ScalarNumber> &solution_values,
3313 const EvaluationFlags::EvaluationFlags &integration_flags,
3314 const bool sum_into_values)
3315{
3316 test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
3317 solution_values.size()),
3318 integration_flags,
3319 sum_into_values);
3320}
3321
3322
3323
3324template <int n_components_, int dim, int spacedim, typename Number>
3325template <bool do_JxW, bool is_linear, std::size_t stride_view>
3326void
3328 const StridedArrayView<ScalarNumber, stride_view> &solution_values,
3329 const EvaluationFlags::EvaluationFlags &integration_flags,
3330 const bool sum_into_values)
3331{
3332 if (!sum_into_values && this->fe->n_components() > n_components)
3333 for (unsigned int i = 0; i < solution_values.size(); ++i)
3334 solution_values[i] = 0;
3335
3336 do_integrate_in_face<do_JxW, is_linear, 1>(this->scratch_data_scalar.begin(),
3337 integration_flags,
3338 false);
3339
3340 ScalarNumber *input = this->scratch_data_scalar.begin();
3341
3342 if (stride_view == 1 && this->component_in_base_element == 0 &&
3343 (is_linear || this->renumber.empty()))
3344 {
3345 if (sum_into_values)
3346 internal::
3347 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3348 template interpolate<false, true>(n_components,
3349 integration_flags,
3350 this->shape_info,
3351 input,
3352 solution_values.data(),
3353 this->current_face_number);
3354 else
3355 internal::
3356 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3357 template interpolate<false, false>(n_components,
3358 integration_flags,
3359 this->shape_info,
3360 input,
3361 solution_values.data(),
3362 this->current_face_number);
3363 }
3364 else
3365 {
3366 const unsigned int dofs_per_comp_face =
3367 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3368
3369 const unsigned int size_input = 3 * dofs_per_comp_face * n_components;
3370 ScalarNumber *output = input + size_input;
3371
3372 internal::
3373 FEFaceNormalEvaluationImpl<dim, is_linear ? 1 : -1, ScalarNumber>::
3374 template interpolate<false, false>(n_components,
3375 integration_flags,
3376 this->shape_info,
3377 input,
3378 output,
3379 this->current_face_number);
3380
3381 const unsigned int dofs_per_comp =
3382 is_linear ? Utilities::pow(2, dim) : this->dofs_per_component;
3383
3384 for (unsigned int comp = 0; comp < n_components; ++comp)
3385 {
3386 const std::size_t offset =
3387 (this->component_in_base_element + comp) * dofs_per_comp;
3388
3389 if (is_linear || this->renumber.empty())
3390 {
3391 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3392 if (sum_into_values)
3393 solution_values[i + offset] +=
3394 output[i + comp * dofs_per_comp];
3395 else
3396 solution_values[i + offset] =
3397 output[i + comp * dofs_per_comp];
3398 }
3399 else
3400 {
3401 const unsigned int *renumber_ptr = this->renumber.data() + offset;
3402 for (unsigned int i = 0; i < dofs_per_comp; ++i)
3403 if (sum_into_values)
3404 solution_values[renumber_ptr[i]] +=
3405 output[i + comp * dofs_per_comp];
3406 else
3407 solution_values[renumber_ptr[i]] =
3408 output[i + comp * dofs_per_comp];
3409 }
3410 }
3411 }
3412}
3413
3414
3415
3416template <int n_components_, int dim, int spacedim, typename Number>
3417template <int stride_face_dof>
3418void
3420 const ScalarNumber *face_dof_values,
3421 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3422{
3423 if (this->use_linear_path)
3424 do_evaluate_in_face<true, stride_face_dof>(face_dof_values,
3425 evaluation_flags);
3426 else
3427 do_evaluate_in_face<false, stride_face_dof>(face_dof_values,
3428 evaluation_flags);
3429}
3430
3431
3432
3433template <int n_components_, int dim, int spacedim, typename Number>
3434template <bool is_linear, int stride_face_dof>
3435inline void
3437 do_evaluate_in_face(const ScalarNumber *face_dof_values,
3438 const EvaluationFlags::EvaluationFlags &evaluation_flags)
3439{
3440 const scalar_value_type *face_dof_values_ptr;
3441 if constexpr (n_components == 1)
3442 face_dof_values_ptr = face_dof_values;
3443 else
3444 {
3445 const unsigned int dofs_per_comp_face =
3446 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3447 for (unsigned int comp = 0; comp < n_components; ++comp)
3448 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3449 ETT::read_value(face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3450 stride_face_dof],
3451 comp,
3452 this->solution_renumbered[i]);
3453
3454 face_dof_values_ptr = this->solution_renumbered.data();
3455 }
3456
3457 constexpr int stride_face_dof_actual =
3458 n_components == 1 ? stride_face_dof : 1;
3459
3460 // loop over quadrature batches qb
3461 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3462
3463 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3464 {
3467
3468 if (evaluation_flags & EvaluationFlags::gradients)
3469 {
3470 const std::array<vectorized_value_type, dim + 1> interpolated_value =
3471 is_linear ?
3473 dim - 1,
3476 2,
3477 stride_face_dof_actual>(face_dof_values_ptr,
3478 this->unit_point_faces_ptr[qb]) :
3480 dim - 1,
3483 2,
3484 false,
3485 stride_face_dof_actual>(this->shapes_faces.data() +
3486 qb * n_shapes,
3487 n_shapes,
3488 face_dof_values_ptr);
3489
3490 value = interpolated_value[dim - 1];
3491 // reorder derivative from tangential/normal derivatives into tensor
3492 // in physical coordinates
3493 if (this->current_face_number / 2 == 0)
3494 {
3495 gradient[0] = interpolated_value[dim];
3496 if (dim > 1)
3497 gradient[1] = interpolated_value[0];
3498 if (dim > 2)
3499 gradient[2] = interpolated_value[1];
3500 }
3501 else if (this->current_face_number / 2 == 1)
3502 {
3503 if (dim > 1)
3504 gradient[1] = interpolated_value[dim];
3505 if (dim == 3)
3506 {
3507 gradient[0] = interpolated_value[1];
3508 gradient[2] = interpolated_value[0];
3509 }
3510 else if (dim == 2)
3511 gradient[0] = interpolated_value[0];
3512 else
3514 }
3515 else if (this->current_face_number / 2 == 2)
3516 {
3517 if (dim > 2)
3518 {
3519 gradient[0] = interpolated_value[0];
3520 gradient[1] = interpolated_value[1];
3521 gradient[2] = interpolated_value[dim];
3522 }
3523 else
3525 }
3526 else
3528 }
3529 else
3530 {
3531 value = is_linear ?
3533 dim - 1,
3536 stride_face_dof_actual>(face_dof_values_ptr,
3537 this->unit_point_faces_ptr[qb]) :
3539 dim - 1,
3542 false,
3543 stride_face_dof_actual>(this->shapes_faces.data() +
3544 qb * n_shapes,
3545 n_shapes,
3546 face_dof_values_ptr);
3547 }
3548
3549 if (evaluation_flags & EvaluationFlags::values)
3550 {
3551 for (unsigned int v = 0, offset = qb * stride;
3552 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3553 ++v, ++offset)
3554 ETT::set_value(value, v, this->values[offset]);
3555 }
3556 if (evaluation_flags & EvaluationFlags::gradients)
3557 {
3558 Assert(this->update_flags & update_gradients ||
3559 this->update_flags & update_inverse_jacobians,
3561
3562 for (unsigned int v = 0, offset = qb * stride;
3563 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3564 ++v, ++offset)
3565 {
3566 unit_gradient_type unit_gradient;
3567 ETT::set_gradient(gradient, v, unit_gradient);
3568 this->gradients[offset] =
3569 this->cell_type <=
3572 this->inverse_jacobian_ptr[0].transpose(), unit_gradient) :
3574 this
3575 ->inverse_jacobian_ptr[this->cell_type <=
3577 GeometryType::affine ?
3578 0 :
3579 offset]
3580 .transpose(),
3581 unit_gradient);
3582 }
3583 }
3584 }
3585}
3586
3587
3588
3589template <int n_components_, int dim, int spacedim, typename Number>
3590template <int stride_face_dof>
3591void
3593 ScalarNumber *face_dof_values,
3594 const EvaluationFlags::EvaluationFlags &integration_flags,
3595 const bool sum_into_values)
3596{
3597 if (this->use_linear_path)
3598 do_integrate_in_face<true, true, stride_face_dof>(face_dof_values,
3599 integration_flags,
3600 sum_into_values);
3601 else
3602 do_integrate_in_face<true, false, stride_face_dof>(face_dof_values,
3603 integration_flags,
3604 sum_into_values);
3605}
3606
3607
3608
3609template <int n_components_, int dim, int spacedim, typename Number>
3610template <bool do_JxW, bool is_linear, int stride_face_dof>
3611inline void
3614 ScalarNumber *face_dof_values,
3615 const EvaluationFlags::EvaluationFlags &integration_flags,
3616 const bool sum_into_values)
3617{
3618 // zero out lanes of incomplete last quadrature point batch
3619 if constexpr (stride == 1)
3620 if (const unsigned int n_filled_lanes =
3621 this->n_q_points_scalar & (n_lanes_internal - 1);
3622 n_filled_lanes > 0)
3623 {
3624 if (integration_flags & EvaluationFlags::values)
3625 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3626 ETT::set_zero_value(this->values.back(), v);
3627 if (integration_flags & EvaluationFlags::gradients)
3628 for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
3629 ETT::set_zero_gradient(this->gradients.back(), v);
3630 }
3631
3632 std::array<vectorized_value_type,
3633 is_linear ? 2 * Utilities::pow(2, dim - 1) : 0>
3634 solution_values_vectorized_linear = {};
3635
3636 // loop over quadrature batches qb
3637 const unsigned int n_shapes = is_linear ? 2 : this->poly.size();
3638
3639 const bool cartesian_cell =
3641 const bool affine_cell =
3643 for (unsigned int qb = 0; qb < this->n_q_batches; ++qb)
3644 {
3645 vectorized_value_type value = {};
3647
3648 if (integration_flags & EvaluationFlags::values)
3649 for (unsigned int v = 0, offset = qb * stride;
3650 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3651 ++v, ++offset)
3652 ETT::get_value(value,
3653 v,
3654 do_JxW ? this->values[offset] * this->JxW_ptr[offset] :
3655 this->values[offset]);
3656
3657 if (integration_flags & EvaluationFlags::gradients)
3658 for (unsigned int v = 0, offset = qb * stride;
3659 v < stride && (stride == 1 || offset < this->n_q_points_scalar);
3660 ++v, ++offset)
3661 {
3662 const auto grad_w =
3663 do_JxW ? this->gradients[offset] * this->JxW_ptr[offset] :
3664 this->gradients[offset];
3665 ETT::get_gradient(
3666 gradient,
3667 v,
3668 cartesian_cell ?
3669 apply_diagonal_transformation(this->inverse_jacobian_ptr[0],
3670 grad_w) :
3672 this->inverse_jacobian_ptr[affine_cell ? 0 : offset],
3673 grad_w));
3674 }
3675
3676 if (integration_flags & EvaluationFlags::gradients)
3677 {
3678 std::array<vectorized_value_type, 2> value_face = {};
3679 Tensor<1, dim - 1, vectorized_value_type> gradient_in_face;
3680
3681 value_face[0] = value;
3682 // fill derivative in physical coordinates into tangential/normal
3683 // derivatives
3684 if (this->current_face_number / 2 == 0)
3685 {
3686 value_face[1] = gradient[0];
3687 if (dim > 1)
3688 gradient_in_face[0] = gradient[1];
3689 if (dim > 2)
3690 gradient_in_face[1] = gradient[2];
3691 }
3692 else if (this->current_face_number / 2 == 1)
3693 {
3694 if (dim > 1)
3695 value_face[1] = gradient[1];
3696 if (dim == 3)
3697 {
3698 gradient_in_face[0] = gradient[2];
3699 gradient_in_face[1] = gradient[0];
3700 }
3701 else if (dim == 2)
3702 gradient_in_face[0] = gradient[0];
3703 else
3705 }
3706 else if (this->current_face_number / 2 == 2)
3707 {
3708 if (dim > 2)
3709 {
3710 value_face[1] = gradient[2];
3711 gradient_in_face[0] = gradient[0];
3712 gradient_in_face[1] = gradient[1];
3713 }
3714 else
3716 }
3717 else
3719
3721 is_linear,
3722 dim - 1,
3725 2>(this->shapes_faces.data() + qb * n_shapes,
3726 n_shapes,
3727 value_face.data(),
3728 gradient_in_face,
3729 is_linear ? solution_values_vectorized_linear.data() :
3730 this->solution_renumbered_vectorized.data(),
3731 this->unit_point_faces_ptr[qb],
3732 qb != 0);
3733 }
3734 else
3736 dim - 1,
3739 this->shapes_faces.data() + qb * n_shapes,
3740 n_shapes,
3741 value,
3742 is_linear ? solution_values_vectorized_linear.data() :
3743 this->solution_renumbered_vectorized.data(),
3744 this->unit_point_faces_ptr[qb],
3745 qb != 0);
3746 }
3747
3748 const unsigned int dofs_per_comp_face =
3749 is_linear ? Utilities::pow(2, dim - 1) : this->dofs_per_component_face;
3750
3751 for (unsigned int comp = 0; comp < n_components; ++comp)
3752 for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
3753 if (sum_into_values)
3754 face_dof_values[(i + comp * 3 * dofs_per_comp_face) *
3755 stride_face_dof] +=
3756 ETT::sum_value(comp,
3757 is_linear ?
3758 *(solution_values_vectorized_linear.data() + i) :
3759 this->solution_renumbered_vectorized[i]);
3760 else
3761 face_dof_values[(i + comp * 3 * dofs_per_comp_face) * stride_face_dof] =
3762 ETT::sum_value(comp,
3763 is_linear ?
3764 *(solution_values_vectorized_linear.data() + i) :
3765 this->solution_renumbered_vectorized[i]);
3766}
3767
3768
3769
3770template <int n_components_, int dim, int spacedim, typename Number>
3773 const unsigned int point_index) const
3774{
3775 AssertIndexRange(point_index, this->n_q_points);
3776 Assert(this->normal_ptr != nullptr,
3778 ExcFEPointEvaluationAccessToUninitializedMappingField(
3779 "update_normal_vectors"));
3780 if (this->cell_type <= ::internal::MatrixFreeFunctions::affine)
3781 {
3783 for (unsigned int d = 0; d < dim; ++d)
3784 normal[d] =
3785 internal::VectorizedArrayTrait<Number>::get(this->normal_ptr[0][d],
3786 0);
3787 if (this->is_interior)
3788 return normal;
3789 else
3790 return -normal;
3791 }
3792 else
3793 {
3794 if (this->is_interior)
3795 return this->normal_ptr[point_index];
3796 else
3797 return -(this->normal_ptr[point_index]);
3798 }
3799}
3800
3802
3803#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
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void do_evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
static constexpr std::size_t stride
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
void do_evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr unsigned int n_components
static constexpr std::size_t n_lanes_user_interface
typename ETT::real_gradient_type gradient_type
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)
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
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
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
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
const unsigned int n_q_points_scalar
const Point< dim, VectorizedArrayType > * unit_point_ptr
FEPointEvaluationBase(FEPointEvaluationBase &&other) noexcept
AlignedVector< ScalarNumber > scratch_data_scalar
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
const value_type & get_value(const unsigned int point_index) const
const Point< dim - 1, VectorizedArrayType > * unit_point_faces_ptr
std::vector< scalar_value_type > solution_renumbered
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
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 ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::vectorized_value_type vectorized_value_type
static constexpr unsigned int n_components
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
FEPointEvaluationBase(FEPointEvaluationBase &other) noexcept
DerivativeForm< 1, dim, spacedim, Number > jacobian(const unsigned int point_index) const
typename ETT::real_gradient_type gradient_type
const Point< spacedim, Number > * real_point_ptr
unsigned int component_in_base_element
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 internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static constexpr std::size_t stride
void integrate_slow(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr unsigned int dimension
typename ETT::value_type value_type
void do_integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
void evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
typename ETT::scalar_value_type scalar_value_type
typename ETT::unit_gradient_type unit_gradient_type
void evaluate(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void prepare_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values)
void test_and_sum(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void integrate(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
typename ETT::vectorized_value_type vectorized_value_type
void compute_integrate_fast(const EvaluationFlags::EvaluationFlags &integration_flags, const unsigned int n_shapes, const unsigned int qb, const vectorized_value_type value, const interface_vectorized_unit_gradient_type gradient, vectorized_value_type *solution_values_vectorized_linear)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
typename ETT::real_gradient_type gradient_type
FEPointEvaluation(const Mapping< dim, spacedim > &mapping, const FiniteElement< dim, spacedim > &fe, const UpdateFlags update_flags, const unsigned int first_selected_component=0)
void evaluate_slow(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
static constexpr std::size_t n_lanes_internal
static constexpr unsigned int n_components
void finish_integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, vectorized_value_type *solution_values_vectorized_linear, const bool sum_into_values)
typename internal::FEPointEvaluation::EvaluatorTypeTraits< dim, spacedim, n_components, Number > ETT
void compute_evaluate_fast(const StridedArrayView< const ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &evaluation_flags, const unsigned int n_shapes, const unsigned int qb, vectorized_value_type &value, interface_vectorized_unit_gradient_type &gradient)
void integrate_fast(const StridedArrayView< ScalarNumber, stride_view > &solution_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values)
static constexpr std::size_t n_lanes_user_interface
typename ETT::interface_vectorized_unit_gradient_type interface_vectorized_unit_gradient_type
Tensor< 1, spacedim, Number > normal_vector(const unsigned int point_index) const
Abstract base class for mapping classes.
Definition mapping.h: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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, dim, ScalarNumber > &shape_gradient)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int, const vectorized_unit_gradient_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int vector_lane, scalar_unit_gradient_type &result)
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static scalar_value_type sum_value(const vectorized_value_type &result)
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
static void set_gradient(const vectorized_unit_gradient_type &value, const unsigned int, vectorized_unit_gradient_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
static void access(value_type &value, const unsigned int vector_lane, const unsigned int, const ScalarNumber &shape_value)
static void get_gradient(vectorized_unit_gradient_type &value, const unsigned int vector_lane, const scalar_unit_gradient_type &result)
static ScalarNumber sum_value(const unsigned int, const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int, scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static Tensor< 1, spacedim, ScalarNumber > access(const real_gradient_type &value, const unsigned int vector_lane, const unsigned int component)
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const DerivativeForm< 1, dim, n_components, Number > &result)
static void get_value(vectorized_value_type &value, const unsigned int, const vectorized_value_type &result)
static void get_value(vectorized_value_type &value, const unsigned int vector_lane, const scalar_value_type &result)
static void set_value(const vectorized_value_type &value, const unsigned int, vectorized_value_type &result)
static scalar_value_type sum_value(const scalar_value_type &result)
typename ::internal::VectorizedArrayTrait< Number >::vectorized_value_type VectorizedArrayType
static scalar_value_type sum_value(const vectorized_value_type &result)
static void read_value(const ScalarNumber vector_entry, const unsigned int component, scalar_value_type &result)
Tensor< 1, n_components, ScalarNumber > scalar_value_type
static void set_value(const vectorized_value_type &value, const unsigned int vector_lane, scalar_value_type &result)
static ScalarNumber access(const value_type &value, const unsigned int vector_lane, const unsigned int component)
static ScalarNumber sum_value(const unsigned int component, const vectorized_value_type &result)
Tensor< 1, n_components, Tensor< 1, dim, VectorizedArrayType > > vectorized_unit_gradient_type
static void set_zero_value(value_type &value, const unsigned int vector_lane)
static void set_gradient(const interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, unit_gradient_type &result)
static void set_zero_gradient(real_gradient_type &value, const unsigned int vector_lane)
typename internal::VectorizedArrayTrait< Number >::value_type ScalarNumber
static void get_gradient(interface_vectorized_unit_gradient_type &value, const unsigned int vector_lane, const unit_gradient_type &result)
std::conditional_t< n_components==spacedim, Tensor< 2, spacedim, Number >, Tensor< 1, n_components, Tensor< 1, spacedim, Number > > > real_gradient_type
Tensor< 1, n_components, VectorizedArrayType > vectorized_value_type
static void access(real_gradient_type &value, const unsigned int vector_lane, const unsigned int component, const Tensor< 1, spacedim, ScalarNumber > &shape_gradient)
Tensor< 1, n_components, Tensor< 1, dim, Number > > unit_gradient_type
static void access(value_type &value, const unsigned int vector_lane, const unsigned int component, const ScalarNumber &shape_value)
static constexpr std::size_t width()
static constexpr std::size_t stride()
static value_type & get(value_type &value, unsigned int c)